$\phantom{\rule{0.3em}{0ex}}$

COLLECTIVE VARIABLES MODULE

Reference manual for VMD

Code version: 2018-12-21

Alejandro Bernardin, Jeﬀrey R. Comer, Giacomo Fiorin, Haohao Fu, Jérôme Hénin, Axel Kohlmeyer, Fabrizio Marinelli, Joshua V. Vermaas, Andrew D. White

### 1 Overview

In molecular dynamics simulations, it is often useful to reduce the large number of degrees of freedom of a physical system into few parameters whose statistical distributions can be analyzed individually, or used to deﬁne biasing potentials to alter the dynamics of the system in a controlled manner. These have been called ‘order parameters’, ‘collective variables’, ‘(surrogate) reaction coordinates’, and many other terms.

Here we use primarily the term ‘collective variable’, often shortened to colvar, to indicate any diﬀerentiable function of atomic Cartesian coordinates, ${\text{}x\text{}}_{i}$, with $i$ between $1$ and $N$, the total number of atoms:

 $\xi \left(t\right)\phantom{\rule{3.04074pt}{0ex}}=\xi \left(\text{}X\text{}\left(t\right)\right)\phantom{\rule{3.04074pt}{0ex}}=\xi \left({\text{}x\text{}}_{i}\left(t\right),{\text{}x\text{}}_{j}\left(t\right),{\text{}x\text{}}_{k}\left(t\right),\dots \right)\phantom{\rule{3.04074pt}{0ex}},\phantom{\rule{3.04074pt}{0ex}}\phantom{\rule{3.04074pt}{0ex}}1\le i,j,k\dots \le N$ (1)

This manual documents the collective variables module (Colvars), a software that provides an implementation for the functions $\xi \left(\text{}X\text{}\right)$ with a focus on ﬂexibility, robustness and high performance. The module is designed to perform multiple tasks concurrently during or after a simulation, the most common of which are:

• apply restraints or biasing potentials to multiple variables, tailored on the system by choosing from a wide set of basis functions, without limitations on their number or on the number of atoms involved;
• calculate potentials of mean force (PMFs) along any set of variables, using diﬀerent enhanced sampling methods, such as Adaptive Biasing Force (ABF), metadynamics, steered MD and umbrella sampling; variants of these methods that make use of an ensemble of replicas are supported as well;
• calculate statistical properties of the variables, such as running averages and standard deviations, correlation functions of pairs of variables, and multidimensional histograms: this can be done either at run-time without the need to save very large trajectory ﬁles, or after a simulation has been completed using VMD and the cv command.

Note: although restraints and PMF algorithms are primarily used during simulations, they are also available in VMD to test a new input for a simulation, or to evaluate the relative free energy of a new structure based on data from a previous calculation. Options that only have an eﬀect during a simulation are also included for compatibility purposes.

Detailed explanations of the design of the Colvars module are provided in reference [1]. Please cite this reference whenever publishing work that makes use of this module.

### 2 A crash course

Suppose that we want to run a steered MD experiment where a small molecule is pulled away from a protein binding site. In Colvars terms, this is done by applying a moving restraint to the distance between the two objects. The conﬁguration will contain two blocks, one deﬁning the distance variable (see section 4 and 4.2.1), and the other the moving harmonic restraint (6.4).

colvar {
name dist
distance {
group1 { atomNumbersRange 42-55 }
group2 {
psfSegID PR
atomNameResidueRange CA 15-30
}
}
}

harmonic {
colvars dist
forceConstant 20.0
centers 4.0         # initial distance
targetCenters 15.0  # final distance
targetNumSteps 500000
}

Reading this input in plain English: the variable here named dist consists in a distance function between the centers of two groups: the ligand (atoms 42 to 55) and the $\alpha$-carbon atoms of residues 15 to 30 in the protein (segment name PR). To the “dist” variable, we apply a harmonic potential of force constant 20 kcal/mol/Å${}^{2}$, initially centered around a value of 4 Å, which will increase to 15 Å over 500,000 simulation steps.

The atom selection keywords is detailed in section 5.

### 3 General parameters

Here, we document the syntax of the commands and parameters used to set up and use the Colvars module in VMD. One of these parameters is the conﬁguration ﬁle or the conﬁguration text for the module itself, whose syntax is described in 3.2 and in the following sections.

