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-11) – Can be set to determine the relative tolerance convergence criteria (defailt 1.0e-11)
verbose (bool) – Should be set to True if verbse debug output is desired (default False)
method (str, optional, defualt='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 statisitcal 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
’dDelta_f’ (float) – Estimated standard deviation of free energy difference
Notes
If you are prodividing 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 statisitcal 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