Action: OPES_METAD

Module	opes
Description	Usage
On-the-fly probability enhanced sampling with metadynamics-like target distribution.

Details and examples

On-the-fly probability enhanced sampling with metadynamics-like target distribution.

This On-the-fly probability enhanced sampling (see OPES) method with metadynamics-like target distribution is described in the first paper cited below.

This OPES_METAD action samples target distributions defined via their marginal $p^{\text{tg}}(\mathbf{s})$ over some collective variables (CVs), $\mathbf{s}=\mathbf{s}(\mathbf{x})$ . By default OPES_METAD targets the well-tempered distribution, $p^{\text{WT}}(\mathbf{s})\propto [P(\mathbf{s})]^{1/\gamma}$ , where $\gamma$ is known as BIASFACTOR. Similarly to METAD, OPES_METAD optimizes the bias on-the-fly, with a given PACE. It does so by reweighting via kernel density estimation the unbiased distribution in the CV space, $P(\mathbf{s})$ . A compression algorithm is used to prevent the number of kernels from growing linearly with the simulation time. The bias at step $n$ is

$V_n(\mathbf{s}) = (1-1/\gamma)\frac{1}{\beta}\log\left(\frac{P_n(\mathbf{s})}{Z_n}+\epsilon\right)\, .$

See the first paper cited below for a complete description of the method.

As an intuitive picture, rather than gradually filling the metastable basins, OPES_METAD quickly tries to get a coarse idea of the full free energy surface (FES), and then slowly refines its details. It has a fast initial exploration phase, and then becomes extremely conservative and does not significantly change the shape of the deposited bias any more, reaching a regime of quasi-static bias. For this reason, it is possible to use standard umbrella sampling reweighting (see REWEIGHT_BIAS) to analyse the trajectory. At this link you can find some python scripts that work in a similar way to sum_hills, but the preferred way to obtain a FES with OPES is via reweighting (this is discussed in the tutorials you can access by clicking the button above). The estimated $c(t)$ is printed for reference only, since it should converge to a fixed value as the bias converges. This $c(t)$ should NOT be used for reweighting. Similarly, the $Z_n$ factor is printed only for reference, and it should converge when no new region of the CV-space is explored.

Notice that OPES_METAD is more sensitive to degenerate CVs than METAD. If the employed CVs map different metastable basins onto the same CV-space region, then OPES_METAD will remain stuck rather than completely reshaping the bias. This can be useful to diagnose problems with your collective variable. If it is not possible to improve the set of CVs and remove this degeneracy, then you might instead consider to use OPES_METAD_EXPLORE or METAD. In this way you will be able to obtain an estimate of the FES, but be aware that you most likely will not reach convergence and thus this estimate will be subjected to systematic errors (see e.g. Fig.3 in the second paper cited below). On the contrary, if your CVs are not degenerate but only suboptimal, you should converge faster by using OPES_METAD instead of METAD as discussed in the first paper cited below.

The parameter BARRIER should be set to be at least equal to the highest free energy barrier you wish to overcome. If it is much lower than that, you will not cross the barrier, if it is much higher, convergence might take a little longer. If the system has a basin that is clearly more stable than the others, it is better to start the simulation from there.

By default, the kernels SIGMA is adaptive, estimated from the fluctuations over ADAPTIVE_SIGMA_STRIDE simulation steps (similar to METAD ADAPTIVE=DIFF, but contrary to that, no artifacts are introduced and the bias will converge to the correct one). However, notice that depending on the system this might not be the optimal choice for SIGMA.

You can target a uniform flat distribution by explicitly setting BIASFACTOR=inf. However, this should be useful only in very specific cases.

It is possible to take into account also of other bias potentials besides the one of OPES_METAD during the internal reweighting for \f $P(\mathbf{s})\f$ estimation. To do so, one has to add those biases with the EXTRA_BIAS keyword, as in the example below. This allows one to define a custom target distribution by adding another bias potential equal to the desired target free energy and setting BIASFACTOR=inf (see example below). Another possible usage of EXTRA_BIAS is to make sure that OPES_METAD does not push against another fixed bias added to restrain the CVs range (e.g. UPPER_WALLS).

