PLUMED: EDS

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

Add a linear bias on a set of observables.

This force is the same as the linear part of the bias in RESTRAINT, but this bias has the ability to compute the prefactors adaptively using the scheme of White and Voth [144] in order to match target observable values for a set of CVs. Further updates to the algorithm are described in [66] and you can read a review on the method and its applications here: [2].

You can see a tutorial on EDS specifically for biasing coordination number at Andrew White's webpage.

The addition to the potential is of the form

\[ \sum_i \frac{\alpha_i}{s_i} x_i \]

where for CV \(x_i\), a coupling constant \({\alpha}_i\) is determined adaptively or set by the user to match a target value for \(x_i\). \(s_i\) is a scale parameter, which by default is set to the target value. It may also be set separately.

Warning: It is not possible to set the target value of the observable to zero with the default value of \(s_i\) as this will cause a divide-by-zero error. Instead, set \(s_i=1\) or modify the CV so the desired target value is no longer zero.

Notice that a similar method is available as MAXENT, although with different features and using a different optimization algorithm.

Virial

The bias forces modify the virial and this can change your simulation density if the bias is used in an NPT simulation. One way to avoid changing the virial contribution from the bias is to add the keyword VIRIAL=1.0, which changes how the bias is computed to minimize its contribution to the virial. This can also lead to smaller magnitude biases that are more robust if transferred to other systems. VIRIAL=1.0 can be a reasonable starting point and increasing the value changes the balance between matching the set-points and minimizing the virial. See [2] for details on the equations. Since the coupling constants are unique with a single CV, VIRIAL is not applicable with a single CV. When used with multiple CVs, the CVs should be correlated which is almost always the case.

Weighting

EDS computes means and variances as part of its algorithm. If you are also using a biasing method like metadynamics, you may wish to remove the effect of this bias in your EDS computations so that EDS works on the canonical values (reweighted to be unbiased). For example, you may be using metadynamics to bias a dihedral angle to enhance sampling and be using EDS to set the average distance between two particular atoms. Specifically:

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

# set-up metadynamics
t: TORSION ATOMSthe four atoms involved in the torsional angle
=1,2,3,4 You cannot view the components that are calculated by each action for this input file. Sorry
md: METAD ARGthe input for this action is the scalar output from one or more other actions.
=d SIGMAcompulsory keyword
the widths of the Gaussian hills
=0.2 HEIGHTthe heights of the Gaussian hills.
=0.3 PACEcompulsory keyword
the frequency for hill addition
=500 TEMPthe system temperature - this is only needed if you are doing well-tempered metadynamics
=300 You cannot view the components that are calculated by each action for this input file. Sorry
# compute bias weights
bias: REWEIGHT_METAD TEMPthe system temperature.
=300 You cannot view the components that are calculated by each action for this input file. Sorry
# now do EDS on distance while removing effect of metadynamics
d: DISTANCE ATOMSthe pair of atom that we are calculating the distance between.
=4,7 You cannot view the components that are calculated by each action for this input file. Sorry
eds: EDS ARGthe input for this action is the scalar output from one or more other actions.
=d CENTERThe desired centers (equilibrium values) which will be sought during the adaptive
linear biasing.
=3.0 PERIODSteps over which to adjust bias for adaptive or ramping
=100 TEMPThe system temperature.
=300 LOGWEIGHTSAdd weights to use for computing statistics.
=bias You cannot view the components that are calculated by each action for this input file. Sorry

This is an approximation though because EDS uses a finite samples while running to get means/variances. At the end of a run, you should ensure this approach worked and indeed your reweighted CV matches the target value.

Examples

The following input for a harmonic oscillator of two beads will adaptively find a linear bias to change the mean and variance to the target values. The PRINT line shows how to access the value of the coupling constants.

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

dist: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=1,2 The DISTANCE action with label dist calculates a single scalar value
# this is the squared of the distance
dist2: COMBINE ARGthe input for this action is the scalar output from one or more other actions. 
=dist POWERScompulsory keyword ( default=1.0 )
the powers to which you are raising each of the arguments in your function 
=2 PERIODICcompulsory keyword 
if the output of your function is periodic then you should specify the periodicity
of the function. 
=NO The COMBINE action with label dist2 calculates a single scalar value
# bias mean and variance
eds: EDS ARGthe input for this action is the scalar output from one or more other actions. 
=dist,dist2 CENTERThe desired centers (equilibrium values) which will be sought during the adaptive
linear biasing. 
=2.0,1.0 PERIODSteps over which to adjust bias for adaptive or ramping 
=100 TEMPThe system temperature. 
=1.0 The EDS action with label eds calculates the following quantities:

 Quantity    Description  
eds.bias the instantaneous value of the bias potential
eds.force2 squared value of force from the bias
eds.dist_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named dist
eds.dist2_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named dist2


PRINT ARGthe input for this action is the scalar output from one or more other actions. 
=dist,dist2,eds.dist_coupling,eds.dist2_coupling,eds.bias,eds.force2 FILEthe name of the file on which to output these quantities 
=colvars.dat STRIDEcompulsory keyword ( default=1 )
the frequency with which the quantities of interest should be output 
=100 The PRINT action with label