#### 3.1 Using the cv command

At any moment during the execution of VMD, several options in the Colvars module can be read or modiﬁed by the command cv with the following syntax:
cv $<$subcommand$>$ [args...]
For example, to record the value of a collective variable named myVar into the Tcl variable value, use the following syntax:
set value [cv colvar myVar value]
All subcommands of cv are documented below.

Note: in VMD, Colvars must be attached to one molecule (system). Therefore, the cv command must be used for the ﬁrst time as cv molid $<$molid$>$ to set up the Colvars module for the molecule identiﬁed by $<$molid$>$. In all following invocations, the cv command will continue operating on the same molecule, regardless of its “top” status. To use the cv command on a diﬀerent molecule, use cv delete ﬁrst and then cv molid $<$molid$>$. Invoking the cv command with no arguments prints a help screen.

##### 3.1.1 Example use of the cv command: analyze a trajectory

By far the most typical use of Colvars in VMD is computing the values of one or more variables along an existing trajectory:
# Activate the module on the current VMD molecule
cv molid top
# Load a Colvars config file
cv configfile test.in
set out [open "test.colvars.traj" "w"]
# Write the labels to the file
puts -nonewline ${out} [cv printframelabels] for { set fr 0 } {${fr} < [molinfo top get numframes] } { incr fr } {
# Point Colvars to this trajectory frame
cv frame ${fr} # Recompute variables and biases cv update # Print variables and biases to the file puts -nonewline${out} [cv printframe]
}

##### 6.9.1 Grid deﬁnition for multidimensional histograms

Like the ABF and metadynamics biases, histogram uses the parameters lowerBoundary, upperBoundary, and width to deﬁne its grid. These values can be overridden if a conﬁguration block histogramGrid { } is provided inside the conﬁguration of histogram. The options supported inside this conﬁguration block are:

• lowerBoundaries $⟨\phantom{\rule{0.3em}{0ex}}$Lower boundaries of the grid$\phantom{\rule{0.3em}{0ex}}⟩$
Context: histogramGrid
Acceptable values: list of space-separated decimals
Description: This option deﬁnes the lower boundaries of the grid, overriding any values deﬁned by the lowerBoundary keyword of each colvar. Note that when gatherVectorColvars is on, each vector variable is automatically treated as a scalar, and a single value should be provided for it.
• upperBoundaries: analogous to lowerBoundaries
• widths: analogous to lowerBoundaries

#### 6.10 Probability distribution-restraints

The histogramRestraint bias implements a continuous potential of many variables (or of a single high-dimensional variable) aimed at reproducing a one-dimensional statistical distribution that is provided by the user. The $M$ variables $\left({\xi }_{1},\dots ,{\xi }_{M}\right)$ are interpreted as multiple observations of a random variable $\xi$ with unknown probability distribution. The potential is minimized when the histogram $h\left(\xi \right)$, estimated as a sum of Gaussian functions centered at $\left({\xi }_{1},\dots ,{\xi }_{M}\right)$, is equal to the reference histogram ${h}_{0}\left(\xi \right)$:

 $V\left({\xi }_{1},\dots ,{\xi }_{M}\right)=\frac{1}{2}k\int {\left(h\left(\xi \right)-{h}_{0}\left(\xi \right)\right)}^{2}\mathrm{d}\xi$ (28)
 $h\left(\xi \right)=\frac{1}{M\sqrt{2\pi {\sigma }^{2}}}\sum _{i=1}^{M}exp\left(-\frac{{\left(\xi -{\xi }_{i}\right)}^{2}}{2{\sigma }^{2}}\right)$ (29)

When used in combination with a distancePairs multi-dimensional variable, this bias implements the reﬁnement algorithm against ESR/DEER experiments published by Shen et al [25].

This bias behaves similarly to the histogram bias with the gatherVectorColvars option, with the important diﬀerence that all variables are gathered, resulting in a one-dimensional histogram. Future versions will include support for multi-dimensional histograms.

The list of options is as follows:

