Multistate Bennett Acceptance Ratio

Free Energy Fundamentals
Theory Definitions Exponential Averaging / Zwanzig Thermodynamic Integration Bennett Acceptance Ratio WHAM MBAR
Methods of Free Energy Simulations Thermodynamic Cycle Intermediate States (Constructing) Soft Core Potentials Running Simulations Extracting Simulation Information Analysis / Bootstrap Sampling
Improving Sampling and Accelerating Simulations Umbrella Sampling Hamiltonian Replica Exchange [math]\displaystyle{ \lambda }[/math]-dynamics / Expanded Ensemble

Free Energy How-to's
Example Calculations Absolute Solvation Free Energy Relative Binding Affinity Absolute Binding Affinity
Tutorials Ethanol Solvation - GROMACS Expanded Ensemble - GROMACS Absolute Binding Free Energy - GROMACS

The Multistate Bennett Acceptance Ratio (MBAR) is a direct extension to BAR as it allows for assessing data from all states, and predicting the free energy at an unsampled state. MBAR reduces to BAR in the limit that only two states are sampled. This equation of free energy calculations can also be seen as a zero-width bin WHAM.

Much like WHAM, the free energies provided by this method are only a statistical estimator, however, MBAR has been shown to have the lowest variance estimator to date.

Derivation

MBAR is derived from a set of [math]\displaystyle{ K \times K }[/math] weighting functions, [math]\displaystyle{ \alpha_{i,j}(\vec{q}) }[/math], that minimized the variance during the reweighting across the board. Starting from our core free energy equation, we have

[math]\displaystyle{ \displaystyle\Delta A_{ij} = -\beta^{-1} \ln\frac{Q_j}{Q_i} }[/math]

We can also manipulate the same identity Bennett used when starting his derivation and then a one extra bit to come up with the relation

[math]\displaystyle{ Q_i\left\langle\alpha_{ij}\exp\left(-\beta U_j\right)\right\rangle_i = Q_j\left\langle\alpha_{ij}\exp\left(-\beta U_i\right)\right\rangle_j }[/math]

Although this may not seem like much, it does allow us to write out the following:

[math]\displaystyle{ \displaystyle \sum\limits_{j=1}^K \frac{\hat{Q_i}}{N_i} \sum\limits_{n=1}^{N_i} \alpha_{ij}\exp\left(-\beta U_j(\vec{q}_{i,n})\right) = \sum\limits_{j=1}^K \frac{\hat{Q_j}}{N_j} \sum\limits_{n=1}^{N_j} \alpha_{ij}\exp\left(-\beta U_i(\vec{q}_{j,n})\right) }[/math]

assuming we use the empirical estimator for the expectation values of [math]\displaystyle{ \left\langle g \right\rangle_i = N_i^{-1}\sum_{n=1}^{N_i}g(\vec{q}_{i,n}) }[/math]

Choosing the optimal [math]\displaystyle{ \alpha_{ij} }[/math] can be done by looking through the literature at extended bridge sampling. We then get an [math]\displaystyle{ \alpha_{ij} }[/math] of:

[math]\displaystyle{ \displaystyle \alpha_{ij} = \frac{N_j \hat{c_j}^{-1}}{\sum\limits_{k=1}^K N_k \hat{c_k}^{-1} \exp\left(-\beta U_k \right)} }[/math].

After making all the necessary substitutions, we can get an expression for an estimated free energy of:

[math]\displaystyle{ \displaystyle \hat{A_i} = -\beta^{-1} \ln \sum\limits_{j=1}^K \sum\limits_{n=1}^{N_j} \frac{\exp\left[-\beta U_i\right]}{\sum_{k=1}^K N_k \exp\left[\beta\hat{A_k} - \beta U_k\right]} }[/math].

One of the first things you should notice is that we have a single free energy, not a difference. This is not a typo, but the free energies for a given set of states is only uniquely determined up to an additive constant. Because of this, one free energy must be taken in reference and thus we are once again calculating free energy differences.

Reduced potential