Rather than trying to find the coupling constants adaptively, one can ramp up to a constant value.

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

dist: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=1,2 The DISTANCE action with label dist calculates a single scalar value
dist2: COMBINE ARGthe input for this action is the scalar output from one or more other actions. 
=dist POWERScompulsory keyword ( default=1.0 )
the powers to which you are raising each of the arguments in your function 
=2 PERIODICcompulsory keyword 
if the output of your function is periodic then you should specify the periodicity
of the function. 
=NO The COMBINE action with label dist2 calculates a single scalar value
# ramp couplings from 0,0 to -1,1 over 50000 steps
eds: EDS ARGthe input for this action is the scalar output from one or more other actions. 
=dist,dist2 CENTERThe desired centers (equilibrium values) which will be sought during the adaptive
linear biasing. 
=2.0,1.0 FIXEDcompulsory keyword ( default=0 )
Fixed target values for coupling constant. 
=-1,1 RAMP( default=off ) Slowly increase bias constant to a fixed value 
 PERIODSteps over which to adjust bias for adaptive or ramping 
=50000 TEMPThe system temperature. 
=1.0 The EDS action with label eds calculates the following quantities:

 Quantity    Description  
eds.bias the instantaneous value of the bias potential
eds.force2 squared value of force from the bias
eds.dist_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named dist
eds.dist2_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named dist2


# same as above, except starting at -0.5,0.5 rather than default of 0,0
eds2: EDS ARGthe input for this action is the scalar output from one or more other actions. 
=dist,dist2 CENTERThe desired centers (equilibrium values) which will be sought during the adaptive
linear biasing. 
=2.0,1.0 FIXEDcompulsory keyword ( default=0 )
Fixed target values for coupling constant. 
=-1,1 INITcompulsory keyword ( default=0 )
Starting value for coupling constant 
=-0.5,0.5 RAMP( default=off ) Slowly increase bias constant to a fixed value 
 PERIODSteps over which to adjust bias for adaptive or ramping 
=50000 TEMPThe system temperature. 
=1.0 The EDS action with label eds2 calculates the following quantities:

 Quantity    Description  
eds2.bias the instantaneous value of the bias potential
eds2.force2 squared value of force from the bias
eds2.dist_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named dist
eds2.dist2_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named dist2

A restart file can be added to dump information needed to restart/continue simulation using these parameters every PERIOD.

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

dist: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=1,2 The DISTANCE action with label dist calculates a single scalar value
dist2: COMBINE ARGthe input for this action is the scalar output from one or more other actions. 
=dist POWERScompulsory keyword ( default=1.0 )
the powers to which you are raising each of the arguments in your function 
=2 PERIODICcompulsory keyword 
if the output of your function is periodic then you should specify the periodicity
of the function. 
=NO The COMBINE action with label dist2 calculates a single scalar value
# add the option to write to a restart file
eds: EDS ARGthe input for this action is the scalar output from one or more other actions. 
=dist,dist2 CENTERThe desired centers (equilibrium values) which will be sought during the adaptive
linear biasing. 
=2.0,1.0 PERIODSteps over which to adjust bias for adaptive or ramping 
=100 TEMPThe system temperature. 
=1.0 OUT_RESTARTOutput file for all information needed to continue EDS simulation. 
=checkpoint.eds The EDS action with label eds calculates the following quantities:

 Quantity    Description  
eds.bias the instantaneous value of the bias potential
eds.force2 squared value of force from the bias
eds.dist_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named dist
eds.dist2_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named dist2

The first few lines of the restart file that is output if we run a calculation with one CV will look something like this:

#! FIELDS time d1_center d1_set d1_target d1_coupling d1_maxrange d1_maxgrad d1_accum d1_mean d1_std
#! SET adaptive  1
#! SET update_period  1
#! SET seed  0
#! SET kbt    2.4943
   0.0000   1.0000   0.0000   0.0000   0.0000   7.4830   0.1497   0.0000   0.0000   0.0000
   1.0000   1.0000   0.0000   0.0000   0.0000   7.4830   0.1497   0.0000   0.0000   0.0000
   2.0000   1.0000  -7.4830   0.0000   0.0000   7.4830   0.1497   0.0224   0.0000   0.0000
   3.0000   1.0000  -7.4830   0.0000  -7.4830   7.4830   0.1497   0.0224   0.0000   0.0000
   4.0000   1.0000  -7.4830   0.0000  -7.4830   7.4830   0.1497   0.0224   0.0000   0.0000

Read in a previous restart file. Adding RESTART flag makes output append

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

