Usage#
Note: The following description of the usage of FastMBAR is provided here to help better undertand when and how to use FastMBAR in general. It might be too general and thus too abstract for firt-time users. If it is the case, we recommand that you read the Examples first , which can help you understand the following general descriptions.
1. When to use MBAR equations?#
A general question in computational chemistry is to calculate relative free energies of different thermodynamic states. Let us assume there are \(M\) thermodynamic states and the \(k\) th state has a reduced potential energy function of \(U_k(x)\), where \(x\) represents a system conformation. The probability distribution of conformations \(x\) in the \(k\) th state is the Boltzmann distribution \(p_k(x) \propto \exp\left(-U_k(x)\right)\). (Note that \(U_k(x)\) is the reduced potential energy function and has included the inverse temperature \(\beta = 1/k_bT\), i.e., if the original potential energy function is \(E_k(x)\), \(U_k(x)\) is defined as \(\beta E_k(x)\). Therefore, \(U_k(x)\) is a unitless quantity.) The reduced free energy of the \(k\) th state is \(F_k = -\ln \int{\exp\left(-U_k(x)\right) dx}\). (Note \(F_k\) is the reduced free energy and is also unitless. To convert it into a free energy with an energy unit, you need to mutiply it by \(k_bT\), i.e., \(k_bT F_k\).) There are many ways to compute the relative free energies of the \(M\) states and one of them is as follows:
Sample conformations from each of the \(M\) states separately using either molecular dynamcis or Monte Carlo sampling.
Reweight the above sampled conformations from all the \(M\) states to compute the relative free energies.
Let us assume that \(N_k\) conformations are sampled from the \(k\) th state in the first step and they are represented as \(\{x^{n}_k, n = 1, ..., N_k\}\) for \(k = 1, ..., M\). MBAR equations are used in the second step to reweight conformations and compute relative free energies \(F_k\) for \(k = 1, ..., M\). After MBAR equations are solved in the second step, they can also be used to calculate relative free energies of perturbed states from which no conformations are sampled. Let us assume there are \(L\) perturbed states and the \(l\) th perturbed state has a reduced potential energy function of \(U^{\prime}_l(x)\).
The situation described above is quite general and it includes alchemical free energy calculation and calculating potential of mean force (PMF) with umbrella sampling or temperature/Hamiltonian replica exchange simulations.
2. How to use FastMBAR?#
The following two inputs are required to initialize FastMBAR: a matrix \(A\) of reduced potential energies in sampled states and a vector \(v\) of numbers of sampled conformations from each state, both of which can be calculated based on sampling.
If the sampling is conducted as described above, the matrix of reduced potential energies in sampled states should be calculated as the matrix \(A_{M,N}\) as shown in the blue bracket in the following figure:
The matrix \(A\) has \(M\) rows and \(N\) columns, where \(M\) is the number of states from which conformations are sampled and \(N\) is the total number of conformations sampled from all states, i.e., \(N = \sum_{k=1}^{M} N_k\). The reduced potential energy of each conformation, no matter which state it is sampled from, is evaluated in all \(M\) states. These \(M\) reduced potential energies constitue one column of the matrix \(A\). The vector \(v\) is constructed as \(v = (N_1, N_2, ..., N_M)\).
In FastMBAR, both the matrix \(A\) and the vector \(v\) have to be numpy arrays: \(A\) is a two dimensional numpy array and \(v\) is a one dimensional numpy array. With the matrix \(A\) and the vector \(v\), we can create a FastMBAR object using the following Python command:
# import both Numpy and FastMBAR package
import numpy as np
from FastMBAR import *
# construct the energy matrix A and the vector v
...
# construct and initialize a FastMBAR object with the energy matrix A and the vector v
# set cuda = True in the following command if you want to run the calcuation on GPUs
fastmbar = FastMBAR(energy = A, num_conf = v, cuda=False, verbose = True)
# after initialization, the relative free energies of the M states is stored in fastmbar.F
print(fastmbar.F)
FastMBAR can also compute the uncertainty of the estimated relative free energies using the following command.
# import both Numpy and FastMBAR package
import numpy as np
from FastMBAR import *
# construct the energy matrix A and the vector v
...
# construct and initialize a FastMBAR object with the energy matrix A and the vector v
# set cuda = True in the following command if you want to run the calcuation on GPUs
fastmbar = FastMBAR(energy = A, num_conf = v, cuda=False, bootstrap = True, verbose = True)
# after initialization, the relative free energies of the M states and standard deviations
# are stored in fastmbar.F and fastmbar.F_std, respectively.
print(fastmbar.F)
print(fastmbar.F_std)
As stated above, the FastMBAR object constructed using the matrix \(A\) and the vector \(v\) can also be used to compute the relative free energies for perturbed states from which no conformations are sampled. In order to do that, a seperate reduced potential energy matrix \(B\) is required. The matrix \(B\) has \(L\) rows and \(N\) columns, where \(L\) is the number of perturbed states and \(N\) is the total number of conformations sampled from the previous \(M\) states, i.e., \(N = \sum_{k=1}^{M} N_k\). As shown in the red bracket in above figure, the energy term \(B_{l,n}\) in the matrix \(B\) is calculated as the reduced potential energy of the \(n\) th conformation in in the \(l\) th perturbed states. With the matrix \(B\) and the FastMBAR object constructed using the matrix \(A\) and the vector \(v\), the relative free energies of the perturbed states can be calculated using the following command:
## calculate relative free enegies of perturbed states
F_perturbed, F_perturbed_std = fastmbar.calculate_free_energies_of_perturbed_states(B)
## relative free enegies of perturbed states
print(F_perturbed)
## standard deviation of relative free enegies of perturbed states
print(F_perturbed_std)