• name: see deﬁnition of name in sec. 6 (biasing and analysis methods)
• colvars: see deﬁnition of colvars in sec. 6 (biasing and analysis methods)
• outputEnergy: see deﬁnition of outputEnergy in sec. 6 (biasing and analysis methods)
• lowerBoundary $⟨\phantom{\rule{0.3em}{0ex}}$Lower boundary of the colvar grid$\phantom{\rule{0.3em}{0ex}}⟩$
Context: histogramRestraint
Acceptable values: decimal
Description: Deﬁnes the lowest end of the interval where the reference distribution ${h}_{0}\left(\xi \right)$ is deﬁned. Exactly one value must be provided, because only one-dimensional histograms are supported by the current version.
• upperBoundary: analogous to lowerBoundary
• width $⟨\phantom{\rule{0.3em}{0ex}}$Width of the colvar grid$\phantom{\rule{0.3em}{0ex}}⟩$
Context: histogramRestraint
Acceptable values: positive decimal
Description: Deﬁnes the spacing of the grid where the reference distribution ${h}_{0}\left(\xi \right)$ is deﬁned.
• gaussianSigma $⟨\phantom{\rule{0.3em}{0ex}}$Standard deviation of the approximating Gaussian$\phantom{\rule{0.3em}{0ex}}⟩$
Context: histogramRestraint
Acceptable values: positive decimal
Default value: 2 $×$ width
Description: Deﬁnes the parameter $\sigma$ in eq. 29.
• forceConstant $⟨\phantom{\rule{0.3em}{0ex}}$Force constant (kcal/mol)$\phantom{\rule{0.3em}{0ex}}⟩$
Context: histogramRestraint
Acceptable values: positive decimal
Default value: 1.0
Description: Deﬁnes the parameter $k$ in eq. 28.
• refHistogram $⟨\phantom{\rule{0.3em}{0ex}}$Reference histogram ${h}_{0}\left(\xi \right)$$\phantom{\rule{0.3em}{0ex}}⟩$
Context: histogramRestraint
Acceptable values: space-separated list of $M$ positive decimals
Description: Provides the values of ${h}_{0}\left(\xi \right)$ consecutively. The mid-point convention is used, i.e. the ﬁrst point that should be included is for $\xi$ = lowerBoundary+width/2. If the integral of ${h}_{0}\left(\xi \right)$ is not normalized to 1, ${h}_{0}\left(\xi \right)$ is rescaled automatically before use.
• refHistogramFile $⟨\phantom{\rule{0.3em}{0ex}}$Reference histogram ${h}_{0}\left(\xi \right)$$\phantom{\rule{0.3em}{0ex}}⟩$
Context: histogramRestraint
Acceptable values: UNIX ﬁle name
Description: Provides the values of ${h}_{0}\left(\xi \right)$ as contents of the corresponding ﬁle (mutually exclusive with refHistogram). The format is that of a text ﬁle, with each line containing the space-separated values of $\xi$ and ${h}_{0}\left(\xi \right)$. The same numerical conventions as refHistogram are used.
• writeHistogram $⟨\phantom{\rule{0.3em}{0ex}}$Periodically write the instantaneous histogram $h\left(\xi \right)$$\phantom{\rule{0.3em}{0ex}}⟩$
Acceptable values: boolean
Default value: off
Description: If on, the histogram $h\left(\xi \right)$ is written every colvarsRestartFrequency steps to a ﬁle with the name outputName.$<$name$>$.hist.dat This is useful to diagnose the convergence of $h\left(\xi \right)$ against ${h}_{0}\left(\xi \right)$.

#### 6.11 Deﬁning scripted biases

Rather than using the biasing methods described above, it is possible to apply biases provided at run time as a Tcl script. This option, also available in NAMD, can be useful to test a new algorithm to be used in a MD simulation.

• scriptedColvarForces $⟨\phantom{\rule{0.3em}{0ex}}$Enable custom, scripted forces on colvars $\phantom{\rule{0.3em}{0ex}}⟩$
Context: global
Acceptable values: boolean
Default value: off
Description: If this ﬂag is enabled, a Tcl procedure named calc_colvar_forces accepting one parameter should be deﬁned by the user. It is executed at each timestep, with the current step number as parameter, between the calculation of colvars and the application of bias forces. This procedure may use the cv command to access the values of colvars and apply forces on them, eﬀectively deﬁning custom collective variable biases.

#### 6.12 Performance of scripted biases

