 |
 |
 |
| 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 [31].
\[ 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 [2] 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.
Another option that is available in plumed is well-tempered metadynamics [4]. 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 [11] 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 [3]. 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:
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 [40] 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 [49]. The expression used to calculate c(t) follows directly from using Eq. 12 in Eq. 3 in [49] 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. [22] . This indeed can be obtained by using the METAD action multiple times in the same input file.
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 |
| ARG | the input for this action is the scalar output from one or more other actions. The particular scalars that you will use are referenced using the label of the action. If the label appears on its own then it is assumed that the Action calculates a single scalar value. The value of this scalar is thus used as the input to this new action. If * or *.* appears the scalars calculated by all the proceding actions in the input file are taken. Some actions have multi-component outputs and each component of the output has a specific label. For example a DISTANCE action labelled dist may have three componets x, y and z. To take just the x component you should use dist.x, if you wish to take all three components then use dist.*.More information on the referencing of Actions can be found in the section of the manual on the PLUMED Getting started. Scalar values can also be referenced using POSIX regular expressions as detailed in the section on Regular Expressions. To use this feature you you must compile PLUMED with the appropriate flag. |
| 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 |
| NUMERICAL_DERIVATIVES | ( default=off ) calculate the derivatives for these quantities numerically |
| GRID_SPARSE | ( default=off ) use a sparse grid to store hills |
| GRID_NOSPLINE | ( default=off ) don't use spline interpolation with grids |
| STORE_GRIDS | ( default=off ) store all the grid files the calculation generates. They will be deleted if this keyword is not present |
| WALKERS_MPI | ( default=off ) Switch on MPI version of multiple walkers - not compatible with other WALKERS_* options |
| ACCELERATION | ( default=off ) Set to TRUE if you want to compute the metadynamics acceleration factor. |
| HEIGHT | the heights of the Gaussian hills. Compulsory unless TAU, TEMP and BIASFACTOR are given |
| FMT | specify format for HILLS files (useful for decrease the number of digits in regtests) |
| BIASFACTOR | use well tempered metadynamics and use this biasfactor. Please note you must also specify temp |
| TEMP | the system temperature - this is only needed if you are doing well-tempered metadynamics |
| TAU | in well tempered metadynamics, sets height to (kb*DeltaT*pace*timestep)/tau |
| GRID_MIN | the lower bounds for the grid |
| GRID_MAX | the upper bounds for the grid |
| GRID_BIN | the number of bins for the grid |
| GRID_SPACING | the approximate grid spacing (to be used as an alternative or together with GRID_BIN) |
| REWEIGHTING_NGRID | calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-c(t)]. Here you should specify the number of grid points required in each dimension. The number of grid points should be equal or larger to the number of grid points given in GRID_BIN.This method is not compatible with metadynamics not on a grid. |
| REWEIGHTING_NHILLS | how many Gaussian hills should be deposited between calculating the c(t) reweighting factor. The default is to do this every 50 hills. |
| GRID_WSTRIDE | write the grid to a file every N steps |
| GRID_WFILE | the file on which to write the grid |
| ADAPTIVE | use a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme. Sigma is one number that has distance units or timestep dimensions |
| WALKERS_ID | walker id |
| WALKERS_N | number of walkers |
| WALKERS_DIR | shared directory with the hills files from all the walkers |
| WALKERS_RSTRIDE | stride for reading hills files |
| INTERVAL | monodimensional lower and upper limits, outside the limits the system will not feel the biasing force. |
| GRID_RFILE | a grid file from which the bias should be read at the initial step of the simulation |
| SIGMA_MAX | the upper bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds |
| SIGMA_MIN | the lower bounds for the sigmas (in CV units) when using adaptive hills. Negative number means no bounds |
| RESTART | allows per-action setting of restart (YES/NO/AUTO) |
| UPDATE_FROM | Only update this action from this time |
| UPDATE_UNTIL | Only update this action until this time |
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(See also DISTANCE PRINT).
DISTANCE ATOMS=3,5 LABEL=d1 DISTANCE ATOMS=2,4 LABEL=d2 METAD ARG=d1,d2 SIGMA=20 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=DIFF PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
DISTANCE ATOMS=3,5 LABEL=d1 DISTANCE ATOMS=2,4 LABEL=d2 METAD ARG=d1,d2 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint ADAPTIVE=GEOM PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
DISTANCE ATOMS=3,5 LABEL=d1 DISTANCE ATOMS=2,4 LABEL=d2 METAD ... 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 ... METAD PRINT ARG=d1,d2,restraint.bias STRIDE=100 FILE=COLVAR
DISTANCE ATOMS=3,5 LABEL=d1 METAD ... ARG=d1 SIGMA=0.05 HEIGHT=0.3 PACE=500 LABEL=restraint WALKERS_N=10 WALKERS_ID=3 WALKERS_DIR=../ WALKERS_RSTRIDE=100 ... METADwhere 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.
METAD ... LABEL=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 ... METADHere 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.
1.8.10 
