Action: METAD

Module bias
Description Usage
Used to performed metadynamics on one or more collective variables. used in 17 tutorialsused in 171 eggs

Output components

This action can calculate the values in the following table when the associated keyword is included in the input for the action. These values can be referenced elsewhere in the input by using this Action's label followed by a dot and the name of the value required from the list below.

Name Type Keyword Description
bias scalar default the instantaneous value of the bias potential
rbias scalar CALC_RCT the instantaneous value of the bias normalized using the c(t) reweighting factor [rbias=bias-rct]
rct scalar CALC_RCT the reweighting factor c(t)
work scalar CALC_WORK accumulator for work
acc scalar ACCELERATION the metadynamics acceleration factor
maxbias scalar CALC_MAX_BIAS the maximum of the metadynamics V(s, t)
transbias scalar CALC_TRANSITION_BIAS the metadynamics transition bias V*(t)
pace scalar FREQUENCY_ADAPTIVE the hill addition frequency when employing frequency adaptive metadynamics
nlker scalar NLIST number of hills in the neighbor list
nlsteps scalar NLIST number of steps from last neighbor list update

Input

The arguments that serve as the input for this action are specified using one or more of the keywords in the following table.

Keyword Type Description
ARG scalar the labels of the scalars on which the bias will act

Further details and examples

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 as discussed in the first paper cited below.

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:

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

The following example provides an example input for a simple metadynamics calculation:

Click on the labels of the actions for more information on what each action computes
tested on2.11
d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,5
d2: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=2,4
restraint: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ARGthe labels of the scalars on which the bias will act=d1,d2 SIGMAthe widths of the Gaussian hills=0.2,0.2 HEIGHTthe heights of the Gaussian hills=0.3 PACEthe frequency for hill addition=500
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d1,d2,restraint.bias STRIDE the frequency with which the quantities of interest should be output=100  FILEthe name of the file on which to output these quantities=COLVAR

In this calculation the distance between atoms 3 and 5 and the distance between atoms 2 and 4 are used as collective variables. Gaussian hills with a heigh of 0.3 kJ/mol are and widths of 0.2 nm in the two CVs are added every 500 steps. The values of the CVs and the metadynamics bias potential are written to the COLVAR file every 100 steps.

Using grids

In the simplest possible implementation of a metadynamics calculation the expense associated with calculating the bias potential 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 the second paper cited below but has the advantage that the grid spacing is independent of the Gaussian width. An example input for a calculation that uses a grid to store the bias is shown below:

Click on the labels of the actions for more information on what each action computes
tested on2.11
t1: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=1,2,3,4
t2: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=5,6,7,8
m: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ARGthe labels of the scalars on which the bias will act=t1,t2 SIGMAthe widths of the Gaussian hills=0.1,0.1 HEIGHTthe heights of the Gaussian hills=1.2 PACEthe frequency for hill addition=500 GRID_MINthe lower bounds for the grid=-pi,-pi GRID_MAXthe upper bounds for the grid=pi,pi GRID_BINthe number of bins for the grid=100,100
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=t1,t2,m.bias STRIDE the frequency with which the quantities of interest should be output=100 FILEthe name of the file on which to output these quantities=COLVAR

Notice that to use the grid you need to provide the grid boundaries (GRID_MIN and GRID_MAX). You can also specify the number of grid points for every collective variable (GRID_BIN) or the desired grid spacing (GRID_SPACING). If you specify the GRID_SPACING and the number of bins PLUMED uses the most conservative choice. In other words, it uses the hoice that gives you the largest number of bins for each dimension.

If you specify neither GRID_BIN nor GRID_SPACING PLUMED uses 1/5 of the Gaussian width (SIGMA) as as grid spacing if the width is fixed or 1/5 of the minimum Gaussian width (SIGMA_MIN - see next section) if the width is variable. This default choice should be reasonable for most applications.

Using neighbour lists

Another way to decrease the cost associated with evaluating the metadynamics bias is to use a neighbor list. If you wish to use this alternative to grids you use the NLIST keyword as illustrated below.