If concurrent computation over multiple threads is available (this is indicated by the message “SMP parallelism is available.” printed at initialization time), it is useful to take advantage of the scripting interface to combine many components, all computed in parallel, into a single variable.

The default SMP schedule is the following:

1. distribute the computation of all components across available threads;
2. on a single thread, collect the results of multi-component variables using polynomial combinations (see 4.12), or scripted functions (see 4.13);
3. distribute the computation of all biases across available threads;
4. compute on a single thread any scripted biases implemented via the keyword scriptedColvarForces.
5. communicate on a single thread forces to VMD.

The following options allow to ﬁne-tune this schedule:

• scriptingAfterBiases $⟨\phantom{\rule{0.3em}{0ex}}$Scripted colvar forces need updated biases?$\phantom{\rule{0.3em}{0ex}}⟩$
Context: global
Acceptable values: boolean
Default value: on
Description: This ﬂag speciﬁes that the calc_colvar_forces procedure (last step in the list above) is executed only after all biases have been updated (next-to-last step) For example, this allows using the energy of a restraint bias, or the force applied on a colvar, to calculate additional scripted forces, such as boundary constraints. When this ﬂag is set to off, it is assumed that only the values of the variables (but not the energy of the biases or applied forces) will be used by calc_colvar_forces: this can be used to schedule the calculation of scripted forces and biases concurrently to increase performance.

### 7 Syntax changes from older versions

The following is a list of syntax changes in Colvars since its ﬁrst release. Many of the older keywords are still recognized by the current code, thanks to speciﬁc compatibility code. This is not a list of new features: its primary purpose is to make you aware of those improvements that aﬀect the use of old conﬁguration ﬁles with new versions of the code.

Note: if you are using any of the NAMD and VMD tutorials:
https://www.ks.uiuc.edu/Training/Tutorials/
please be aware that several of these tutorials are not actively maintained: for those cases, this list will help you reconcile any inconsistencies.

• Colvars version 2016-06-09 or later (VMD version 1.9.3 or later).
The legacy keyword refPositionsGroup has been renamed fittingGroup for clarity (the legacy version is still supported).
• Colvars version 2016-08-10 or later (VMD version 1.9.3 or later).
“System forces” have been replaced by “total forces” (see for example outputTotalForce). See the following page for more information:
• Colvars version 2017-01-09 or later (VMD version 1.9.4 or later).
A new type of restraint, harmonicWalls (see 6.6), replaces and improves upon the legacy keywords lowerWall and upperWall: these are still supported as short-hands.
• Colvars version 2018-11-15 or later (VMD version 1.9.4 or later).
The global analysis keyword has been discontinued: speciﬁc analysis tasks are controlled directly by the keywords corrFunc and runAve, which continue to remain off by default.
• Deprecation warning for calculations including wall potentials.
The legacy keywords lowerWall and upperWall will stop having default values and will need to be set explicitly (preferably as part of the harmonicWalls restraint). When using an ABF bias, it is recommended to set the two walls equal to lowerBoundary and upperBoundary, respectively. When using a metadynamics bias, it is recommended to set the two walls within lowerBoundary and upperBoundary. This guarantees that the tails of each Gaussian hill are accounted in the region between the grid boundaries and the wall potentials. See also expandBoundaries for an automatic deﬁnition of the PMF grid boundaries.

Up-to-date documentation can always be accessed at:
https://colvars.github.io/colvars-refman-vmd/colvars-refman-vmd.html

### References

[1]   G. Fiorin, M. L. Klein, and J. Hénin. Using collective variables to drive molecular dynamics simulations. Mol. Phys., 111(22-23):3345–3362, 2013.

[2]   M. Iannuzzi, A. Laio, and M. Parrinello. Eﬃcient exploration of reactive potential energy surfaces using car-parrinello molecular dynamics. Phys. Rev. Lett., 90(23):238302, 2003.

[3]   E A Coutsias, C Seok, and K A Dill. Using quaternions to calculate RMSD. J. Comput. Chem., 25(15):1849–1857, 2004.

[4]   Davide Branduardi, Francesco Luigi Gervasio, and Michele Parrinello. From a to b in free energy space. J Chem Phys, 126(5):054103, 2007.