Through the EXCLUDED_REGION keywork, it is possible to specify a region of CV space where no kernels will be deposited. This can be useful for example for making sure the bias does not modify the transition region, thus allowing for rate calculation. See below for an example of how to use this keyword.

Restart can be done from a KERNELS file, but it might be not perfect (due to limited precision when printing kernels to file, or if adaptive SIGMA is used). For an exact restart you must use STATE_RFILE to read a checkpoint with all the needed info. To save such checkpoints, define a STATE_WFILE and choose how often to print them with STATE_WSTRIDE. By default this file is overwritten, but you can instead append to it using the flag STORE_STATES.

Multiple walkers are supported only with MPI communication, via the keyword WALKERS_MPI.

Examples

Several examples can be found on the PLUMED-NEST website, by searching for the OPES keyword. The tutorials on this method that can be accessed by clicking the tutorial button at the top of this page can also be useful to get started with the method.

The following is a minimal working example:

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

cvThe DISTANCE action with label cv calculates the following quantities: Quantity    Type    Description  
cv scalar the DISTANCE between this pair of atoms: DISTANCECalculate the distance/s between pairs of atoms. More details ATOMSthe pair of atom that we are calculating the distance between=1,2
opesThe OPES_METAD action with label opes calculates the following quantities: Quantity    Type    Description  
opes.bias scalar the instantaneous value of the bias potential
opes.rct scalar estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
opes.zed scalar estimate of Z_n. should become flat once no new CV-space region is explored
opes.neff scalar effective sample size
opes.nker scalar total number of compressed kernels used to represent the bias: OPES_METADOn-the-fly probability enhanced sampling with metadynamics-like target distribution. This action has hidden defaults. More details ARGthe labels of the scalars on which the bias will act=cv PACEthe frequency for kernel deposition=200 BARRIERthe free energy barrier to be overcome=40
opes: OPES_METADOn-the-fly probability enhanced sampling with metadynamics-like target distribution. This action uses the defaults shown here. More details ARGthe labels of the scalars on which the bias will act=cv PACEthe frequency for kernel deposition=200 BARRIERthe free energy barrier to be overcome=40  TEMP temperature=-1 SIGMA the initial widths of the kernels=ADAPTIVE COMPRESSION_THRESHOLD merge kernels if closer than this threshold, in units of sigma=1 FILE a file in which the list of all deposited kernels is stored=KERNELS
PRINTPrint quantities to a file. More details STRIDE the frequency with which the quantities of interest should be output=200 FILEthe name of the file on which to output these quantities=COLVAR ARGthe labels of the values that you would like to print to the file=*

Another more articulated one:

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

phiThe TORSION action with label phi calculates the following quantities: Quantity    Type    Description  
phi scalar the TORSION involving these atoms: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=5,7,9,15
psiThe TORSION action with label psi calculates the following quantities: Quantity    Type    Description  
psi scalar the TORSION involving these atoms: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=7,9,15,17
opesThe OPES_METAD action with label opes calculates the following quantities: Quantity    Type    Description  
opes.bias scalar the instantaneous value of the bias potential
opes.rct scalar estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
opes.zed scalar estimate of Z_n. should become flat once no new CV-space region is explored
opes.neff scalar effective sample size
opes.nker scalar total number of compressed kernels used to represent the bias
opes.nlker scalar number of kernels in the neighbor list
opes.nlsteps scalar number of steps from last neighbor list update: OPES_METADOn-the-fly probability enhanced sampling with metadynamics-like target distribution. This action has hidden defaults. More details ...
  FILE a file in which the list of all deposited kernels is stored=Kernels.data
  TEMP temperature=300
  ARGthe labels of the scalars on which the bias will act=phi,psi
  PACEthe frequency for kernel deposition=500
  BARRIERthe free energy barrier to be overcome=50
  SIGMA the initial widths of the kernels=0.15,0.15
  SIGMA_MINnever reduce SIGMA below this value=0.01,0.01
  STATE_RFILEread from this file the compressed kernels and all the info needed to RESTART the simulation=Restart.data
  STATE_WFILEwrite to this file the compressed kernels and all the info needed to RESTART the simulation=State.data
  STATE_WSTRIDEnumber of MD steps between writing the STATE_WFILE=500*100
  STORE_STATES append to STATE_WFILE instead of ovewriting it each time
  WALKERS_MPI switch on MPI version of multiple walkers
  NLIST use neighbor list for kernels summation, faster but experimental
