next up previous contents index
Next: Alchemical Free Energy Methods1 Up: Collective Variable-based Calculations1 Previous: Declaring and using collective   Contents   Index


Biasing and analysis methods

All of the biasing and analysis methods implemented (abf, harmonic, histogram and metadynamics) recognize the following options:

Adaptive Biasing Force

For a full description of the Adaptive Biasing Force method, see reference [20]. For details about this implementation, see references [32] and [33]. When publishing research that makes use of this functionality, please cite references [20] and [33].

An alternate usage of this feature is the application of custom tabulated biasing potentials to one or more colvars. See inputPrefix and updateBias below.

ABF is based on the thermodynamic integration (TI) scheme for computing free energy profiles. The free energy as a function of a set of collective variables $ {\mbox{\boldmath {$\xi$}}}=(\xi_{i})_{i\in[1,n]}$ is defined from the canonical distribution of $ {\mbox{\boldmath {$\xi$}}}$, $ {\mathcal P}({\mbox{\boldmath {$\xi$}}})$:

$\displaystyle A({\mbox{\boldmath {$\xi$}}}) = -\frac{1}{\beta} \ln {\mathcal P}({\mbox{\boldmath {$\xi$}}}) + A_0$ (45)

In the TI formalism, the free energy is obtained from its gradient, which is generally calculated in the form of the average of a force $ {\mbox{\boldmath {$F$}}}_\xi$ exerted on $ {\mbox{\boldmath {$\xi$}}}$, taken over an iso- $ {\mbox{\boldmath {$\xi$}}}$ surface:

$\displaystyle {\mbox{\boldmath {$\nabla$}}}_\xi A({\mbox{\boldmath {$\xi$}}}) =...
...t\langle -{\mbox{\boldmath {$F$}}}_\xi \right\rangle_{\mbox{\boldmath {$\xi$}}}$ (46)

Several formulae that take the form of (47) have been proposed. This implementation relies partly on the classic formulation [14], and partly on a more versatile scheme originating in a work by Ruiz-Montero et al. [58], generalized by den Otter [21] and extended to multiple variables by Ciccotti et al. [17]. Consider a system subject to constraints of the form $ \sigma_{k}({\mbox{\boldmath {$x$}}}) = 0$. Let ( $ {\mbox{\boldmath {$v$}}}_{i})_{i\in[1,n]}$ be arbitrarily chosen vector fields ( $ \mathbb{R}^{3N}\rightarrow\mathbb{R}^{3N}$) verifying, for all $ i$, $ j$, and $ k$:

$\displaystyle {\mbox{\boldmath {$v$}}}_{i} \cdot \mbox{\boldmath$\nabla_{\!\!x}\,$}\xi_{j}$ $\displaystyle =$ $\displaystyle \delta_{ij}$ (47)
$\displaystyle {\mbox{\boldmath {$v$}}}_{i} \cdot \mbox{\boldmath$\nabla_{\!\!x}\,$}\sigma_{k}$ $\displaystyle =$ 0 (48)

then the following holds [17]:

