This is part of the bias module

Used to performed MetaDynamics on one or more collective variables.

In a metadynamics simulations a history dependent bias composed of intermittently added Gaussian functions is added to the potential [57].

$V(\vec{s},t) = \sum_{ k \tau < t} W(k \tau) \exp\left( -\sum_{i=1}^{d} \frac{(s_i-s_i^{(0)}(k \tau))^2}{2\sigma_i^2} \right).$

This potential forces the system away from the kinetic traps in the potential energy surface and out into the unexplored parts of the energy landscape. Information on the Gaussian functions from which this potential is composed is output to a file called HILLS, which is used both the restart the calculation and to reconstruct the free energy as a function of the CVs. The free energy can be reconstructed from a metadynamics calculation because the final bias is given by:

$V(\vec{s}) = -F(\vec(s))$

During post processing the free energy can be calculated in this way using the sum_hills utility.

In the simplest possible implementation of a metadynamics calculation the expense of a metadynamics calculation increases with the length of the simulation as one has to, at every step, evaluate the values of a larger and larger number of Gaussians. To avoid this issue you can store the bias on a grid. This approach is similar to that proposed in [4] but has the advantage that the grid spacing is independent on the Gaussian width. Notice that you should provide either the number of bins for every collective variable (GRID_BIN) or the desired grid spacing (GRID_SPACING). In case you provide both PLUMED will use the most conservative choice (highest number of bins) for each dimension. In case you do not provide any information about bin size (neither GRID_BIN nor GRID_SPACING) and if Gaussian width is fixed PLUMED will use 1/5 of the Gaussian width as grid spacing. This default choice should be reasonable for most applications.

Metadynamics can be restarted either from a HILLS file as well as from a GRID, in this second case one can first save a GRID using GRID_WFILE (and GRID_WSTRIDE) and at a later stage read it using GRID_RFILE.

Another option that is available in plumed is well-tempered metadynamics [7]. In this varient of metadynamics the heights of the Gaussian hills are rescaled at each step so the bias is now given by:

