	This is part of the opes module
	It is only available if you configure PLUMED with ./configure –enable-modules=opes . Furthermore, this feature is still being developed so take care when using it and report any problems on the mailing list.

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

This On-the-fly probability enhanced sampling (OPES) method with metadynamics-like target distribution is described in [58].

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 Ref.[58] 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 (see opes-metad). 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 [90]). On the contrary, if your CVs are not degenerate but only suboptimal, you should converge faster by using OPES_METAD instead of METAD [58].

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 \(P(\mathbf{s})\) 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 opes-metad 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

cv: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=1,2 The DISTANCE action with label cv calculates a single scalar value
opes: OPES_METAD ARGcompulsory keyword 
the labels of the scalars on which the bias will act 
=cv PACEcompulsory keyword 
the frequency for kernel deposition 
=200 BARRIERcompulsory keyword 
the free energy barrier to be overcome. 
=40 The OPES_METAD action with label opes calculates the following quantities:

 Quantity    Description  
opes.bias the instantaneous value of the bias potential
opes.rct estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
opes.zed estimate of Z_n. should become flat once no new CV-space region is explored
opes.neff effective sample size
opes.nker total number of compressed kernels used to represent the bias


PRINT STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output 
=200 FILEthe name of the file on which to output these quantities 
=COLVAR ARGcompulsory keyword 
the labels of the values that you would like to print to the file 
=* The PRINT action with label

Another more articulated one:

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

phi: TORSION ATOMSthe four atoms involved in the torsional angle 
=5,7,9,15 The TORSION action with label phi calculates a single scalar value
psi: TORSION ATOMSthe four atoms involved in the torsional angle 
=7,9,15,17 The TORSION action with label psi calculates a single scalar value
opes: OPES_METAD ...
   FILEcompulsory keyword ( default=KERNELS )
a file in which the list of all deposited kernels is stored =Kernels.data 
   TEMPcompulsory keyword ( default=-1 )
temperature. =300 
   ARGcompulsory keyword 
the labels of the scalars on which the bias will act 
=phi,psi 
   PACEcompulsory keyword 
the frequency for kernel deposition =500 
   BARRIERcompulsory keyword 
the free energy barrier to be overcome. =50 
   SIGMAcompulsory keyword ( default=ADAPTIVE )
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( default=off ) append to STATE_WFILE instead of ovewriting it each time  
   WALKERS_MPI( default=off ) switch on MPI version of multiple walkers  
   NLIST( default=off ) use neighbor list for kernels summation, faster but experimental
 
...The OPES_METAD action with label opes calculates the following quantities:

 Quantity    Description  
opes.bias the instantaneous value of the bias potential
opes.rct estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
opes.zed estimate of Z_n. should become flat once no new CV-space region is explored
opes.neff effective sample size
opes.nker total number of compressed kernels used to represent the bias
opes.nlker number of kernels in the neighbor list
opes.nlsteps number of steps from last neighbor list update


PRINT FMTthe format that should be used to output real numbers 
=%g STRIDEcompulsory keyword ( default=1 )
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 ARGcompulsory keyword 
the labels of the values that you would like to print to the file 
=phi,psi,opes.* The PRINT action with label

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

phi: TORSION ATOMSthe four atoms involved in the torsional angle 
=5,7,9,15 The TORSION action with label phi calculates a single scalar value
FtgValue: CUSTOM ARGcompulsory keyword 
the values input to this function 
=phi PERIODICcompulsory keyword 
if the output of your function is periodic then you should specify the periodicity
of the function. 
=NO FUNCcompulsory keyword 
the function you wish to evaluate 
=(x/0.4)^2 The CUSTOM action with label FtgValue calculates a single scalar value
Ftg: BIASVALUE ARGcompulsory keyword 
the labels of the scalar/vector arguments whose values will be used as a bias on
the system 
=FtgValue The BIASVALUE action with label Ftg calculates the following quantities:

 Quantity    Description  
Ftg.bias the instantaneous value of the bias potential
Ftg.FtgValue_bias 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


opes: OPES_METAD ...
   ARGcompulsory keyword 
the labels of the scalars on which the bias will act 
=phi 
   PACEcompulsory keyword 
the frequency for kernel deposition =500 
   BARRIERcompulsory keyword 
the free energy barrier to be overcome. =50 
   SIGMAcompulsory keyword ( default=ADAPTIVE )
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 
...The OPES_METAD action with label opes calculates the following quantities:

 Quantity    Description  
opes.bias the instantaneous value of the bias potential
opes.rct estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
opes.zed estimate of Z_n. should become flat once no new CV-space region is explored
opes.neff effective sample size
opes.nker total number of compressed kernels used to represent the bias


PRINT FMTthe format that should be used to output real numbers 
=%g STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output 
=500 FILEthe name of the file on which to output these quantities 
=COLVAR ARGcompulsory keyword 
the labels of the values that you would like to print to the file 
=phi,Ftg.bias,opes.bias The PRINT action with label