[5]   Yuguang Mu, Phuong H. Nguyen, and Gerhard Stock. Energy landscape of a small peptide revealed by dihedral angle principal component analysis. Proteins, 58(1):45–52, 2005.

[6]   Alexandros Altis, Phuong H. Nguyen, Rainer Hegger, and Gerhard Stock. Dihedral angle principal component analysis of molecular dynamics simulations. J. Chem. Phys., 126(24):244111, 2007.

[7]   Nicholas M Glykos. Carma: a molecular dynamics analysis program. J. Comput. Chem., 27(14):1765–1768, 2006.

[8]   Eric Darve, David Rodríguez-Gómez, and Andrew Pohorille. Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys., 128(14):144120, 2008.

[9]   J. Hénin and C. Chipot. Overcoming free energy barriers using unconstrained molecular dynamics simulations. J. Chem. Phys., 121:2904–2914, 2004.

[10]   J. Hénin, G. Fiorin, C. Chipot, and M. L. Klein. Exploring multidimensional free energy landscapes using time-dependent biases on collective variables. J. Chem. Theory Comput., 6(1):35–47, 2010.

[11]   A. Carter, E, G. Ciccotti, J. T. Hynes, and R. Kapral. Constrained reaction coordinate dynamics for the simulation of rare events. Chem. Phys. Lett., 156:472–477, 1989.

[12]   M. J. Ruiz-Montero, D. Frenkel, and J. J. Brey. Eﬃcient schemes to compute diﬀusive barrier crossing rates. Mol. Phys., 90:925–941, 1997.

[13]   W. K. den Otter. Thermodynamic integration of the free energy along a reaction coordinate in cartesian coordinates. J. Chem. Phys., 112:7283–7292, 2000.

[14]   Giovanni Ciccotti, Raymond Kapral, and Eric Vanden-Eijnden. Blue moon sampling, vectorial reaction coordinates, and unbiased constrained dynamics. ChemPhysChem, 6(9):1809–1814, 2005.

[15]   Adrien Lesage, Tony Lelièvre, Gabriel Stoltz, and Jérôme Hénin. Smoothed biasing forces yield unbiased free energies with the extended-system adaptive biasing force method. J. Phys. Chem. B, 121(15):3676–3685, 2017.

[16]   A. Laio and M. Parrinello. Escaping free-energy minima. Proc. Natl. Acad. Sci. USA, 99(20):12562–12566, 2002.

[17]   Helmut Grubmüller. Predicting slow structural transitions in macromolecular systems: Conformational ﬂooding. Phys. Rev. E, 52(3):2893–2906, Sep 1995.

[18]   T. Huber, A. E. Torda, and W.F. van Gunsteren. Local elevation - A method for improving the searching properties of molecular-dynamics simulation. Journal of Computer-Aided Molecular Design, 8(6):695–708, DEC 1994.

[19]   G. Bussi, A. Laio, and M. Parrinello. Equilibrium free energies from nonequilibrium metadynamics. Phys. Rev. Lett., 96(9):090601, 2006.

[20]   Alessandro Barducci, Giovanni Bussi, and Michele Parrinello. Well-tempered metadynamics: A smoothly converging and tunable free-energy method. Phys. Rev. Lett., 100:020603, 2008.

[21]   P. Raiteri, A. Laio, F. L. Gervasio, C. Micheletti, and M. Parrinello. Eﬃcient reconstruction of complex free energy landscapes by multiple walkers metadynamics. J. Phys. Chem. B, 110(8):3533–9, 2005.

[22]   Yuqing Deng and Benoît Roux. Computations of standard binding free energies with molecular dynamics simulations. J. Phys. Chem. B, 113(8):2234–2246, 2009.

[23]   Jed W. Pitera and John D. Chodera. On the use of experimental observations to bias simulated ensembles. J. Chem. Theory Comput., 8:3445–3451, 2012.

[24]   A. D. White and G. A. Voth. Eﬃcient and minimal method to bias molecular simulations with experimental data. J. Chem. Theory Comput., ASAP, 2014.

[25]   Rong Shen, Wei Han, Giacomo Fiorin, Shahidul M Islam, Klaus Schulten, and Benoît Roux. Structural reﬁnement of proteins by restrained molecular dynamics simulations with non-interacting molecular fragments. PLoS Comput. Biol., 11(10):e1004368, 2015.