...
opes: OPES_METADOn-the-fly probability enhanced sampling with metadynamics-like target distribution. This action uses the defaults shown here. More details ...
  FILE a file in which the list of all deposited kernels is stored=Kernels.data
  TEMP temperature=300
  ARGthe labels of the scalars on which the bias will act=phi,psi
  PACEthe frequency for kernel deposition=500
  BARRIERthe free energy barrier to be overcome=50
  SIGMA the initial widths of the kernels=0.15,0.15
  SIGMA_MINnever reduce SIGMA below this value=0.01,0.01
  STATE_RFILEread from this file the compressed kernels and all the info needed to RESTART the simulation=Restart.data
  STATE_WFILEwrite to this file the compressed kernels and all the info needed to RESTART the simulation=State.data
  STATE_WSTRIDEnumber of MD steps between writing the STATE_WFILE=500*100
  STORE_STATES append to STATE_WFILE instead of ovewriting it each time
  WALKERS_MPI switch on MPI version of multiple walkers
  NLIST use neighbor list for kernels summation, faster but experimental
 COMPRESSION_THRESHOLD merge kernels if closer than this threshold, in units of sigma=1
...
PRINTPrint quantities to a file. More details FMT the format that should be used to output real numbers=%g STRIDE the frequency with which the quantities of interest should be output=500 FILEthe name of the file on which to output these quantities=Colvar.data ARGthe labels of the values that you would like to print to the file=phi,psi,opes.*

Next is an example of how to define a custom target distribution different from the well-tempered one. Here we chose to focus more on the transition state, that is around $\phi=0$ . Our target distribution is a Gaussian centered there, thus the target free energy we want to sample is a parabola, $F^{\text{tg}}(\mathbf{s})=-\frac{1}{\beta} \log [p^{\text{tg}}(\mathbf{s})]$ .

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

phiThe TORSION action with label phi calculates the following quantities: Quantity    Type    Description  
phi scalar the TORSION involving these atoms: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=5,7,9,15
FtgValueThe CUSTOM action with label FtgValue calculates the following quantities: Quantity    Type    Description  
FtgValue scalar an arbitrary function: CUSTOMCalculate a combination of variables using a custom expression. More details ARGthe values input to this function=phi PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO FUNCthe function you wish to evaluate=(x/0.4)^2
FtgThe BIASVALUE action with label Ftg calculates the following quantities: Quantity    Type    Description  
Ftg.bias scalar the instantaneous value of the bias potential
Ftg.FtgValue_bias scalar one or multiple instances of this quantity can be referenced elsewhere in the input file. these quantities will named with  the arguments of the bias followed by the character string _bias. These quantities tell the user how much the bias is due to each of the colvars. This particular component measures this quantity for the input CV named FtgValue: BIASVALUETakes the value of one variable and use it as a bias More details ARGthe labels of the scalar/vector arguments whose values will be used as a bias on the system=FtgValue
opesThe OPES_METAD action with label opes calculates the following quantities: Quantity    Type    Description  
opes.bias scalar the instantaneous value of the bias potential
opes.rct scalar estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
opes.zed scalar estimate of Z_n. should become flat once no new CV-space region is explored
opes.neff scalar effective sample size
opes.nker scalar total number of compressed kernels used to represent the bias: OPES_METADOn-the-fly probability enhanced sampling with metadynamics-like target distribution. This action has hidden defaults. More details ...
  ARGthe labels of the scalars on which the bias will act=phi
  PACEthe frequency for kernel deposition=500
  BARRIERthe free energy barrier to be overcome=50
  SIGMA the initial widths of the kernels=0.2
  BIASFACTORthe gamma bias factor used for the well-tempered target distribution=inf
  EXTRA_BIASconsider also these other bias potentials for the internal reweighting=Ftg.bias
