Adaptive Biasing Force Method

Navigation: Documentation / Methods / Adaptive Biasing Force Method

---

Chapters

• ABF:Description
• ABF:Controls
• ABF:Collective variables
• ABF:Post-processing
• ABF:Examples
• ABF:Multiple walkers approach

Description

Free Energy from Unconstrained MD Simulations

The Adaptive Biasing Force (ABF) method calculates the free energy as a function of selected collective variables from unconstrained molecular dynamics (MD) simulations. The method does not provide the free energy directly, but instead, it provides the free energy gradient $\frac{\partial G(\boldsymbol \xi)}{\partial \boldsymbol \xi}$, which must be integrated to get the free energy:

$\Delta G = G(\boldsymbol \xi_2) - G(\boldsymbol \xi_1) = \int_{\boldsymbol \xi_1}^{\boldsymbol \xi_2} \frac{\partial G(\boldsymbol \xi)} {\partial \boldsymbol \xi} \boldsymbol d \xi$ ... (1)

The free energy gradient is calculated as a mean of instantaneous collective force $\bold F_{ICF}$:

$\frac{\partial G(\boldsymbol \xi^{*})}{\partial \boldsymbol \xi} = - {\langle \bold F_{ICF} \rangle}_{\boldsymbol \xi=\boldsymbol \xi^{*}}$ ... (2)

with the instantaneous collective force calculated from the time evolution of the collective variable:

$\bold F_{ICF} = \frac {d} {dt} \left( \bold Z^{-1} \frac{d \boldsymbol \xi} {dt} \right)$ ... (3)

where $\bold Z$ is the matrix in the form:

$[Z]_{ij}=\sum_{k=1}^{N_{atoms}} \frac{1}{m_k} \frac{\partial \xi_i}{\partial \bold x_k} \frac{\partial \xi_j}{\partial \bold x_k}$ ... (4)

The analytical calculation of instantaneous collective force by Equation 3 requires the second derivatives of collective variables with respect to Cartesian coordinates. Since this can be prohibitive for complex collective variables such as the simple base-pair parameters, Equation 3 is evaluated numerically by a finite-difference approach.

Sampling Space Discretization

Due to numerical reasons, mean forces are collected on a regular grid. The averaging of instantaneous collective force is then done in small intervals centered at discrete CV values:

$\frac{\partial G(\boldsymbol \xi^{*})}{\partial \boldsymbol \xi} = - {\langle \bold F_{ICF} \rangle}_{\boldsymbol \xi=\boldsymbol \xi^{*} \pm \bold h/2}$ ... (5)

with the standard error:

$\sigma_e \left( \frac{\partial G(\boldsymbol \xi^{*})}{\partial \boldsymbol \xi} \right) = \sqrt{ \frac {N_{corr}} {N_b} } \sigma(\bold F_{ICF})_ {\boldsymbol \xi=\boldsymbol \xi^{*} \pm \bold h/2}$ ... (6)

where the standard deviation is given by:

$\sigma(\bold F_{ICF})_ {\boldsymbol \xi=\boldsymbol \xi^{*} \pm \bold h/2} = \sqrt { \frac{1}{N_b} \sum_{i=1}^{N_b} \left( F_{ICF}(\boldsymbol \xi_i) - \overline {F_{ICF}} \right)^2 \vert _{\boldsymbol \xi=\boldsymbol \xi^{*} \pm \bold h/2} }$ ... (7)

where $\bold h$ is the interval size also called a bin, $N_b$ is the number of samples collected in a bin centered at $\boldsymbol \xi^{*}$, and $N_{corr}$ is a statistical inefficiency due to correlation in time series.

Therefore, each CV involved in ABF simulations must be discretized by specifying an interval in which the sampling is performed and the number of intervals (bins) for discretization, for further details, see ABF:Collective variables. The increasing number of bins improves the accuracy of Equation 5 and subsequently the quality of the integrated free energy (Equation 1) but it also increases the noise because of a smaller number of samples collected in a bin. A reasonable compromise is the number of bins, which leads to 0.1 Å or 1-2° bin widths.

Adaptive Bias

ABF introduces an adaptive bias, which improves the sampling of rare events. The bias removes barriers or higher free energy regions in the space described by predefined collective variables $\boldsymbol \xi$. As a result, the system evolves by free diffusion alongside these collective variables. The bias is derived from the free energy, which is projected in the form of biasing force $\bold F_{bias}$ to the Cartesian space and removed from force $\bold F_{pot}$ originated from inter-atomic interaction potential $V(\bold R)$. The application of the bias thus leads to the modified equations of motions:

