OPES_METAD
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 (OPES) with metadynamics-like target distribution [64].

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{\tilde{P}_n(\mathbf{s})}{Z_n}+\epsilon\right)\, . \]

See Ref.[64] 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 Tutorial: Running and post-processing). 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 [97]). On the contrary, if your CVs are not degenerate but only suboptimal, you should converge faster by using OPES_METAD instead of METAD [64].

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.

Restart can be done from a KERNELS file, but it might not be perfect (due to limited precision when printing kernels to file, or usage of adaptive SIGMA). 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 following OPES_METAD Tutorial: Running and post-processing 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
tested on master
cv: DISTANCE 
ATOMS
the pair of atom that we are calculating the distance between.
=1,2 opes: OPES_METAD
ARG
the input for this action is the scalar output from one or more other actions.
=cv
PACE
compulsory keyword the frequency for kernel deposition
=200
BARRIER
compulsory keyword the free energy barrier to be overcome.
=40 PRINT
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=200
FILE
the name of the file on which to output these quantities
=COLVAR
ARG
the input for this action is the scalar output from one or more other actions.
=*

Another more articulated one:

Click on the labels of the actions for more information on what each action computes
tested on master
phi: TORSION 
ATOMS
the four atoms involved in the torsional angle
=5,7,9,15 psi: TORSION
ATOMS
the four atoms involved in the torsional angle
=7,9,15,17 opes: OPES_METAD ...
FILE
compulsory keyword ( default=KERNELS ) a file in which the list of all deposited kernels is stored
=Kernels.data
TEMP
compulsory keyword ( default=-1 ) temperature.
=300
ARG
the input for this action is the scalar output from one or more other actions.
=phi,psi
SIGMA
compulsory keyword ( default=ADAPTIVE ) the initial widths of the kernels.
=0.15,0.15
PACE
compulsory keyword the frequency for kernel deposition
=500
BARRIER
compulsory keyword the free energy barrier to be overcome.
=50
BIASFACTOR
the \f$\gamma\f$ bias factor used for the well-tempered target \f$p(\mathbfs)\f$.
=inf
STATE_RFILE
read from this file the compressed kernels and all the info needed to RESTART the simulation
=Restart.data
STATE_WFILE
write to this file the compressed kernels and all the info needed to RESTART the simulation
=State.data
STATE_WSTRIDE
number 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
... PRINT
FMT
the format that should be used to output real numbers
=%g
STRIDE
compulsory keyword ( default=1 ) the frequency with which the quantities of interest should be output
=500
FILE
the name of the file on which to output these quantities
=Colvar.data
ARG
the input for this action is the scalar output from one or more other actions.
=phi,psi,opes.*
Glossary 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)\): \(\frac{1}{\beta}\log \langle e^{\beta V} \rangle\), should become flat as the simulation converges. Do NOT use for reweighting
zed estimate of \(Z_n\), should become flat as 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
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, adaptive sigma will be used and the ADAPTIVE_SIGMA_STRIDE should be set
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 simulation goes on. 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

ARG the input for this action is the scalar output from one or more other actions. The particular scalars that you will use are referenced using the label of the action. If the label appears on its own then it is assumed that the Action calculates a single scalar value. The value of this scalar is thus used as the input to this new action. If * or *.* appears the scalars calculated by all the proceeding actions in the input file are taken. Some actions have multi-component outputs and each component of the output has a specific label. For example a DISTANCE action labelled dist may have three components x, y and z. To take just the x component you should use dist.x, if you wish to take all three components then use dist.*.More information on the referencing of Actions can be found in the section of the manual on the PLUMED Getting Started. Scalar values can also be referenced using POSIX regular expressions as detailed in the section on Regular Expressions. To use this feature you you must compile PLUMED with the appropriate flag. You can use multiple instances of this keyword i.e. ARG1, ARG2, ARG3...
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 \(p(\mathbf{s})\). 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)
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