...
opes: OPES_METADOn-the-fly probability enhanced sampling with metadynamics-like target distribution. This action uses the defaults shown here. More details ...
  ARGthe labels of the scalars on which the bias will act=phi
  PACEthe frequency for kernel deposition=500
  BARRIERthe free energy barrier to be overcome=50
  SIGMA the initial widths of the kernels=0.2
  BIASFACTORthe gamma bias factor used for the well-tempered target distribution=inf
  EXTRA_BIASconsider also these other bias potentials for the internal reweighting=Ftg.bias
 TEMP temperature=-1 COMPRESSION_THRESHOLD merge kernels if closer than this threshold, in units of sigma=1 FILE a file in which the list of all deposited kernels is stored=KERNELS
...
PRINTPrint quantities to a file. More details FMT the format that should be used to output real numbers=%g STRIDE the frequency with which the quantities of interest should be output=500 FILEthe name of the file on which to output these quantities=COLVAR ARGthe labels of the values that you would like to print to the file=phi,Ftg.bias,opes.bias

Notice that in order to reweight for the unbiased $P(\mathbf{s})$ during postprocessing, the total bias Ftg.bias+opes.bias must be used.

Finally, an example of how to use the EXCLUDED_REGION keyword. It expects a characteristic function that is different from zero in the region to be excluded. You can use CUSTOM and a combination of the step function to define it. With the following input no kernel is deposited in the transition state region of alanine dipeptide, defined by the interval $\phi \in [-0.6, 0.7]$ :

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

phiThe TORSION action with label phi calculates the following quantities: Quantity    Type    Description  
phi scalar the TORSION involving these atoms: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=5,7,9,15
psiThe TORSION action with label psi calculates the following quantities: Quantity    Type    Description  
psi scalar the TORSION involving these atoms: TORSIONCalculate one or multiple torsional angles. More details ATOMSthe four atoms involved in the torsional angle=7,9,15,17
xxThe CUSTOM action with label xx calculates the following quantities: Quantity    Type    Description  
xx scalar an arbitrary function: CUSTOMCalculate a combination of variables using a custom expression. More details PERIODICif the output of your function is periodic then you should specify the periodicity of the function=NO ARGthe values input to this function=phi FUNCthe function you wish to evaluate=step(x+0.6)-step(x-0.7)
opesThe OPES_METAD action with label opes calculates the following quantities: Quantity    Type    Description  
opes.bias scalar the instantaneous value of the bias potential
opes.rct scalar estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
opes.zed scalar estimate of Z_n. should become flat once no new CV-space region is explored
opes.neff scalar effective sample size
opes.nker scalar total number of compressed kernels used to represent the bias
opes.nlker scalar number of kernels in the neighbor list
opes.nlsteps scalar number of steps from last neighbor list update: OPES_METADOn-the-fly probability enhanced sampling with metadynamics-like target distribution. This action has hidden defaults. More details ...
  ARGthe labels of the scalars on which the bias will act=phi,psi
  PACEthe frequency for kernel deposition=500
  BARRIERthe free energy barrier to be overcome=30
  EXCLUDED_REGIONkernels are not deposited when the action provided here has a nonzero value, see example above=xx
  NLIST use neighbor list for kernels summation, faster but experimental
...
opes: OPES_METADOn-the-fly probability enhanced sampling with metadynamics-like target distribution. This action uses the defaults shown here. More details ...
  ARGthe labels of the scalars on which the bias will act=phi,psi
  PACEthe frequency for kernel deposition=500
  BARRIERthe free energy barrier to be overcome=30
  EXCLUDED_REGIONkernels are not deposited when the action provided here has a nonzero value, see example above=xx
  NLIST use neighbor list for kernels summation, faster but experimental
 TEMP temperature=-1 SIGMA the initial widths of the kernels=ADAPTIVE COMPRESSION_THRESHOLD merge kernels if closer than this threshold, in units of sigma=1 FILE a file in which the list of all deposited kernels is stored=KERNELS