phi: TORSION ATOMSthe four atoms involved in the torsional angle 
=5,7,9,15 The TORSION action with label phi calculates a single scalar value
Ftg: BIASVALUE ARGcompulsory keyword 
the labels of the scalar/vector arguments whose values will be used as a bias on
the system 
=FtgValue The BIASVALUE action with label FtgNotice 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 




phi: TORSION ATOMSthe four atoms involved in the torsional angle 
=5,7,9,15 The TORSION action with label phi calculates a single scalar value
psi: TORSION ATOMSthe four atoms involved in the torsional angle 
=7,9,15,17 The TORSION action with label psi calculates a single scalar value
xx: CUSTOM PERIODICcompulsory keyword 
if the output of your function is periodic then you should specify the periodicity
of the function. 
=NO ARGcompulsory keyword 
the values input to this function 
=phi FUNCcompulsory keyword 
the function you wish to evaluate 
=step(x+0.6)-step(x-0.7) The CUSTOM action with label xx calculates a single scalar value
opes: OPES_METAD ...
   ARGcompulsory keyword 
the labels of the scalars on which the bias will act 
=phi,psi 
   PACEcompulsory keyword 
the frequency for kernel deposition =500 
   BARRIERcompulsory keyword 
the 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( default=off ) use neighbor list for kernels summation, faster but experimental
 
...The OPES_METAD action with label opes calculates the following quantities:

 Quantity    Description  
opes.bias the instantaneous value of the bias potential
opes.rct estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
opes.zed estimate of Z_n. should become flat once no new CV-space region is explored
opes.neff effective sample size
opes.nker total number of compressed kernels used to represent the bias
opes.nlker number of kernels in the neighbor list
opes.nlsteps number of steps from last neighbor list update


PRINT FMTthe format that should be used to output real numbers 
=%g STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output 
=500 FILEthe name of the file on which to output these quantities 
=COLVAR ARGcompulsory keyword 
the labels of the values that you would like to print to the file 
=phi,psi,xx,opes.* The PRINT action with label 


phi: TORSION ATOMSthe four atoms involved in the torsional angle 
=5,7,9,15 The TORSION action with label phi calculates a single scalar value
psi: TORSION ATOMSthe four atoms involved in the torsional angle 
=7,9,15,17 The TORSION action with label psi calculates a single scalar value
opes: OPES_METAD ...
   ARGcompulsory keyword 
the labels of the scalars on which the bias will act 
=phi,psi 
   PACEcompulsory keyword 
the frequency for kernel deposition =500 
   BARRIERcompulsory keyword 
the 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( default=off ) use neighbor list for kernels summation, faster but experimental
 
...The OPES_METAD action with label opesGlossary of keywords and components

Description of components

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   

 rct  estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting   

 zed  estimate of Z_n. should become flat once no new CV-space region is explored   

 neff  effective sample size   

 nker  total number of compressed kernels used to represent the bias   

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


 Quantity    Keyword    Description    

 work   CALC_WORK   total accumulated work done by the bias   

 nlker   NLIST   number of kernels in the neighbor list   

 nlsteps   NLIST   number of steps from last neighbor list update   

Compulsory keywords



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

 TEMP  ( default=-1 ) temperature. If not set, it is taken from MD engine, but not all MD codes provide it   

 PACE  the frequency for kernel deposition   

 SIGMA  ( default=ADAPTIVE ) the initial widths of the kernels. If not set, an adaptive sigma will be used with the given ADAPTIVE_SIGMA_STRIDE   

 BARRIER  the free energy barrier to be overcome. It is used to set BIASFACTOR, EPSILON, and KERNEL_CUTOFF to reasonable values   

 COMPRESSION_THRESHOLD  ( default=1 ) merge kernels if closer than this threshold, in units of sigma   

 FILE  ( default=KERNELS ) a file in which the list of all deposited kernels is stored   