$\displaystyle \frac{\partial A}{\partial \xi_{i}} = \left\langle {\mbox{\boldma...
...\,$}\cdot {\mbox{\boldmath {$v$}}}_{i} \right\rangle_{\mbox{\boldmath {$\xi$}}}$ (49)

where $ V$ is the potential energy function. $ {\mbox{\boldmath {$v$}}}_{i}$ can be interpreted as the direction along which the force acting on variable $ \xi_{i}$ is measured, whereas the second term in the average corresponds to the geometric entropy contribution that appears as a Jacobian correction in the classic formalism [14]. Condition (48) states that the direction along which the system force on $ \xi_{i}$ is measured is orthogonal to the gradient of $ \xi_{j}$, which means that the force measured on $ \xi_{i}$ does not act on $ \xi_{j}$.

Equation (49) implies that constraint forces are orthogonal to the directions along which the free energy gradient is measured, so that the measurement is effectively performed on unconstrained degrees of freedom. In NAMD, constraints are typically applied to the lengths of bonds involving hydrogen atoms, for example in TIP3P water molecules (parameter rigidBonds, section 5.6.1).

In the framework of ABF, $ {\bf F}_\xi$ is accumulated in bins of finite size, $ \delta \xi$, thereby providing an estimate of the free energy gradient according to equation (47). The biasing force applied along the colective variables to overcome free energy barriers is calculated as:

$\displaystyle {\bf F}^{\rm ABF} =$   $\displaystyle \mbox{\boldmath$\nabla_{\!\!x}\,$}$$\displaystyle \widetilde A({\mbox{\boldmath {$\xi$}}})$ (50)

where $ \nabla_{\!\!x}\,$$ \widetilde A$ denotes the current estimate of the free energy gradient at the current point $ {\mbox{\boldmath {$\xi$}}}$ in the collective variable subspace.

As sampling of the phase space proceeds, the estimate $ \nabla_{\!\!x}\,$$ \widetilde A$ is progressively refined. The biasing force introduced in the equations of motion guarantees that in the bin centered around $ {\mbox{\boldmath {$\xi$}}}$, the forces acting along the selected collective variables average to zero over time. Eventually, as the undelying free energy surface is canceled by the adaptive bias, evolution of the system along $ {\mbox{\boldmath {$\xi$}}}$ is governed mainly by diffusion. Although this implementation of ABF can in principle be used in arbitrary dimension, a higher-dimension collective variable space is likely to result in sampling difficulties. Most commonly, the number of variables is one or two.

ABF requirements on collective variables

  1. Only linear combinations of colvar components can be used in ABF calculations.
  2. Availability of system forces is necessary. The following colvar components can be used in ABF calculations: distance, distance_xy, distance_z, angle, dihedral, gyration, rmsd and eigenvector. Atom groups may not be replaced by dummy atoms, unless they are excluded from the force measurement by specifying oneSiteSystemForce, if available.
  3. Mutual orthogonality of colvars. In a multidimensional ABF calculation, equation (48) must be satisfied for any two colvars $ \xi_{i}$ and $ \xi_{j}$. Various cases fulfill this orthogonality condition:
  4. Mutual orthogonality of components: when several components are combined into a colvar, it is assumed that their vectors $ {\mbox{\boldmath {$v$}}}_{i}$ (equation (50)) are mutually orthogonal. The cases described for colvars in the previous paragraph apply.
  5. Orthogonality of colvars and constraints: equation 49 can be satisfied in two simple ways, if either no constrained atoms are involved in the force measurement (see point 3 above) or pairs of atoms joined by a constraint bond are part of an atom group which only intervenes through its center (center of mass or geometric center) in the force measurement. In the latter case, the contributions of the two atoms to the left-hand side of equation 49 cancel out. For example, all atoms of a rigid TIP3P water molecule can safely be included in an atom group used in a distance component.

Parameters for ABF

The following parameters can be set in the ABF configuration block (in addition to generic bias parameters such as colvars):

ABF also depends on parameters from collective variables to define the grid on which free energy gradients are computed. In the direction of each colvar, the grid ranges from lowerBoundary to upperBoundary, and the bin width (grid spacing) is set by the width parameter.

Output files

The ABF bias produces the following files, all in multicolumn ASCII format:

If several ABF biases are defined concurrently, their name is inserted to produce unique filenames for output, as in $ <$outputName$ >$.abf1.grad. This should not be done routinely and could lead to meaningless results: only do it if you know what you are doing!

If the colvar space has been partitioned into sections (windows) in which independent ABF simulations have been run, the resulting data can be merged using the inputPrefix option described above (a NAMD run of 0 steps is enough).

Reconstructing a multidimensional free energy surface

If a one-dimensional calculation is performed, the estimated free energy gradient is automatically integrated and a potential of mean force is written under the file name <outputName>.pmf, in a plain text format that can be read by most data plotting and analysis programs (e.g. gnuplot).

In dimension 2 or greater, integrating the discretized gradient becomes non-trivial. The standalone utility abf_integrate is provided to perform that task. abf_integrate reads the gradient data and uses it to perform a Monte-Carlo (M-C) simulation in discretized collective variable space (specifically, on the same grid used by ABF to discretize the free energy gradient). By default, a history-dependent bias (similar in spirit to metadynamics) is used: at each M-C step, the bias at the current position is incremented by a preset amount (the hill height). Upon convergence, this bias counteracts optimally the underlying gradient; it is negated to obtain the estimate of the free energy surface.

abf_integrate is invoked using the command-line:

integrate <gradient_file> [-n <nsteps>] [-t <temp>] [-m (0|1)]
                          [-h <hill_height>] [-f <factor>]

The gradient file name is provided first, followed by other parameters in any order. They are described below, with their default value in square brackets:

Using the default values of all parameters should give reasonable results in most cases.

abf_integrate produces the following output files:

Note: Typically, the ``deviation'' vector field does not vanish as the integration converges. This happens because the numerical estimate of the gradient does not exactly derive from a potential, due to numerical approximations used to obtain it (finite sampling and discretization on a grid).


Many methods have been introduced in the past that make use of an artificial energy term, that changes and adapts over time, to reconstruct a potential of mean force from a conventional molecular dynamics simulation [34,27,71,19,42,35]. One of the most recent, metadynamics, was first designed as a stepwise algorithm, which may be roughly described as an ``adaptive umbrella sampling'' [42], and was later made continuous over time [36]. This implementation provides only he latter version, which is the most commonly used.

In metadynamics, the external potential on the colvars $ {\mbox{\boldmath {$\xi$}}} =
(\xi_{1}, \xi_{2}, \ldots, \xi_{N_{\mathrm{cv}}})$ is:

$\displaystyle V_{\mathrm{meta}}({\mbox{\boldmath {$\xi$}}}) \; = \; { \sum_{t' ...
...t(-\frac{(\xi_{i}-\xi_{i}(t'))^{2}}{2\delta_{\xi_{i}}^{2}}\right) } }\mathrm{,}$ (51)

that is, $ V_{\mathrm{meta}}$ is a history-dependent potential, which acts on the current values of the colvars $ {\mbox{\boldmath {$\xi$}}}$ and depends parametrically on the previous values of the colvars. It is constructed as a sum of $ N_{\mathrm{cv}}$-dimensional repulsive Gaussian ``hills'' with a height $ W$: their centers are located at the previously explored configurations $ \left({\mbox{\boldmath {$\xi$}}}(\delta{}t),
{\mbox{\boldmath {$\xi$}}}(2\delta{}t), \ldots\right)$, and they extend by approximately $ 2\delta_{\xi_{i}}$ in the direction of the $ i$-th colvar.

As the system evolves according to the underlying potential of mean force $ A({\mbox{\boldmath {$\xi$}}})$ incremented by the metadynamics potential $ V_{\mathrm{meta}}({\mbox{\boldmath {$\xi$}}})$, new hills will tend to accumulate in the regions with a lower effective free energy $ \tilde{A}({\mbox{\boldmath {$\xi$}}}) =
A({\mbox{\boldmath {$\xi$}}})+V_{\mathrm{meta}}({\mbox{\boldmath {$\xi$}}})$. That is, the probability of having a given system configuration $ {\mbox{\boldmath {$\xi^{*}$}}}$ being explored (and thus, a hill being added there) is proportional to $ \exp\left(-\tilde{A}({\mbox{\boldmath {$\xi^{*}$}}})/\kappa_{\mathrm{B}}T\right)$, which tends to a nearly flat histogram when the simulation is continued until the system has deposited hills across the whole free energy landscape. In this situation, $ -V_{\mathrm{meta}}({\mbox{\boldmath {$\xi$}}})$ is a good approximant of the free energy $ A({\mbox{\boldmath {$\xi$}}})$, and the only dependence on the specific conformational history $ {\mbox{\boldmath {$\xi$}}}(\delta{}t), {\mbox{\boldmath {$\xi$}}}(2\delta{}t), \ldots$ is by an irrelevant additive constant:

$\displaystyle A({\mbox{\boldmath {$\xi$}}}) \; \simeq \; { -V_{\mathrm{meta}}({\mbox{\boldmath {$\xi$}}}) + K }$ (52)

Provided that the set of collective variables fully describes the relevant degrees of freedom, the accuracy of the reconstructed profile is a function of the ratio between $ W$ and $ \delta{}t$ [13]. For the optimal choice of $ \delta_{\xi_{i}}$ and $ D_{\xi_{i}}$, the diffusion constant of the variable $ \xi_{i}$, see reference [13]. As a rule of thumb, the very upper limit for the ratio $ W/\delta{}t$ is given by $ \kappa_{\mathrm{B}}T/\tau_{{\mbox{\boldmath {$\xi$}}}}$, where $ \tau_{{\mbox{\boldmath {$\xi$}}}}$ is the longest among $ {\mbox{\boldmath {$\xi$}}}$'s correlation times. In the most typical conditions, to achieve a good statistical convergence the user would prefer to keep $ W/\delta{}t$ much smaller than $ \kappa_{\mathrm{B}}T/\tau_{{\mbox{\boldmath {$\xi$}}}}$.

Given $ \Delta\xi$ the extension of the free energy profile along the colvar $ \xi$, and $ A^{*}=A({\mbox{\boldmath {$\xi$}}}^{*})$ the highest free energy that needs to be sampled (e.g. that of a transition state), the upper bound for the required simulation time is of the order of $ N_{\mathrm{s}}(\xi) = (A^{*}\Delta\xi)/(W2\delta_{\xi})$ multiples of $ \delta{}t$. When several colvars $ {\mbox{\boldmath {$\xi$}}}$ are used, the upper bound amounts to $ N_{\mathrm{s}}(\xi_{1}) \times
N_{\mathrm{s}}(\xi_{2}) \times \ldots \times
N_{\mathrm{s}}(\xi_{N_{\mathrm{cv}}}) \times \delta{}t$.

In metadynamics runs performed with this module, the parameter $ \delta_{\xi_{i}}$ for each hill (eq. 52) is chosen as approximately half the width of the corresponding colvar $ \xi_{i}$, while all the other parameters must be provided within the metadynamics {...} block.

The only mandatory parameter is the colvars option, listing all the variables to which this bias is applied. Note: multidimensional PMFs are obtained with one metadynamics instance applied to all the colvars, and not with multiple instances applied to individual colvars.

The following two options have default values that are reasonable in typical situations, but it is strongly recommended that the user chooses them according to the above discussion on the diffusion times of the variables, $ \tau_{{\mbox{\boldmath {$\xi$}}}}$:

The following options define the configuration for the ``well-tempered'' metadynamics approach [4]:

The following options control the performance of metadynamics calculations, but do not affect the results:

The following options define metadynamics calculations with more than one replica:

The following options may be useful for applications that go beyond the direct application of metadynamics for a calculation of a PMF.

Harmonic restraints and Steered Molecular Dynamics

The harmonic biasing method may be used to enforce fixed or moving restraints, including variants of Steered and Targeted MD. Within energy minimization runs, it allows for restrained minimization, e.g. to calculate relaxed potential energy surfaces. In the context of the colvars module, harmonic potentials are meant according to their textbook definition: $ \displaystyle V({\mbox{\boldmath {$x$}}}) = \frac{1}{2} k ({\mbox{\boldmath {$x$}}} - {\mbox{\boldmath {$x_0$}}})^2$. Note that this differs from harmonic bond and angle potentials in common force fields, where the factor of one half is typically omitted, resulting in a non-standard definition of the force constant. The restraint energy is reported by NAMD under the MISC title. A harmonic restraint is set up by a harmonic {...} block, which may contain (in addition to the standard option colvars) the following keywords:

Tip: A complex set of restraints can be applied to a system, by defining several colvars, and applying one or more harmonic restraints to different groups of colvars. In some cases, dozens of colvars can be defined, but their value may not be relevant: to limit the size of the colvars trajectory file, it may be wise to disable outputValue for such ``ancillary'' variables, and leave it enabled only for ``relevant'' ones.

Multidimensional histograms

The histogram feature is used to record the distribution of a set of collective variables in the form of a N-dimensional histogram. It functions as a ``collective variable bias'', and is invoked by adding a histogram block to the colvars configuration file.

In addition to the common parameters name and colvars described above, a histogram block may define the following parameter:

Like the ABF and metadynamics biases, histogram uses parameters from the colvars to define its grid. The grid ranges from lowerBoundary to upperBoundary, and the bin width is set by the width parameter.

next up previous contents index
Next: Alchemical Free Energy Methods1 Up: Collective Variable-based Calculations1 Previous: Declaring and using collective   Contents   Index