Umbrella sampling is a sampling method used in computational physics and chemistry. This sampling can sample the rare states which normal molecular dynamic sampling ignored. Therefore, umbrella sampling can improve free energy calculation when a system is undergoing a systematic change.

Biased molecular dynamics simulations

Normal MD simulations samples system in equilibrium. In an MD simulation of the time series of the C-C-C-C dihedral angle of n-butane(aq), only gauche states and trans states are sampled. Because this simulation is only performed in 2 ns, states with high free energy (e.g. cis state) are less likely to happen. These configurations are ignored and it is impossible to calculate the free energy of these states from this simulation. An artificial bias potential is needed in this case to help the molecule cross the energy barrier. With bias potential, rare states can be effectively sampled.

The time series of the C-C-C-C dihedral angle of n-butane(aq) from a 2 ns molecular dynamics simulation.

In this case, a harmonic bias potential ${\displaystyle w(\phi )=-V(\phi )}$ is needed to counteract the dihedral barrier.

${\displaystyle w(\phi )=-k_{\phi }(1+\cos(n\phi -\delta ))}$

The time series of the C-C-C-C dihedral angle of n-butane(aq) from a 2 ns molecular dynamics simulation with a bias potential to stabilize the gauche conformations.

High free energy states were captured by biased simulation. In order to calculate the free energy profile of these states, biased probability distribution has to be converted to an unbiased probability distribution.

Acquire free energy profile from biased simulations

The potential energy ${\displaystyle V'(r)}$ includes the bias potential ${\displaystyle w(s)}$at the reaction coordinate ${\displaystyle s}$ is ${\displaystyle V'(r)=V(r)+w(s)}$

The probability distribution of this potential is

${\displaystyle P'(s)={\frac {\int {\exp \left({\frac {-V'(r)}{k_{B}T}}\right)\delta (f_{\alpha }(r_{1},...,r_{N})-s)d^{N}r}}{\int {\exp \left({\frac {-V'(r)}{k_{B}T}}\right)}d^{N}r}}}$

The probability distribution of unbiased potential is

${\displaystyle P(s)={\frac {\int {\exp \left({\frac {-V(r)}{k_{B}T}}\right)\delta (f_{\alpha }(r_{1},...,r_{N})-s)d^{N}r}}{\int {\exp \left({\frac {-V(r)}{k_{B}T}}\right)}d^{N}r}}}$

From this equation, we can derive,

${\displaystyle P'(s)=\exp \left({\frac {w(s)}{k_{B}T}}\right)P(s)\left\langle \exp \left({\frac {-w(s)}{k_{B}T}}\right)\right\rangle ^{-1}}$

Free energy profile can be calculated from probability distribution by,

${\displaystyle A(s)=k_{B}T\ln P(s)}$

Using this relation, the PMF of the biased simulation can be converted to unbiased PMF by:

${\displaystyle A(s)=A'(s)-w(s)-k_{B}T\ln \left\langle \exp \left({\frac {-w(s)}{k_{B}T}}\right)\right\rangle }$

${\displaystyle -k_{B}T\ln \left\langle \exp \left({\frac {-w(s)}{k_{B}T}}\right)\right\rangle }$ term is denoted as ${\displaystyle F_{i}}$. It is generally a constant and in some cases does not affect the relative energy and no needed to calculate. It can be calculated by ^[1]:

${\displaystyle \exp \left({\frac {-F_{i}}{k_{B}T}}\right)=\int {P(s)\exp \left({\frac {-w(s)}{k_{B}T}}\right)}ds=\int {\exp \left({\frac {-w(s)-A(s)}{k_{B}T}}\right)}ds}$

This method provides free energy profile of all possible states. In umbrella sampling of n-butane(aq), the chosen bias potential covered all reaction coordinates. General cases are more complex, which leads to a more complex determination of bias potential.

Choice of Bias Potential

The previous section discussed the biased molecular dynamic simulation of n-butane(aq). The reaction coordinate is one-dimensional and periodic, and the bias potential was chosen to be the negative the dihedral potential of n-butane^[2]. The optimum bias potential is the opposite of the free energy ${\displaystyle A(s)}$ ^[1]. However, ${\displaystyle A(s)}$ is unknown for most cases. For general cases, the bias potential needs to be adjusted along the reaction coordinate. Thus, a harmonic bias potential restrained on a reference point ${\displaystyle s_{0,i}}$ with repect to a window ${\displaystyle i}$ on the reaction coordinate is introduced^[2]:

${\displaystyle w_{i}(s)={\frac {1}{2}}k_{i}(s-s_{0,i})^{2}}$

Therefore, a full umbrella sampling can be obtained by a series of biased MD simulation on different reference points on the reaction coordinate.

Calculation of the Potential of Mean Force from Umbrella Sampling Data

The Weighted Histogram Analysis Method (WHAM)^[3] transferred a series of biased distribution functions to a single distribution function. The value of ${\displaystyle F_{i}}$ needs to be estimated to give the correct value of ${\displaystyle A(s)}$: ${\displaystyle A(s)=A'(s)-w(s)+F_{i}}$

The true distribution P(s) is the weight average of each step^[1]:

${\displaystyle P(s)=\sum _{i=1}^{window}p_{i}(s)P_{i}(s)}$

And ${\displaystyle p_{i}=N_{i}\exp \left({\frac {-w(s)+F_{i}}{k_{B}T}}\right)}$, where ${\displaystyle N_{i}}$ is the total number of steps sampled for window ${\displaystyle i}$^[3].

Combined with ${\displaystyle \exp \left({\frac {-F_{i}}{k_{B}T}}\right)=\int {P(s)\exp \left({\frac {-w(s)}{k_{B}T}}\right)}ds=\int {\exp \left({\frac {-w(s)-A(s)}{k_{B}T}}\right)}ds}$, both ${\displaystyle P(s)}$ and ${\displaystyle F_{i}}$ can be obtained.

The other way to analyze umbrella sampling is Umbrella Integration, see^[1].

See Also

Wikipedia:Umbrella sampling

For more information about umbrella sampling, see^[4]

References