It is important to note that the [math]\displaystyle{ U }[/math] that appears in the MBAR derivation and equation is not only valid for potential energies, but any generalized/reduced potential as a function of pressure, volume, chemical potential, and number of particles. For example, in a general form, we can take some subset of the additive terms in the following to define the reduced potential [math]\displaystyle{ u_i(x) }[/math] for thermodynamic state [math]\displaystyle{ i }[/math]:

[math]\displaystyle{ \displaystyle u_i(x) \equiv \beta_i [ U_i(x) + p_i V(x) + {\bf \mu}_i^\mathrm{T} {\bf N}(x) ] }[/math]

The reduced potential function is uniquely defined by some combination of thermodynamic parameters [math]\displaystyle{ \beta_i }[/math] denoting the inverse temperature, [math]\displaystyle{ U_i(x) }[/math] denoting the potential energy function, [math]\displaystyle{ p_i }[/math] denoting the pressure, and [math]\displaystyle{ {\bf \mu}_i }[/math] denoting the chemical potential of one or more components of the system. These latter two thermodynamic variables are conjugate to the box volume [math]\displaystyle{ V(x) }[/math] and particle numbers [math]\displaystyle{ {\bf N}(x) }[/math]. Use of the reduced potential simplifies the MBAR equations and generalizes them to the computation of arbitrary reduced free energy differences [math]\displaystyle{ \Delta f_{ij} \equiv \beta_j A_j - \beta_i A_i }[/math] among states.

[math]\displaystyle{ \displaystyle \hat{f_i} = -\ln \sum\limits_{j=1}^K \sum\limits_{n=1}^{N_j} \frac{\exp\left[-u_i(x_{nj})\right]}{\sum_{k=1}^K N_k \exp\left[\hat{f_k} - u_k(x_{nj})\right]} }[/math].

In the sum above [math]\displaystyle{ x_{nj} }[/math] indicates the [math]\displaystyle{ n }[/math]th sample from the [math]\displaystyle{ j }[/math]th state. Interestingly, in this formula, it actually doesn't matter which sample comes from which state, so it can be rewritten as a sum over all [math]\displaystyle{ N = \sum_{j=1}^K N_j }[/math], so that we have: [math]\displaystyle{ \displaystyle \hat{f_i} = -\ln \sum\limits_{n=1}^{N} \frac{\exp\left[-u_i(x_{n})\right]}{\sum_{k=1}^K N_k \exp\left[\hat{f_k} - u_k(x_{n})\right]} }[/math]

Finally, one can rewrite this as: [math]\displaystyle{ \displaystyle \hat{f_i} = -\ln \left\langle \frac{\exp\left[-u_i(x_{n})\right]}{\sum_{k=1}^K \frac{N_k}{N} \exp\left[\hat{f_k} - u_k(x_{n})\right]}\right\rangle }[/math]

Where the ensemble averaged over is the mixture ensemble that consist of samples taken from each state [math]\displaystyle{ i }[/math] [math]\displaystyle{ \frac{N_i}{N} }[/math]th of the time.

Remember always the free energy will change depending on which reduced potential you use, so please take this into account when working with MBAR.

Estimating Free Energies with MBAR

MBAR provides the direct equation to find free energies as show above assuming you have the energies. This can be rather difficult to implement for beginners since this does require iterative solutions like BAR (note: see the free, Python implementation below). Despite this, MBAR still has the lowest variance of all the other methods listed under the theory section of the fundamentals. Further, MBAR has a direct way to calculate errors (see paper for derivation).

You will need a complete set of [math]\displaystyle{ \Delta U_{k,j} }[/math] for all [math]\displaystyle{ K }[/math] states for MBAR to work, just like you do in WHAM. Fortunately, you can do most of this in post processing if need be. Once again, you do not need to calculate [math]\displaystyle{ \Delta U_{k,k} }[/math], but it is still a good self check to ensure the reweighting method is working correctly; all the remaining checks for WHAM hold as well and should be followed when analyzing data.

Download MBAR

MBAR may seem like a daunting set of equations to program yourself, so the authors have provided a Python implementation of MBAR for anyone to use, free of charge at http://github.com/choderalab/pymbar, with a number of examples for a variety of situations at http://github.com/choderalab/pymbar-examplesThe software comes with examples and uses cases. Also bundled with it are the tools to compute expectation values and an implementation of BAR.