d1: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=1,2 The DISTANCE action with label d1 calculates a single scalar value
eds: EDS ARGthe input for this action is the scalar output from one or more other actions. 
=d1 CENTERThe desired centers (equilibrium values) which will be sought during the adaptive
linear biasing. 
=2.0 PERIODSteps over which to adjust bias for adaptive or ramping 
=100 TEMPThe system temperature. 
=1.0 IN_RESTARTRead this file to continue an EDS simulation. 
=restart.eds RESTARTallows per-action setting of restart (YES/NO/AUTO) 
=YES The EDS action with label eds calculates the following quantities:

 Quantity    Description  
eds.bias the instantaneous value of the bias potential
eds.force2 squared value of force from the bias
eds.d1_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named d1

Read in a previous restart file and freeze the bias at the final level from the previous simulation

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

d1: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=1,2 The DISTANCE action with label d1 calculates a single scalar value
eds: EDS ARGthe input for this action is the scalar output from one or more other actions. 
=d1 CENTERThe desired centers (equilibrium values) which will be sought during the adaptive
linear biasing. 
=2.0 TEMPThe system temperature. 
=1.0 IN_RESTARTRead this file to continue an EDS simulation. 
=restart.eds FREEZE( default=off ) Fix bias at current level (only used for restarting). 
 The EDS action with label eds calculates the following quantities:

 Quantity    Description  
eds.bias the instantaneous value of the bias potential
eds.force2 squared value of force from the bias
eds.d1_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named d1

Read in a previous restart file and freeze the bias at the mean from the previous simulation

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

d1: DISTANCE ATOMSthe pair of atom that we are calculating the distance between. 
=1,2 The DISTANCE action with label d1 calculates a single scalar value
eds: EDS ARGthe input for this action is the scalar output from one or more other actions. 
=d1 CENTERThe desired centers (equilibrium values) which will be sought during the adaptive
linear biasing. 
=2.0 TEMPThe system temperature. 
=1.0 IN_RESTARTRead this file to continue an EDS simulation. 
=restart.eds FREEZE( default=off ) Fix bias at current level (only used for restarting). 
 MEAN( default=off ) Instead of using final bias level from restart, use average. 
 The EDS action with label eds calculates the following quantities:

 Quantity    Description  
eds.bias the instantaneous value of the bias potential
eds.force2 squared value of force from the bias
eds.d1_coupling For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor. This particular component measures this quantity for the input CV named d1

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
force2	squared value of force from the bias
pressure	If using virial keyword, this is the current sum of virial terms. It is in units of pressure (energy / vol^3)
_coupling	For each named CV biased, there will be a corresponding output CV_coupling storing the current linear bias prefactor.

Compulsory keywords

RANGE	( default=25.0 ) The (starting) maximum increase in coupling constant per PERIOD (in \(k_B T\)/[BIAS_SCALE unit]) for each CV biased
SEED	( default=0 ) Seed for random order of changing bias
INIT	( default=0 ) Starting value for coupling constant
FIXED	( default=0 ) Fixed target values for coupling constant. Non-adaptive.
LM_MIXING	( default=1 ) Initial mixing parameter when using Levenberg-Marquadt minimization.

Options

NUMERICAL_DERIVATIVES	( default=off ) calculate the derivatives for these quantities numerically
LM	( default=off ) Use Levenberg-Marquadt algorithm along with simultaneous keyword. Otherwise use gradient descent.
RAMP	( default=off ) Slowly increase bias constant to a fixed value
COVAR	( default=off ) Utilize the covariance matrix when updating the bias. Default Off, but may be enabled due to other options
FREEZE	( default=off ) Fix bias at current level (only used for restarting).
MEAN	( default=off ) Instead of using final bias level from restart, use average. Can only be used in conjunction with FREEZE
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...
CENTER	The desired centers (equilibrium values) which will be sought during the adaptive linear biasing. This is for fixed centers
CENTER_ARG	The desired centers (equilibrium values) which will be sought during the adaptive linear biasing. CENTER_ARG is for calculated centers, e.g. from a CV or analysis.
PERIOD	Steps over which to adjust bias for adaptive or ramping
BIAS_SCALE	A divisor to set the units of the bias. If not set, this will be the CENTER value by default (as is done in White and Voth 2014).
TEMP	The system temperature. If not provided will be taken from MD code (if available)
MULTI_PROP	What proportion of dimensions to update at each step. Must be in interval [1,0), where 1 indicates all and any other indicates a stochastic update. If not set, default is 1 / N, where N is the number of CVs.
VIRIAL	Add an update penalty for having non-zero virial contributions. Only makes sense with multiple correlated CVs.
LOGWEIGHTS	Add weights to use for computing statistics. For example, if biasing with metadynamics.
RESTART_FMT	the format that should be used to output real numbers in EDS restarts
OUT_RESTART	Output file for all information needed to continue EDS simulation. If you have the RESTART directive set (global or for EDS), this file will be appended to. Note that the header will be printed again if appending.
IN_RESTART	Read this file to continue an EDS simulation. If same as OUT_RESTART and you have not set the RESTART directive, the file will be backed-up and overwritten with new output. If you do have the RESTART flag set and it is the same name as OUT_RESTART, this file will be appended.
RESTART	allows per-action setting of restart (YES/NO/AUTO)

1.8.17