$V({s},t)= \sum_{t'=0,\tau_G,2\tau_G,\dots}^{t'<t} W e^{-V({s}({q}(t'),t')/\Delta T} \exp\left( -\sum_{i=1}^{d} \frac{(s_i({q})-s_i({q}(t'))^2}{2\sigma_i^2} \right),$

This method ensures that the bias converges more smoothly. It should be noted that, in the case of well-tempered metadynamics, in the output printed the Gaussian height is re-scaled using the bias factor. Also notice that with well-tempered metadynamics the HILLS file does not contain the bias, but the negative of the free-energy estimate. This choice has the advantage that one can restart a simulation using a different value for the $$\Delta T$$. The applied bias will be scaled accordingly.

Note that you can use here also the flexible gaussian approach [22] in which you can adapt the gaussian to the extent of Cartesian space covered by a variable or to the space in collective variable covered in a given time. In this case the width of the deposited gaussian potential is denoted by one value only that is a Cartesian space (ADAPTIVE=GEOM) or a time (ADAPTIVE=DIFF). Note that a specific integration technique for the deposited gaussians should be used in this case. Check the documentation for utility sum_hills.

With the keyword INTERVAL one changes the metadynamics algorithm setting the bias force equal to zero outside boundary [6]. If, for example, metadynamics is performed on a CV s and one is interested only to the free energy for s > sw, the history dependent potential is still updated according to the above equations but the metadynamics force is set to zero for s < sw. Notice that Gaussians are added also if s < sw, as the tails of these Gaussians influence VG in the relevant region s > sw. In this way, the force on the system in the region s > sw comes from both metadynamics and the force field, in the region s < sw only from the latter. This approach allows obtaining a history-dependent bias potential VG that fluctuates around a stable estimator, equal to the negative of the free energy far enough from the boundaries. Note that:

• It works only for one-dimensional biases;
• It works both with and without GRID;
• The interval limit sw in a region where the free energy derivative is not large;
• If in the region outside the limit sw the system has a free energy minimum, the INTERVAL keyword should be used together with a UPPER_WALLS or LOWER_WALLS at sw.

As a final note, since version 2.0.2 when the system is outside of the selected interval the force is set to zero and the bias value to the value at the corresponding boundary. This allows acceptances for replica exchange methods to be computed correctly.

Multiple walkers [75] can also be used. See below the examples.

The c(t) reweighting factor can also be calculated on the fly using the equations presented in [85]. The expression used to calculate c(t) follows directly from using Eq. 12 in Eq. 3 in [85] and gives smoother results than equivalent Eqs. 13 and Eqs. 14 in that paper. The c(t) is given by the rct component while the bias normalized by c(t) is given by the rbias component (rbias=bias-ct) which can be used to obtain a reweighted histogram. The calculation of c(t) is enabled by using the keyword REWEIGHTING_NGRID where the grid used for the calculation is specified. This grid should have a size that is equal or larger than the grid given in GRID_BIN./ By default c(t) is updated every 50 Gaussian hills but this can be changed by using the REWEIGHTING_NHILLS keyword. This option can only be employed together with Well-Tempered Metadynamics and requires that a grid is used.

Additional material and examples can be also found in the tutorials:

Notice that at variance with PLUMED 1.3 it is now straightforward to apply concurrent metadynamics as done e.g. in Ref. [45] . This indeed can be obtained by using the METAD action multiple times in the same input file.

Description of components

By default this Action calculates the following quantities. These quanties can be referenced elsewhere in the input by using this Action's label followed by a dot and the name of the quantity required from the list below.

 Quantity Description bias the instantaneous value of the bias potential work accumulator for work

In addition the following quantities can be calculated by employing the keywords listed below

 Quantity Keyword Description rbias REWEIGHTING_NGRID the instantaneous value of the bias normalized using the $$c(t)$$ reweighting factor [rbias=bias-c(t)].This component can be used to obtain a reweighted histogram. rct REWEIGHTING_NGRID the reweighting factor $$c(t)$$. acc ACCELERATION the metadynamics acceleration factor maxbias CALC_MAX_BIAS the maximum of the metadynamics V(s, t) transbias CALC_TRANSITION_BIAS the metadynamics transition bias V*(t)
Compulsory keywords
 SIGMA the widths of the Gaussian hills PACE the frequency for hill addition FILE ( default=HILLS ) a file in which the list of added hills is stored
Options
Examples

The following input is for a standard metadynamics calculation using as collective variables the distance between atoms 3 and 5 and the distance between atoms 2 and 4. The value of the CVs and the metadynamics bias potential are written to the COLVAR file every 100 steps.

DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
METAD ARG=d1,d2 SIGMA=0.2,0.2 HEIGHT=0.3 PACE=500 LABEL=restraint
PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR


If you use adaptive Gaussians, with diffusion scheme where you use a Gaussian that should cover the space of 20 timesteps in collective variables. Note that in this case the histogram correction is needed when summing up hills.
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR

If you use adaptive Gaussians, with geometrical scheme where you use a Gaussian that should cover the space of 0.05 nm in Cartesian space. Note that in this case the histogram correction is needed when summing up hills.
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR

When using adaptive Gaussians you might want to limit how the hills width can change. You can use SIGMA_MIN and SIGMA_MAX keywords. The sigmas should specified in terms of CV so you should use the CV units. Note that if you use a negative number, this means that the limit is not set. Note also that in this case the histogram correction is needed when summing up hills.
DISTANCE ATOMS=3,5 LABEL=d1
DISTANCE ATOMS=2,4 LABEL=d2
ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM
SIGMA_MIN=0.2,0.1 SIGMA_MAX=0.5,1.0
PRINT ARG=d1,d2,restraint.bias STRIDE=100  FILE=COLVAR

Multiple walkers can be also use as in [75] These are enabled by setting the number of walker used, the id of the current walker which interprets the input file, the directory where the hills containing files resides, and the frequency to read the other walkers. Here is an example
DISTANCE ATOMS=3,5 LABEL=d1
ARG=d1 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint
WALKERS_N=10
WALKERS_ID=3
WALKERS_DIR=../
WALKERS_RSTRIDE=100

where WALKERS_N is the total number of walkers, WALKERS_ID is the id of the present walker (starting from 0 ) and the WALKERS_DIR is the directory where all the walkers are located. WALKERS_RSTRIDE is the number of step between one update and the other. Since version 2.2.5, hills files are automatically flushed every WALKERS_RSTRIDE steps.
The c(t) reweighting factor can be calculated on the fly using the equations presented in [85] as described above. This is enabled by using the keyword REWEIGHTING_NGRID where the grid used for the calculation is set. The number of grid points given in REWEIGHTING_NGRID should be equal or larger than the number of grid points given in GRID_BIN.
METAD ...
ARG=phi,psi SIGMA=0.20,0.20 HEIGHT=1.20 BIASFACTOR=5 TEMP=300.0 PACE=500
GRID_MIN=-pi,-pi GRID_MAX=pi,pi GRID_BIN=150,150
REWEIGHTING_NGRID=150,150
REWEIGHTING_NHILLS=20

Here we have asked that the calculation is performed every 20 hills by using REWEIGHTING_NHILLS keyword. If this keyword is not given the calculation will by default be performed every 50 hills. The c(t) reweighting factor will be given in the rct component while the instantaneous value of the bias potential normalized using the c(t) reweighting factor is given in the rbias component [rbias=bias-c(t)] which can be used to obtain a reweighted histogram or free energy surface using the HISTOGRAM analysis.
The kinetics of the transitions between basins can also be analysed on the fly as in [84]. The flag ACCELERATION turn on accumulation of the acceleration factor that can then be used to determine the rate. This method can be used together with COMMITTOR analysis to stop the simulation when the system get to the target basin. It must be used together with Well-Tempered Metadynamics.
You can also provide a target distribution using the keyword TARGET [96] [65] [46] The TARGET should be a grid containing a free-energy (i.e. the -kbT*log of the desired target distribution). Gaussians will then be scaled by a factor

$e^{\beta(\tilde{F}(s)-\tilde{F}_{max})}$

Here $$\tilde{F}(s)$$ is the free energy defined on the grid and $$\tilde{F}_{max}$$ its maximum value. Notice that we here used the maximum value as in ref [46] This choice allows to avoid exceedingly large Gaussians to be added. However, it could make the Gaussian too small. You should always choose carefully the HEIGHT parameter in this case. The grid file should be similar to other PLUMED grid files in that it should contain both the target free-energy and its derivatives.

Notice that if you wish your simulation to converge to the target free energy you should use the DAMPFACTOR command to provide a global tempering [37] Alternatively, if you use a BIASFACTOR yout simulation will converge to a free energy that is a linear combination of the target free energy and of the intrinsic free energy determined by the original force field.

DISTANCE ATOMS=3,5 LABEL=d1
LABEL=t1
ARG=d1 SIGMA=0.05 TAU=200 DAMPFACTOR=100 PACE=250
GRID_MIN=0 GRID_MAX=2 GRID_BIN=200
TARGET=dist.dat

PRINT ARG=d1,t1.bias STRIDE=100 FILE=COLVAR


#! FIELDS d1 t1.target der_d1
#! SET min_d1 0
#! SET max_d1 2
#! SET nbins_d1  200
#! SET periodic_d1 false


Notice that BIASFACTOR can also be chosen as equal to 1. In this case one will perform unbiased sampling. Instead of using HEIGHT, one should provide the TAU parameter.

d: DISTANCE ATOMS=3,5
METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 BIASFACTOR=1.0


The HILLS file obtained will still work with plumed sum_hills so as to plot a free-energy. The case where this makes sense is probably that of RECT simulations.

Regarding RECT simulations, you can also use the RECT keyword so as to avoid using multiple input files. For instance, a single input file will be

d: DISTANCE ATOMS=3,5
METAD ARG=d SIGMA=0.1 TAU=4.0 TEMP=300 PACE=100 RECT=1.0,1.5,2.0,3.0


The number of elements in the RECT array should be equal to the number of replicas.