next up previous contents index
Next: Harmonic restraints Up: Biasing and analysis methods Previous: Extended-system Adaptive Biasing Force   Contents   Index



The metadynamics method uses a history-dependent potential [67] that generalizes to any type of colvars the conformational flooding [68] and local elevation [69] methods, originally formulated to use as colvars the principal components of a covariance matrix or a set of dihedral angles, respectively. The metadynamics potential on the colvars $ {\mbox{\boldmath {$\xi$}}} = (\xi_{1}, \xi_{2}, \ldots, \xi_{N_{\mathrm{cv}}})$ is defined as:

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

where $ V_{\mathrm{meta}}$ is the history-dependent potential acting on the current values of the colvars $ {\mbox{\boldmath {$\xi$}}}$ , and depends only parametrically on the previous values of the colvars. $ V_{\mathrm{meta}}$ is constructed as a sum of $ N_{\mathrm{cv}}$ -dimensional repulsive Gaussian ``hills'', whose height is a chosen energy constant $ W$ , and whose centers are the previously explored configurations $ \left({\mbox{\boldmath {$\xi$}}}(\delta{}t), {\mbox{\boldmath {$\xi$}}}(2\delta{}t), \ldots\right)$ .

During the simulation, the system evolves towards the nearest minimum of the ``effective'' potential of mean force $ \tilde{A}({\mbox{\boldmath {$\xi$}}})$ , which is the sum of the ``real'' underlying potential of mean force $ A({\mbox{\boldmath {$\xi$}}})$ and the the metadynamics potential, $ V_{\mathrm{meta}}({\mbox{\boldmath {$\xi$}}})$ . Therefore, at any given time the probability of observing the configuration $ {\mbox{\boldmath {$\xi^{*}$}}}$ is proportional to $ \exp\left(-\tilde{A}({\mbox{\boldmath {$\xi^{*}$}}})/\kappa_{\mathrm{B}}T\right)$ : this is also the probability that a new Gaussian ``hill'' is added at that configuration. If the simulation is run for a sufficiently long time, each local minimum is canceled out by the sum of the Gaussian ``hills''. At that stage the ``effective'' potential of mean force $ \tilde{A}({\mbox{\boldmath {$\xi$}}})$ is constant, and $ -V_{\mathrm{meta}}({\mbox{\boldmath {$\xi$}}})$ is an estimator of the ``real'' potential of mean force $ A({\mbox{\boldmath {$\xi$}}})$ , save for an additive constant:

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

Such estimate of the free energy can be provided by enabling writeFreeEnergyFile. Assuming that the set of collective variables includes all relevant degrees of freedom, the predicted error of the estimate is a simple function of the correlation times of the colvars $ \tau_{\xi_{i}}$ , and of the user-defined parameters $ W$ , $ \sigma_{\xi_{i}}$ and $ \delta{}t$ [70]. In typical applications, a good rule of thumb can be to choose the ratio $ W/\delta{}t$ much smaller than $ \kappa_{\mathrm{B}}T/\tau_{{\mbox{\boldmath {$\xi$}}}}$ , where $ \tau_{{\mbox{\boldmath {$\xi$}}}}$ is the longest among $ {\mbox{\boldmath {$\xi$}}}$ 's correlation times: $ \sigma_{\xi_{i}}$ then dictates the resolution of the calculated PMF.

If the metadynamics parameters are chosen correctly, after an equilibration time, $ t_{e}$ , the estimator provided by eq. 13.31 oscillates on time around the ``real'' free energy, thereby a better estimate of the latter can be obtained as the time average of the bias potential after $ t_{e}$ [71,72]:

$\displaystyle A({\mbox{\boldmath {$\xi$}}}) \; = \; {-\frac{1}{t_{tot}-t_{e}} \int_{t_{e}}^{t_{tot}} { V_{\mathrm{meta}}({\mbox{\boldmath {$\xi$}}},t)dt} }$ (13.31)

where $ t_{e}$ is the time after which the bias potential grows (approximately) evenly during the simulation and $ t_{tot}$ is the total simulation time. The free energy calculated according to eq. 13.32 can thus be obtained averaging on time mutiple time-dependent free energy estimates, that can be printed out through the keyword keepFreeEnergyFiles. An alternative is to obtain the free energy profiles by summing the hills added during the simulation; the hills trajectory can be printed out by enabling the option writeHillsTrajectory.

Treatment of the PMF boundaries

In typical scenarios the Gaussian hills of a metadynamics potential are interpolated and summed together onto a grid, which is much more efficient than computing each hill independently at every step (the keyword useGrids is on by default). This numerical approximation typically yields neglibile errors in the resulting PMF [48]. However, due to the finite thickness of the Gaussian function, the metadynamics potential would suddenly vanish each time a variable exceeds its grid boundaries.

