pymbar
Python implementation of the multistate Bennett acceptance ratio (MBAR) method for estimating expectations and free energy differences from equilibrium samples from multiple probability densities.
Installation
The easiest way to install the pymbar release is via conda:
conda install -c conda-forge pymbar
You can also install pymbar from the Python package index using pip:
pip install pymbar
The development version can be installed directly from github via pip:
pip install git+https://github.com/choderalab/pymbar.git
Usage
Basic usage involves importing pymbar and constructing an MBAR object from the reduced potential of simulation or experimental data.
Suppose we sample a 1D harmonic oscillator from a few thermodynamic states:
>>> from pymbar import testsystems
>>> [x_n, u_kn, N_k, s_n] = testsystems.HarmonicOscillatorsTestCase().sample()
We have the nsamples sampled oscillator positions x_n (with samples from all states concatenated), reduced potentials in the (nstates,nsamples) matrix u_kn, number of samples per state in the nsamples array N_k, and indices s_n denoting which thermodynamic state each sample was drawn from.
To analyze this data, we first initialize the MBAR object:
>>> mbar = MBAR(u_kn, N_k)
Estimating dimensionless free energy differences between the sampled thermodynamic states and their associated uncertainties (standard errors) simply requires a call to getFreeEnergyDifferences():
>>> results = mbar.getFreeEnergyDifferences(return_dict=True)
Here results is a dictionary with keys Deltaf_ij, dDeltaf, and Theta. Deltaf_ij[i,j] is the matrix of dimensionless free energy differences f_j - f_i, dDeltaf_ij[i,j] is the matrix of standard errors in this matrices estimate, and Theta is a covariance matrix that can be used to propagate error into quantities derived from the free energies.
By default, return_dict is False and the return will be a tuple.
Expectations and associated uncertainties can easily be estimated for observables A(x) for all states:
>>> A_kn = x_kn # use position of harmonic oscillator as observable
>>> results = mbar.computeExpectations(A_kn, return_dict=True)
where results is a dictionary with keys mu, sigma, and Theta, where mu[i] is the array of the estimate for the average of the observable for in state i, sigma[i] is the estimated standard deviation of the mu estimates, and Theta[i,j] is the covarinace matrix of the log weights.
See the docstring help for these individual methods for more information on exact usage; in Python or IPython, you can view the docstrings with help().
Authors
- Kyle A. Beauchamp [email protected]
- John D. Chodera [email protected]
- Levi N. Naden [email protected]
- Michael R. Shirts [email protected]
