Other estimators

pymbar implements other reweighting estimators, specifically the Bennett Acceptance Ratio, exponential averaging, and a Gaussian approximation to exponential averaging.

Please reference the following if you use this code in your research:

[1] Shirts MR and Chodera JD. Statistically optimal analysis of samples from multiple equilibrium states. J. Chem. Phys. 129:124105, 2008. http://dx.doi.org/10.1063/1.2978177

This module contains implementations of

  • bar - bidirectional estimator for free energy differences / Bennett acceptance ratio estimator

  • exp - unidirectional estimator for free energy differences based on Zwanzig relation / exponential averaging

  • exp_gauss - unidirectional estimator for free energy differences based on Zwanzig relation / exponential averaging, assuming the distribution is Gaussian.

pymbar.other_estimators.bar(w_F, w_R, DeltaF=0.0, compute_uncertainty=True, uncertainty_method='BAR', maximum_iterations=500, relative_tolerance=1e-12, verbose=False, method='false-position', iterated_solution=True)

Compute free energy difference using the Bennett acceptance ratio (BAR) method.

Parameters:
  • w_F (np.ndarray) – w_F[t] is the forward work value from snapshot t. t = 0…(T_F-1) Length T_F is deduced from vector.

  • w_R (np.ndarray) – w_R[t] is the reverse work value from snapshot t. t = 0…(T_R-1) Length T_R is deduced from vector.

  • DeltaF (float, optional, default=0.0) – DeltaF can be set to initialize the free energy difference with a guess

  • compute_uncertainty (bool, optional, default=True) – if False, only the free energy is returned

  • uncertainty_method (string, optional, default=''BAR'') – There are two possible uncertainty estimates for BAR. One agrees with MBAR for two states exactly, and is indicated by “MBAR”. The other estimator, which is the one originally derived for BAR, only agrees with MBAR in the limit of good overlap, and is designated ‘BAR’ See code comments below for derivations of the two methods.

  • maximum_iterations (int, optional, default=500) – Can be set to limit the maximum number of iterations performed

  • relative_tolerance (float, optional, default=1E-12) – Can be set to determine the relative tolerance convergence criteria (default 1.0e-12)

  • verbose (bool) – Should be set to True if verbse debug output is desired (default False)

  • method (str, optional, default='false-position') – Choice of method to solve bar nonlinear equations: one of ‘bisection’, ‘self-consistent-iteration’ or ‘false-position’ (default : ‘false-position’).

  • iterated_solution (bool, optional, default=True) – whether to fully solve the optimized bar equation to consistency, or to stop after one step, to be equivalent to transition matrix sampling.

Returns:

‘Delta_f’float

Free energy difference

’dDelta_f’float

Estimated standard deviation of free energy difference

Return type:

dict

References

[1] Shirts MR, Bair E, Hooker G, and Pande VS. Equilibrium free energies from nonequilibrium measurements using maximum-likelihood methods. PRL 91(14):140601, 2003.

Notes

The false position method is used to solve the implicit equation.

Examples

Compute free energy difference between two specified samples of work values.

>>> from pymbar import testsystems
>>> [w_F, w_R] = testsystems.gaussian_work_example(mu_F=None, DeltaF=1.0, seed=0)
>>> results = bar(w_F, w_R)
>>> print('Free energy difference is {:.3f} +- {:.3f} kT'.format(results['Delta_f'], results['dDelta_f']))
Free energy difference is 1.088 +- 0.050 kT

Test completion of various other schemes.

>>> results = bar(w_F, w_R, method='self-consistent-iteration')
>>> results = bar(w_F, w_R, method='false-position')
>>> results = bar(w_F, w_R, method='bisection')
pymbar.other_estimators.bar_overlap(w_F, w_R)

Compute overlap between forward and backward ensembles (using MBAR definition of overlap)

Parameters:
  • w_F (np.ndarray) – w_F[t] is the forward work value from snapshot t. t = 0…(T_F-1) Length T_F is deduced from vector.

  • w_R (np.ndarray) – w_R[t] is the reverse work value from snapshot t. t = 0…(T_R-1) Length T_R is deduced from vector.

Returns:

overlap – The overlap: 0 denotes no overlap, 1 denotes complete overlap

Return type:

float

pymbar.other_estimators.bar_zero(w_F, w_R, DeltaF)

A function that when zeroed is equivalent to the solution of the Bennett acceptance ratio.

from http://journals.aps.org/prl/pdf/10.1103/PhysRevLett.91.140601

D_F = M + w_F - Delta F D_R = M + w_R - Delta F

