 |
 |
 |
| |
|
|
| 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 [64].
\[ 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 Gaussian kernels. To avoid this issue you can store the bias on a grid. This approach is similar to that proposed in [8] but has the advantage that the grid spacing is independent on the Gaussian width. Notice that you should provide the grid boundaries (GRID_MIN and GRID_MAX) and 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) PLUMED will use 1/5 of the Gaussian width (SIGMA) as grid spacing if the width is fixed or 1/5 of the minimum Gaussian width (SIGMA_MIN) if the width is variable. This default choice should be reasonable for most applications.
Alternatively to the use of grids, it is possible to use a neighbor list to decrease the cost of evaluating the bias, this can be enabled using NLIST. NLIST can be beneficial with more than 2 collective variables, where GRID becomes expensive and memory consuming. The neighbor list will be updated everytime the CVs go farther than a cut-off value from the position they were at last neighbor list update. Gaussians are added to the neigbhor list if their center is within 6.*DP2CUTOFF*sigma*sigma. While the list is updated if the CVs are farther from the center than 0.5 of the standard deviation of the Gaussian center distribution of the list. These parameters (6 and 0.5) can be modified using NLIST_PARAMETERS. Note that the use of neighbor list does not provide the exact bias.
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.
The work performed by the METAD bias can be calculated using CALC_WORK, note that this is expensive when not using grids.
Another option that is available in plumed is well-tempered metadynamics [11]. In this variant of metadynamics the heights of the Gaussian hills are scaled 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 [25] 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 Gaussian kernels 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 [10]. If, for example, metadynamics is performed on a CV s and one is interested only to the free energy for s > boundary, the history dependent potential is still updated according to the above equations but the metadynamics force is set to zero for s < boundary. Notice that Gaussian kernels are added also if s < boundary, as the tails of these Gaussian kernels influence VG in the relevant region s > boundary. In this way, the force on the system in the region s > boundary comes from both metadynamics and the force field, in the region s < boundary 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 [95] 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 [104]. The expression used to calculate \(c(t)\) follows directly from Eq. 3 in [104], where \(F(\vec{s})=-\gamma/(\gamma-1) V(\vec{s})\). This 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-rct) which can be used to obtain a reweighted histogram. The calculation of \(c(t)\) is enabled by using the keyword CALC_RCT. By default \(c(t)\) is updated every time the bias changes, but if this slows down the simulation the keyword RCT_USTRIDE can be set to a value higher than 1. This option requires that a grid is used.
Additional material and examples can be also found in the tutorials:
Concurrent metadynamics as done e.g. in Ref. [43] . This indeed can be obtained by using the METAD action multiple times in the same input file.
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.
d1: DISTANCEATOMS=3,5 d2: DISTANCEthe pair of atom that we are calculating the distance between.ATOMS=2,4 restraint: METADthe pair of atom that we are calculating the distance between.ARG=d1,d2compulsory keyword the labels of the scalars on which the bias will actSIGMA=0.2,0.2compulsory keyword the widths of the Gaussian hillsHEIGHT=0.3the heights of the Gaussian hills.PACE=500 PRINTcompulsory keyword the frequency for hill additionARG=d1,d2,restraint.biascompulsory keyword the labels of the values that you would like to print to the fileSTRIDE=100compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be outputFILE=COLVARthe name of the file on which to output these quantities
d1: DISTANCEATOMS=3,5 d2: DISTANCEthe pair of atom that we are calculating the distance between.ATOMS=2,4 restraint: METADthe pair of atom that we are calculating the distance between.ARG=d1,d2compulsory keyword the labels of the scalars on which the bias will actSIGMA=20compulsory keyword the widths of the Gaussian hillsHEIGHT=0.3the heights of the Gaussian hills.PACE=500compulsory keyword the frequency for hill additionADAPTIVE=DIFF PRINTuse a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme.ARG=d1,d2,restraint.biascompulsory keyword the labels of the values that you would like to print to the fileSTRIDE=100compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be outputFILE=COLVARthe name of the file on which to output these quantities
d1: DISTANCEATOMS=3,5 d2: DISTANCEthe pair of atom that we are calculating the distance between.ATOMS=2,4 restraint: METADthe pair of atom that we are calculating the distance between.ARG=d1,d2compulsory keyword the labels of the scalars on which the bias will actSIGMA=0.05compulsory keyword the widths of the Gaussian hillsHEIGHT=0.3the heights of the Gaussian hills.PACE=500compulsory keyword the frequency for hill additionADAPTIVE=GEOM PRINTuse a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme.ARG=d1,d2,restraint.biascompulsory keyword the labels of the values that you would like to print to the fileSTRIDE=100compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be outputFILE=COLVARthe name of the file on which to output these quantities
d1: DISTANCEATOMS=3,5 d2: DISTANCEthe pair of atom that we are calculating the distance between.ATOMS=2,4 restraint: METAD ...the pair of atom that we are calculating the distance between.ARG=d1,d2compulsory keyword the labels of the scalars on which the bias will actSIGMA=0.05compulsory keyword the widths of the Gaussian hillsHEIGHT=0.3the heights of the Gaussian hills.PACE=500compulsory keyword the frequency for hill additionADAPTIVE=GEOMuse a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme.SIGMA_MIN=0.2,0.1the lower bounds for the sigmas (in CV units) when using adaptive hills.SIGMA_MAX=0.5,1.0 ... PRINTthe upper bounds for the sigmas (in CV units) when using adaptive hills.ARG=d1,d2,restraint.biascompulsory keyword the labels of the values that you would like to print to the fileSTRIDE=100compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be outputFILE=COLVARthe name of the file on which to output these quantities
d1: DISTANCEwhere 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.ATOMS=3,5 restraint: METAD ...the pair of atom that we are calculating the distance between.ARG=d1compulsory keyword the labels of the scalars on which the bias will actSIGMA=0.05compulsory keyword the widths of the Gaussian hillsHEIGHT=0.3the heights of the Gaussian hills.PACE=500compulsory keyword the frequency for hill additionWALKERS_N=10number of walkersWALKERS_ID=3walker idWALKERS_DIR=../shared directory with the hills files from all the walkersWALKERS_RSTRIDE=100 ...stride for reading hills files
phi: TORSIONHere we have asked that the calculation is performed every 10 hills deposition by using RCT_USTRIDE keyword. If this keyword is not given, the calculation will by default be performed every time the bias changes. 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-rct] which can be used to obtain a reweighted histogram or free energy surface using the HISTOGRAM analysis.ATOMS=1,2,3,4 psi: TORSIONthe four atoms involved in the torsional angleATOMS=5,6,7,8 metad: METAD ...the four atoms involved in the torsional angleARG=phi,psicompulsory keyword the labels of the scalars on which the bias will actSIGMA=0.20,0.20compulsory keyword the widths of the Gaussian hillsHEIGHT=1.20the heights of the Gaussian hills.BIASFACTOR=5use well tempered metadynamics and use this bias factor.TEMP=300.0the system temperature - this is only needed if you are doing well-tempered metadynamicsPACE=500compulsory keyword the frequency for hill additionGRID_MIN=-pi,-pithe lower bounds for the gridGRID_MAX=pi,pithe upper bounds for the gridGRID_BIN=150,150the number of bins for the gridCALC_RCT( default=off ) calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-rct].ThisRCT_USTRIDE=10 ...the update stride for calculating the c(t) reweighting factor.The
\[ \tau_{\mathrm{dep}}(t) = \min\left[ \tau_0 \cdot \max\left[\frac{\alpha(t)}{\theta},1\right] ,\tau_{c} \right] \]
where \(\tau_0\) is the initial hill addition frequency given by the PACE keyword, \(\tau_{c}\) is the maximum allowed frequency given by the FA_MAX_PACE keyword, \(\alpha(t)\) is the instantaneous acceleration factor at time \(t\), and \(\theta\) is a threshold value that acceleration factor has to reach before triggering a change in the hill addition frequency given by the FA_MIN_ACCELERATION keyword. The frequency for updating the hill addition frequency according to this equation is given by the FA_UPDATE_FREQUENCY keyword, by default it is the same as the value given in PACE. The hill hill addition frequency increase monotonously such that if the instantaneous acceleration factor is lower than in the previous updating step the previous \(\tau_{\mathrm{dep}}\) is kept rather than updating it to a lower value. The instantaneous hill addition frequency \(\tau_{\mathrm{dep}}(t)\) is outputted to pace component. Note that if restarting from a previous metadynamics run you need to use the ACCELERATION_RFILE keyword to read in the acceleration factors from the previous run, otherwise the hill addition frequency will start from the initial frequency.\[ 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 [44] This choice allows to avoid exceedingly large Gaussian kernels 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 [35] Alternatively, if you use a BIASFACTOR your 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.
d1: DISTANCEATOMS=3,5 t1: METAD ...the pair of atom that we are calculating the distance between.ARG=d1compulsory keyword the labels of the scalars on which the bias will actSIGMA=0.05compulsory keyword the widths of the Gaussian hillsTAU=200in well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tauDAMPFACTOR=100damp hills with exp(-max(V)/(kT*DAMPFACTOR)PACE=250compulsory keyword the frequency for hill additionGRID_MIN=1.14the lower bounds for the gridGRID_MAX=1.32the upper bounds for the gridGRID_BIN=6the number of bins for the gridTARGET=dist.grid ... PRINTtarget to a predefined distributionARG=d1,t1.biascompulsory keyword the labels of the values that you would like to print to the fileSTRIDE=100compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be outputFILE=COLVARthe name of the file on which to output these quantities
The file dist.dat for this calculation would read:
#! FIELDS d1 t1.target der_d1 #! SET min_d1 1.14 #! SET max_d1 1.32 #! SET nbins_d1 6 #! SET periodic_d1 false 1.1400 0.0031 0.1101 1.1700 0.0086 0.2842 1.2000 0.0222 0.6648 1.2300 0.0521 1.4068 1.2600 0.1120 2.6873 1.2900 0.2199 4.6183 1.3200 0.3948 7.1055
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: DISTANCEATOMS=3,5 METADthe pair of atom that we are calculating the distance between.ARG=dcompulsory keyword the labels of the scalars on which the bias will actSIGMA=0.1compulsory keyword the widths of the Gaussian hillsTAU=4.0in well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tauTEMP=300the system temperature - this is only needed if you are doing well-tempered metadynamicsPACE=100compulsory keyword the frequency for hill additionBIASFACTOR=1.0use well tempered metadynamics and use this bias factor.
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: DISTANCEATOMS=3,5 METADthe pair of atom that we are calculating the distance between.ARG=dcompulsory keyword the labels of the scalars on which the bias will actSIGMA=0.1compulsory keyword the widths of the Gaussian hillsTAU=4.0in well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tauTEMP=300the system temperature - this is only needed if you are doing well-tempered metadynamicsPACE=100compulsory keyword the frequency for hill additionRECT=1.0,1.5,2.0,3.0list of bias factors for all the replicas
The number of elements in the RECT array should be equal to the number of replicas.
By default this Action calculates the following quantities. These quantities 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 |
In addition the following quantities can be calculated by employing the keywords listed below
| Quantity | Keyword | Description |
| rbias | CALC_RCT | the instantaneous value of the bias normalized using the c(t) reweighting factor [rbias=bias-rct].This component can be used to obtain a reweighted histogram. |
| rct | CALC_RCT | the reweighting factor c(t). |
| work | CALC_WORK | accumulator for work |
| 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) |
| pace | FREQUENCY_ADAPTIVE | the hill addition frequency when employing frequency adaptive metadynamics |
| nlker | NLIST | number of hills in the neighbor list |
| nlsteps | NLIST | number of steps from last neighbor list update |
| ARG | the labels of the scalars on which the bias will act |
| 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 |
| CALC_WORK | ( default=off ) calculate the total accumulated work done by the bias since last restart |
| CALC_RCT | ( default=off ) calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-rct].This method is not compatible with metadynamics not on a grid. |
| 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 |
| NLIST | ( default=off ) Use neighbor list for kernels summation, faster but experimental |
| WALKERS_MPI | ( default=off ) Switch on MPI version of multiple walkers - not compatible with WALKERS_* options other than WALKERS_DIR |
| FLYING_GAUSSIAN | ( default=off ) Switch on flying Gaussian method, must be used with WALKERS_MPI |
| ACCELERATION | ( default=off ) Set to TRUE if you want to compute the metadynamics acceleration factor. |
| CALC_MAX_BIAS | ( default=off ) Set to TRUE if you want to compute the maximum of the metadynamics V(s, t) |
| CALC_TRANSITION_BIAS | ( default=off ) Set to TRUE if you want to compute a metadynamics transition bias V*(t) |
| FREQUENCY_ADAPTIVE | ( default=off ) Set to TRUE if you want to enable frequency adaptive metadynamics such that the frequency for hill addition to change dynamically based on the acceleration factor. |
| HEIGHT | the heights of the Gaussian hills. Compulsory unless TAU and either BIASFACTOR or DAMPFACTOR 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 bias factor. Please note you must also specify temp |
| RECT | list of bias factors for all the replicas |
| DAMPFACTOR | damp hills with exp(-max(V)/(kT*DAMPFACTOR) |
| TTBIASFACTOR | use transition tempered metadynamics with this bias factor. Please note you must also specify temp |
| TTBIASTHRESHOLD | use transition tempered metadynamics with this bias threshold. Please note you must also specify TTBIASFACTOR |
| TTALPHA | use transition tempered metadynamics with this hill size decay exponent parameter. Please note you must also specify TTBIASFACTOR |
| TARGET | target to a predefined distribution |
| TEMP | the system temperature - this is only needed if you are doing well-tempered metadynamics |
| TAU | in well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tau |
| RCT_USTRIDE | the update stride for calculating the c(t) reweighting factor.The default 1, so c(t) is updated every time the bias is updated. |
| 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) |
| GRID_WSTRIDE | write the grid to a file every N steps |
| GRID_WFILE | the file on which to write the grid |
| GRID_RFILE | a grid file from which the bias should be read at the initial step of the simulation |
| NLIST_PARAMETERS | (default=6.,0.5) the two cutoff parameters for the Gaussians neighbor list |
| ADAPTIVE | use a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme. Sigma is one number that has distance units or time step dimensions |
| 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 |
| 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 | one dimensional lower and upper limits, outside the limits the system will not feel the biasing force. |
| ACCELERATION_RFILE | a data file from which the acceleration should be read at the initial step of the simulation |
| TRANSITIONWELL | This keyword appears multiple times as TRANSITIONWELL followed by an integer. Each specifies the coordinates for one well as in transition-tempered metadynamics. At least one must be provided.. You can use multiple instances of this keyword i.e. TRANSITIONWELL1, TRANSITIONWELL2, TRANSITIONWELL3... |
| FA_UPDATE_FREQUENCY | the frequency for updating the hill addition pace in frequency adaptive metadynamics, by default this is equal to the value given in PACE |
| FA_MAX_PACE | the maximum hill addition frequency allowed in frequency adaptive metadynamics. By default there is no maximum value. |
| FA_MIN_ACCELERATION | only update the hill addition pace in frequency adaptive metadynamics after reaching the minimum acceleration factor given here. By default it is 1.0. |
| 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 |
1.8.17