Click on the labels of the actions for more information on what each action computes
tested on2.11
d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2
d2: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,4
d3: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=5,6
d4: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=7,8
m: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ARGthe labels of the scalars on which the bias will act=d1,d2,d3,d4 SIGMAthe widths of the Gaussian hills=0.1,0.1,0.1,0.1 HEIGHTthe heights of the Gaussian hills=1.2 PACEthe frequency for hill addition=500 NLIST Use neighbor list for kernels summation, faster but experimental
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d1,d2,d3,d4,m.bias,m.nlker FILEthe name of the file on which to output these quantities=colvar

Using NLIST is beneficial if you are using more than 2 collective variables as if you are using large numbers of CVs GRID becomes expensive and memory consuming. Furthermore, as indicated above, if you are using this option you can have access to a value called nlker, which tells you how many Gaussians were evaluated in each step.

The neighbor list is updated everytime the CVs move 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.DP2CUTOFFsigma*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.

Restrating a metadynamics simulation

Metadynamics can be restarted from a HILLS file as illustrated in the following input:

Click on the labels of the actions for more information on what each action computes
tested on2.11
#SETTINGS INPUTFILES=regtest/basic/rt-mpi6/HILLS_multi

RESTARTActivate restart. More details
d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,10
c: COORDINATIONCalculate coordination numbers. This action has hidden defaults. More details GROUPAFirst list of atoms=1-108 GROUPBSecond list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)=1-108 R_0The r_0 parameter of the switching function=0.5
m: METADUsed to performed metadynamics on one or more collective variables. More details ARGthe labels of the scalars on which the bias will act=c,d SIGMAthe widths of the Gaussian hills=0.1,0.2 HEIGHTthe heights of the Gaussian hills=0.1 PACEthe frequency for hill addition=500 FILE a file in which the list of added hills is stored=
regtest/basic/rt-mpi6/HILLS_multi
 Click here to see an extract from this file.
  
  ×  
FILE: regtest/basic/rt-mpi6/HILLS_multi
    
  

PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d,c FILEthe name of the file on which to output these quantities=colvar STRIDE the frequency with which the quantities of interest should be output=10



If you don't specify the location of the HILLS file then PLUMED looks for a file called HILLS to read in.
Metadynamics can also be restrated from a GRID, as in this second example.




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


#SETTINGS INPUTFILES=extras/1d_bias.grid

d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2
m: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ...
  ARGthe labels of the scalars on which the bias will act=d1 SIGMAthe widths of the Gaussian hills=0.1
  HEIGHTthe heights of the Gaussian hills=0.1 PACEthe frequency for hill addition=500
  GRID_MINthe lower bounds for the grid=1.14 GRID_MAXthe upper bounds for the grid=1.32 GRID_BINthe number of bins for the grid=5
  GRID_RFILEa grid file from which the bias should be read at the initial step of the simulation=
extras/1d_bias.grid
 Click here to see an extract from this file.
  
  ×  
FILE: extras/1d_bias.grid
    
#! FIELDS d1 m.bias 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

  

  RESTARTallows per-action setting of restart (YES/NO/AUTO)=YES
...
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d1,m.bias FILEthe name of the file on which to output these quantities=colvar STRIDE the frequency with which the quantities of interest should be output=10



Notice that if you want to use this option you first need to save the GRID using GRID_WFILE (and GRID_WSTRIDE).  Using this
option ensures that the grid file that is read in with GRID_RFILE is output.


Notice, also that in the first input above the
RESTART action was used, while the second input used METAD's RESTART keyword.  If you use the first input
every action in the input is restarted so the PRINT command will append to the output file from the earlier calculation.
In this second input, however, only the METAD action is restarted so the PRINT file will output data to a new file and
back up any old files it finds with the same name.


Calculating the work done by the bias


The work performed by the METAD bias can be calculated by using the CALC_WORK option. You can then output the work as shown below:




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


mu1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,10
vol: VOLUMECalculate the volume the simulation box. More details