we want:

sum_N_F (1+exp(D_F))^-1 = sum N_R N_R <(1+exp(-D_R))^-1> ln sum N_F (1+exp(D_F))^-1>_F = ln sum N_R exp((1+exp(-D_R))^(-1)>_R ln sum N_F (1+exp(D_F))^-1>_F - ln sum N_R exp((1+exp(-D_R))^(-1)>_R = 0

Parameters:
  • w_F (np.ndarray) – w_F[t] is the forward work value from snapshot t. t = 0…(T_F-1) Length T_F is deduced from vector.

  • w_R (np.ndarray) – w_R[t] is the reverse work value from snapshot t. t = 0…(T_R-1) Length T_R is deduced from vector.

  • DeltaF (float) – Our current guess

Returns:

fzero – a variable that is zeroed when DeltaF satisfies bar.

Return type:

float

Examples

Compute free energy difference between two specified samples of work values.

>>> from pymbar import testsystems
>>> [w_F, w_R] = testsystems.gaussian_work_example(mu_F=None, DeltaF=1.0, seed=0)
>>> DeltaF = bar_zero(w_F, w_R, 0.0)
pymbar.other_estimators.exp(w_F, compute_uncertainty=True, is_timeseries=False)

Estimate free energy difference using one-sided (unidirectional) exponential averaging (EXP).

Parameters:
  • w_F (np.ndarray, float) – w_F[t] is the forward work value from snapshot t. t = 0…(T-1) Length T is deduced from vector.

  • compute_uncertainty (bool, optional, default=True) – if False, will disable computation of the statistical uncertainty (default: True)

  • is_timeseries (bool, default=False) – if True, correlation in data is corrected for by estimation of statistical inefficiency (default: False) Use this option if you are providing correlated timeseries data and have not subsampled the data to produce uncorrelated samples.

Returns:

dict_vals – Dictionary with keys Delta_f and dDelta_f for the free energy difference and its estimated deviation, respectively.

Return type:

dict[float]

Notes

If you are providing correlated timeseries data, be sure to set the ‘timeseries’ flag to True

Examples

Compute the free energy difference given a sample of forward work values.

>>> from pymbar import testsystems
>>> [w_F, w_R] = testsystems.gaussian_work_example(mu_F=None, DeltaF=1.0, seed=0)
>>> results = exp(w_F)
>>> print('Forward free energy difference is {:.3f} +- {:.3f} kT'.format(results['Delta_f'], results['dDelta_f']))
Forward free energy difference is 1.088 +- 0.076 kT
>>> results = exp(w_R)
>>> print('Reverse free energy difference is {:.3f} +- {:.3f} kT'.format(results['Delta_f'], results['dDelta_f']))
Reverse free energy difference is -1.073 +- 0.082 kT
pymbar.other_estimators.exp_gauss(w_F, compute_uncertainty=True, is_timeseries=False)

Estimate free energy difference using gaussian approximation to one-sided (unidirectional) exponential averaging.

Parameters:
  • w_F (np.ndarray, float) – w_F[t] is the forward work value from snapshot t. t = 0…(T-1) Length T is deduced from vector.

  • compute_uncertainty (bool, optional, default=True) – if False, will disable computation of the statistical uncertainty (default: True)

  • is_timeseries (bool, default=False) – if True, correlation in data is corrected for by estimation of statistical inefficiency (default: False) Use this option if you are providing correlated timeseries data and have not subsampled the data to produce uncorrelated samples.

Returns:

‘Delta_f’float

Free energy difference between the two states

’dDelta_f’float

Estimated standard deviation of free energy difference between the two states

Return type:

Results dictionary with keys

Notes

If you are providing correlated timeseries data, be sure to set the ‘timeseries’ flag to True

Examples

Compute the free energy difference given a sample of forward work values.

>>> from pymbar import testsystems
>>> [w_F, w_R] = testsystems.gaussian_work_example(mu_F=None, DeltaF=1.0, seed=0)
>>> results = exp_gauss(w_F)
>>> print('Forward Gaussian approximated free energy difference is {:.3f} +- {:.3f} kT'.format(results['Delta_f'], results['dDelta_f']))
Forward Gaussian approximated free energy difference is 1.049 +- 0.089 kT
>>> results = exp_gauss(w_R)
>>> print('Reverse Gaussian approximated free energy difference is {:.3f} +- {:.3f} kT'.format(results['Delta_f'], results['dDelta_f']))
Reverse Gaussian approximated free energy difference is -1.073 +- 0.080 kT