Options



 NUMERICAL_DERIVATIVES  ( default=off ) calculate the derivatives for these quantities numerically   

 NLIST  ( default=off ) use neighbor list for kernels summation, faster but experimental   

 NLIST_PACE_RESET  ( default=off ) force the reset of the neighbor list at each PACE. Can be useful with WALKERS_MPI   

 FIXED_SIGMA  ( default=off ) do not decrease sigma as the simulation proceeds. Can be added in a RESTART, to keep in check the number of compressed kernels   

 RECURSIVE_MERGE_OFF  ( default=off ) do not recursively attempt kernel merging when a new one is added   

 NO_ZED  ( default=off ) do not normalize over the explored CV space, Z_n=1   

 STORE_STATES  ( default=off ) append to STATE_WFILE instead of ovewriting it each time   

 CALC_WORK  ( default=off ) calculate the total accumulated work done by the bias since last restart   

 WALKERS_MPI  ( default=off ) switch on MPI version of multiple walkers   

 SERIAL  ( default=off ) perform calculations in serial  



 ADAPTIVE_SIGMA_STRIDE  number of steps for measuring adaptive sigma. Default is 10xPACE   

 SIGMA_MIN  never reduce SIGMA below this value   

 BIASFACTOR  the gamma bias factor used for the well-tempered target distribution. Set to 'inf' for uniform flat target   

 EPSILON  the value of the regularization constant for the probability   

 KERNEL_CUTOFF  truncate kernels at this distance, in units of sigma   

 NLIST_PARAMETERS  ( default=3.0,0.5 ) the two cutoff parameters for the kernels neighbor list   

 FMT  specify format for KERNELS file   

 STATE_RFILE  read from this file the compressed kernels and all the info needed to RESTART the simulation   

 STATE_WFILE  write to this file the compressed kernels and all the info needed to RESTART the simulation   

 STATE_WSTRIDE  number of MD steps between writing the STATE_WFILE. Default is only on CPT events (but not all MD codes set them)   

 EXCLUDED_REGION  kernels are not deposited when the action provided here has a nonzero value, see example above   

 EXTRA_BIAS  consider also these other bias potentials for the internal reweighting. This can be used e.g. for sampling a custom target distribution (see example above)   

 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

Quantity	Description
bias	the instantaneous value of the bias potential
rct	estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
zed	estimate of Z_n. should become flat once no new CV-space region is explored
neff	effective sample size
nker	total number of compressed kernels used to represent the bias

Quantity	Keyword	Description
work	CALC_WORK	total accumulated work done by the bias
nlker	NLIST	number of kernels in the neighbor list
nlsteps	NLIST	number of steps from last neighbor list update

Quantity	Description
opes.bias	the instantaneous value of the bias potential
opes.rct	estimate of c(t). log(exp(beta V)/beta, should become flat as the simulation converges. Do NOT use for reweighting
opes.zed	estimate of Z_n. should become flat once no new CV-space region is explored
opes.neff	effective sample size
opes.nker	total number of compressed kernels used to represent the bias

Quantity	Description
Ftg.bias	the instantaneous value of the bias potential
Ftg.FtgValue_bias	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

ARG	the labels of the scalars on which the bias will act
TEMP	( default=-1 ) temperature. If not set, it is taken from MD engine, but not all MD codes provide it
PACE	the frequency for kernel deposition
SIGMA	( default=ADAPTIVE ) the initial widths of the kernels. If not set, an adaptive sigma will be used with the given ADAPTIVE_SIGMA_STRIDE
BARRIER	the free energy barrier to be overcome. It is used to set BIASFACTOR, EPSILON, and KERNEL_CUTOFF to reasonable values
COMPRESSION_THRESHOLD	( default=1 ) merge kernels if closer than this threshold, in units of sigma
FILE	( default=KERNELS ) a file in which the list of all deposited kernels is stored

NUMERICAL_DERIVATIVES	( default=off ) calculate the derivatives for these quantities numerically
NLIST	( default=off ) use neighbor list for kernels summation, faster but experimental
NLIST_PACE_RESET	( default=off ) force the reset of the neighbor list at each PACE. Can be useful with WALKERS_MPI
FIXED_SIGMA	( default=off ) do not decrease sigma as the simulation proceeds. Can be added in a RESTART, to keep in check the number of compressed kernels
RECURSIVE_MERGE_OFF	( default=off ) do not recursively attempt kernel merging when a new one is added
NO_ZED	( default=off ) do not normalize over the explored CV space, Z_n=1
STORE_STATES	( default=off ) append to STATE_WFILE instead of ovewriting it each time
CALC_WORK	( default=off ) calculate the total accumulated work done by the bias since last restart
WALKERS_MPI	( default=off ) switch on MPI version of multiple walkers
SERIAL	( default=off ) perform calculations in serial
ADAPTIVE_SIGMA_STRIDE	number of steps for measuring adaptive sigma. Default is 10xPACE
SIGMA_MIN	never reduce SIGMA below this value
BIASFACTOR	the gamma bias factor used for the well-tempered target distribution. Set to 'inf' for uniform flat target
EPSILON	the value of the regularization constant for the probability
KERNEL_CUTOFF	truncate kernels at this distance, in units of sigma
NLIST_PARAMETERS	( default=3.0,0.5 ) the two cutoff parameters for the kernels neighbor list
FMT	specify format for KERNELS file
STATE_RFILE	read from this file the compressed kernels and all the info needed to RESTART the simulation
STATE_WFILE	write to this file the compressed kernels and all the info needed to RESTART the simulation
STATE_WSTRIDE	number of MD steps between writing the STATE_WFILE. Default is only on CPT events (but not all MD codes set them)
EXCLUDED_REGION	kernels are not deposited when the action provided here has a nonzero value, see example above
EXTRA_BIAS	consider also these other bias potentials for the internal reweighting. This can be used e.g. for sampling a custom target distribution (see example above)
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