md: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ARGthe labels of the scalars on which the bias will act=mu1,vol SIGMAthe widths of the Gaussian hills=0.1,0.2 HEIGHTthe heights of the Gaussian hills=1.0 PACEthe frequency for hill addition=10 GRID_MINthe lower bounds for the grid=-5,20 GRID_MAXthe upper bounds for the grid=5,40 GRID_BINthe number of bins for the grid=100,100 CALC_WORK calculate the total accumulated work done by the bias since last restart

PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=mu1,vol,md.bias,md.work FILEthe name of the file on which to output these quantities=colvar



Note that this calculation is expensive when not using grids.


Well-tempered metadynamics


Another option that is available in plumed is well-tempered metadynamics that is introduced in the third paper cited below. In this
variant of metadynamics the heights of the Gaussian hills are scaled at each step so the bias is now
given by:







This method ensures that the bias converges more smoothly. The following example shows an example input for a well-tempered metadynamics
calculation.




 Click on the labels of the actions for more information on what each action computes 


tested on2.11




d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,5
d2: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=2,4
restraint: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ARGthe labels of the scalars on which the bias will act=d1,d2 SIGMAthe widths of the Gaussian hills=0.2,0.2 HEIGHTthe heights of the Gaussian hills=0.3 BIASFACTORuse well tempered metadynamics and use this bias factor=10 PACEthe frequency for hill addition=500
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d1,d2,restraint.bias STRIDE the frequency with which the quantities of interest should be output=100  FILEthe name of the file on which to output these quantities=COLVAR



It should be noted that, if you do well-tempered metadynamics,
the hills file output contains the value of Gaussian height re-scaled by the bias factor.
This rescaling is done to ensure that the HILLS file contains
negative of the free-energy estimate and not the bias potential. An added bonus of this choice is that
one can restart a simulation using a different value for the . The applied bias will be scaled accordingly.


Metadynamics with flexible Gaussians


You can also use the flexible Gaussian approach that is discussed in the fourth paper cited below.
This method adapts the Gaussian to the extent of Cartesian space covered by a variable or
to the space in collective variable covered in a given time. When you use these methods 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 that is discussed in the documenation for sum_hills
mist be used when calculating the free energy from the deposited Gaussian kernels
should be used in this case.


The following example shows how to use adaptive Gaussian kernels, with the diffusion scheme.  In this example
Gaussians that cover the space of 20 time steps in collective variables are used.




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,5
d2: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=2,4
restraint: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ARGthe labels of the scalars on which the bias will act=d1,d2 SIGMAthe widths of the Gaussian hills=20 HEIGHTthe heights of the Gaussian hills=0.3 PACEthe frequency for hill addition=500 ADAPTIVEuse a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme=DIFF
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d1,d2,restraint.bias STRIDE the frequency with which the quantities of interest should be output=100  FILEthe name of the file on which to output these quantities=COLVAR



This second example shows you how to use adaptive Gaussian kernels, with the geometrical scheme.  In this example
Gaussians that cover the space of 0.05 nm in Cartesian space are used




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,5
d2: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=2,4
restraint: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ARGthe labels of the scalars on which the bias will act=d1,d2 SIGMAthe widths of the Gaussian hills=0.05 HEIGHTthe heights of the Gaussian hills=0.3 PACEthe frequency for hill addition=500 ADAPTIVEuse a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme=GEOM
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d1,d2,restraint.bias STRIDE the frequency with which the quantities of interest should be output=100  FILEthe name of the file on which to output these quantities=COLVAR



With both these methods you can limit how much the hills width can change by using the SIGMA_MIN and SIGMA_MAX keywords
as shown below.




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,5
d2: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=2,4
restraint: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ...
  ARGthe labels of the scalars on which the bias will act=d1,d2 SIGMAthe widths of the Gaussian hills=0.05 HEIGHTthe heights of the Gaussian hills=0.3 PACEthe frequency for hill addition=500 ADAPTIVEuse a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme=GEOM
  SIGMA_MINthe lower bounds for the sigmas (in CV units) when using adaptive hills=0.2,0.1 SIGMA_MAXthe upper bounds for the sigmas (in CV units) when using adaptive hills=0.5,1.0
