pymbar: Multistate Bennett Acceptance Ratio

This dataset contains parallel tempering simulation data for alanine dipeptide in explicit solvent, based on the protocol described in Chodera et al. (2006).

Dataset Specifications

• Simulation Scale: 500 production iterations, each with 20 ps of Hamiltonian dynamics between exchange attempts.
• Total Time: 10 ns per replica.
• Temperature Range: 40 temperatures spanning 273–400 K.
• Sampling Frequency: Data is collected every 1 ps per temperature.

Data Manifest

The following files and structures are available in the dataset:

The pymbar-examples repository provides various implementations of MBAR for computational and experimental chemistry applications. Key use cases include:

• Alchemical Free Energy: Calculating free energy differences from simulations with alchemical intermediates. (Note: This specific example has moved to https://github.com/MobleyLab/alchemical-analysis)
• Constant-Force Optical Trap: Calculating the Potential of Mean Force (PMF) of dsDNA extension from constant-force optical trap experiments.
• Harmonic Oscillators: Demonstrates full pymbar functionality, including validation of error estimate distributions.
• Heat Capacity: Calculating heat capacities using the fluctuation theorem, including $C_v$ at interpolated temperatures.
• Parallel Tempering: Computing 2D PMF from parallel tempering data.
• Umbrella Sampling: Computing 1D PMF from umbrella sampling data.

The heat-capacity.py script implements three distinct methods for computing heat capacity ($C_V$):

1. Fluctuations Method: Uses the formula $\frac{\langle E^2 \rangle - \langle E \rangle^2}{R T^2}$. Uncertainty is estimated via error propagation, which requires an estimate of the number of effective samples.
2. Numerical Differentiation (First Derivative): Computes $C_V$ using finite differences of the average energy. This can be done using either temperature ($C_V = \frac{d\langle E \rangle}{dT}$) or beta ($C_V = -k_B \beta^2 \frac{d\langle E \rangle}{d\beta}$) differences. The method is selected via the -t flag.
3. Second Derivative of Free Energy: Computes $C_V$ based on the second derivative of the reduced free energy $f$ (where $f = -k_B \ln Q$).
   • Temperature differences: $C_V = 2T \frac{df}{dT} + T^2 \frac{d^2f}{dT^2}$
   • Beta differences: $C_V = -k_B \beta^2 \frac{d^2f}{d^2\beta}$
   • Note: Due to central difference computation complications at endpoints, these points are omitted. Error propagation for this method is currently not correctly implemented.

MBAR solves for the free energies of $K$ states. Because solutions are unique only up to an overall constant, pymbar sets the first free energy to zero, leaving $K-1$ free energies.

Built-in Solvers

• SciPy Solvers: Uses scipy.optimize.minimize() and scipy.optimize.roots(). Supported methods include L-BFGS-B, dogleg, CG, BFGS, Newton-CG, TNC, trust-ncg, trust-krylov, trust-exact, and SLSQP for minimization, and hybr and lm for root finding. Methods requiring a Hessian receive the analytical Hessian.
• Adaptive Solver: A built-in solver that alternates between a self-consistent iterative formula and a Newton-Raphson step. It selects the step that minimizes the gradient magnitude. It uses the self-consistent iteration for the first min_sc_iter iterations (default is 2) to improve stability when overlap is poor.

Solver Protocols

To handle different convergence properties of various algorithms, pymbar uses solver protocols. A protocol is a sequence of optimization operations applied one after another.

• Predefined Protocols: Use solver_protocol='default' or solver_protocol='robust'.
• Custom Protocols: Defined as a tuple of dictionaries. Each dictionary specifies a method and an optional continuation boolean (whether to start from the previous method's result or restart from initial values).

# Example of a custom solver protocol
solver_protocol = (
    dict(method="hybr", continuation=True),
    dict(method="adaptive", options=dict(min_sc_iter=0)),
)

You can construct a custom solver_protocol as a tuple of dictionaries to chain different optimization methods.

• method: The name of the solver (e.g., 'hybr', 'L-BFGS-B', 'adaptive').
• continuation: A boolean. If True, the next method starts from the free energies resulting from the previous method. If False, it restarts from the initialized free energies.
• options: A dictionary of arguments passed directly to the underlying solver (e.g., maxiter, min_sc_iter).
• tol: The tolerance level. Note that tol is a direct option to the solver and is passed directly in the protocol, not inside the options dictionary.

If a method in the sequence fails, pymbar will attempt to continue with the next method using the last known free energies (if continuation=True) and will issue a warning.

solver_protocol = (
        dict(method="hybr"),
        dict(method="L-BFGS-B", tol=1.0e-5, continuation=True, options=dict(maxiter=1000)),
        dict(method="adaptive", tol=1.0e-12, options=dict(maxiter=1000, min_sc_iter=5))
)