$m_{i} \frac { d^2 \bold r_i }{dt^2} = \bold F_{pot,i}(\bold r_i) - \bold F_{bias,i}(\bold r_i) = - \frac{\partial V(\bold R)} {\partial \bold r_i} + \frac{\partial G(\boldsymbol \xi)}{\partial \boldsymbol \xi} \frac{\partial \boldsymbol \xi}{\partial \bold r_i}$ ... (8)

where $m_i$ is mass of atom i, $\bold r_i$ is atom position, and $t$ is time.

Since the biasing force is not known prior to the simulation, it is calculated during the simulation and adaptively applied. To accelerate sampling, the biasing force is applied even if an inadequate number of samples is collected in a bin. In this case, the biasing force is scaled in the early stages to avoid artifacts from applications of overestimated biasing forces. The biasing force can also be smoothed to decrease noise in collected data. For further details, see feimode in ABF:Controls.

Instantaneous Collective Forces

PMFLib implements two algorithms to calculate Equation 3: simplified and original ABF algorithms. The original ABF algorithm is theoretically more accurate than the simplified one but it works only for MD without constraints. Therefore, the simplified algorithm is the default one as it is numerically more robust at all MD setups.

Simplified ABF Algorithm

The simplified ABF algorithm (fmode=1) uses the result of the product rule for derivatives applied on Equation 3 leading into the two distinct contributions into the instantaneous collective forces:

$\frac {d} {dt} \left( \bold Z^{-1} \frac{d \boldsymbol \xi} {dt} \right) \big|_{(t)} = + \left( \bold Z^{-1}(t) \nabla \boldsymbol \xi (t) \right) \cdot \bold a(t) + \frac{d}{dt} \left( \bold Z^{-1}(t) \nabla \boldsymbol \xi (t)) \right) \cdot \bold v(t) = \bold F_{ICF,pot}(t) + \bold F_{ICF,kin}(t)$ ... (9)

The potential contribution $\bold F_{ICF,pot}(t)$ employs the accelerations $\bold a(t)$ calculated back from the velocities. This ensures that effect of SHAKE and other constraints on acceleration is properly counted (see the leap-frog integration algorithm for further details):

$\bold a(t) = \frac{\bold v(t+dt/2) - \bold v(t-dt/2)}{dt}$ ... (10)

The kinetic contribution $\bold F_{ICF,kin}(t-dt/2)$ is calculated from numerical differentiation:

$\bold F_{ICF,kin}(t-dt/2) = \frac {\bold Z^{-1}(t) \nabla \boldsymbol \xi (t) - \bold Z^{-1}(t-dt) \nabla \boldsymbol \xi (t-dt)}{dt} \bold v(t-dt/2)$ ... (11)

Finally, to get $\bold F_{ICF,kin}$ at the same time as $\bold F_{ICF,kin}$, two values are averaged:

$\bold F_{ICF,kin}(t) = \frac {\bold F_{ICF,kin}(t+dt/2) + \bold F_{ICF,kin}(t-dt/2) }{2}$ ... (12)

The algorithm uses a history of values collected in two consecutive time steps. But, the first result is available from the fourth time step to be compatible with the original algorithm.

Original ABF Algorithm

The original ABF algorithm (fmode=2) implements the algorithm presented by Darve at al.. This algorithm should not be used when SHAKE is enabled because the accelerations entering into the calculations are not corrected for constraints. The algorithm uses a history of values collected in four consecutive time steps. Thus, the first result is available from the fourth time step.

$\frac {d} {dt} \left( \bold Z^{-1} \frac{d \boldsymbol \xi} {dt} \right) \big|_{(t)} = \frac{1}{2} \left( \frac{\bold p^+_{\xi}(t+dt) - \bold p^+_{\xi}(t) } {dt} + \frac{\bold p^-_{\xi}(t+dt) - \bold p^-_{\xi}(t) } {dt} \right)$ ... (11)

where

$\bold p^+_{\xi}(t) = \bold Z^{-1}(t) \nabla \boldsymbol \xi(t) \cdot \left[ \bold v(t+dt/2) - \frac{dt}{6} \bold a(t+dt) \right]$ ... (12) $\bold p^-_{\xi}(t) = \bold Z^{-1}(t) \nabla \boldsymbol \xi(t) \cdot \left[ \bold v(t-dt/2) + \frac{dt}{6} \bold a(t-dt) \right]$ ... (13)