...
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d1,d2,restraint.bias STRIDE the frequency with which the quantities of interest should be output=100  FILEthe name of the file on which to output these quantities=COLVAR



The sigma values that are provided to these two keywords are specified in the units of the CV.
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.


Metadynamics with the interval keyword


When you use the keyword INTERVAL the metadynamics bias force equal is set equal to zero
outside of a boundary as discussed in the fifth paper cited below. This feature is useful if, for example, metadynamics is performed on one or more CVs and if one is only interested
in the free energy for s > boundary.  When you use this option the history dependent potential is still updated in accordance with the forumulas above.
The difference is simply that the metadynamics force is set to zero for s < boundary. This means 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, while in the region
s < boundary the only force the system experiences if from the forcefield. This approach allows one to obtain 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. Notice that:


    

  • This method only works for one-dimensional biases;
    • 

  • It works both with and without GRID;
    • 

  • The interval limit boundary should be set in a region where the derivative of the free energy is not large;
    • 

  • If in the region outside the limit boundary the system has a free energy minimum, the INTERVAL keyword should
  be used together with a UPPER_WALLS or LOWER_WALLS that act when the system is beyond the boundary as illustrated in the following example.
    • 





 Click on the labels of the actions for more information on what each action computes 


tested on2.11


r2g: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2
m2g: METADUsed to performed metadynamics on one or more collective variables. More details ARGthe labels of the scalars on which the bias will act=r2g HEIGHTthe heights of the Gaussian hills=0.1 PACEthe frequency for hill addition=500 SIGMAthe widths of the Gaussian hills=1.0 INTERVALone dimensional lower and upper limits, outside the limits the system will not feel the biasing force=4,8 FILE a file in which the list of added hills is stored=H2G GRID_MINthe lower bounds for the grid=0 GRID_MAXthe upper bounds for the grid=20 GRID_BINthe number of bins for the grid=100
lw: LOWER_WALLSDefines a wall for the value of one or more collective variables, More details ARGthe arguments on which the bias is acting=r2g ATthe positions of the wall=4 KAPPAthe force constant for the wall=100
uw: UPPER_WALLSDefines a wall for the value of one or more collective variables, More details ARGthe arguments on which the bias is acting=r2g ATthe positions of the wall=8 KAPPAthe force constant for the wall=100
PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=r2g,*.bias FILEthe name of the file on which to output these quantities=colvar STRIDE the frequency with which the quantities of interest should be output=10



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


Multiple walkers metadynamics


The multiple walkers metadynamics method that is described in the sixth paper cited below can be used by using an input
like the one shown below:




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,5
restraint: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ...
   ARGthe labels of the scalars on which the bias will act=d1 SIGMAthe widths of the Gaussian hills=0.05 HEIGHTthe heights of the Gaussian hills=0.3 PACEthe frequency for hill addition=500
   WALKERS_Nnumber of walkers=10
   WALKERS_IDwalker id=3
   WALKERS_DIRshared directory with the hills files from all the walkers=../
   WALKERS_RSTRIDEstride for reading hills files=100
...



To use this method you need to set the number of walkers (WALKERS_N), the ID of the
current walker (WALKERS-ID), the directory where the files containing the hills files reside
(WALKERS_DIR) and the frequency to read the other walkers (WALKERS_RSTRIDE).
The bias potential experienced by all the walkers is synched between walkers every
WALKERS_RSTRIDE steps.  In addition, since version 2.2.5, hills files are automatically
flushed every WALKERS_RSTRIDE steps.


Reweighting a metadynamics simulation


The  reweighting factor can also be calculated on the fly using the equations
presented in the seventh paper in the references section by using the keyword CALC_RCT as illustrated in the example input below:




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


phi: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=1,2,3,4
psi: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=5,6,7,8