...
PRINTPrint quantities to a file. More details FMT the format that should be used to output real numbers=%g STRIDE the frequency with which the quantities of interest should be output=500 FILEthe name of the file on which to output these quantities=COLVAR ARGthe labels of the values that you would like to print to the file=phi,psi,xx,opes.*

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
EXCLUDED_REGION	scalar	kernels are not deposited when the action provided here has a nonzero value, see example above
EXTRA_BIAS	scalar	consider also these other bias potentials for the internal reweighting

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
rct	scalar	default	estimate of c(t)
zed	scalar	default	estimate of Z_n
neff	scalar	default	effective sample size
nker	scalar	default	total number of compressed kernels used to represent the bias
work	scalar	CALC_WORK	total accumulated work done by the bias
nlker	scalar	NLIST	number of kernels in the neighbor list
nlsteps	scalar	NLIST	number of steps from last neighbor list update

Full list of keywords

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
EXCLUDED_REGION	input	none	kernels are not deposited when the action provided here has a nonzero value, see example above
EXTRA_BIAS	input	none	consider also these other bias potentials for the internal reweighting
TEMP	compulsory	-1	temperature
PACE	compulsory	none	the frequency for kernel deposition
SIGMA	compulsory	ADAPTIVE	the initial widths of the kernels
BARRIER	compulsory	none	the free energy barrier to be overcome
COMPRESSION_THRESHOLD	compulsory	1	merge kernels if closer than this threshold, in units of sigma
FILE	compulsory	KERNELS	a file in which the list of all deposited kernels is stored
NUMERICAL_DERIVATIVESThis keyword do not have examples	optional	false	calculate the derivatives for these quantities numerically
ADAPTIVE_SIGMA_STRIDEThis keyword do not have examples	optional	not used	number of steps for measuring adaptive sigma
SIGMA_MIN	optional	not used	never reduce SIGMA below this value
BIASFACTOR	optional	not used	the gamma bias factor used for the well-tempered target distribution
EPSILONThis keyword do not have examples	optional	not used	the value of the regularization constant for the probability
KERNEL_CUTOFFThis keyword do not have examples	optional	not used	truncate kernels at this distance, in units of sigma
NLIST_PARAMETERSThis keyword do not have examples	optional	not used	the two cutoff parameters for the kernels neighbor list
NLIST	optional	false	use neighbor list for kernels summation, faster but experimental
NLIST_PACE_RESETThis keyword do not have examples	optional	false	force the reset of the neighbor list at each PACE
FIXED_SIGMAThis keyword do not have examples	optional	false	do not decrease sigma as the simulation proceeds
RECURSIVE_MERGE_OFFThis keyword do not have examples	optional	false	do not recursively attempt kernel merging when a new one is added
NO_ZEDThis keyword do not have examples	optional	false	do not normalize over the explored CV space, Z_n=1
FMTThis keyword do not have examples	optional	not used	specify format for KERNELS file
STATE_RFILE	optional	not used	read from this file the compressed kernels and all the info needed to RESTART the simulation
STATE_WFILE	optional	not used	write to this file the compressed kernels and all the info needed to RESTART the simulation
STATE_WSTRIDE	optional	not used	number of MD steps between writing the STATE_WFILE
STORE_STATES	optional	false	append to STATE_WFILE instead of ovewriting it each time
CALC_WORKThis keyword do not have examples	optional	false	calculate the total accumulated work done by the bias since last restart
WALKERS_MPI	optional	false	switch on MPI version of multiple walkers
SERIALThis keyword do not have examples	optional	false	perform calculations in serial
RESTARTThis keyword do not have examples	optional	not used	allows per-action setting of restart (YES/NO/AUTO)
UPDATE_FROMThis keyword do not have examples	optional	not used	Only update this action from this time
UPDATE_UNTILThis keyword do not have examples	optional	not used	Only update this action until this time

References

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

Quantity	Type	Description
cv	scalar	the DISTANCE between this pair of atoms