To avoid such discontinuity the Colvars metadynamics code will keep an explicit copy of each hill that straddles a grid's boundary, and will use it to compute metadynamics forces outside the grid. This measure is taken to protect the accuracy and stability of a metadynamics simulation, except in cases of ``natural'' boundaries (for example, the $ [0:180]$ interval of an angle colvar) or when the flags hardLowerBoundary and hardUpperBoundary are explicitly set by the user. Unfortunately, processing explicit hills alongside the potential and force grids could easily become inefficient, slowing down the simulation and increasing the state file's size.

In general, it is a good idea to define a repulsive potential to avoid hills from coming too close to the grid's boundaries, for example as a harmonicWalls restraint (see [*]).

Example: Using harmonic walls to protect the grid's boundaries.
colvar {
  name r
  distance { ... }
  upperBoundary 15.0
  width 0.2

metadynamics {
  name meta_r
  colvars r
  hillWeight 0.001
  hillWidth 2.0

harmonicWalls {
  name wall_r
  colvars r
  upperWalls 13.0
  upperWallConstant 2.0

In the colvar r, the distance function used has a lowerBoundary automatically set to 0 Å by default, thus the keyword lowerBoundary itself is not mandatory and hardLowerBoundary is set to yes internally. However, upperBoundary does not have such a ``natural'' choice of value. The metadynamics potential meta_r will individually process any hill whose center is too close to the upperBoundary, more precisely within fewer grid points than 6 times the Gaussian $ \sigma$ parameter plus one. It goes without saying that if the colvar r represents a distance between two freely-moving molecules, it will cross this ``threshold'' rather frequently.

In this example, where the value of hillWidth ($ 2\sigma$ ) amounts to 2 grid points, the threshold is 6+1 = 7 grid points away from upperBoundary. In explicit units, the width of $ r$ is $ w_r =$ 0.2 Å, and the threshold is 15.0 - 7$ \times$ 0.2 = 13.6 Å.

The wall_r restraint included in the example prevents this: the position of its upperWall is 13 Å, i.e. 3 grid points below the buffer's threshold (13.6 Å). For the chosen value of upperWallConstant, the energy of the wall_r bias at r = $ r_{\mathrm{upper}}$ = 13.6 Å is:

$\displaystyle E^* = \frac{1}{2} \- k \left(\frac{r - r_{\mathrm{upper}}}{w_r}\right)^2 = \frac{1}{2} \- 2.0 \left(-3\right)^2 = 9~\mathrm{kcal/mol}$    

which results in a relative probability $ \exp(-E^*/\kappa_{\mathrm{B}}T) \simeq$ $ 3\times{}10^{-7}$ that r crosses the threshold. The probability that r exceeds upperBoundary, which is further away, has also become vanishingly small. At that point, you may want to set hardUpperBoundary to yes for r, and let meta_r know that no special treatment near the grid's boundaries will be needed.

What is the impact of the wall restraint onto the PMF? Not a very complicated one: the PMF reconstructed by metadynamics will simply show a sharp increase in free-energy where the wall potential kicks in (r $ >$ 13 Å). You may then choose between using the PMF only up until that point and discard the rest, or subtracting the energy of the harmonicWalls restraint from the PMF itself. Keep in mind, however, that the statistical convergence of metadynamics may be less accurate where the wall potential is strong.

In summary, although it would be simpler to set the wall's position upperWall and the grid's boundary upperBoundary to the same number, the finite width of the Gaussian hills calls for setting the former strictly within the latter.

Basic configuration keywords

To enable a metadynamics calculation, a metadynamics {...} block must be defined in the Colvars configuration file. Its mandatory keywords are colvars, the variables involved, hillWeight, the weight parameter $ W$ , and the widths $ 2\sigma$ of the Gaussian hills in each dimension given by the single dimensionless parameter hillWidth, or more explicitly by the gaussianSigmas.

Output files

When interpolating grids are enabled (default behavior), the PMF is written by default every colvarsRestartFrequency steps to the file outputName.pmf in multicolumn text format ([*]). The following two options allow to disable or control this behavior and to track statistical convergence:

Performance optimization

The following options control the computational cost of metadynamics calculations, but do not affect results. Default values are chosen to minimize such cost with no loss of accuracy.

Ensemble-Biased Metadynamics

The ensemble-biased metadynamics (EBMetaD) approach [73] is designed to reproduce a target probability distribution along selected collective variables. Standard metadynamics can be seen as a special case of EBMetaD with a flat distribution as target. This is achieved by weighing the Gaussian functions used in the metadynamics approach by the inverse of the target probability distribution:

$\displaystyle V_{\mathrm{EBmetaD}}({\mbox{\boldmath {$\xi$}}}(t)) \; = \; { \su...
...\frac{(\xi_{i}(t)-\xi_{i}(t'))^{2}}{2\sigma_{\xi_{i}}^{2}}\right) } }\mathrm{,}$ (13.32)

where $ \rho_{exp}({\mbox{\boldmath {$\xi$}}})$ is the target probability distribution and $ S_{\rho} = - \int \rho_{exp}({\mbox{\boldmath {$\xi$}}}) \log \rho_{exp}({\mbox{\boldmath {$\xi$}}}) \, \mathrm{d}{\mbox{\boldmath {$\xi$}}}$ its corresponding differential entropy. The method is designed so that during the simulation the resulting distribution of the collective variable $ {\mbox{\boldmath {$\xi$}}}$ converges to $ \rho_{exp}({\mbox{\boldmath {$\xi$}}})$ . A practical application of EBMetaD is to reproduce an ``experimental'' probability distribution, for example the distance distribution between spectroscopic labels inferred from Förster resonance energy transfer (FRET) or double electron-electron resonance (DEER) experiments [73].

The PMF along $ \xi$ can be estimated from the bias potential and the target ditribution [73]:

$\displaystyle A({\mbox{\boldmath {$\xi$}}}) \; \simeq \; { -V_{\mathrm{EBmetaD}...
...{$\xi$}}}) - \kappa_{\mathrm{B}}T \log \rho_{exp}({\mbox{\boldmath {$\xi$}}}) }$ (13.33)

and obtained by enabling writeFreeEnergyFile. Similarly to eq. 13.32, a more accurate estimate of the free energy can be obtained by averaging (after an equilibration time) multiple time-dependent free energy estimates (see keepFreeEnergyFiles).

The following additional options define the configuration for the ensemble-biased metadynamics approach:

As with standard metadynamics, multidimensional probability distributions can be targeted using a single metadynamics block using multiple colvars and a multidimensional target distribution file (see [*]). Instead, multiple probability distributions on different variables can be targeted separately in the same simulation by introducing multiple metadynamics blocks with the ebMeta option.

Example: EBmetaD configuration for a single variable.
colvar {
  name r
  distance {
    group1 { atomNumbers 991 992 }
    group2 { atomNumbers 1762 1763 }
  upperBoundary 100.0
  width 0.1

metadynamics {
  name ebmeta
  colvars r
  hillWeight 0.01
  hillWidth 3.0
  ebMeta on
  targetDistFile targetdist1.dat
  ebMetaEquilSteps 500000

where targetdist1.dat is a text file in ``multicolumn'' format ([*]) with the same width as the variable r (0.1 in this case):
# 1        
# 0.0 0.1 1000 0  
  0.05 0.0012      
  0.15 0.0014      
  ... ...      
  99.95 0.0010      

Tip: Besides setting a meaninful value for targetDistMinVal, the exploration of unphysically low values of the target distribution (which would lead to very large hills and possibly numerical instabilities) can be also prevented by restricting sampling to a given interval, using e.g. harmonicWalls restraint ([*]).

Well-tempered metadynamics

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

Multiple-walker metadynamics

Metadynamics calculations can be performed concurrently by multiple replicas that share a common history. This variant of the method is called multiple-walker metadynamics [75]: the Gaussian hills of all replicas are periodically combined into a single biasing potential, intended to converge to a single PMF.

In the implementation here described [48], replicas communicate through files. This arrangement allows launching the replicas either (1) as a bundle (i.e. a single job in a cluster's queueing system) or (2) as fully independent runs (i.e. as separate jobs for the queueing system). One advantage of the use case (1) is that an identical Colvars configuration can be used for all replicas (otherwise, replicaID needs to be manually set to a different string for each replica). However, the use case (2) is less demanding in terms of high-performance computing resources: a typical scenario would be a computer cluster (including virtual servers from a cloud provider) where not all nodes are connected to each other at high speed, and thus each replica runs on a small group of nodes or a single node.

Whichever way the replicas are started (coupled or not), a shared filesystem is needed so that each replica can read the files created by the others: paths to these files are stored in the shared file replicasRegistry. This file, and those listed in it, are read every replicaUpdateFrequency steps. Each time the Colvars state file is written (for example, colvarsRestartFrequency steps), the file named:
is written as well; this file contains only the state of the metadynamics bias, which the other replicas will read in turn. In between the times when this file is modified/replaced, new hills are also temporarily written to the file named:
Both files are only used for communication, and may be deleted after the replica begins writing files with a new outputName.

Example: Multiple-walker metadynamics with file-based communication.
metadynamics {
  name mymtd
  colvars x
  hillWeight 0.001
  newHillFrequency 1000
  hillWidth 3.0
  multipleReplicas on
  replicasRegistry /shared-folder/mymtd-replicas.txt
  replicaUpdateFrequency 50000 # Best if larger than newHillFrequency

The following are the multiple-walkers related options:

next up previous contents index
Next: Harmonic restraints Up: Biasing and analysis methods Previous: Extended-system Adaptive Biasing Force   Contents   Index