m: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ...
 ARGthe labels of the scalars on which the bias will act=phi,psi SIGMAthe widths of the Gaussian hills=0.20,0.20 HEIGHTthe heights of the Gaussian hills=1.20 BIASFACTORuse well tempered metadynamics and use this bias factor=5 TEMPthe system temperature - this is only needed if you are doing well-tempered metadynamics=300.0 PACEthe frequency for hill addition=500
 GRID_MINthe lower bounds for the grid=-pi,-pi GRID_MAXthe upper bounds for the grid=pi,pi GRID_BINthe number of bins for the grid=150,150
 CALC_RCT calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-rct]
 RCT_USTRIDEthe update stride for calculating the c(t) reweighting factor=10
...

rw: REWEIGHT_METADCalculate the weights configurations should contribute to the histogram in a simulation in which a metadynamics bias acts upon the system. This action has hidden defaults. More details TEMPthe system temperature=300

h: HISTOGRAMAccumulate the average probability density along a few CVs from a trajectory. This action is a shortcut and it has hidden defaults. More details ...
   ARGthe quantities that are being used to construct the histogram=phi,psi
   GRID_MIN the lower bounds for the grid=-pi,-pi GRID_MAX the upper bounds for the grid=pi,pi GRID_BINthe number of bins for the grid=150,150
   BANDWIDTHthe bandwidths for kernel density esimtation=0.1,0.1
   LOGWEIGHTSthe logarithm of the quantity to use as the weights when calculating averages=rw
...

ff: CONVERT_TO_FESConvert a histogram to a free energy surface. This action is a shortcut. More details ARGthe histogram that you would like to convert into a free energy surface=h TEMPthe temperature at which you are operating=300
DUMPGRIDOutput the function on the grid to a file with the PLUMED grid format. More details ARGthe label for the grid that you would like to output=ff FILE the file on which to write the grid=fes.dat




The expression used to calculate  follows directly from Eq. 3 in the seventh paper cited below,
where .
This gives smoother results than Eqs. 13 and Eqs. 14 in that paper.


Notice that you can only use this option if you are accumulating the metadynamics bias on a grid.
 is given by the rct component of the METAD action above while the bias
normalized by  is given by the rbias component (rbias=bias-rct).  As illustrated above this
value  can be used to obtain a reweighted histogram.


By default  is updated every time the bias changes, but if this slows down the simulation.
If you use the keyword RCT_USTRIDE and set to a value higher than 1 the calculation is sped up.
In the example input above we are calcluating  after every Gaussian deposition event.


Concurrent metadynamics


To implement the method of concurrent metadynamics that is discussed in the eigth paper cited below you insert multiple
METAD actions in your plumed input file.


Infrequent metadynamics


The kinetics of the transitions between basins can also be analyzed on the fly as
in discussed in the 9th paper cited below. The flag ACCELERATION turn on accumulation of the acceleration
factor that can then be used to determine the rate as illustrated by the following input.




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


psi: ANGLECalculate one or multiple angle/s. More details ATOMSthe list of atoms involved in this collective variable (either 3 or 4 atoms)=7,9,15

md1: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ...
 ARGthe labels of the scalars on which the bias will act=psi SIGMAthe widths of the Gaussian hills=0.20
 HEIGHTthe heights of the Gaussian hills=1.20 TEMPthe system temperature - this is only needed if you are doing well-tempered metadynamics=300.0 BIASFACTORuse well tempered metadynamics and use this bias factor=10
 PACEthe frequency for hill addition=500 GRID_MINthe lower bounds for the grid=0 GRID_MAXthe upper bounds for the grid=pi GRID_BINthe number of bins for the grid=10
 ACCELERATION Set to TRUE if you want to compute the metadynamics acceleration factor
...

PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=psi,md1.bias,md1.acc FILEthe name of the file on which to output these quantities=colvar STRIDE the frequency with which the quantities of interest should be output=10



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. If restarting from a previous
metadynamics you need to use the ACCELERATION_RFILE keyword to give the name of the
data file from which the previous value of the acceleration factor should be read, otherwise the
calculation of the acceleration factor will be wrong.


Changing the deposition frequency adaptively


By using the flag FREQUENCY_ADAPTIVE the frequency adaptive scheme introduced in the 10th paper cited below
is turned on. The frequency for hill addition then changes dynamically based on the acceleration factor
according to the following equation







where  is the initial hill addition frequency given by the PACE keyword,
 is the maximum allowed frequency given by the FA_MAX_PACE keyword,
 is the instantaneous acceleration factor at time ,
and  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  is kept rather than updating it to a lower value.
The instantaneous hill addition frequency  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.


Targetted metadynamics


You can also provide a target distribution for your metadynamics simulations by using the keyword TARGET.
This idea is discussed in the last four papers cited in the references section.


The TARGET should be a grid containing a free-energy (i.e. the -T*log for the desired target distribution).
Gaussian kernels will then be scaled by a factor







Here  is the free energy defined on the grid and  its maximum value.
Notice that we here used the maximum value as was done in the penultimate paper in the reference section.
This choice allows one to avoid the addition of exceedingly large Gaussian kernels. However,
it can also make the Gaussian too small. You should always choose the HEIGHT parameter carefully
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 as is discussed in the last paper from the reference section.
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.




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


#SETTINGS INPUTFILES=extras/target.fes

d1: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,5
t1: METADUsed to performed metadynamics on one or more collective variables. This action has hidden defaults. More details ...
 ARGthe labels of the scalars on which the bias will act=d1 SIGMAthe widths of the Gaussian hills=0.05 TAUin well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tau=200 DAMPFACTORdamp hills with exp(-max(V)/(kT*DAMPFACTOR)=100 PACEthe frequency for hill addition=250
 GRID_MINthe lower bounds for the grid=1.14 GRID_MAXthe upper bounds for the grid=1.32 GRID_BINthe number of bins for the grid=6
 TARGETtarget to a predefined distribution=
extras/target.fes
 Click here to see an extract from this file.
  
  ×  
FILE: extras/target.fes
    
#! 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

  

...

PRINTPrint quantities to a file. More details ARGthe labels of the values that you would like to print to the file=d1,t1.bias STRIDE the frequency with which the quantities of interest should be output=100 FILEthe name of the file on which to output these quantities=COLVAR



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




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,5
METADUsed to performed metadynamics on one or more collective variables. More details ARGthe labels of the scalars on which the bias will act=d SIGMAthe widths of the Gaussian hills=0.1 TAUin well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tau=4.0 TEMPthe system temperature - this is only needed if you are doing well-tempered metadynamics=300 PACEthe frequency for hill addition=100 BIASFACTORuse well tempered metadynamics and use this bias factor=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




 Click on the labels of the actions for more information on what each action computes 


tested on2.11


d: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=3,5
METADUsed to performed metadynamics on one or more collective variables. More details ARGthe labels of the scalars on which the bias will act=d SIGMAthe widths of the Gaussian hills=0.1 TAUin well tempered metadynamics, sets height to (k_B Delta T*pace*timestep)/tau=4.0 TEMPthe system temperature - this is only needed if you are doing well-tempered metadynamics=300 PACEthe frequency for hill addition=100 RECTlist of bias factors for all the replicas=1.0,1.5,2.0,3.0



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


References


More information about how this action can be used is available in the following articles:



Syntax


The following table describes the keywords and options that can be used with this action 






Keyword
Type
Default
Description






ARG
input
none
the labels of the scalars on which the bias will act




SIGMA
compulsory
none
the widths of the Gaussian hills




PACE
compulsory
none
the frequency for hill addition




FILE
compulsory
HILLS
a file in which the list of added hills is stored




NUMERICAL_DERIVATIVES
optional
false
calculate the derivatives for these quantities numerically




HEIGHT
optional
not used
the heights of the Gaussian hills




FMT
optional
not used
specify format for HILLS files (useful for decrease the number of digits in regtests)




BIASFACTOR
optional
not used
use well tempered metadynamics and use this bias factor




CALC_WORK
optional
false
calculate the total accumulated work done by the bias since last restart




RECT
optional
not used
list of bias factors for all the replicas




DAMPFACTOR
optional
not used
damp hills with exp(-max(V)/(kT*DAMPFACTOR)




TTBIASFACTOR
optional
not used
use transition tempered metadynamics with this bias factor




TTBIASTHRESHOLD
optional
not used
use transition tempered metadynamics with this bias threshold




TTALPHA
optional
not used
use transition tempered metadynamics with this hill size decay exponent parameter




TARGET
optional
not used
target to a predefined distribution




TEMP
optional
not used
the system temperature - this is only needed if you are doing well-tempered metadynamics




TAU
optional
not used
in well tempered metadynamics, sets height to (k_B Delta Tpacetimestep)/tau




CALC_RCT
optional
false
calculate the c(t) reweighting factor and use that to obtain the normalized bias [rbias=bias-rct]




RCT_USTRIDE
optional
not used
the update stride for calculating the c(t) reweighting factor




GRID_MIN
optional
not used
the lower bounds for the grid




GRID_MAX
optional
not used
the upper bounds for the grid




GRID_BIN
optional
not used
the number of bins for the grid




GRID_SPACING
optional
not used
the approximate grid spacing (to be used as an alternative or together with GRID_BIN)




GRID_SPARSE
optional
false
use a sparse grid to store hills




GRID_NOSPLINE
optional
false
don't use spline interpolation with grids




GRID_WSTRIDE
optional
not used
write the grid to a file every N steps




GRID_WFILE
optional
not used
the file on which to write the grid




GRID_RFILE
optional
not used
a grid file from which the bias should be read at the initial step of the simulation




STORE_GRIDS
optional
false
store all the grid files the calculation generates




NLIST
optional
false
Use neighbor list for kernels summation, faster but experimental




NLIST_PARAMETERS
optional
not used
the two cutoff parameters for the Gaussians neighbor list




ADAPTIVE
optional
not used
use a geometric (=GEOM) or diffusion (=DIFF) based hills width scheme




SIGMA_MAX
optional
not used
the upper bounds for the sigmas (in CV units) when using adaptive hills




SIGMA_MIN
optional
not used
the lower bounds for the sigmas (in CV units) when using adaptive hills




WALKERS_ID
optional
not used
walker id




WALKERS_N
optional
not used
number of walkers




WALKERS_DIR
optional
not used
shared directory with the hills files from all the walkers




WALKERS_RSTRIDE
optional
not used
stride for reading hills files




WALKERS_MPI
optional
false
Switch on MPI version of multiple walkers - not compatible with WALKERS_* options other than WALKERS_DIR




INTERVAL
optional
not used
one dimensional lower and upper limits, outside the limits the system will not feel the biasing force




FLYING_GAUSSIAN
optional
false
Switch on flying Gaussian method, must be used with WALKERS_MPI




ACCELERATION
optional
false
Set to TRUE if you want to compute the metadynamics acceleration factor




ACCELERATION_RFILE
optional
not used
a data file from which the acceleration should be read at the initial step of the simulation




CALC_MAX_BIAS
optional
false
Set to TRUE if you want to compute the maximum of the metadynamics V(s, t)




CALC_TRANSITION_BIAS
optional
false
Set to TRUE if you want to compute a metadynamics transition bias V*(t)




TRANSITIONWELL
optional
not used
This keyword appears multiple times as TRANSITIONWELL followed by an integer




FREQUENCY_ADAPTIVE
optional
false
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




FA_UPDATE_FREQUENCY
optional
not used
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
optional
not used
the maximum hill addition frequency allowed in frequency adaptive metadynamics




FA_MIN_ACCELERATION
optional
not used
only update the hill addition pace in frequency adaptive metadynamics after reaching the minimum acceleration factor given here




RESTART
optional
not used
allows per-action setting of restart (YES/NO/AUTO)




UPDATE_FROM
optional
not used
Only update this action from this time




UPDATE_UNTIL
optional
not used
Only update this action until this time