Update Colvars to version 2018-12-14
[namd.git] / ug / ug_colvars.tex
1 % This file is part of the Collective Variables module (Colvars).
2 % The original version of Colvars and its updates are located at:
3 % https://github.com/colvars/colvars
4 % Please update all Colvars source files before making any changes.
5 % If you wish to distribute your changes, please submit them to the
6 % Colvars repository at GitHub.
7
8 \cvnamdugonly{
9 \section{Collective Variable-based Calculations (Colvars)}
10 \label{section:colvars}
11
12 The features described in this section were originally contributed to NAMD by Giacomo Fiorin (NIH) and J\'er\^ome H\'enin (CNRS, France) and are currently developed at this external repository:\\
13 \cvurl{https://github.com/Colvars/colvars}\\
14 An updated version of this section can also be downloaded as a separate manual: \\
15 HTML: \cvurl{https://colvars.github.io/colvars-refman-namd/colvars-refman-namd.html} \\
16 PDF: \cvurl{https://colvars.github.io/pdf/colvars-refman-namd.pdf}\\
17
18 See section~\ref{sec:colvars_config_changes} for specific changes that affect compatibility between versions.
19 Please ask any usage questions through the NAMD mailing list, and development questions through GitHub.
20
21 \subsection*{Overview}
22 }
23
24 \cvvmdugonly{
25 \chapter{Collective Variables Interface (Colvars)}
26 \label{section:colvars}
27
28 The features described in this section were originally contributed to VMD by Giacomo Fiorin (NIH) and J\'er\^ome H\'enin (CNRS, France) and are currently developed at this external repository:\\
29 \cvurl{https://github.com/Colvars/colvars}\\
30 An updated version of this section can also be downloaded as a separate manual: \\
31 HTML: \cvurl{https://colvars.github.io/colvars-refman-vmd/colvars-refman-vmd.html} \\
32 PDF: \cvurl{https://colvars.github.io/pdf/colvars-refman-vmd.pdf}\\
33
34 See section~\ref{sec:colvars_config_changes} for specific changes that affect compatibility between versions.
35 Please ask any usage questions through the VMD mailing list, and development questions through GitHub.
36
37 \subsection*{Overview}
38 }
39 \cvrefmanonly{
40 \cvsec{Overview}{sec:colvars_intro}
41 }
42
43 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 define biasing potentials to alter the dynamics of the system in a controlled manner.
44 These have been called `order parameters', `collective variables', `(surrogate) reaction coordinates', and many other terms.
45
46 Here we use primarily the term `collective variable', often shortened to \textit{colvar}, to indicate any differentiable function of atomic Cartesian coordinates, $\bm{x}_{i}$, with $i$ between $1$ and $N$, the total
47 number of atoms:
48 \begin{equation} 
49   \label{eq:colvar_basic}
50   \xi(t) \; = \xi(\bm{X}(t)) \; = \xi\left(\bm{x}_{i}(t), \bm{x}_{j}(t), \bm{x}_{k}(t),
51   \ldots \right)\;, \;\; 1 \leq i,j,k\ldots \leq N
52 \end{equation}
53 \cvvmdugonly{%
54 The Colvars module in VMD may be used to calculate these functions over a molecular structure, and to analyze the results of previous simulations.}
55 \cvrefmanonly{%
56 This manual documents the collective variables module (\textbf{Colvars}), a software that provides an implementation for the functions $\xi(\bm{X})$ with a focus on flexibility, robustness and high performance.}
57 The module is designed to perform multiple tasks concurrently during or after a simulation, the most common of which are:
58 \begin{itemize}
59
60 \item 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; \cvnamdonly{while this can in principle be done through a TclForces script, using the Colvars module is both easier and computationally more efficient;}
61
62 \item calculate potentials of mean force (PMFs) along any set of variables, using different 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;
63
64 \item 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 files, or after a simulation has been completed using VMD and the \texttt{cv} command\cvnamdugonly{ or NAMD and the \texttt{coorfile read} command as illustrated in \ref{section:sample}}.
65
66 \end{itemize}
67
68 \cvvmdonly{\textbf{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.  \emph{Options that only have an effect during a simulation are also included for compatibility purposes.}} 
69
70 Detailed explanations of the design of the Colvars module are provided in reference~\cite{Fiorin2013}. Please cite this reference whenever publishing work that makes use of this module.
71
72
73 \cvsec{A crash course}{sec:colvars_crash_course}
74
75 Suppose that we want to run a steered MD experiment where a small molecule is pulled away from a protein binding site.
76 In Colvars terms, this is done by applying a moving restraint to the distance between the two objects.
77 The configuration will contain two blocks, one defining the distance variable (see section \ref{sec:colvar} and \ref{sec:cvc_distance}), and the other the moving harmonic restraint (\ref{sec:colvarbias_harmonic}).
78
79 \bigskip
80 % verbatim can't appear within commands
81 {\noindent\ttfamily
82 \cvlammpsonly{indexFile index.ndx\\}colvar \{\\
83 \-~~name dist\\
84 \-~~distance \{\\
85 \-~~~~group1 \{ atomNumbersRange 42-55 \}\\
86 \cvnamebasedonly{\-~~~~group2 \{\\
87 \-~~~~~~psfSegID  PR\\
88 \-~~~~~~atomNameResidueRange  CA 15-30\\
89 \-~~~~\}\\}
90 \cvlammpsonly{\-~~~~group2 \{ indexGroup C-alpha\_15-30 \}\\}
91 \-~~\}\\
92 \}\\
93 \\
94 harmonic \{\\
95 \-~~colvars dist\\
96 \-~~forceConstant 20.0\\
97 \-~~centers 4.0~~~~~~~~~\# initial distance\\
98 \-~~targetCenters 15.0~~\# final distance\\
99 \-~~targetNumSteps 500000\\
100 \}\\
101 }
102
103 Reading this input in plain English: the variable here named \emph{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 \cvnamebasedonly{(segment name PR)}\cvlammpsonly{(selected from the index group ``\texttt{C-alpha\_15-30}'')}.
104 To the ``\emph{dist}'' variable, we apply a \texttt{harmonic} potential of force constant 20~\cvnamdonly{kcal/mol}\cvvmdonly{kcal/mol}\cvlammpsonly{(energy unit)}/\cvnamdonly{\AA}\cvvmdonly{\AA}\cvlammpsonly{(length unit)}$^2$, initially centered around a value of 4~\cvnamdonly{\AA}\cvvmdonly{\AA}\cvlammpsonly{(length unit)}, which will increase to 15~\cvnamdonly{\AA}\cvvmdonly{\AA}\cvlammpsonly{(length unit)} over 500,000 simulation steps.
105
106 The atom selection keywords is detailed in section \ref{sec:colvar_atom_groups}.
107 \cvlammpsonly{If the selection is too complex to implement only via internal keywords, an external index file may be created following the NDX format used in GROMACS:\\
108 \url{http://manual.gromacs.org/documentation/current/user-guide/file-formats.html}\\
109 or by using the \texttt{group2ndx} LAMMPS command.}
110
111
112
113 \cvsec{General parameters}{sec:colvarmodule}
114
115 Here, we document the syntax of the commands and parameters used to set up and use the Colvars module in \MDENGINE{}.
116 One of these parameters is the configuration file or the configuration text for the module itself, whose syntax is described in \ref{sec:colvars_config_syntax} and in the following sections.
117
118
119 \cvnamdonly{
120 \cvsubsec{NAMD parameters}{sec:colvars_mdengine_params}
121
122 To enable a Colvars-based calculation, two parameters must be added to the NAMD configuration file, \texttt{colvars} and \texttt{colvarsConfig}.
123 An optional third parameter, \texttt{colvarsInput}, can be used to continue a previous simulation.
124
125 \begin{itemize}
126 %  \setlength{\itemsep}{0.4cm}
127
128 \item %
129   \keydef
130     {colvars}{%
131     NAMD configuration file}{%
132     Enable the Colvars module}{%
133     boolean}{%
134     \texttt{off}}{%
135     If this flag is on, the Colvars module within
136     NAMD is enabled; the module requires a separate configuration
137     file, to be provided with \texttt{colvarsConfig}.}
138
139 \item %
140   \key
141     {colvarsConfig}{%
142     NAMD configuration file}{%
143     Configuration file for the collective variables}{%
144     UNIX filename}{%
145     This file contains the definition of all collective variables and
146     their biasing or analysis methods.
147     This file can also be provided by the Tcl command \texttt{cv configfile}; alternatively, the contents of the file itself can be given as an argument to the command \texttt{cv config}.
148   }
149
150 \item %
151   \key
152     {colvarsInput}{%
153     NAMD configuration file}{%
154     Input state file for the collective variables}{%
155     UNIX filename}{%
156     When continuing a previous simulation run, this file contains the current state of all collective variables and of their associated algorithms.
157     It is written automatically at the end of any simulation with collective variables.
158     The step number contained by this file will be used internally by Colvars to control time-dependent biases, unless \texttt{firstTimestep} is given, in which case that value will be used.
159     This file can also be loaded with the Tcl command \texttt{cv load}.
160   }
161
162 \end{itemize}
163 }
164
165
166 \cvlammpsonly{
167 \subsection{LAMMPS keywords}
168 \label{sec:colvars_mdengine_parameters}
169
170 To enable a Colvars-based calculation, the following line must be added to the LAMMPS configuration file:\\
171 \\
172 \texttt{fix } \emph{ID } \texttt{all }  \texttt{colvars } \emph{configfile } \emph{keyword value pairs ...}\\
173 \\
174 where \emph{ID} is a string that uniquely identifies this fix command inside a LAMMPS script, \emph{configfile} is the name of the configuration file for the Colvars module, followed by one or more of the following optional keywords with their corresponding arguments:
175
176 \begin{itemize}
177
178 \item %
179   \key
180     {input}{%
181     Keyword of the \texttt{fix colvars} command}{%
182     Name or prefix of the input state file}{%
183     string}{%
184     If a value is provided, it is interpreted as either the name of the input state file, or as the prefix of the file named \emph{input}\texttt{.colvars.state}.
185     This file contains information needed to continue a previous collective variables-based calculation, including the number of the last computed step (useful for time-dependent biases).
186     The same information is also stored in the binary restart files written by LAMMPS, so this option is not needed when continuing a calculation from a LAMMPS restart.
187   }
188
189 \item %
190   \keydef
191     {output}{%
192     Keyword of the \texttt{fix colvars} command}{%
193     Prefix of the output state file}{%
194     string}{%
195     ``out''}{%
196     If a value is provided, it is interpreted as the prefix to all output files that will be written by the Colvars module (see \ref{sec:colvars_output}).}
197
198 \item %
199   \keydef
200     {unwrap}{%
201     keyword of the \texttt{fix colvars} command}{%
202     Whether to unwrap coordinates passed to the Colvars module}{%
203     ``yes'' or ``no''}{%
204     ``yes''}{%
205     This keyword controls whether wrapped or unwrapped coordinates are passed to the Colvars module for calculation of the collective variables and of the resulting forces. The default is to use the image flags to reconstruct the absolute atom positions: under this convention, centers of mass and centers of geometry are calculated as a weighted vector sum (see \ref{sec:colvar_atom_groups_wrapping}).
206 Setting this to \emph{no} will use the current local coordinates that are wrapped back into the simulation cell at each re-neighboring instead.}
207
208 \item %
209   \keydef
210     {seed}{%
211     Keyword of the \texttt{fix colvars} command}{%
212     Seed for the random number generator}{%
213     positive integer}{%
214     1966}{%
215     If defined, the value of this keyword is provided as seed to the random number generator.
216     This is only meaningful when the \texttt{extendedLangevinDamping} keyword is used (see \ref{sec:colvar_extended}).}
217
218 \item %
219   \keydef
220     {tstat}{%
221     Keyword of the \texttt{fix colvars} command}{%
222     Thermostating fix}{%
223     string}{%
224     NULL}{%
225     This keyword provides the \emph{ID} of an applicable thermostating fix command. This will be used to provide the Colvars module with the current thermostat target temperature when using a method that needs this information.}
226
227 \end{itemize}
228
229 }
230
231
232 \cvscriptonly{
233 \cvsubsec{Using the \texttt{cv} command}{sec:cv_scripting}
234 \cvvmdugonly{\index{Colvars!\texttt{cv} command}}
235
236 At any moment during the execution of \MDENGINE{}, several options in the Colvars module can be read or modified by the command \texttt{cv} with the following syntax:\\
237 \noindent\texttt{cv~$<$subcommand$>$ [args...]}\\
238 \noindent{}For example, to record the value of a collective variable named \texttt{myVar} into the Tcl variable \texttt{value}, use the following syntax:\\
239 \noindent\texttt{set value [cv colvar myVar value]}\\
240 \noindent{}All subcommands of \texttt{cv} are documented below.
241
242 \cvvmdonly{
243 \textbf{Note:} in VMD, Colvars must be attached to one molecule (system).
244 Therefore, the \texttt{cv} command must be used for the first time as \texttt{cv molid }\emph{$<$molid$>$} to set up the Colvars module for the molecule identified by \emph{$<$molid$>$}.
245 In all following invocations, the \texttt{cv} command will continue operating on the same molecule, regardless of its ``top'' status.
246 To use the \texttt{cv} command on a different molecule, use \texttt{cv delete} first and then \texttt{cv molid }\emph{$<$molid$>$}.
247 Invoking the \texttt{cv} command with no arguments prints a help screen.
248 }
249
250 \cvvmdonly{
251 \cvsubsubsec{Example use of the \texttt{cv} command: analyze a trajectory}{sec:cv_command_example}
252
253 By far the most typical use of Colvars in VMD is computing the values of one or more variables along an existing trajectory:
254 \bigskip
255 {\noindent\ttfamily\\
256 \-{\bfseries\# Activate the module on the current VMD molecule}\\
257 \-cv molid top\\
258 \-{\bfseries\# Load a Colvars config file}\\
259 \-cv configfile test.in\\
260 \-set out [open "test.colvars.traj" "w"]\\
261 \-{\bfseries\# Write the labels to the file}\\
262 \-puts -nonewline \$\{out\} [cv printframelabels]\\
263 \-for \{ set fr 0 \} \{ \$\{fr\} < [molinfo top get numframes] \} \{ incr fr \} \{\\
264 \-~~{\bfseries\# Point Colvars to this trajectory frame}\\
265 \-~~cv frame \$\{fr\}\\
266 \-~~{\bfseries\# Recompute variables and biases}\\
267 \-~~cv update\\
268 \-~~{\bfseries\# Print variables and biases to the file}\\
269 \-~~puts -nonewline \$\{out\} [cv printframe]\\
270 \-\}\\
271 \-close \$\{out\}\\
272 }
273
274 The following sections explain the syntax of the various sub-commands of \texttt{cv}.
275 }
276
277
278 \cvsubsubsec{Managing the Colvars module}{sec:cv_command_general}
279
280 \begin{itemize}
281 \cvvmdonly{
282 \item \texttt{molid} \emph{$<$molid$>$}: setup the Colvars module for the given molecule; \emph{$<$molid$>$} is the identifier of a valid VMD molecule, either an integer or \texttt{top} for the top molecule;
283 \item \texttt{delete}: delete the Colvars module; after this, the module can be setup again using \texttt{cv molid};
284 }
285 \item \texttt{configfile} \emph{$<$file name$>$}: read configuration from a file;
286 \item \texttt{config} \emph{$<$string$>$}: read configuration from the given string; both \texttt{config} and \texttt{configfile} subcommands may be invoked multiple times;
287 \cvvmdonly{see \ref{sec:colvars_config_syntax} for the configuration syntax;}
288 \item \texttt{reset}: delete all internal configuration and atom selections of the Colvars module;
289 \item \texttt{version}: return the version of the Colvars module (see \ref{sec:colvars_config_changes} for any changes to older keywords).
290 \end{itemize}
291
292 \cvsubsubsec{Input and output commands}{sec:cv_command_io}
293
294 \begin{itemize}
295 \item \texttt{list}: return a list of all currently defined variables;
296 \cvvmdonly{See \ref{sec:colvar} and \ref{sec:colvar_atom_groups} for how to configure a collective variable;}
297 \item \texttt{list biases}: return a list of all currently defined biases (i.e.~sampling and analysis algorithms);
298 \item \texttt{load} \emph{$<$file name$>$}: load a collective variables state file, typically produced during a previous simulation;
299   \cvvmdonly{this parameter requires that the corresponding configuration has already been loaded by \texttt{configfile} or \texttt{config}; see \ref{sec:colvars_output} for a brief description of this file; the contents of this file are not required for as long as the VMD molecule has valid coordinates and \texttt{cv update} is used;}
300   \cvnamdonly{the step number contained by this file will be used internally by Colvars to control time-dependent biases, unless \texttt{firstTimestep} is given, in which case that value will be used;}
301 \item \texttt{load} \emph{$<$prefix$>$}: same as above, but without the \texttt{.colvars.state} suffix;
302 \item \texttt{save} \emph{$<$prefix$>$}: save the current state in a file whose name begins with the given argument; if any of the biases have additional output files defined, those are saved as well using the same prefix;
303 \item \texttt{update}: recalculate all variables and biases based on the current atomic coordinates\cvnamdonly{ (this is typically not needed in NAMD, as these functions are called already during a timestep)};
304 \item \texttt{addenergy} \emph{$<$E$>$}: add value \emph{E} to the total bias energy; this can be used within \texttt{calc\_colvar\_forces};
305 \item \texttt{printframe}: return a summary of the current frame, in a format equivalent to a line of the collective variables trajectory file;
306 \item \texttt{printframelabels}: return text labels for the columns of \texttt{printframe}'s output;
307 \cvvmdonly{
308 \item \texttt{frame}: return the current frame number, from which colvar values are calculated;
309 \item \texttt{frame} \emph{$<$new frame$>$}: set the frame number; returns 0 if successful, nonzero if the requested frame does not exist.
310 }
311 \end{itemize}
312
313 \cvsubsubsec{Managing collective variables}{sec:cv_command_colvar}
314
315 \begin{itemize}
316 \item \texttt{colvar} \emph{$<$name$>$} \texttt{value}: return the current value of colvar \emph{$<$name$>$};
317 \item \texttt{colvar} \emph{$<$name$>$} \texttt{update}: recalculate colvar \emph{$<$name$>$};
318 \item \texttt{colvar} \emph{$<$name$>$} \texttt{type}: return the type of colvar \emph{$<$name$>$};
319 \cvvmdonly{
320 \item \texttt{colvar} \emph{$<$name$>$} \texttt{delete}: delete colvar \emph{$<$name$>$};
321 }
322 \cvnamdonly{
323 \item \texttt{colvar} \emph{$<$name$>$} \texttt{addforce} \emph{$<$F$>$}: apply given force on colvar \emph{$<$name$>$}; this can be used within \texttt{calc\_colvar\_forces}
324 \item \texttt{colvar} \emph{$<$name$>$} \texttt{getappliedforce}: get the force currently being applied on colvar \emph{$<$name$>$};
325 \item \texttt{colvar} \emph{$<$name$>$} \texttt{gettotalforce}: get the total force acting on colvar \emph{$<$name$>$} at the previous step (see \ref{sec:cvc_sys_forces});
326 }
327 \item \texttt{colvar} \emph{$<$name$>$} \texttt{getconfig}: return config string of colvar \emph{$<$name$>$}.
328 \item \texttt{colvar} \emph{$<$name$>$} \texttt{cvcflags} \emph{$<$flags$>$}: for a colvar with several CVCs (numbered according to their name string order), set which CVCs are enabled or disabled in subsequent evaluations according to a list of 0/1 flags (one per CVC).
329 \item \texttt{colvar} \emph{$<$name$>$} \texttt{modifycvcs} \emph{$<$strings$>$}: for a colvar with one or more CVCs (numbered according to their name string order), pass a list of new configuration strings to modify each CVC without needing to delete the colvar.
330 This option is currently limited to changing the values of \refkey{componentCoeff}{sec:cvc_superp} and \refkey{componentExp}{sec:cvc_superp} (e.g.{} to update the polynomial superposition parameters on the fly), of \refkey{period}{sec:cvc_periodic} and \refkey{wrapAround}{sec:cvc_periodic} (e.g.{} to update the period of variables such as \texttt{distanceZ}), and of the \texttt{forceNoPBC} option for those CVCs that support it.
331 Changes in configuration done by this command are not saved to state files, and are lost when restarting a simulation, deleting the parent colvar, or resetting the module with \cvvmdonly{\texttt{cv delete} or} \texttt{cv reset}.
332
333 \end{itemize}
334
335 \cvsubsubsec{Managing biases}{sec:cv_command_bias}
336
337 \begin{itemize}
338 \item \texttt{bias} \emph{$<$name$>$} \texttt{energy}: return the current energy of the bias \emph{$<$name$>$};
339 \item \texttt{bias} \emph{$<$name$>$} \texttt{update}: recalculate the bias \emph{$<$name$>$};
340 \item \texttt{bias} \emph{$<$name$>$} \texttt{delete}: delete the bias \emph{$<$name$>$};
341 \item \texttt{bias} \emph{$<$name$>$} \texttt{getconfig}: return config string of bias \emph{$<$name$>$}.
342 \end{itemize}
343
344 }
345
346
347 \cvsubsec{Configuration syntax}{sec:colvars_config_syntax}
348
349 \cvnamdonly{All the parameters defining variables and their biasing or analysis algorithms are read from the file specified by the configuration option \texttt{colvarsConfig}, or by the Tcl commands \texttt{cv config} and \texttt{cv configfile}.
350 \emph{None of the keywords described in the remainder of this manual are recognized directly in the NAMD configuration file, unless as arguments of \texttt{cv config}.}}
351 \cvvmdonly{The Colvars configuration is usually read using the commands \texttt{cv configfile} or \texttt{cv config}.}
352 The syntax of the Colvars configuration is ``\texttt{keyword value}'', where the keyword and its value are separated by any white space.
353 The following rules apply:
354
355 \begin{itemize}
356
357 \item keywords are case-insensitive (\texttt{upperBoundary} is the same as \texttt{upperboundary} and \texttt{UPPERBOUNDARY}): their string values are however case-sensitive (e.g.~file names);
358
359 \item a long value or a list of multiple values can be distributed across multiple lines by using curly braces, ``\texttt{\{}'' and ``\texttt{\}}'': the opening brace ``\texttt{\{}'' must occur on the same line as the keyword, following a space character or other white space; the closing brace ``\texttt{\}}'' can be at any position after that;
360
361 \item many keywords are nested, and are only meaningful within a specific context: for every keyword documented in the following, the ``parent'' keyword that defines such context is also indicated\cvnamdugonly{ in parentheses};
362
363 \cvnamdonly{%
364 \item the `\texttt{=}' sign between a keyword and its value, deprecated in the NAMD main configuration file, is not allowed;
365
366 \item Tcl syntax is generally not available, but it is possible to use Tcl variables or bracket expansion of commands within a configuration string, when this is passed via the command \texttt{cv config \ldots}; this is particularly useful when combined with parameter introspection\cvnamdugonly{ (see \ref{section:tclscripting})}, e.g.{} \texttt{cv config "colvarsTrajFrequency [DCDFreq]"};
367 }
368
369 \cvvmdonly{%
370 \item Tcl syntax is generally not available, but it is possible to use Tcl variables or bracket expansion of commands within a configuration string, when this is passed via the command \texttt{cv config \ldots}: for example, it is possible to convert the atom selection \$\emph{sel} into an atom group (see \ref{sec:colvar_atom_groups_sel}) using \texttt{cv config "atomNumbers \{ [\$sel get serial] \}"};}
371
372 \item if a keyword requiring a boolean value (\texttt{yes|on|true} or \texttt{no|off|false}) is provided without an explicit value, it defaults to `\texttt{yes|on|true}'; for example, `\texttt{outputAppliedForce}' may be used as shorthand for `\texttt{outputAppliedForce on}';
373
374 \item the hash character \texttt{\#} indicates a comment: all text in the same line following this character will be ignored.
375
376 \end{itemize}
377
378
379 \cvsubsec{Global keywords}{sec:colvars_global}
380
381 The following keywords are available in the global context of the Colvars configuration, i.e.~they are not nested inside other keywords:
382 \begin{itemize}
383
384 \item %
385   \labelkey{Colvars-global|colvarsTrajFrequency}
386   \keydef
387     {colvarsTrajFrequency}{%
388     global}{%
389   Colvar value trajectory frequency}{%
390     positive integer}{%
391     \texttt{100}}{%
392     The values of each colvar (and of other related quantities, if requested) are written to the file \outputName\texttt{.colvars.traj} every these many steps throughout the simulation.
393     If the value is \texttt{0}, such trajectory file is not written.
394     For optimization the output is buffered, and synchronized with the disk only when the restart file is being written.}
395
396 \item %
397   \labelkey{Colvars-global|colvarsRestartFrequency}
398   \keydef
399     {colvarsRestartFrequency}{%
400     global}{%
401     Colvar module restart frequency}{%
402     positive integer}{%
403     \texttt{restartFreq}}{%
404     Allows to choose a different restart frequency for the Colvars module.
405     Redefining it may be useful to trace the time
406     evolution of those few properties which are not written to the
407     trajectory file for reasons of disk space.}
408
409 \item %
410   \labelkey{Colvars-global|indexFile}
411   \key
412     {indexFile}{%
413     global}{%
414     Index file for atom selection (GROMACS ``ndx'' format)}{%
415     UNIX filename}{%
416     This option reads an index file (usually with a \texttt{.ndx}
417     extension) as produced by the \texttt{make\_ndx} tool of GROMACS.
418     This keyword may be repeated to load multiple index files: the same group name cannot appear in multiple index files.
419     \cvlammpsonly{In LAMMPS, the \texttt{group2ndx} command can be used to generate such file from existing groups.
420     Note that the Colvars module reads the indices of atoms from the index file: therefore, the LAMMPS groups do not need to remain active during the simulation, and could be deleted right after issuing \texttt{group2ndx}.
421 }
422     The names of index groups contained in this file can then be used to define
423     atom groups with the \texttt{indexGroup} keyword.
424     Other supported methods to select atoms are described in \ref{sec:colvar_atom_groups}.
425   }
426
427 \item %
428   \labelkey{Colvars-global|smp}
429   \keydef
430     {smp}{%
431     global}{%
432     Whether SMP parallelism should be used}{%
433     boolean}{%
434     \texttt{on}}{%
435     If this flag is enabled (default), SMP parallelism over threads will be used to compute variables and biases, provided that this is supported by the \MDENGINE{} build in use.}
436     
437 \end{itemize}
438
439
440 To illustrate the flexibility of the Colvars module, a non-trivial setup is represented in Figure~\ref{fig:colvars_diagram}.
441 The corresponding configuration is given below. The options within the \texttt{colvar} blocks are described in \ref{sec:colvar}, those within the \texttt{harmonic} and \texttt{histogram} blocks in \ref{sec:colvarbias}.
442 \textbf{Note:} \emph{except }\texttt{colvar}\emph{, none of the keywords shown is mandatory}.\\
443
444 \begin{figure}[!ht]
445   \centering
446 \cvnamdugonly{\includegraphics[width=12cm]{figures/colvars_diagram}}
447 \cvvmdugonly{\includegraphics[width=12cm]{pictures/colvars_diagram}}
448 \cvrefmanonly{\ifdefined\HCode{\HCode{<img class="diagram" src="colvars_diagram.png" alt="Colvars diagram">}}\else\includegraphics[width=12cm]{colvars_diagram}\fi}
449   \caption[Graphical representation of a Colvars configuration.]{Graphical representation of a Colvars configuration\cvlammpsonly{ \textbf{(note:} \emph{currently, the $\alpha$-helical content colvar is unavailable in LAMMPS)}}.
450     The colvar called ``$d$'' is defined as the difference between two distances: the first distance ($d_{1}$) is taken between the center of mass of atoms 1 and 2 and that of atoms 3 to 5, the second ($d_{2}$) between atom 7 and the center of mass of atoms 8 to 10.
451 The difference $d = d_{1} - d_{2}$ is obtained by multiplying the two by a coefficient $C = +1$ or $C = -1$, respectively.
452 The colvar called ``$c$'' is the coordination number calculated between atoms 1 to 10 and atoms 11 to 20.  A harmonic restraint is applied to both $d$ and $c$: to allow using the same force constant $K$, both $d$ and $c$ are scaled by their respective fluctuation widths $w_d$ and $w_c$.
453 \cvnamebasedonly{A third colvar ``alpha'' is defined as the $\alpha$-helical content of residues 1 to 10.}
454 The values of ``$c$''\cvnamebasedonly{ and ``alpha''} are also recorded throughout the simulation as a joint 2-dimensional histogram.
455 }
456   \label{fig:colvars_diagram}
457 \end{figure}
458
459 {%
460 % verbatim can't appear within commands
461 \noindent\ttfamily
462 colvar \{\\
463 \-~~\# difference of two distances\\
464 \-~~name d \\
465 \-~~width 0.2  \# 0.2 \AA{} of estimated fluctuation width \\
466 \-~~distance \{\\
467 \-~~~~componentCoeff  1.0\\
468 \-~~~~group1 \{ atomNumbers 1 2 \}\\
469 \-~~~~group2 \{ atomNumbers 3 4 5 \}\\
470 \-~~\}\\
471 \-~~distance \{\\
472 \-~~~~componentCoeff -1.0\\
473 \-~~~~group1 \{ atomNumbers 7 \}\\
474 \-~~~~group2 \{ atomNumbers 8 9 10 \}\\
475 \-~~\}\\
476 \}\\
477 \\
478 colvar \{\\
479 \-~~name c\\
480 \-~~coordNum \{\\
481 \-~~~~cutoff 6.0\\
482 \-~~~~group1 \{ atomNumbersRange  1-10 \}\\
483 \-~~~~group2 \{ atomNumbersRange 11-20 \}\\
484 \-~~\}\\
485 \}\\}
486 \cvnamebasedonly{{%
487 \noindent\ttfamily\\
488 colvar \{\\
489 \-~~name alpha\\
490 \-~~alpha \{\\
491 \-~~~~psfSegID PROT\\
492 \-~~~~residueRange 1-10\\
493 \-~~\}\\
494 \}}}
495 {%
496 \noindent\ttfamily\\
497 \\
498 harmonic \{\\
499 \-~~colvars d c\\
500 \-~~centers 3.0 4.0\\
501 \-~~forceConstant 5.0\\
502 \}\\
503
504 \noindent histogram \{\\
505 \-~~colvars c\cvnamebasedonly{ alpha}\\
506 \}\\}
507
508 %\cvlammpsonly{\textbf{Note:} \emph{currently, the $\alpha$-helical content colvar is unavailable in LAMMPS, as it requires a name-based topology; future releases will overcome this limitation.}}
509
510 Section \ref{sec:colvar} explains how to define a colvar and its behavior, regardless of its specific functional form.
511 To define colvars that are appropriate to a specific physical system, Section \ref{sec:colvar_atom_groups} documents how to select atoms, and section \ref{sec:colvar} lists all of the available functional forms, which we call ``colvar components''.
512 Finally, section \ref{sec:colvarbias} lists the available methods and algorithms to perform biased simulations and multidimensional analysis of colvars.
513
514
515 \cvsubsec{Output files}{sec:colvars_output}
516
517 During a simulation with collective variables defined, the following three output files are written:
518
519 \begin{itemize}
520
521 \item A \emph{state file}, named \outputName\texttt{.colvars.state}; this file is in ASCII (plain text) format\cvnamdonly{, regardless of the value of \texttt{binaryOutput} in the NAMD configuration}.  This file is written at the end of the specified run\cvscriptonly{, but can also be written at any time with the command \texttt{cv save} (\ref{sec:cv_command_io})}.\\
522   \emph{This is the only Colvars output file needed to continue a simulation.}
523
524 \item If the parameter \refkey{colvarsRestartFrequency}{Colvars-global|colvarsRestartFrequency} is larger than zero, a \emph{restart file} is written every that many steps: this file is fully equivalent to the final state file.
525   The name of this file is \restartName\texttt{.colvars.state}.
526
527 \item If the parameter \refkey{colvarsTrajFrequency}{Colvars-global|colvarsTrajFrequency} is greater than 0 (default: 100), a \emph{trajectory file} is written during the simulation: its name is \outputName\texttt{.colvars.traj}; unlike the state file, it is not needed to restart a simulation, but can be used later for post-processing and analysis.
528
529 \end{itemize}
530
531 Other output files may also be written by specific methods, e.g.{} the ABF or metadynamics methods (\ref{sec:colvarbias_abf}, \ref{sec:colvarbias_meta}).
532 Like the trajectory file, they are needed only for analyzing, not continuing a simulation.
533 All such files' names also begin with the prefix \outputName.
534
535 \cvnamdonly{Lastly, the total energy of all biases or restraints applied to the colvars appears under the NAMD standard output, under the MISC column.}
536
537
538 \cvsec{Defining collective variables}{sec:colvar}
539
540 A collective variable is defined by the keyword \texttt{colvar} followed by its configuration options contained within curly braces:
541
542 {\bigskip\noindent\ttfamily
543 colvar \{\\
544 \-~~name xi\\
545 \-~~$<$other options$>$\\
546 \-~~function\_name \{\\
547 \-~~~~$<$parameters$>$\\
548 \-~~~~$<$atom selection$>$\\
549 \-~~\}\\
550 \}\\
551 }
552
553 \noindent{}There are multiple ways of defining a variable:
554 \begin{itemize}
555 \item The \emph{simplest and most common way} way is using one of the precompiled functions (here called ``components''), which are listed in section~\ref{sec:cvc_list}.  For example, using the keyword \texttt{rmsd} (section \ref{sec:cvc_rmsd}) defines the variable as the root mean squared deviation (RMSD) of the selected atoms.
556 \item A new variable may also be constructed as a linear or polynomial combination of the components listed in section~\ref{sec:cvc_list} (see \ref{sec:cvc_superp} for details).
557 \cvleptononly{
558 \item A user-defined mathematical function of the existing components (see list in section~\ref{sec:cvc_list}), or of the atomic coordinates directly (see the \texttt{cartesian} keyword in \ref{sec:cvc_cartesian}).
559 The function is defined through the keyword \refkey{customFunction}{colvar|customFunction} (see \ref{sec:colvar_custom_function} for details).
560 }
561 \cvscriptonly{
562 \item A user-defined Tcl function of the existing components (see list in section~\ref{sec:cvc_list}), or of the atomic coordinates directly (see the \texttt{cartesian} keyword in \ref{sec:cvc_cartesian}).
563 The function is provided by a separate Tcl script, and referenced through the keyword \refkey{scriptedFunction}{colvar|scriptedFunction} (see \ref{sec:colvar_scripted} for details).
564 }
565 \end{itemize}
566 Choosing a component (function) is the only parameter strictly required to define a collective variable.
567 It is also highly recommended to specify a name for the variable:
568 \begin{itemize}
569 \label{sec:colvar_general}
570 \item %
571   \labelkey{colvar|name}
572   \keydef
573     {name}{%
574     \texttt{colvar}}{%
575     Name of this colvar}{%
576     string}{%
577     ``\texttt{colvar}'' + numeric id}{%
578     The name is an unique case-sensitive string which allows the
579     Colvars module to identify this colvar unambiguously; it is also
580     used in the trajectory file to label to the columns corresponding
581     to this colvar.}
582
583 \end{itemize}
584
585
586 \cvsubsec{Choosing a function}{sec:cvc_list}
587
588 In this context, the function that computes a colvar is called a \emph{component}.
589 A component's choice and definition consists of including in the variable's configuration a keyword indicating the type of function (e.g.{} \texttt{rmsd}), followed by a definition block specifying the atoms involved (see \ref{sec:colvar_atom_groups}) and any additional parameters (cutoffs, ``reference'' values, \ldots).
590 \emph{At least one component must be chosen to define a variable:} if none of the keywords listed below is found, an error is raised.
591
592 The following components implement functions with a scalar value (i.e.{} a real number):
593 \begin{itemize}
594 \item \refkey{distance}{colvar|distance}: distance between two groups;
595 \item \refkey{distanceZ}{colvar|distanceZ}: projection of a distance vector on an axis;
596 \item \refkey{distanceXY}{colvar|distanceXY}: projection of a distance vector on a plane;
597 \item \refkey{distanceInv}{colvar|distanceInv}: mean distance between two groups of atoms (e.g.~NOE-based distance);
598 \item \refkey{angle}{colvar|angle}: angle between three groups;
599 \item \refkey{dihedral}{colvar|dihedral}: torsional (dihedral) angle between four groups;
600 \item \refkey{dipoleAngle}{colvar|dipoleAngle}: angle between two groups and dipole of a third group;
601 \item \refkey{polarTheta}{colvar|polarTheta}: polar angle of a group in spherical coordinates;
602 \item \refkey{polarPhi}{colvar|polarPhi}: azimuthal angle of a group in spherical coordinates;
603 \item \refkey{coordNum}{colvar|coordNum}: coordination number between two groups;
604 \item \refkey{selfCoordNum}{colvar|selfCoordNum}: coordination number of atoms within a
605   group;
606 \item \refkey{hBond}{colvar|hBond}: hydrogen bond between two atoms;
607 \item \refkey{rmsd}{colvar|rmsd}: root mean square deviation (RMSD) from a set of
608   reference coordinates;
609 \item \refkey{eigenvector}{colvar|eigenvector}: projection of the atomic coordinates on a
610   vector;
611 \item \refkey{orientationAngle}{colvar|orientationAngle}: angle of the best-fit rotation from
612   a set of reference coordinates;
613 \item \refkey{orientationProj}{colvar|orientationProj}: cosine of \refkey{orientationProj}{colvar|orientationProj};
614 \item \refkey{spinAngle}{colvar|spinAngle}: projection orthogonal to an axis of the best-fit rotation
615   from a set of reference coordinates;
616 \item \refkey{tilt}{colvar|tilt}: projection on an axis of the best-fit rotation
617   from a set of reference coordinates;
618 \item \refkey{gyration}{colvar|gyration}: radius of gyration of a group of atoms;
619 \item \refkey{inertia}{colvar|inertia}: moment of inertia of a group of atoms;
620 \item \refkey{inertiaZ}{colvar|inertiaZ}: moment of inertia of a group of atoms around a chosen axis;
621 \cvnamebasedonly{
622 \item \refkey{alpha}{colvar|alpha}: $\alpha$-helix content of a protein segment.
623 \item \refkey{dihedralPC}{colvar|dihedralPC}: projection of protein backbone dihedrals onto a dihedral principal component.
624 }
625 \end{itemize}
626
627 Some components do not return scalar, but vector values:
628 \begin{itemize}
629 \item \refkey{distanceVec}{colvar|distanceVec}: distance vector between two groups (length: 3);
630 \item \refkey{distanceDir}{colvar|distanceDir}: unit vector parallel to distanceVec (length: 3);
631 \item \refkey{cartesian}{colvar|cartesian}: vector of atomic Cartesian coordinates (length: $N$ times the number of Cartesian components requested, X, Y or Z);
632 \item \refkey{distancePairs}{colvar|distancePairs}: vector of mutual distances (length: $N_{\mathrm{1}}\times{}N_{\mathrm{2}}$);
633 \item \refkey{orientation}{colvar|orientation}: best-fit rotation, expressed as a unit quaternion (length: 4).
634 \end{itemize}
635
636 The types of components used in a colvar (scalar or not) determine the
637 properties of that colvar, and particularly which biasing or analysis methods
638 can be applied.
639
640 \textbf{What if ``X'' is not listed?} If a function type is not available on this list, it may be possible to define it as a polynomial superposition of existing ones (see \ref{sec:cvc_superp})\cvleptononly{, a custom function (see \ref{sec:colvar_custom_function})}\cvscriptonly{, or a scripted function (see \ref{sec:colvar_scripted})}.
641
642 In the rest of this section, all available component types are listed, along with their physical units and the ranges of values, if limited.
643 Such limiting values can be used to define \refkey{lowerBoundary}{colvar|lowerBoundary} and \refkey{upperBoundary}{colvar|upperBoundary} in the parent colvar.
644
645 For each type of component, the available configurations keywords are listed:
646 when two components share certain keywords, the second component references to
647 the documentation of the first one that uses that keyword.
648 The very few keywords that are available for all types of components are listed in a separate section \ref{sec:cvc_common}.
649
650 \newenvironment{cvcoptions}%
651   {\noindent\textbf{List of keywords} (see also \ref{sec:cvc_superp} for additional options):
652   \begin{itemize}}
653   {
654   \end{itemize}
655 }
656
657 \cvsubsubsec{\texttt{distance}: center-of-mass distance between two groups.}{sec:cvc_distance}
658 \labelkey{colvar|distance}
659
660 The \texttt{distance \{...\}} block defines a distance component between the two atom groups, \texttt{group1} and \texttt{group2}.
661
662 \begin{cvcoptions}
663 \item %
664   \labelkey{colvar|distance|group1}
665   \key
666     {group1}{%
667     \texttt{distance}}{%
668     First group of atoms}{%
669     Block \texttt{group1 \{...\}}}{%
670     First group of atoms.}
671
672 \item %
673   \labelkey{colvar|distance|group2}
674   \simkey{group2}{\texttt{distance}}{group1}
675
676 \item %
677   \labelkey{colvar|distance|forceNoPBC}
678   \keydef
679     {forceNoPBC}{%
680     \texttt{distance}}{%
681     Calculate absolute rather than minimum-image distance?}{%
682     boolean}{%
683     \texttt{no}}{%
684     By default, in calculations with periodic boundary conditions, the
685     \texttt{distance} component returns the distance according to the
686     minimum-image convention. If this parameter is set to \texttt{yes},
687     PBC will be ignored and the distance between the coordinates as maintained
688     internally will be used. This is only useful in a limited number of
689     special cases, e.g. to describe the distance between remote points
690     of a single macromolecule, which cannot be split across periodic cell
691     boundaries, and for which the minimum-image distance might give the
692     wrong result because of a relatively small periodic cell.}
693
694 \item %
695   \labelkey{colvar|distance|oneSiteTotalForce}
696   \keydef
697     {oneSiteTotalForce}{%
698     \texttt{angle}, \texttt{dipoleAngle}, \texttt{dihedral}}{%
699     Measure total force on group 1 only?}{%
700     boolean}{%
701     \texttt{no}}{%
702     If this is set to \texttt{yes}, the total force is measured along
703     a vector field (see equation~(\ref{eq:gradient_vector}) in
704     section~\ref{sec:colvarbias_abf}) that only involves atoms of
705     \texttt{group1}.  This option is only useful for ABF, or custom
706     biases that compute total forces.  See
707     section~\ref{sec:colvarbias_abf} for details.}
708
709 \end{cvcoptions}
710
711 The value returned is a positive number (in \AA), ranging from $0$
712 to the largest possible interatomic distance within the chosen
713 boundary conditions (with PBCs, the minimum image convention is used
714 unless the \texttt{forceNoPBC} option is set).
715
716
717 \cvsubsubsec{\texttt{distanceZ}: projection of a distance vector on an axis.}{sec:cvc_distanceZ}
718 \labelkey{colvar|distanceZ}
719
720 The \texttt{distanceZ~\{...\}} block defines a distance projection
721 component, which can be seen as measuring the distance between two
722 groups projected onto an axis, or the position of a group along such
723 an axis.  The axis can be defined using either one reference group and
724 a constant vector, or dynamically based on two reference groups.
725 One of the groups can be set to a dummy atom to allow the use of an absolute Cartesian coordinate.
726
727 \begin{cvcoptions}
728 \item %
729   \labelkey{colvar|distanceZ|main}
730   \key
731     {main}{%
732     \texttt{distanceZ}}{%
733     Main group of atoms}{%
734     Block \texttt{main \{...\}}}{%
735     Group of atoms whose position $\bm{r}$ is measured.}
736
737 \item %
738   \labelkey{colvar|distanceZ|ref}
739   \key
740     {ref}{%
741     \texttt{distanceZ}}{%
742     Reference group of
743     atoms}{%
744     Block \texttt{ref \{...\}}}{%
745     Reference group of atoms.  The position of its center of mass is
746     noted $\bm{r}_1$ below.}
747
748 \item %
749   \labelkey{colvar|distanceZ|ref2}
750   \keydef
751     {ref2}{%
752     \texttt{distanceZ}}{%
753     Secondary reference
754     group}{%
755     Block \texttt{ref2 \{...\}}}{%
756     none}{%
757     Optional group of reference atoms, whose position $\bm{r}_2$ can
758     be used to define a dynamic projection axis: $\bm{e}=(\| \bm{r}_2
759     - \bm{r}_1\|)^{-1} \times (\bm{r}_2 - \bm{r}_1)$.  In this case,
760     the origin is $\bm{r}_m = 1/2 (\bm{r}_1+\bm{r}_2)$, and the value
761     of the component is $\bm{e} \cdot (\bm{r}-\bm{r}_m)$.}
762
763 \item %
764   \labelkey{colvar|distanceZ|axis}
765   \keydef
766     {axis}{%
767     \texttt{distanceZ}}{%
768     Projection axis (\AA{})}{%
769     \texttt{(x, y, z)} triplet}{%
770     \texttt{(0.0, 0.0, 1.0)}}{%
771     The three components of this vector define a
772     projection axis $\bm{e}$ for the distance vector $\bm{r} -
773     \bm{r}_1$ joining the centers of groups \texttt{ref} and
774     \texttt{main}. The value of the component is then $\bm{e} \cdot
775     (\bm{r}-\bm{r}_1)$.  The vector should be written as three
776     components separated by commas and enclosed in parentheses.}
777
778 \item %
779   \dupkey{forceNoPBC}{\texttt{distanceZ}}{colvar|distance|forceNoPBC}{\texttt{distance} component}
780
781 \item %
782   \dupkey{oneSiteTotalForce}{\texttt{distanceZ}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component}
783 \end{cvcoptions}
784 This component returns a number (in \AA{}) whose range is determined
785 by the chosen boundary conditions.  For instance, if the $z$ axis is
786 used in a simulation with periodic boundaries, the returned value ranges
787 between $-b_{z}/2$ and $b_{z}/2$, where $b_{z}$ is the box length
788 along $z$ (this behavior is disabled if \texttt{forceNoPBC} is set).
789
790
791 \cvsubsubsec{\texttt{distanceXY}: modulus of the projection of a distance vector on a plane.}{sec:cvc_distanceXY}
792 \labelkey{colvar|distanceXY}
793
794 The \texttt{distanceXY~\{...\}} block defines a distance projected on
795 a plane, and accepts the same keywords as the component \texttt{distanceZ}, i.e.
796 \texttt{main}, \texttt{ref}, either \texttt{ref2} or \texttt{axis},
797 and \texttt{oneSiteTotalForce}.  It returns the norm of the
798 projection of the distance vector between \texttt{main} and
799 \texttt{ref} onto the plane orthogonal to the axis.  The axis is
800 defined using the \texttt{axis} parameter or as the vector joining
801 \texttt{ref} and \texttt{ref2} (see \texttt{distanceZ} above).
802
803 \begin{cvcoptions}
804 \item %
805   \dupkey{main}{\texttt{distanceXY}}{colvar|distanceZ|main}{\texttt{distanceZ} component}
806 \item %
807   \dupkey{ref}{\texttt{distanceXY}}{colvar|distanceZ|ref}{\texttt{distanceZ} component}
808 \item %
809   \dupkey{ref2}{\texttt{distanceXY}}{colvar|distanceZ|ref2}{\texttt{distanceZ} component}
810 \item %
811   \dupkey{axis}{\texttt{distanceXY}}{colvar|distanceZ|axis}{\texttt{distanceZ} component}
812 \item %
813   \dupkey{forceNoPBC}{\texttt{distanceXY}}{colvar|distance|forceNoPBC}{\texttt{distance} component}
814 \item %
815   \dupkey{oneSiteTotalForce}{\texttt{distanceZ}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component}
816 \end{cvcoptions}
817
818
819
820 \cvsubsubsec{\texttt{distanceVec}: distance vector  between two groups.}{sec:cvc_distanceVec}
821 \labelkey{colvar|distanceVec}
822
823 The \texttt{distanceVec~\{...\}} block defines
824 a distance vector component, which accepts the same keywords as
825 the component \texttt{distance}: \texttt{group1}, \texttt{group2}, and
826 \texttt{forceNoPBC}. Its value is the 3-vector joining the centers
827 of mass of \texttt{group1} and \texttt{group2}.
828
829 \begin{cvcoptions}
830 \item %
831   \dupkey{group1}{\texttt{distanceVec}}{colvar|distance|group1}{\texttt{distance} component}
832 \item %
833   \simkey{group2}{\texttt{distanceVec}}{group1}
834 \item %
835   \dupkey{forceNoPBC}{\texttt{distanceVec}}{colvar|distance|forceNoPBC}{\texttt{distance} component}
836 \item %
837   \dupkey{oneSiteTotalForce}{\texttt{distanceVec}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component}
838 \end{cvcoptions}
839
840
841
842 \cvsubsubsec{\texttt{distanceDir}: distance unit vector between two groups.}{sec:cvc_distanceDir}
843 \labelkey{colvar|distanceDir}
844
845 The \texttt{distanceDir~\{...\}} block defines
846 a distance unit vector component, which accepts the same keywords as
847 the component \texttt{distance}: \texttt{group1}, \texttt{group2}, and
848 \texttt{forceNoPBC}.  It returns a
849 3-dimensional unit vector $\mathbf{d} = (d_{x}, d_{y}, d_{z})$, with
850 $|\mathbf{d}| = 1$.
851
852 \begin{cvcoptions}
853 \item %
854   \dupkey{group1}{\texttt{distanceDir}}{colvar|distance|group1}{\texttt{distance} component}
855 \item %
856   \simkey{group2}{\texttt{distanceDir}}{group1}
857 \item %
858   \dupkey{forceNoPBC}{\texttt{distanceDir}}{colvar|distance|forceNoPBC}{\texttt{distance} component}
859 \item %
860   \dupkey{oneSiteTotalForce}{\texttt{distanceDir}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component}
861 \end{cvcoptions}
862
863
864 \cvsubsubsec{\texttt{distanceInv}: mean distance between two groups of atoms.}{sec:cvc_distanceInv}
865 \labelkey{colvar|distanceInv}
866
867 The \texttt{distanceInv~\{...\}} block defines a generalized mean distance between two groups of atoms 1 and 2, weighted with exponent $1/n$:
868 \begin{equation}
869   \label{eq:distanceInv}
870   d_{\mathrm{1,2}}^{[n]} \; = \;   \left(\frac{1}{N_{\mathrm{1}}N_{\mathrm{2}}}\sum_{i,j} \left(\frac{1}{\Vert\mathbf{d}^{ij}\Vert}\right)^{n} \right)^{-1/n}
871 \end{equation}
872 where $\Vert\mathbf{d}^{ij}\Vert$ is the distance between atoms $i$ and $j$ in groups 1 and 2 respectively, and $n$ is an even integer.
873
874 \begin{cvcoptions}
875 \item %
876   \dupkey{group1}{\texttt{distanceInv}}{colvar|distance|group1}{\texttt{distance} component}
877 \item %
878   \simkey{group2}{\texttt{distanceInv}}{group1}
879 \item %
880   \dupkey{oneSiteTotalForce}{\texttt{distanceInv}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component}
881 \item %
882   \keydef
883     {exponent}{%
884     \texttt{distanceInv}}{%
885     Exponent $n$ in equation~\ref{eq:distanceInv}}{%
886     positive even integer}{%
887     6}{Defines the exponent to which the individual distances are elevated before averaging.  The default value of 6 is useful for example to applying restraints based on NOE-measured distances.}
888 \end{cvcoptions}
889 This component returns a number in \AA{}, ranging from $0$ to the largest possible distance within the chosen boundary conditions.
890
891
892 \cvsubsubsec{\texttt{distancePairs}: set of pairwise distances between two groups.}{sec:cvc_distancePairs}
893 \labelkey{colvar|distancePairs}
894
895 The \texttt{distancePairs~\{...\}} block defines a $N_{\mathrm{1}}\times{}N_{\mathrm{2}}$-dimensional variable that includes all mutual distances between the atoms of two groups.
896 This can be useful, for example, to develop a new variable defined over two groups, by using the \texttt{scriptedFunction} feature.
897
898 \begin{cvcoptions}
899 \item %
900   \dupkey{group1}{\texttt{distancePairs}}{colvar|distance|group1}{\texttt{distance} component}
901 \item %
902   \simkey{group2}{\texttt{distancePairs}}{group1}
903 \item %
904   \dupkey{forceNoPBC}{\texttt{distancePairs}}{colvar|distance|forceNoPBC}{\texttt{distance} component}
905 \end{cvcoptions}
906 This component returns a $N_{\mathrm{1}}\times{}N_{\mathrm{2}}$-dimensional vector of numbers, each ranging from $0$ to the largest possible distance within the chosen boundary conditions.
907
908
909 \cvsubsubsec{\texttt{cartesian}: vector of atomic Cartesian coordinates.}{sec:cvc_cartesian}
910 \labelkey{colvar|cartesian}
911
912 The \texttt{cartesian~\{...\}} block defines a component returning a flat vector containing
913 the Cartesian coordinates of all participating atoms, in the order
914 $(x_1, y_1, z_1, \cdots, x_n, y_n, z_n)$.
915
916 \begin{cvcoptions}
917 \item %
918   \key
919     {atoms}{%
920     \texttt{cartesian}}{%
921     Group of atoms}{%
922     Block \texttt{atoms \{...\}}}{%
923     Defines the atoms whose coordinates make up the value of the component.
924     If \texttt{rotateReference} or \texttt{centerReference} are defined, coordinates
925     are evaluated within the moving frame of reference.}
926 \end{cvcoptions}
927
928
929 \cvsubsubsec{\texttt{angle}: angle between three groups.}{sec:cvc_angle}
930 \labelkey{colvar|angle}
931
932 The \texttt{angle~\{...\}} block defines an angle, and contains the
933 three blocks \texttt{group1}, \texttt{group2} and \texttt{group3}, defining
934 the three groups.  It returns an angle (in degrees) within the
935 interval $[0:180]$.
936
937 \begin{cvcoptions}
938 \item %
939   \dupkey{group1}{\texttt{angle}}{colvar|distance|group1}{\texttt{distance} component}
940 \item %
941   \simkey{group2}{\texttt{angle}}{group1}
942 \item %
943   \simkey{group3}{\texttt{angle}}{group1}
944 \item %
945   \dupkey{forceNoPBC}{\texttt{angle}}{colvar|distance|forceNoPBC}{\texttt{distance} component}
946 \item %
947   \dupkey{oneSiteTotalForce}{\texttt{angle}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component}
948 \end{cvcoptions}
949
950
951
952 \cvsubsubsec{\texttt{dipoleAngle}: angle between two groups and dipole of a third group.}{sec:cvc_dipoleAngle}
953 \labelkey{colvar|dipoleAngle}
954
955 The \texttt{dipoleAngle~\{...\}} block defines an angle, and contains the
956 three blocks \texttt{group1}, \texttt{group2} and \texttt{group3}, defining
957 the three groups, being \texttt{group1} the group where dipole is calculated. 
958 It returns an angle (in degrees) within the interval $[0:180]$.
959
960 \begin{cvcoptions}
961 \item %
962   \dupkey{group1}{\texttt{dipoleAngle}}{colvar|distance|group1}{\texttt{distance} component}
963 \item %
964   \simkey{group2}{\texttt{dipoleAngle}}{group1}
965 \item %
966   \simkey{group3}{\texttt{dipoleAngle}}{group1}
967 \item %
968   \dupkey{forceNoPBC}{\texttt{dipoleAngle}}{colvar|distance|forceNoPBC}{\texttt{distance} component}
969 \item %
970   \dupkey{oneSiteTotalForce}{\texttt{dipoleAngle}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component}
971 \end{cvcoptions}
972
973
974 \cvsubsubsec{\texttt{dihedral}: torsional angle between four groups.}{sec:cvc_dihedral}
975 \labelkey{colvar|dihedral}
976
977 The \texttt{dihedral~\{...\}} block defines a torsional angle, and
978 contains the blocks \texttt{group1}, \texttt{group2}, \texttt{group3}
979 and \texttt{group4}, defining the four groups.  It returns an angle
980 (in degrees) within the interval $[-180:180]$.  The Colvars module
981 calculates all the distances between two angles taking into account
982 periodicity.  For instance, reference values for restraints or range
983 boundaries can be defined by using any real number of choice.
984
985 \begin{cvcoptions}
986 \item %
987   \dupkey{group1}{\texttt{dihedral}}{colvar|distance|group1}{\texttt{distance} component}
988 \item %
989   \simkey{group2}{\texttt{dihedral}}{group1}
990 \item %
991   \simkey{group3}{\texttt{dihedral}}{group1}
992 \item %
993   \simkey{group4}{\texttt{dihedral}}{group1}
994 \item %
995   \dupkey{forceNoPBC}{\texttt{dihedral}}{colvar|distance|forceNoPBC}{\texttt{distance} component}
996 \item %
997   \dupkey{oneSiteTotalForce}{\texttt{dihedral}}{colvar|distance|oneSiteTotalForce}{\texttt{distance} component}
998 \end{cvcoptions}
999
1000
1001 \cvsubsubsec{\texttt{polarTheta}: polar angle in spherical coordinates.}{sec:cvc_polarTheta}
1002 \labelkey{colvar|polarTheta}
1003
1004 The \texttt{polarTheta~\{...\}} block defines the polar angle in
1005 spherical coordinates, for the center of mass of a group of atoms 
1006 described by the block \texttt{atoms}.  It returns an angle
1007 (in degrees) within the interval $[0:180]$.
1008 To obtain spherical coordinates in a frame of reference tied to
1009 another group of atoms, use the \texttt{fittingGroup} (\ref{sec:colvar_atom_groups_ref_frame}) option
1010 within the \texttt{atoms} block.
1011 An example is provided in file \texttt{examples/11\_polar\_angles.in} of the Colvars public repository.
1012
1013 \begin{cvcoptions}
1014 \item %
1015   \labelkey{colvar|polarTheta|atoms}
1016   \key
1017     {atoms}{%
1018     \texttt{polarPhi}}{%
1019     Atom group}{%
1020     \texttt{atoms~\{...\}} block}{%
1021     Defines the group of atoms for the COM of which the angle should be calculated.
1022     }
1023 \end{cvcoptions}
1024
1025
1026 \cvsubsubsec{\texttt{polarPhi}: azimuthal angle in spherical coordinates.}{sec:cvc_polarPhi}
1027 \labelkey{colvar|polarPhi}
1028
1029 The \texttt{polarPhi~\{...\}} block defines the azimuthal angle in
1030 spherical coordinates, for the center of mass of a group of atoms 
1031 described by the block \texttt{atoms}. It returns an angle
1032 (in degrees) within the interval $[-180:180]$.  The Colvars module
1033 calculates all the distances between two angles taking into account
1034 periodicity.  For instance, reference values for restraints or range
1035 boundaries can be defined by using any real number of choice.
1036 To obtain spherical coordinates in a frame of reference tied to
1037 another group of atoms, use the \texttt{fittingGroup} (\ref{sec:colvar_atom_groups_ref_frame}) option
1038 within the \texttt{atoms} block.
1039 An example is provided in file \texttt{examples/11\_polar\_angles.in} of the Colvars public repository.
1040
1041
1042 \begin{cvcoptions}
1043 \item %
1044   \labelkey{colvar|polarPhi|atoms}
1045   \key
1046     {atoms}{%
1047     \texttt{polarPhi}}{%
1048     Atom group}{%
1049     \texttt{atoms~\{...\}} block}{%
1050     Defines the group of atoms for the COM of which the angle should be calculated.
1051     }
1052 \end{cvcoptions}
1053
1054
1055 \cvsubsubsec{\texttt{coordNum}: coordination number between two groups.}{sec:cvc_coordNum}
1056 \labelkey{colvar|coordNum}
1057
1058 The \texttt{coordNum \{...\}} block defines
1059 a coordination number (or number of contacts), which calculates the
1060 function $(1-(d/d_0)^{n})/(1-(d/d_0)^{m})$, where $d_0$ is the
1061 ``cutoff'' distance, and $n$ and $m$ are exponents that can control
1062 its long range behavior and stiffness \cite{Iannuzzi2003}.  This
1063 function is summed over all pairs of atoms in \texttt{group1} and
1064 \texttt{group2}:
1065 \begin{equation}
1066   \label{eq:cvc_coordNum}
1067   C (\mathtt{group1}, \mathtt{group2}) \; = \; 
1068   \sum_{i\in\mathtt{group1}}\sum_{j\in\mathtt{group2}} {
1069     \frac{1 - (|\mathbf{x}_{i}-\mathbf{x}_{j}|/d_{0})^{n}}{
1070       1 - (|\mathbf{x}_{i}-\mathbf{x}_{j}|/d_{0})^{m} }
1071   }
1072 \end{equation}
1073
1074 \begin{cvcoptions}
1075
1076 \item %
1077   \labelkey{colvar|coordNum|group1}
1078   \dupkey{group1}{\texttt{coordNum}}{colvar|distance|group1}{\texttt{distance} component}
1079
1080 \item %
1081   \labelkey{colvar|coordNum|group2}
1082   \simkey{group2}{\texttt{coordNum}}{group1}
1083
1084 \item %
1085   \labelkey{colvar|coordNum|cutoff}
1086   \keydef
1087     {cutoff}{%
1088     \texttt{coordNum}}{%
1089     ``Interaction'' distance (\AA)}{%
1090     positive decimal}{%
1091     4.0}{%
1092     This number defines the switching distance to define an
1093     interatomic contact: for $d \ll d_0$, the switching function
1094     $(1-(d/d_0)^{n})/(1-(d/d_0)^{m})$ is close to 1, at $d = d_0$ it
1095     has a value of $n/m$ ($1/2$ with the default $n$ and $m$), and at
1096     $d \gg d_0$ it goes to zero approximately like $d^{m-n}$.  Hence,
1097     for a proper behavior, $m$ must be larger than $n$.}
1098
1099 \item %
1100   \labelkey{colvar|coordNum|cutoff3}
1101   \keydef
1102     {cutoff3}{%
1103     \texttt{coordNum}}{%
1104     Reference distance vector (\AA)}{%
1105     ``\texttt{(x, y, z)}'' triplet of positive decimals}{%
1106     \texttt{(4.0, 4.0, 4.0)}}{%
1107     The three components of this vector define three different cutoffs
1108     $d_{0}$ for each direction.  This option is mutually exclusive with
1109     \texttt{cutoff}.}
1110
1111 \item %
1112   \labelkey{colvar|coordNum|expNumer}
1113   \keydef
1114     {expNumer}{%
1115     \texttt{coordNum}}{%
1116     Numerator exponent}{%
1117     positive even integer}{%
1118     6}{%
1119     This number defines the $n$ exponent for the switching function.}
1120
1121 \item %
1122   \labelkey{colvar|coordNum|expDenom}
1123   \keydef
1124     {expDenom}{%
1125     \texttt{coordNum}}{%
1126     Denominator exponent}{%
1127     positive even integer}{%
1128     12}{%
1129     This number defines the $m$ exponent for the switching function.}
1130
1131 \item %
1132   \labelkey{colvar|coordNum|group2CenterOnly}
1133   \keydef
1134     {group2CenterOnly}{%
1135     \texttt{coordNum}}{%
1136     Use only \texttt{group2}'s center of
1137     mass}{%
1138     boolean}{%
1139     \texttt{off}}{%
1140     If this option is \texttt{on}, only contacts between each atoms in \texttt{group1} and the center of mass of     \texttt{group2} are calculated (by default, the sum extends over all pairs of atoms in \texttt{group1} and \texttt{group2}).
1141 If \texttt{group2} is a \texttt{dummyAtom}, this option is set to \texttt{yes} by default.
1142 }
1143
1144 \item %
1145     \labelkey{colvar|coordNum|tolerance}
1146     \keydef
1147      {tolerance}{%
1148      \texttt{coordNum}}{%
1149      Pairlist control}{%
1150     decimal}{%
1151     0.0}{This controls the pairlist feature, dictating the minimum value for each summation element in Eq.~\ref{eq:cvc_coordNum} such that the pair that contributed the summation element is included in subsequent simulation timesteps until the next pairlist recalculation. For most applications, this value should be small (eg. 0.001) to avoid missing important contributions to the overall sum. Higher values will improve performance by reducing the number of pairs that contribute to the sum. Values above 1 will exclude all possible pair interactions. Similarly, values below 0 will never exclude a pair from consideration. To ensure continuous forces, Eq.~\ref{eq:cvc_coordNum} is further modified by subtracting the tolerance and then rescaling so that each pair covers the range $\left[0, 1\right]$.
1152   }
1153
1154 \item %
1155     \labelkey{colvar|coordNum|pairListFrequency}
1156     \keydef
1157      {pairListFrequency}{%
1158      \texttt{coordNum}}{%
1159      Pairlist regeneration frequency}{%
1160     positive integer}{%
1161     100}{This controls the pairlist feature, dictating how many steps are taken between regenerating pairlists if the tolerance is greater than 0.
1162   }
1163 \end{cvcoptions}
1164
1165 This component returns a dimensionless number, which ranges from
1166 approximately 0 (all interatomic distances are much larger than the
1167 cutoff) to $N_{\mathtt{group1}} \times N_{\mathtt{group2}}$ (all distances
1168 are less than the cutoff), or $N_{\mathtt{group1}}$ if
1169 \texttt{group2CenterOnly} is used.  For performance reasons, at least
1170 one of \texttt{group1} and \texttt{group2} should be of limited size or \texttt{group2CenterOnly} should be used: the cost of the loop over all pairs grows as $N_{\mathtt{group1}} \times N_{\mathtt{group2}}$.
1171 Setting $\mathtt{tolerance} > 0$ ameliorates this to some degree, although every pair is still checked to regenerate the pairlist.
1172
1173
1174
1175 \cvsubsubsec{\texttt{selfCoordNum}: coordination number between atoms within a group.}{sec:cvc_selfCoordNum}
1176 \labelkey{colvar|selfCoordNum}
1177
1178 The \texttt{selfCoordNum \{...\}} block defines
1179 a coordination number similarly to the component \texttt{coordNum},
1180 but the function is summed over atom pairs within \texttt{group1}:
1181 \begin{equation}
1182   \label{eq:cvc_selfCoordNum}
1183   C (\mathtt{group1}) \; = \; 
1184   \sum_{i\in\mathtt{group1}}\sum_{j > i} {
1185     \frac{1 - (|\mathbf{x}_{i}-\mathbf{x}_{j}|/d_{0})^{n}}{
1186       1 - (|\mathbf{x}_{i}-\mathbf{x}_{j}|/d_{0})^{m} }
1187   }
1188 \end{equation}
1189 The keywords accepted by \texttt{selfCoordNum} are a subset of
1190 those accepted by \texttt{coordNum}, namely \texttt{group1}
1191 (here defining \emph{all} of the atoms to be considered),
1192 \texttt{cutoff}, \texttt{expNumer}, and \texttt{expDenom}.
1193
1194 \begin{cvcoptions}
1195 \item %
1196   \dupkey{group1}{\texttt{selfCoordNum}}{colvar|coordNum|group1}{\texttt{coordNum} component}
1197 \item %
1198   \dupkey{cutoff}{\texttt{selfCoordNum}}{colvar|coordNum|cutoff}{\texttt{coordNum} component}
1199 \item %
1200   \dupkey{cutoff3}{\texttt{selfCoordNum}}{colvar|coordNum|cutoff3}{\texttt{coordNum} component}
1201 \item %
1202   \dupkey{expNumer}{\texttt{selfCoordNum}}{colvar|coordNum|expNumer}{\texttt{coordNum} component}
1203 \item %
1204   \dupkey{expDenom}{\texttt{selfCoordNum}}{colvar|coordNum|expDenom}{\texttt{coordNum} component}
1205 \item %
1206   \dupkey{tolerance}{\texttt{selfCoordNum}}{colvar|coordNum|tolerance}{\texttt{coordNum} component}
1207 \item %
1208   \dupkey{pairListFrequency}{\texttt{selfCoordNum}}{colvar|coordNum|pairListFrequency}{\texttt{coordNum} component}
1209 \end{cvcoptions}
1210
1211 This component returns a dimensionless number, which ranges from
1212 approximately 0 (all interatomic distances much larger than the
1213 cutoff) to $N_{\mathtt{group1}} \times (N_{\mathtt{group1}} - 1) / 2$ (all
1214 distances within the cutoff).  For performance reasons,
1215 \texttt{group1} should be of limited size, because the cost of the
1216 loop over all pairs grows as $N_{\mathtt{group1}}^2$.
1217
1218
1219
1220 \cvsubsubsec{\texttt{hBond}: hydrogen bond between two atoms.}{sec:cvc_hBond}
1221 \labelkey{colvar|hBond}
1222
1223 The \texttt{hBond \{...\}} block defines a hydrogen
1224 bond, implemented as a coordination number (eq.~\ref{eq:cvc_coordNum})
1225 between the donor and the acceptor atoms.  Therefore, it accepts the
1226 same options \texttt{cutoff} (with a different default value of
1227 3.3~\AA{}), \texttt{expNumer} (with a default value of 6) and
1228 \texttt{expDenom} (with a default value of 8).  Unlike
1229 \texttt{coordNum}, it requires two atom numbers, \texttt{acceptor} and
1230 \texttt{donor}, to be defined.  It returns an adimensional number,
1231 with values between 0 (acceptor and donor far outside the cutoff
1232 distance) and 1 (acceptor and donor much closer than the cutoff).
1233
1234 \begin{cvcoptions}
1235 \item %
1236   \key
1237     {acceptor}{%
1238     \texttt{hBond}}{%
1239     Number of the acceptor atom}{%
1240     positive integer}{%
1241     Number that uses the same convention as \texttt{atomNumbers}.}
1242 \item %
1243   \simkey{donor}{\texttt{hBond}}{acceptor}
1244 \item %
1245   \dupkey{cutoff}{\texttt{hBond}}{colvar|coordNum|cutoff}{\texttt{coordNum} component}\\
1246   \textbf{Note:} default value is 3.3~\AA.
1247 \item %
1248   \dupkey{expNumer}{\texttt{hBond}}{colvar|coordNum|expNumer}{\texttt{coordNum} component}\\
1249   \textbf{Note:} default value is 6.
1250 \item %
1251   \dupkey{expDenom}{\texttt{hBond}}{colvar|coordNum|expDenom}{\texttt{coordNum} component}\\
1252   \textbf{Note:} default value is 8.
1253 \end{cvcoptions}
1254
1255
1256 \cvsubsubsec{\texttt{rmsd}: root mean square displacement (RMSD) from reference positions.}{sec:cvc_rmsd}
1257 \labelkey{colvar|rmsd}
1258
1259 The block \texttt{rmsd~\{...\}} defines the root mean square replacement
1260 (RMSD) of a group of atoms with respect to a reference structure.  For
1261 each set of coordinates $\{ \mathbf{x}_1(t), \mathbf{x}_2(t), \ldots
1262 \mathbf{x}_N(t) \}$, the colvar component \texttt{rmsd} calculates the
1263 optimal rotation
1264 $U^{\{\mathbf{x}_{i}(t)\}\rightarrow\{\mathbf{x}_{i}^{\mathrm{(ref)}}\}}$
1265 that best superimposes the coordinates $\{\mathbf{x}_{i}(t)\}$ onto a
1266 set of reference coordinates $\{\mathbf{x}_{i}^{\mathrm{(ref)}}\}$.
1267 Both the current and the reference coordinates are centered on their
1268 centers of geometry, $\mathbf{x}_{\mathrm{cog}}(t)$ and
1269 $\mathbf{x}_{\mathrm{cog}}^{\mathrm{(ref)}}$.  The root mean square
1270 displacement is then defined as:
1271 \begin{equation}
1272   \label{eq:cvc_rmsd}
1273   { \mathrm{RMSD}(\{\mathbf{x}_{i}(t)\},
1274     \{\mathbf{x}_{i}^{\mathrm{(ref)}}\}) } \; = \; \sqrt{
1275     \frac{1}{N} \sum_{i=1}^{N} \left|
1276       U
1277       \left(\mathbf{x}_{i}(t) - \mathbf{x}_{\mathrm{cog}}(t)\right) -
1278       \left(\mathbf{x}_{i}^{\mathrm{(ref)}} -
1279         \mathbf{x}_{\mathrm{cog}}^{\mathrm{(ref)}} \right) \right|^{2} }
1280 \end{equation}
1281 The optimal rotation
1282 $U^{\{\mathbf{x}_{i}(t)\}\rightarrow\{\mathbf{x}_{i}^{\mathrm{(ref)}}\}}$
1283 is calculated within the formalism developed in
1284 reference~\cite{Coutsias2004}, which guarantees a continuous
1285 dependence of
1286 $U^{\{\mathbf{x}_{i}(t)\}\rightarrow\{\mathbf{x}_{i}^{\mathrm{(ref)}}\}}$
1287 with respect to $\{\mathbf{x}_{i}(t)\}$.
1288
1289 \begin{cvcoptions}
1290
1291 \item %
1292   \labelkey{colvar|rmsd|atoms}
1293   \key
1294     {atoms}{%
1295     \texttt{rmsd}}{%
1296     Atom group}{%
1297     \texttt{atoms~\{...\}} block}{%
1298     Defines the group of atoms of which the RMSD should be calculated.
1299     Optimal fit options (such as \texttt{refPositions} and
1300     \texttt{rotateReference}) should typically NOT be set within this
1301     block. Exceptions to this rule are the special cases discussed in
1302     the \emph{Advanced usage} paragraph below.
1303     }
1304
1305 \item %
1306   \labelkey{colvar|rmsd|refPositions}
1307   \key
1308     {refPositions}{%
1309     \texttt{rmsd}}{%
1310     Reference coordinates}{%
1311     space-separated list of \texttt{(x, y, z)} triplets}{%
1312     This option (mutually exclusive with \texttt{refPositionsFile}) sets the reference coordinates for RMSD calculation, and uses these to compute the roto-translational fit.  
1313     It is functionally equivalent to the option \refkey{refPositions}{atom-group|refPositions} in the atom group definition, which also supports more advanced fitting options.
1314     }
1315
1316 \item %
1317   \labelkey{colvar|rmsd|refPositionsFile}
1318   \key
1319     {refPositionsFile}{%
1320     \texttt{rmsd}}{%
1321     Reference coordinates file}{%
1322     UNIX filename}{%
1323     This option (mutually exclusive with \texttt{refPositions}) sets the reference coordinates for RMSD calculation, and uses these to compute the roto-translational fit.  
1324     It is functionally equivalent to the option \refkey{refPositionsFile}{atom-group|refPositionsFile} in the atom group definition, which also supports more advanced fitting options.
1325     }
1326
1327 \cvnamebasedonly{
1328 \item %
1329   \labelkey{colvar|rmsd|refPositionsCol}
1330   \key
1331     {refPositionsCol}{%
1332     \texttt{rmsd}}{%
1333     PDB column containing atom flags}{%
1334     \texttt{O}, \texttt{B}, \texttt{X}, \texttt{Y}, or \texttt{Z}}{%
1335     If \texttt{refPositionsFile} is a PDB file that contains all the atoms in the topology, this option may be provided to set which PDB field is used to flag the reference coordinates for \texttt{atoms}.
1336   }
1337
1338 \item %
1339   \labelkey{colvar|rmsd|refPositionsColValue}
1340   \key
1341     {refPositionsColValue}{%
1342     \texttt{rmsd}}{%
1343     Atom selection flag in the PDB column}{%
1344     positive decimal}{%
1345     If defined, this value identifies in the PDB column
1346     \texttt{refPositionsCol} of the file \texttt{refPositionsFile}
1347     which atom positions are to be read.  Otherwise, all positions
1348     with a non-zero value are read.
1349   }
1350 }
1351 \end{cvcoptions}
1352 This component returns a positive real number (in \AA).
1353
1354
1355 \cvsubsubsec{Advanced usage of the \texttt{rmsd} component.}{sec:cvc_rmsd_advanced}
1356 In the standard usage as described above, the \texttt{rmsd} component
1357 calculates a minimum RMSD, that is, current coordinates are optimally
1358 fitted onto the same reference coordinates that are used to 
1359 compute the RMSD value. The fit itself is handled by the atom group
1360 object, whose parameters are automatically set by the \texttt{rmsd}
1361 component.
1362 For very specific applications, however, it may be
1363 useful to control the fitting process separately from the definition
1364 of the reference coordinates, to evaluate various types of
1365 non-minimal RMSD values. This can be achieved by setting the
1366 related options (\texttt{refPositions}, etc.) explicitly in the
1367 atom group block. This allows for the following non-standard cases:
1368
1369 \begin{enumerate}
1370 \item applying the optimal translation, but no rotation
1371 (\texttt{rotateReference off}), to bias or restrain the shape and
1372 orientation, but not the position of the atom group;
1373 \item applying the optimal rotation, but no translation
1374 (\texttt{translateReference off}), to bias or restrain the shape and
1375 position, but not the orientation of the atom group;
1376 \item disabling the application of optimal roto-translations, which
1377 lets the RMSD component decribe the deviation of atoms
1378 from fixed positions in the laboratory frame: this allows for custom
1379 positional restraints within the Colvars module;
1380 \item fitting the atomic positions to different reference coordinates
1381 than those used in the RMSD calculation itself;
1382 \item applying the optimal rotation and/or translation from a separate
1383 atom group, defined through \texttt{fittingGroup}: the RMSD then
1384 reflects the deviation from reference coordinates in a separate, moving
1385 reference frame.
1386 \end{enumerate}
1387
1388
1389 \cvscriptonly{
1390 \cvsubsubsec{Path collective variables}{sec:pathcv}
1391
1392 An application of the \texttt{rmsd} component is "path collective variables",\cite{Branduardi2007}
1393 which are implemented as Tcl-scripted combinations or RMSDs.
1394 The implementation is available as file \texttt{colvartools/pathCV.tcl}, and
1395 an example is provided in file \texttt{examples/10\_pathCV.namd} of the Colvars public repository.}
1396
1397
1398 \cvsubsubsec{\texttt{eigenvector}: projection of the atomic  coordinates on a vector.}{sec:cvc_eigenvector}
1399 \labelkey{colvar|eigenvector}
1400
1401 The block \texttt{eigenvector~\{...\}} defines the projection of the coordinates
1402 of a group of atoms (or more precisely, their deviations from the
1403 reference coordinates) onto a vector in $\mathbb{R}^{3n}$, where $n$ is the
1404 number of atoms in the group. The computed quantity is the
1405 total projection:
1406 \begin{equation}
1407   \label{eq:cvc_eigenvector}
1408   { p(\{\mathbf{x}_{i}(t)\},
1409     \{\mathbf{x}_{i}^{\mathrm{(ref)}}\}) } \; = \; {
1410     \sum_{i=1}^{n}  \mathbf{v}_{i} \cdot
1411     \left(U(\mathbf{x}_{i}(t) - \mathbf{x}_{\mathrm{cog}}(t)) -
1412       (\mathbf{x}_{i}^{\mathrm{(ref)}} -
1413       \mathbf{x}_{\mathrm{cog}}^{\mathrm{(ref)}}) \right)\mathrm{,} }
1414 \end{equation}
1415 where, as in the \texttt{rmsd} component, $U$ is the optimal rotation
1416 matrix, $\mathbf{x}_{\mathrm{cog}}(t)$ and
1417 $\mathbf{x}_{\mathrm{cog}}^{\mathrm{(ref)}}$ are the centers of
1418 geometry of the current and reference positions respectively, and
1419 $\mathbf{v}_{i}$ are the components of the vector for each atom.
1420 Example choices for $(\mathbf{v}_{i})$ are an eigenvector
1421 of the covariance matrix (essential mode), or a normal
1422 mode of the system.  It is assumed that $\sum_{i}\mathbf{v}_{i} = 0$:
1423 otherwise, the Colvars module centers the $\mathbf{v}_{i}$
1424 automatically when reading them from the configuration.
1425
1426 \begin{cvcoptions}
1427 \item %
1428   \dupkey{atoms}{\texttt{eigenvector}}{colvar|rmsd|atoms}{\texttt{rmsd} component}
1429 \item %
1430   \dupkey{refPositions}{\texttt{eigenvector}}{colvar|rmsd|refPositions}{\texttt{rmsd} component}
1431 \item %
1432   \dupkey{refPositionsFile}{\texttt{eigenvector}}{colvar|rmsd|refPositionsFile}{\texttt{rmsd} component}
1433 \cvnamebasedonly{
1434 \item %
1435   \dupkey{refPositionsCol}{\texttt{eigenvector}}{colvar|rmsd|refPositionsCol}{\texttt{rmsd} component}
1436 \item %
1437   \dupkey{refPositionsColValue}{\texttt{eigenvector}}{colvar|rmsd|refPositionsColValue}{\texttt{rmsd} component}
1438 }
1439
1440 \item %
1441   \key
1442     {vector}{%
1443     \texttt{eigenvector}}{%
1444     Vector components}{%
1445     space-separated list of \texttt{(x, y, z)} triplets}{%
1446     This option (mutually exclusive with \texttt{vectorFile}) sets the values of the vector components.}
1447
1448 \item %
1449   \key
1450     {vectorFile}{%
1451     \texttt{eigenvector}}{%
1452     file containing vector components}{%
1453     UNIX filename}{%
1454     This option (mutually exclusive with \texttt{vector}) sets the name of a coordinate file containing the vector components; the file is read according to the same format used for \texttt{refPositionsFile}.
1455     \cvnamebasedonly{For a PDB file specifically, the components are read from the X, Y and Z fields.
1456       \textbf{Note:} \emph{The PDB file has limited precision and fixed-point numbers: in some cases, the vector components may not be accurately represented; a XYZ file should be used instead, containing floating-point numbers.}}
1457   }
1458
1459 \cvnamebasedonly{
1460 \item %
1461   \key
1462     {vectorCol}{%
1463     \texttt{eigenvector}}{%
1464     PDB column used to flag participating atoms}{%
1465     \texttt{O} or \texttt{B}}{%
1466     Analogous to \texttt{atomsCol}.}
1467
1468 \item %
1469   \key
1470     {vectorColValue}{%
1471     \texttt{eigenvector}}{%
1472     Value used to flag participating atoms in the PDB file}{%
1473     positive decimal}{%
1474     Analogous to \texttt{atomsColValue}.}
1475 }
1476
1477 \item %
1478   \keydef
1479     {differenceVector}{%
1480     \texttt{eigenvector}}{%
1481     The $3n$-dimensional vector is the difference between \texttt{vector} and \texttt{refPositions}}{%
1482     boolean}{%
1483     \texttt{off}}{%
1484     If this option is \texttt{on}, the numbers provided by \texttt{vector}\cvnamebasedonly{ or \texttt{vectorFile}} are interpreted as another set of positions, $\mathbf{x}_{i}'$: the vector $\mathbf{v}_{i}$ is then defined as $\mathbf{v}_{i} = \left(\mathbf{x}_{i}' - \mathbf{x}_{i}^{\mathrm{(ref)}}\right)$.
1485 This allows to conveniently define a colvar $\xi$ as a projection on the linear transformation between two sets of positions, ``A'' and ``B''.
1486 For convenience, the vector is also normalized so that $\xi = 0$ when the atoms are at the set of positions ``A'' and $\xi = 1$ at the set of positions ``B''.
1487 }
1488 \end{cvcoptions}
1489 This component returns a number (in \AA), whose value ranges between
1490 the smallest and largest absolute positions in the unit cell during
1491 the simulations (see also \texttt{distanceZ}).  Due to the
1492 normalization in eq.~\ref{eq:cvc_eigenvector}, this range does not
1493 depend on the number of atoms involved.
1494
1495
1496 \cvsubsubsec{\texttt{gyration}: radius of gyration of a group  of atoms.}{sec:cvc_gyration}
1497 \labelkey{colvar|gyration}
1498
1499 The block \texttt{gyration~\{...\}} defines the
1500 parameters for calculating the radius of gyration of a group of atomic
1501 positions $\{ \mathbf{x}_1(t), \mathbf{x}_2(t), \ldots \mathbf{x}_N(t)
1502 \}$ with respect to their center of geometry,
1503 $\mathbf{x}_{\mathrm{cog}}(t)$:
1504 \begin{equation}
1505   \label{eq:colvar_gyration}
1506   R_{\mathrm{gyr}} \; = \; \sqrt{ \frac{1}{N}
1507     \sum_{i=1}^{N} \left|\mathbf{x}_{i}(t) -
1508       \mathbf{x}_{\mathrm{cog}}(t)\right|^{2} }
1509 \end{equation}
1510 This component must contain one \texttt{atoms~\{...\}} block to
1511 define the atom group, and returns a positive number, expressed in
1512 \AA{}.
1513
1514 \begin{cvcoptions}
1515 \item %
1516   \dupkey{atoms}{\texttt{gyration}}{colvar|rmsd|atoms}{\texttt{rmsd} component}
1517 \end{cvcoptions}
1518
1519
1520 \cvsubsubsec{\texttt{inertia}: total moment of inertia of a group  of atoms.}{sec:cvc_inertia}
1521 \labelkey{colvar|inertia}
1522
1523 The block \texttt{inertia~\{...\}} defines the
1524 parameters for calculating the total moment of inertia of a group of atomic
1525 positions $\{ \mathbf{x}_1(t), \mathbf{x}_2(t), \ldots \mathbf{x}_N(t)
1526 \}$ with respect to their center of geometry,
1527 $\mathbf{x}_{\mathrm{cog}}(t)$:
1528 \begin{equation}
1529   \label{eq:colvar_inertia}
1530   I \; = \; \sum_{i=1}^{N} \left|\mathbf{x}_{i}(t) -
1531       \mathbf{x}_{\mathrm{cog}}(t)\right|^{2}
1532 \end{equation}
1533 \emph{Note that all atomic masses are set to 1 for simplicity.}
1534 This component must contain one \texttt{atoms~\{...\}} block to
1535 define the atom group, and returns a positive number, expressed in
1536 \AA{}$^{2}$.
1537
1538 \begin{cvcoptions}
1539 \item %
1540   \dupkey{atoms}{\texttt{inertia}}{colvar|rmsd|atoms}{\texttt{rmsd} component}
1541 \end{cvcoptions}
1542
1543
1544 \cvsubsubsec{\texttt{inertiaZ}: total moment of inertia of a group of atoms around a chosen axis.}{sec:cvc_inertiaZ}
1545 \labelkey{colvar|inertiaZ}
1546
1547 The block \texttt{inertiaZ~\{...\}} defines the
1548 parameters for calculating the component along the axis $\mathbf{e}$ of the moment of inertia of a group of atomic
1549 positions $\{ \mathbf{x}_1(t), \mathbf{x}_2(t), \ldots \mathbf{x}_N(t)
1550 \}$ with respect to their center of geometry,
1551 $\mathbf{x}_{\mathrm{cog}}(t)$:
1552 \begin{equation}
1553   \label{eq:colvar_inertia_z}
1554   I_{\mathbf{e}} \; = \; \sum_{i=1}^{N} \left(\left(\mathbf{x}_{i}(t) -
1555       \mathbf{x}_{\mathrm{cog}}(t)\right)\cdot\mathbf{e}\right)^{2}
1556 \end{equation}
1557 \emph{Note that all atomic masses are set to 1 for simplicity.}
1558 This component must contain one \texttt{atoms~\{...\}} block to
1559 define the atom group, and returns a positive number, expressed in
1560 \AA{}$^{2}$.
1561
1562 \begin{cvcoptions}
1563 \item %
1564   \dupkey{atoms}{\texttt{inertiaZ}}{colvar|rmsd|atoms}{\texttt{rmsd} component}
1565 \item %
1566   \keydef
1567     {axis}{%
1568     \texttt{inertiaZ}}{%
1569     Projection axis (\AA{})}{%
1570     \texttt{(x, y, z)} triplet}{%
1571     \texttt{(0.0, 0.0, 1.0)}}{%
1572     The three components of this vector define (when normalized) the
1573     projection axis $\mathbf{e}$.}
1574 \end{cvcoptions}
1575
1576
1577 \cvsubsubsec{\texttt{orientation}: orientation from reference coordinates.}{sec:cvc_orientation}
1578 \labelkey{colvar|orientation}
1579
1580 The block \texttt{orientation~\{...\}} returns the
1581 same optimal rotation used in the \texttt{rmsd} component to
1582 superimpose the coordinates $\{\mathbf{x}_{i}(t)\}$ onto a set of
1583 reference coordinates $\{\mathbf{x}_{i}^{\mathrm{(ref)}}\}$.  Such
1584 component returns a four dimensional vector $\mathsf{q} = (q_0, q_1,
1585 q_2, q_3)$, with $\sum_{i} q_{i}^{2} = 1$; this \emph{quaternion}
1586 expresses the optimal rotation $\{\mathbf{x}_{i}(t)\} \rightarrow
1587 \{\mathbf{x}_{i}^{\mathrm{(ref)}}\}$ according to the formalism in
1588 reference~\cite{Coutsias2004}.  The quaternion $(q_0, q_1, q_2, q_3)$
1589 can also be written as $\left(\cos(\theta/2), \,
1590   \sin(\theta/2)\mathbf{u}\right)$, where $\theta$ is the angle and
1591 $\mathbf{u}$ the normalized axis of rotation; for example, a rotation
1592 of 90$^{\circ}$ around the $z$ axis is expressed as
1593 ``\texttt{(0.707, 0.0, 0.0, 0.707)}''.  The script
1594 \texttt{quaternion2rmatrix.tcl} provides Tcl functions for converting
1595 to and from a $4\times{}4$ rotation matrix in a format suitable for
1596 usage in VMD.
1597
1598 As for the component \texttt{rmsd}, the available options are \texttt{atoms}, \texttt{refPositionsFile}\cvnamebasedonly{, \texttt{refPositionsCol} and \texttt{refPositionsColValue}, } and \texttt{refPositions}.
1599
1600 \textbf{Note:} \texttt{refPositions}and \texttt{refPositionsFile} define the set of positions \emph{from which} the optimal rotation is calculated, but this rotation is not applied to the coordinates of the atoms involved: it is used instead to define the variable itself.
1601
1602 \begin{cvcoptions}
1603 \item %
1604   \dupkey{atoms}{\texttt{orientation}}{colvar|rmsd|atoms}{\texttt{rmsd} component}
1605 \item %
1606   \dupkey{refPositions}{\texttt{orientation}}{colvar|rmsd|refPositions}{\texttt{rmsd} component}
1607 \item %
1608   \dupkey{refPositionsFile}{\texttt{orientation}}{colvar|rmsd|refPositionsFile}{\texttt{rmsd} component}
1609
1610 \cvnamebasedonly{
1611 \item %
1612   \dupkey{refPositionsCol}{\texttt{orientation}}{colvar|rmsd|refPositionsCol}{\texttt{rmsd} component}
1613 \item %
1614   \dupkey{refPositionsColValue}{\texttt{orientation}}{colvar|rmsd|refPositionsColValue}{\texttt{rmsd} component}
1615 }
1616
1617 \item %
1618   \keydef
1619     {closestToQuaternion}{%
1620     \texttt{orientation}}{%
1621     Reference rotation}{%
1622     ``\texttt{(q0, q1, q2, q3)}'' quadruplet}{%
1623     \texttt{(1.0, 0.0, 0.0, 0.0)} (``null'' rotation)}{%
1624     Between the two equivalent quaternions $(q_0, q_1, q_2, q_3)$ and
1625     $(-q_0, -q_1, -q_2, -q_3)$, the closer to \texttt{(1.0, 0.0, 0.0,
1626       0.0)} is chosen.  This simplifies the visualization of the
1627     colvar trajectory when sampled values are a smaller subset of all
1628     possible rotations.  \textbf{Note:} \emph {this only affects the
1629       output, never the dynamics}.}
1630
1631 \end{cvcoptions}
1632
1633 \textbf{Tip: stopping the rotation of a protein.}  To stop the
1634 rotation of an elongated macromolecule in solution (and use an
1635 anisotropic box to save water molecules), it is possible to define a
1636 colvar with an \texttt{orientation} component, and restrain it throuh
1637 the \texttt{harmonic} bias around the identity rotation, \texttt{(1.0,
1638   0.0, 0.0, 0.0)}.  Only the overall orientation of the macromolecule
1639 is affected, and \emph{not} its internal degrees of freedom.  The user
1640 should also take care that the macromolecule is composed by a single
1641 chain, or disable \texttt{wrapAll} otherwise.
1642
1643
1644
1645 \cvsubsubsec{\texttt{orientationAngle}: angle of rotation from reference coordinates.}{sec:cvc_orientationAngle}
1646 \labelkey{colvar|orientationAngle}
1647
1648 The block \texttt{orientationAngle~\{...\}} accepts the same base options as
1649 the component \texttt{orientation}: \texttt{atoms}, \texttt{refPositions}, \texttt{refPositionsFile}\cvnamebasedonly{, \texttt{refPositionsCol} and \texttt{refPositionsColValue}}.
1650 The returned value is the angle of rotation $\theta$ between the current and the reference positions.
1651 This angle is expressed in degrees within the range [0$^{\circ}$:180$^{\circ}$].
1652
1653 \begin{cvcoptions}
1654 \item %
1655   \dupkey{atoms}{\texttt{orientationAngle}}{colvar|rmsd|atoms}{\texttt{rmsd} component}
1656 \item %
1657   \dupkey{refPositions}{\texttt{orientationAngle}}{colvar|rmsd|refPositions}{\texttt{rmsd} component}
1658 \item %
1659   \dupkey{refPositionsFile}{\texttt{orientationAngle}}{colvar|rmsd|refPositionsFile}{\texttt{rmsd} component}
1660
1661 \cvnamebasedonly{
1662 \item %
1663   \dupkey{refPositionsCol}{\texttt{orientationAngle}}{colvar|rmsd|refPositionsCol}{\texttt{rmsd} component}
1664 \item %
1665   \dupkey{refPositionsColValue}{\texttt{orientationAngle}}{colvar|rmsd|refPositionsColValue}{\texttt{rmsd} component}
1666 }
1667 \end{cvcoptions}
1668
1669
1670 \cvsubsubsec{\texttt{orientationProj}: cosine of the angle of rotation from reference coordinates.}  {sec:cvc_orientationProj}
1671 \labelkey{colvar|orientationProj}
1672
1673 The block \texttt{orientationProj~\{...\}} accepts the same base options as
1674 the component \texttt{orientation}: \texttt{atoms}, \texttt{refPositions}, \texttt{refPositionsFile}\cvnamebasedonly{, \texttt{refPositionsCol} and \texttt{refPositionsColValue}}.
1675 The returned value is the cosine of the angle of rotation $\theta$ between the current and the reference positions.
1676 The range of values is [-1:1].
1677
1678 \begin{cvcoptions}
1679 \item %
1680   \dupkey{atoms}{\texttt{orientationProj}}{colvar|rmsd|atoms}{\texttt{rmsd} component}
1681 \item %
1682   \dupkey{refPositions}{\texttt{orientationProj}}{colvar|rmsd|refPositions}{\texttt{rmsd} component}
1683 \item %
1684   \dupkey{refPositionsFile}{\texttt{orientationProj}}{colvar|rmsd|refPositionsFile}{\texttt{rmsd} component}
1685 \cvnamebasedonly{
1686 \item %
1687   \dupkey{refPositionsCol}{\texttt{orientationProj}}{colvar|rmsd|refPositionsCol}{\texttt{rmsd} component}
1688 \item %
1689   \dupkey{refPositionsColValue}{\texttt{orientationProj}}{colvar|rmsd|refPositionsColValue}{\texttt{rmsd} component}
1690 }
1691 \end{cvcoptions}
1692
1693
1694 \cvsubsubsec{\texttt{spinAngle}: angle of rotation around a given axis.}{sec:cvc_spinAngle}
1695 \labelkey{colvar|spinAngle}
1696
1697 The complete rotation described by \texttt{orientation} can optionally be decomposed into two sub-rotations: one is a ``\emph{spin}'' rotation around \textbf{e}, and the other a ``\emph{tilt}'' rotation around an axis orthogonal to \textbf{e}.
1698 The component \texttt{spinAngle} measures the angle of the ``spin'' sub-rotation around \textbf{e}.
1699
1700 \begin{cvcoptions}
1701 \item %
1702   \dupkey{atoms}{\texttt{spinAngle}}{colvar|rmsd|atoms}{\texttt{rmsd} component}
1703 \item %
1704   \dupkey{refPositions}{\texttt{spinAngle}}{colvar|rmsd|refPositions}{\texttt{rmsd} component}
1705 \item %
1706   \dupkey{refPositionsFile}{\texttt{spinAngle}}{colvar|rmsd|refPositionsFile}{\texttt{rmsd} component}
1707 \cvnamebasedonly{
1708 \item %
1709   \dupkey{refPositionsCol}{\texttt{spinAngle}}{colvar|rmsd|refPositionsCol}{\texttt{rmsd} component}
1710 \item %
1711   \dupkey{refPositionsColValue}{\texttt{spinAngle}}{colvar|rmsd|refPositionsColValue}{\texttt{rmsd} component}
1712 }
1713 \item %
1714   \labelkey{colvar|spinAngle|axis}
1715   \keydef
1716     {axis}{%
1717     \texttt{tilt}}{%
1718     Special rotation axis (\AA{})}{%
1719     \texttt{(x, y, z)} triplet}{%
1720     \texttt{(0.0, 0.0, 1.0)}}{%
1721     The three components of this vector define (when normalized) the special rotation axis used to calculate the \texttt{tilt} and \texttt{spinAngle} components.}
1722 \end{cvcoptions}
1723 The component \texttt{spinAngle} returns an angle (in degrees) within the periodic interval $[-180:180]$.  
1724
1725 \textbf{Note:} the value of \texttt{spinAngle} is a continuous function almost everywhere, with the exception of configurations with the corresponding ``tilt'' angle equal to 180$^\circ$ (i.e.~the \texttt{tilt} component is equal to $-1$): in those cases, \texttt{spinAngle} is undefined.  If such configurations are expected, consider defining a \texttt{tilt} colvar using the same axis \textbf{e}, and restraining it with a lower wall away from $-1$.
1726
1727
1728 \cvsubsubsec{\texttt{tilt}: cosine of the rotation orthogonal to a given axis.}{sec:cvc_tilt}
1729 \labelkey{colvar|tilt}
1730
1731 The component \texttt{tilt} measures the cosine of the angle of the ``tilt'' sub-rotation, which combined with the ``spin'' sub-rotation provides the complete rotation of a group of atoms.
1732 The cosine of the tilt angle rather than the tilt angle itself is implemented, because the latter is unevenly distributed even for an isotropic system: consider as an analogy the angle $\theta$ in the spherical coordinate system.
1733 The component \texttt{tilt} relies on the same options as \texttt{spinAngle}, including the definition of the axis \textbf{e}.
1734 The values of \texttt{tilt} are real numbers in the interval $[-1:1]$: the value $1$ represents an orientation fully parallel to \textbf{e} (tilt angle = 0$^\circ$), and the value $-1$ represents an anti-parallel orientation.
1735
1736 \begin{cvcoptions}
1737 \item %
1738   \dupkey{atoms}{\texttt{tilt}}{colvar|rmsd|atoms}{\texttt{rmsd} component}
1739 \item %
1740   \dupkey{refPositions}{\texttt{tilt}}{colvar|rmsd|refPositions}{\texttt{rmsd} component}
1741 \item %
1742   \dupkey{refPositionsFile}{\texttt{tilt}}{colvar|rmsd|refPositionsFile}{\texttt{rmsd} component}
1743 \cvnamebasedonly{
1744 \item %
1745   \dupkey{refPositionsCol}{\texttt{tilt}}{colvar|rmsd|refPositionsCol}{\texttt{rmsd} component}
1746 \item %
1747   \dupkey{refPositionsColValue}{\texttt{tilt}}{colvar|rmsd|refPositionsColValue}{\texttt{rmsd} component}
1748 }
1749 \item %
1750   \dupkey{axis}{\texttt{tilt}}{colvar|spinAngle|axis}{\texttt{spinAngle} component}
1751 \end{cvcoptions}
1752  
1753
1754 \cvnamebasedonly{
1755
1756 \cvsubsubsec{\texttt{alpha}: $\alpha$-helix content of a protein segment.}{sec:cvc_alpha}
1757 \labelkey{colvar|alpha}
1758
1759 The block \texttt{alpha~\{...\}} defines the
1760 parameters to calculate the helical content of a segment of protein
1761 residues.  The $\alpha$-helical content across the $N+1$ residues
1762 $N_{0}$ to $N_{0}+N$ is calculated by the formula:
1763 \begin{eqnarray}
1764   \label{eq:colvars_alpha}
1765   { 
1766     \alpha\left(
1767       \mathrm{C}_{\alpha}^{(N_{0})},
1768       \mathrm{O}^{(N_{0})},
1769       \mathrm{C}_{\alpha}^{(N_{0}+1)},
1770       \mathrm{O}^{(N_{0}+1)},
1771       \ldots
1772       \mathrm{N}^{(N_{0}+5)},
1773       \mathrm{C}_{\alpha}^{(N_{0}+5)},
1774       \mathrm{O}^{(N_{0}+5)},
1775       \ldots
1776       \mathrm{N}^{(N_{0}+N)},
1777       \mathrm{C}_{\alpha}^{(N_{0}+N)}
1778     \right)
1779   } \; = \; \; \; \; \\ \; \; \; \; {
1780     \nonumber
1781     \frac{1}{2(N-2)} 
1782     \sum_{n=N_{0}}^{N_{0}+N-2}
1783     \mathrm{angf}\left(
1784         \mathrm{C}_{\alpha}^{(n)},
1785         \mathrm{C}_{\alpha}^{(n+1)},
1786         \mathrm{C}_{\alpha}^{(n+2)}\right)
1787   } \; + \; {
1788     \frac{1}{2(N-4)} 
1789     \sum_{n=N_{0}}^{N_{0}+N-4}
1790     \mathrm{hbf}\left(
1791       \mathrm{O}^{(n)},
1792       \mathrm{N}^{(n+4)}\right) \mathrm{,}
1793   } \\
1794 \end{eqnarray}
1795 where the score function for the $\mathrm{C}_{\alpha} -
1796 \mathrm{C}_{\alpha} - \mathrm{C}_{\alpha}$ angle is defined as: 
1797 \begin{equation}
1798   \label{eq:colvars_alpha_Calpha}
1799   {
1800     \mathrm{angf}\left(
1801       \mathrm{C}_{\alpha}^{(n)},
1802       \mathrm{C}_{\alpha}^{(n+1)},
1803       \mathrm{C}_{\alpha}^{(n+2)}\right)
1804   } \; = \; {
1805     \frac{1 - \left(\theta(
1806         \mathrm{C}_{\alpha}^{(n)},
1807         \mathrm{C}_{\alpha}^{(n+1)},
1808         \mathrm{C}_{\alpha}^{(n+2)}) -
1809         \theta_{0}\right)^{2} /
1810       \left(\Delta\theta_{\mathrm{tol}}\right)^{2}}{
1811       1 - \left(\theta(
1812         \mathrm{C}_{\alpha}^{(n)},
1813         \mathrm{C}_{\alpha}^{(n+1)},
1814         \mathrm{C}_{\alpha}^{(n+2)}) -
1815         \theta_{0}\right)^{4} /
1816       \left(\Delta\theta_{\mathrm{tol}}\right)^{4}} \mathrm{,}
1817   }
1818 \end{equation}
1819 and the score function for the $\mathrm{O}^{(n)} \leftrightarrow
1820 \mathrm{N}^{(n+4)}$ hydrogen bond is defined through a \texttt{hBond}
1821 colvar component on the same atoms.
1822
1823 \begin{cvcoptions}
1824
1825 \item %
1826   \labelkey{colvar|alpha|residueRange}
1827   \key
1828     {residueRange}{%
1829       \texttt{alpha}}{%
1830       Potential $\alpha$-helical residues}{%
1831     ``$<$Initial residue number$>$-$<$Final residue number$>$''}{%
1832     This option specifies the range of residues on which this
1833     component should be defined.  The Colvars module looks for the
1834     atoms within these residues named ``\texttt{CA}'', ``\texttt{N}''
1835     and ``\texttt{O}'', and raises an error if any of those atoms is
1836     not found.}
1837
1838 \item %
1839   \labelkey{colvar|alpha|psfSegID}
1840   \key
1841     {psfSegID}{%
1842     \texttt{alpha}}{%
1843     PSF segment identifier}{%
1844     string (max 4 characters)}{%
1845     This option sets the PSF segment identifier for the residues
1846     specified in \texttt{residueRange}.  This option is only required
1847     when PSF topologies are used.}
1848
1849
1850 \item %
1851   \labelkey{colvar|alpha|hBondCoeff}
1852   \keydef
1853     {hBondCoeff}{%
1854     \texttt{alpha}}{%
1855     Coefficient for the hydrogen bond term}{%
1856     positive between 0 and 1}{%
1857     0.5}{%
1858     This number specifies the contribution to the total value from the
1859     hydrogen bond terms.  0 disables the hydrogen bond terms, 1
1860     disables the angle terms.}
1861
1862 \item %
1863   \labelkey{colvar|alpha|angleRef}
1864   \keydef
1865     {angleRef}{%
1866     \texttt{alpha}}{%
1867     Reference $\mathrm{C}_{\alpha} -
1868     \mathrm{C}_{\alpha} - \mathrm{C}_{\alpha}$ angle}{%
1869     positive decimal}{%
1870     88$^{\circ}$}{%
1871     This option sets the reference angle used in the score function
1872     (\ref{eq:colvars_alpha_Calpha}).}
1873
1874 \item %
1875   \labelkey{colvar|alpha|angleTol}
1876   \keydef
1877     {angleTol}{%
1878     \texttt{alpha}}{%
1879     Tolerance in the $\mathrm{C}_{\alpha} -
1880     \mathrm{C}_{\alpha} - \mathrm{C}_{\alpha}$ angle}{%
1881     positive decimal}{%
1882     15$^{\circ}$}{%
1883     This option sets the angle tolerance used in the score function
1884     (\ref{eq:colvars_alpha_Calpha}).}
1885
1886 \item %
1887   \labelkey{colvar|alpha|hBondCutoff}
1888   \keydef
1889     {hBondCutoff}{%
1890     \texttt{alpha}}{%
1891     Hydrogen bond cutoff}{%
1892     positive decimal}{%
1893     3.3~\AA{}}{%
1894     Equivalent to the \texttt{cutoff} option in the \texttt{hBond}
1895     component.}
1896
1897 \item %
1898   \labelkey{colvar|alpha|hBondExpNumer}
1899   \keydef
1900     {hBondExpNumer}{%
1901     \texttt{alpha}}{%
1902     Hydrogen bond numerator exponent}{%
1903     positive integer}{%
1904     6}{%
1905     Equivalent to the \texttt{expNumer} option in the \texttt{hBond}
1906     component.}
1907
1908 \item %
1909   \labelkey{colvar|alpha|hBondExpDenom}
1910   \keydef
1911     {hBondExpDenom}{%
1912     \texttt{alpha}}{%
1913     Hydrogen bond denominator exponent}{%
1914     positive integer}{%
1915     8}{%
1916     Equivalent to the \texttt{expDenom} option in the \texttt{hBond}
1917     component.}
1918
1919 \end{cvcoptions}
1920
1921 This component returns positive values, always comprised between 0
1922 (lowest $\alpha$-helical score) and 1 (highest $\alpha$-helical
1923 score).
1924
1925
1926 \cvsubsubsec{\texttt{dihedralPC}: protein dihedral pricipal component}{sec:cvc_dihedralPC}
1927 \labelkey{colvar|dihedralPC}
1928
1929 The block \texttt{dihedralPC~\{...\}} defines the
1930 parameters to calculate the projection of backbone dihedral angles within
1931 a protein segment onto a \emph{dihedral principal component}, following
1932 the formalism of dihedral principal component analysis (dPCA) proposed by
1933 Mu et al.\cite{Mu2005} and documented in detail by Altis et
1934 al.\cite{Altis2007}.
1935 Given a peptide or protein segment of $N$ residues, each with Ramachandran
1936 angles $\phi_i$ and $\psi_i$, dPCA rests on a variance/covariance analysis
1937 of the $4(N-1)$ variables $\cos(\psi_1), \sin(\psi_1), \cos(\phi_2), \sin(\phi_2)
1938 \cdots \cos(\phi_N), \sin(\phi_N)$. Note that angles $\phi_1$ and $\psi_N$
1939 have little impact on chain conformation, and are therefore discarded,
1940 following the implementation of dPCA in the analysis software Carma.\cite{Glykos2006}
1941
1942 For a given principal component (eigenvector) of coefficients
1943 $(k_i)_{1 \leq i \leq 4(N-1)}$,
1944 the projection of the current backbone conformation is:
1945 \begin{equation}
1946 \xi = \sum_{n=1}^{N-1} k_{4n-3} \cos(\psi_n) + k_{4n-2} \sin (\psi_n)
1947 + k_{4n-1} \cos (\phi_{n+1}) + k_{4n} \sin(\phi_{n+1})
1948 \end{equation}
1949
1950 \texttt{dihedralPC} expects the same parameters as the \texttt{alpha}
1951 component for defining the relevant residues (\texttt{residueRange}
1952 and \texttt{psfSegID}) in addition to the following:
1953
1954 \begin{cvcoptions}
1955
1956 \item %
1957   \dupkey{residueRange}{\texttt{dihedralPC}}{colvar|alpha|residueRange}{\texttt{alpha} component}
1958
1959 \item %
1960   \dupkey{psfSegID}{\texttt{dihedralPC}}{colvar|alpha|psfSegID}{\texttt{alpha} component}
1961
1962 \item %
1963   \key
1964     {vectorFile}{%
1965     \texttt{dihedralPC}}{%
1966     File containing dihedral PCA eigenvector(s)}{%
1967     file name}{%
1968     A text file containing the coefficients of dihedral PCA eigenvectors on the
1969     cosine and sine coordinates. The vectors should be arranged in columns,
1970     as in the files output by Carma.\cite{Glykos2006}}
1971
1972 \item %
1973   \key
1974     {vectorNumber}{%
1975     \texttt{dihedralPC}}{%
1976     File containing dihedralPCA eigenvector(s)}{%
1977     positive integer}{%
1978     Number of the eigenvector to be used for this component.}
1979 \end{cvcoptions}
1980
1981 } % end of \cvnamebasedonly
1982
1983
1984 \cvsubsec{Shared keywords for all components}{sec:cvc_common}
1985
1986 The following options can be used for any of the above colvar components in order to obtain a polynomial combination\cvscriptonly{ or any user-supplied function provided by \refkey{scriptedFunction}{sec:cvc_superp}}.
1987 \begin{itemize}
1988 \item %
1989   \keydef
1990     {name}{%
1991     any component}{%
1992     Name of this component}{%
1993     string}{%
1994     type of component + numeric id}{%
1995     The name is an unique case-sensitive string which allows the
1996     Colvars module to identify this component. This is useful, for example,
1997     when combining multiple components via a \texttt{scriptedFunction}.
1998     \cvnamdonly{It also defines the variable name representing the component's value in a \refkey{customFunction}{colvar|customFunction} expression.}}
1999
2000 \item %
2001   \keydef
2002     {scalable}{%
2003     any component}{%
2004     Attempt to calculate this component in parallel?}{%
2005     boolean}{%
2006     \texttt{on}, if available}{%
2007     If set to \texttt{on} (default), the Colvars module will attempt to calculate this component in parallel to reduce overhead.
2008     Whether this option is available depends on the type of component: currently supported are \texttt{distance}, \texttt{distanceZ}, \texttt{distanceXY}, \texttt{distanceVec}, \texttt{distanceDir}, \texttt{angle} and \texttt{dihedral}.
2009     This flag influences computational cost, but does not affect numerical results: therefore, it should only be turned off for debugging or testing purposes.
2010   }
2011 \end{itemize}
2012
2013
2014 \cvsubsec{Periodic components}{sec:cvc_periodic}
2015 The following components returns
2016 real numbers that lie in a periodic interval:
2017 \begin{itemize}
2018 \item \texttt{dihedral}: torsional angle between four groups;
2019 \item \texttt{spinAngle}: angle of rotation around a predefined axis
2020   in the best-fit from a set of reference coordinates.
2021 \end{itemize}
2022 In certain conditions, \texttt{distanceZ} can also be periodic, namely
2023 when periodic boundary conditions (PBCs) are defined in the simulation
2024 and \texttt{distanceZ}'s axis is parallel to a unit cell vector.
2025
2026 In addition, a custom \cvscriptonly{or scripted} scalar colvar may be periodic
2027 depending on its user-defined expression. It will only be treated as such by
2028 the Colvars module if the period is specified using the \texttt{period} keyword,
2029 while \texttt{wrapAround} is optional.
2030
2031 The following keywords can be used within periodic components\cvleptononly{, or within custom variables (\ref{sec:colvar_custom_function})}\cvscriptonly{, or wthin scripted variables \ref{sec:colvar_scripted}}).
2032
2033 \begin{itemize}
2034 \item %
2035   \keydef
2036     {period}{%
2037     \texttt{distanceZ}, custom colvars}{%
2038     Period of the component}{%
2039     positive decimal}{%
2040     0.0}{%
2041     Setting this number enables the treatment of \texttt{distanceZ} as
2042     a periodic component: by default, \texttt{distanceZ} is not
2043     considered periodic.  The keyword is supported, but irrelevant
2044     within \texttt{dihedral} or \texttt{spinAngle}, because their
2045     period is always 360~degrees.}
2046
2047 \item %
2048   \keydef
2049     {wrapAround}{%
2050     \texttt{distanceZ}, \texttt{dihedral}, \texttt{spinAngle}, custom colvars}{%
2051     Center of the wrapping interval for periodic variables}{%
2052     decimal}{%
2053     0.0}{%
2054     By default, values of the periodic components are centered around zero, ranging from $-P/2$ to $P/2$, where $P$ is the period.
2055     Setting this number centers the interval around this value.
2056     This can be useful for convenience of output, or to set the walls for a \texttt{harmonicWalls} in an order that would not otherwise be allowed.}
2057 \end{itemize}
2058
2059 Internally, all differences between two values of a periodic colvar
2060 follow the minimum image convention: they are calculated based on
2061 the two periodic images that are closest to each other.
2062
2063 \emph{Note: linear or polynomial combinations of periodic components (see \ref{sec:cvc_superp}) may become meaningless when components cross the periodic boundary.  Use such combinations carefully: estimate the range of possible values of each component in a given simulation, and make use of \texttt{wrapAround} to limit this problem whenever possible.}
2064
2065
2066 \cvsubsec{Non-scalar components}{sec:cvc_non_scalar}
2067
2068 When one of the following components are used, the defined colvar returns a value that is not a scalar number:
2069 \begin{itemize}
2070 \item \texttt{distanceVec}: 3-dimensional vector of the distance
2071   between two groups;
2072 \item \texttt{distanceDir}: 3-dimensional unit vector of the distance
2073   between two groups;
2074 \item \texttt{orientation}: 4-dimensional unit quaternion representing
2075   the best-fit rotation from a set of reference coordinates.
2076 \end{itemize}
2077 The distance between two 3-dimensional unit vectors is computed as the
2078 angle between them.  The distance between two quaternions is computed
2079 as the angle between the two 4-dimensional unit vectors: because the
2080 orientation represented by $\mathsf{q}$ is the same as the one
2081 represented by $-\mathsf{q}$, distances between two quaternions are
2082 computed considering the closest of the two symmetric images.
2083
2084 Non-scalar components carry the following restrictions:
2085 \begin{itemize}
2086 \item Calculation of total forces (\texttt{outputTotalForce} option)
2087   is currently not implemented.
2088 \item Each colvar can only contain one non-scalar component.
2089 \item Binning on a grid (\texttt{abf}, \texttt{histogram} and
2090   \texttt{metadynamics} with \texttt{useGrids} enabled) is currently
2091   not implemented for colvars based on such components.
2092 \end{itemize}
2093
2094 \emph{Note: while these restrictions apply to individual colvars based
2095   on non-scalar components, no limit is set to the number of scalar
2096   colvars.  To compute multi-dimensional histograms and PMFs, use sets
2097   of scalar colvars of arbitrary size.}
2098
2099
2100 \cvsubsubsec{Calculating total forces}{sec:cvc_sys_forces}
2101 In addition to the restrictions due to the type of value computed (scalar or non-scalar),
2102 a final restriction can arise when calculating total force
2103 (\texttt{outputTotalForce} option or application of a \texttt{abf}
2104 bias).  total forces are available currently only for the following
2105 components: \texttt{distance}, \texttt{distanceZ},
2106 \texttt{distanceXY}, \texttt{angle}, \texttt{dihedral}, \texttt{rmsd},
2107 \texttt{eigenvector} and \texttt{gyration}.
2108
2109
2110
2111 \cvsubsec{Linear and polynomial combinations of components}{sec:cvc_superp}
2112
2113 To extend the set of possible definitions of colvars $\xi(\mathbf{r})$, multiple components
2114 $q_i(\mathbf{r})$ can be summed with the formula:
2115 \begin{equation}
2116   \label{eq:colvar_combination}
2117   \xi(\mathbf{r}) = \sum_i c_i [q_i(\mathbf{r})]^{n_i}
2118 \end{equation}
2119 where each component appears with a unique coefficient $c_i$ (1.0 by
2120 default) the positive integer exponent $n_i$ (1 by default).
2121
2122 Any set of components can be combined within a colvar, provided that
2123 they return the same type of values (scalar, unit vector, vector, or
2124 quaternion).  By default, the colvar is the sum of its components.
2125 Linear or polynomial combinations (following
2126 equation~(\ref{eq:colvar_combination})) can be obtained by setting the
2127 following parameters, which are common to all components:
2128 \begin{itemize}
2129 \item %
2130   \keydef
2131     {componentCoeff}{%
2132     any component}{%
2133     Coefficient of this component in the colvar}{%
2134     decimal}{%
2135     \texttt{1.0}}{%
2136     Defines the coefficient by which this component is multiplied
2137     (after being raised to \texttt{componentExp}) before being added
2138     to the sum.}
2139
2140 \item %
2141   \keydef
2142     {componentExp}{%
2143     any component}{%
2144     Exponent of this component in the colvar}{%
2145     integer}{%
2146     \texttt{1}}{%
2147     Defines the power at which the value of this component is raised
2148     before being added to the sum.  When this exponent is
2149     different than 1 (non-linear sum), total forces and the Jacobian
2150     force are not available, making the colvar unsuitable for ABF calculations.}
2151 \end{itemize}
2152
2153 \textbf{Example:} To define the \emph{average} of a colvar across
2154 different parts of the system, simply define within the same colvar
2155 block a series of components of the same type (applied to different
2156 atom groups), and assign to each component a \texttt{componentCoeff}
2157 of $1/N$.
2158
2159
2160 \cvleptononly{
2161 \cvsubsec{Custom functions}{sec:colvar_custom_function}
2162
2163 Collective variables may be defined by specifying a custom function as an analytical
2164 expression such as \texttt{cos(x) + y\^{}2}.
2165 The expression is parsed by the Lepton expression parser (written by Peter Eastman),
2166 which produces efficient evaluation routines for the function itself as well as its derivatives.
2167 The expression may use the collective variable components as variables, refered to as their \texttt{name} string.
2168 Scalar elements of vector components may be accessed by appending a 1-based index to their \texttt{name}.
2169 When implementing generic functions of Cartesian coordinates rather
2170 than functions of existing components, the \texttt{cartesian} component
2171 may be particularly useful.
2172 A scalar-valued custom variable may be manually defined as periodic by providing
2173 the keyword \texttt{period}, and the optional keyword \texttt{wrapAround}, with the
2174 same meaning as in periodic components (see \ref{sec:cvc_periodic} for details).
2175 A vector variable may be defined by specifying the \texttt{customFunction} parameter several times: each expression defines one scalar element of the vector colvar.
2176 This is illustrated in the example below.
2177
2178 \bigskip
2179 {%
2180 % verbatim can't appear within commands
2181 \noindent\ttfamily colvar \{\\
2182 \-~~name custom\\
2183 \\
2184 \-~~\# A 2-dimensional vector function of a scalar x and a 3-vector r\\
2185 \-~~customFunction cos(x) * (r1 + r2 + r3)\\
2186 \-~~customFunction sqrt(r1 * r2)\\
2187 \\
2188 \-~~distance \{\\
2189 \-~~~~name x\\
2190 \-~~~~group1 \{ atomNumbers 1 \}\\
2191 \-~~~~group2 \{ atomNumbers 50 \}\\
2192 \-~~\}\\
2193 \-~~distanceVec \{\\
2194 \-~~~~name r\\
2195 \-~~~~group1 \{ atomNumbers 10 11 12 \}\\
2196 \-~~~~group2 \{ atomNumbers  20 21 22 \}\\
2197 \-~~\}\\
2198 \}}
2199
2200 \begin{itemize}
2201 \item %
2202    \labelkey{colvar|customFunction}
2203    \key 
2204     {customFunction}{%
2205     \texttt{colvar}}{%
2206     Compute colvar as a custom function of its components}{%
2207     string}{%
2208     Defines the colvar as a scalar expression of its colvar components. Multiple mentions can
2209     be used to define a vector variable  (as in the example above).}
2210
2211 \item %
2212   \keydef
2213     {customFunctionType}{%
2214     \texttt{colvar}}{%
2215     Type of value returned by the scripted colvar}{%
2216     string}{%
2217     \texttt{scalar}}{%
2218     With this flag, the user may specify whether the
2219     colvar is a scalar or one of the following vector types: \texttt{vector3}
2220     (a 3D vector), \texttt{unit\_vector3} (a normalized 3D vector), or
2221     \texttt{unit\_quaternion} (a normalized quaternion), or \texttt{vector}.
2222     Note that the scalar and vector cases are not necessary, as they are detected automatically.}
2223 \end{itemize}
2224 }
2225
2226
2227 \cvscriptonly{
2228 \cvsubsec{Scripted functions}{sec:colvar_scripted}
2229 When scripting is supported\cvnamdonly{ (default in NAMD)}\cvvmdonly{ (default in VMD)},
2230 a colvar may be defined as a scripted function of its components,
2231 rather than a linear or polynomial combination.
2232 When implementing generic functions of Cartesian coordinates rather
2233 than functions of existing components, the \texttt{cartesian} component
2234 may be particularly useful.
2235 A scalar-valued scripted variable may be manually defined as periodic by providing
2236 the keyword \texttt{period}, and the optional keyword \texttt{wrapAround}, with the
2237 same meaning as in periodic components (see \ref{sec:cvc_periodic} for details).
2238
2239 An example of elaborate scripted colvar is given in example 10, in the
2240 form of path-based collective variables as defined by Branduardi et al\cite{Branduardi2007}
2241 (\ref{sec:pathcv}).
2242
2243 \begin{itemize}
2244  \item %
2245    \labelkey{colvar|scriptedFunction}
2246    \key 
2247     {scriptedFunction}{%
2248     \texttt{colvar}}{%
2249     Compute colvar as a scripted function of its components}{%
2250     string}{%
2251     If this option is specified, the colvar will be computed as a
2252     scripted function of the values of its components.
2253     To that effect, the user should define two Tcl procedures:
2254     \texttt{calc\_$<$scriptedFunction$>$} and \texttt{calc\_$<$scriptedFunction$>$\_gradient},
2255     both accepting as many parameters as the colvar has components.
2256     Values of the components will be passed to those procedures in the
2257     order defined by their sorted \texttt{name} strings. Note that if all
2258     components are of the same type, their default names are sorted in the
2259     order in which they are defined, so that names need only be specified
2260     for combinations of components of different types.
2261     \texttt{calc\_$<$scriptedFunction$>$} should return one value of 
2262     type $<$scriptedFunctionType$>$, corresponding to the colvar value.
2263     \texttt{calc\_$<$scriptedFunction$>$\_gradient} should return a Tcl list
2264     containing the derivatives of the function with respect to each
2265     component. 
2266     If both the function and some of the components are vectors, the gradient
2267     is really a Jacobian matrix that should be passed as a linear vector in
2268     row-major order, i.e. for a function $f_i(x_j)$: $\nabla_x f_1 \nabla_x f_2 \cdots$.
2269   }
2270
2271 \item%
2272   \keydef
2273     {scriptedFunctionType}{%
2274     \texttt{colvar}}{%
2275     Type of value returned by the scripted colvar}{%
2276     string}{%
2277     \texttt{scalar}}{%
2278     If a colvar is defined as a scripted function, its type is not constrained by
2279     the types of its components. With this flag, the user may specify whether the
2280     colvar is a scalar or one of the following vector types: \texttt{vector3}
2281     (a 3D vector), \texttt{unit\_vector3} (a normalized 3D vector), or
2282     \texttt{unit\_quaternion} (a normalized quaternion), or \texttt{vector}
2283     (a vector whose size is specified by \texttt{scriptedFunctionVectorSize}).
2284     Non-scalar values should be passed as space-separated lists.}
2285
2286  \item %
2287   \key
2288     {scriptedFunctionVectorSize}{%
2289     \texttt{colvar}}{%
2290     Dimension of the vector value of a scripted colvar}{%
2291     positive integer}{%
2292     This parameter is only valid when \texttt{scriptedFunctionType} is
2293     set to \texttt{vector}. It defines the vector length of the colvar value
2294     returned by the function.}
2295 \end{itemize}
2296 } % \cvscriptonly
2297
2298
2299
2300 \cvsubsec{Defining grid parameters}{sec:colvar_grid_params}
2301
2302 Many algorithms require the definition of boundaries and/or characteristic spacings that can be used to define discrete ``states'' in the collective variable, or to combine variables with very different units.
2303 The parameters described below offer a way to specify these parameters only once for each variable, while using them multiple times in restraints, time-dependent biases or analysis methods.
2304
2305 \begin{itemize}
2306
2307 \item %
2308   \labelkey{colvar|width}
2309   \keydef
2310     {width}{%
2311     \texttt{colvar}}{%
2312     Unit of the variable, or grid spacing}{%
2313     positive decimal}{%
2314     1.0}{%
2315     This number defines the effective unit of measurement for the collective variable, and is used by the biasing methods for the following purposes.
2316     Harmonic (\ref{sec:colvarbias_harmonic}), harmonic walls (\ref{sec:colvarbias_harmonic_walls}) and linear restraints (\ref{sec:colvarbias_linear}) use it to set the physical unit of the force constant, which is useful for multidimensional restraints involving multiple variables with very different units (for examples, $\AA$ or degrees $^\circ$) with a single, scaled force constant.
2317     The values of the scaled force constant in the units of each variable are printed at initialization time. 
2318     Histograms (\ref{sec:colvarbias_histogram}), ABF (\ref{sec:colvarbias_abf}) and metadynamics (\ref{sec:colvarbias_meta}) all use this number as the initial choice for the grid spacing along this variable: for this reason, \texttt{width} should generally be no larger than the standard deviation of the colvar in an unbiased simulation.
2319     Unless it is required to control the spacing, it is usually simplest to keep the default value of 1, so that restraint force constants are provided with their full physical unit.
2320   }
2321
2322 \item %
2323   \labelkey{colvar|lowerBoundary}
2324   \key
2325     {lowerBoundary}{%
2326     \texttt{colvar}}{%
2327     Lower boundary of the colvar}{%
2328     decimal}{%
2329     Defines the lowest end of the interval of ``relevant'' values for the colvar.
2330     This number can be either a true physical boundary, or a user-defined number.  
2331     Together with \texttt{upperBoundary} and \texttt{width}, it is used to define a grid of values along the variable (not available for variables with vector values, \ref{sec:cvc_non_scalar}).
2332     \emph{This option does not affect dynamics: to confine a colvar within a certain interval, use a \texttt{harmonicWalls} bias.}
2333 }
2334
2335 \item %
2336   \labelkey{colvar|upperBoundary} %
2337   \key
2338     {upperBoundary}{%
2339     \texttt{colvar}}{%
2340     Upper boundary of the colvar}{%
2341     decimal}{%
2342     Similarly to \texttt{lowerBoundary}, defines the highest possible or allowed value.}
2343
2344 \item %
2345   \keydef
2346     {hardLowerBoundary}{%
2347     \texttt{colvar}}{%
2348     Whether the lower boundary is the physical lower limit}{%
2349     boolean}{%
2350     \texttt{off}}{%
2351     This option does not affect simulation results, but enables some internal optimizations.
2352     Depending on its mathematical definition, a colvar may have ``natural'' boundaries: for example, a \texttt{distance} colvar has a ``natural'' lower boundary at 0.  Setting this option instructs the Colvars module that the user-defined lower boundary is ``natural''.
2353 See Section~\ref{sec:cvc_list} for the physical ranges of values of each component.}
2354
2355 \item %
2356   \keydef
2357     {hardUpperBoundary}{%
2358     \texttt{colvar}}{%
2359     Whether the upper boundary is the physical upper limit of the colvar's values}{%
2360     boolean}{%
2361     \texttt{off}}{%
2362     Analogous to \texttt{hardLowerBoundary}.}
2363
2364 \item %
2365   \labelkey{colvar|expandBoundaries} %
2366   \keydef
2367     {expandBoundaries}{%
2368     \texttt{colvar}}{%
2369     Allow to expand the two boundaries if needed}{%
2370     boolean}{%
2371     \texttt{off}}{%
2372     If defined, biasing and analysis methods may keep their own copies
2373     of \texttt{lowerBoundary} and \texttt{upperBoundary}, and expand
2374     them to accommodate values that do not fit in the initial range.
2375     Currently, this option is used by the metadynamics bias
2376     (\ref{sec:colvarbias_meta}) to keep all of its hills fully within
2377     the grid.  This option cannot be used when
2378       the initial boundaries already span the full period of a periodic
2379       colvar.}
2380
2381 \end{itemize}
2382
2383
2384
2385 \cvsubsec{Trajectory output}{sec:colvar_traj_output}
2386
2387 \begin{itemize}
2388 \item %
2389   \keydef
2390     {outputValue}{%
2391     \texttt{colvar}}{%
2392     Output a trajectory for this colvar}{%
2393     boolean}{%
2394     \texttt{on}}{%
2395     If \texttt{colvarsTrajFrequency} is non-zero, the value of this
2396     colvar is written to the trajectory file every
2397     \texttt{colvarsTrajFrequency} steps in the column labeled
2398     ``$<$\texttt{name}$>$''.}
2399
2400 \item %
2401   \keydef
2402     {outputVelocity}{%
2403     \texttt{colvar}}{%
2404     Output a velocity trajectory for this colvar}{%
2405     boolean}{%
2406     \texttt{off}}{%
2407     If \texttt{colvarsTrajFrequency} is defined, the
2408     finite-difference calculated velocity of this colvar are written
2409     to the trajectory file under the label
2410     ``\texttt{v\_}$<$\texttt{name}$>$''.}
2411
2412 \item %
2413   \keydef
2414     {outputEnergy}{%
2415     \texttt{colvar}}{%
2416     Output an energy trajectory for this colvar}{%
2417     boolean}{%
2418     \texttt{off}}{%
2419     This option applies only to extended Lagrangian colvars. If
2420     \texttt{colvarsTrajFrequency} is defined, the kinetic energy of
2421     the extended degree and freedom and the potential energy of the
2422     restraining spring are are written to the trajectory file under
2423     the labels ``\texttt{Ek\_}$<$\texttt{name}$>$'' and
2424     ``\texttt{Ep\_}$<$\texttt{name}$>$''.}
2425
2426 \item %
2427   \labelkey{colvar|outputTotalForce}
2428   \keydef
2429     {outputTotalForce}{%
2430     \texttt{colvar}}{%
2431     Output a total force trajectory for this
2432     colvar}{%
2433     boolean}{%
2434     \texttt{off}}{%
2435     If \texttt{colvarsTrajFrequency} is defined, the total force on this
2436     colvar (i.e.~the projection of all atomic total forces
2437     onto this colvar --- see
2438     equation~(\ref{eq:gradient_vector}) in
2439     section~\ref{sec:colvarbias_abf}) are written to the trajectory
2440     file under the label ``\texttt{fs\_}$<$\texttt{name}$>$''.
2441     For extended Lagrangian colvars, the ``total force'' felt by the extended degree of freedom
2442     is simply the force from the harmonic spring.
2443     \textbf{Note:} not all components support this option.  The
2444     physical unit for this force is \cvnamdonly{kcal/mol}\cvvmdonly{kcal/mol}\cvlammpsonly{the unit of energy specified by \texttt{units}}, divided by the colvar unit U.}
2445
2446 \item %
2447   \keydef
2448     {outputAppliedForce}{%
2449     \texttt{colvar}}{%
2450     Output an applied force trajectory for this
2451     colvar}{%
2452     boolean}{%
2453     \texttt{off}}{%
2454     If \texttt{colvarsTrajFrequency} is defined, the total force
2455     applied on this colvar by Colvars biases are
2456     written to the trajectory under the label
2457     ``\texttt{fa\_}$<$\texttt{name}$>$''. 
2458     For extended Lagrangian colvars, this force is actually applied to the
2459     extended degree of freedom rather than the geometric colvar itself.
2460     The physical unit for this
2461     force is \cvnamdonly{kcal/mol}\cvvmdonly{kcal/mol}\cvlammpsonly{the unit of energy specified by \texttt{units}} divided by the colvar unit.}
2462
2463 \end{itemize}
2464
2465
2466 \cvsubsec{Extended Lagrangian}{sec:colvar_extended}
2467
2468 The following options enable extended-system
2469 dynamics, where a colvar is coupled to an additional degree of freedom 
2470 (fictitious particle) by a harmonic spring.
2471 All biasing and confining forces are then applied to the extended degree
2472 of freedom. The ``actual'' geometric colvar (function of Cartesian 
2473 coordinates) only feels the force from the harmonic spring.
2474 This is particularly useful when combined with an ABF bias (\ref{sec:colvarbias_abf})
2475 to perform eABF simulations (\ref{sec:eABF}).
2476
2477 \begin{itemize}
2478 \item %
2479   \keydef
2480     {extendedLagrangian}{%
2481     \texttt{colvar}}{%
2482     Add extended degree of freedom}{%
2483     boolean}{%
2484     \texttt{off}}{%
2485     Adds a fictitious particle to be coupled to the colvar by a harmonic
2486     spring. The fictitious mass and the force constant of the coupling
2487     potential are derived from the parameters \texttt{extendedTimeConstant}
2488     and \texttt{extendedFluctuation}, described below. Biasing forces on the
2489     colvar are applied to this fictitious particle, rather than to the
2490     atoms directly.  This implements the extended Lagrangian formalism
2491     used in some metadynamics simulations~\cite{Iannuzzi2003}.
2492     \cvnamdonly{The energy associated with the extended degree of freedom is reported
2493     under the MISC title in NAMD's energy output.}
2494     }
2495
2496 \item %
2497   \key
2498     {extendedFluctuation}{%
2499     \texttt{colvar}}{%
2500     Standard deviation between the colvar and the fictitious
2501     particle (colvar unit)}{%
2502     positive decimal}{%
2503     Defines the spring stiffness for the \texttt{extendedLagrangian}
2504     mode, by setting the typical deviation between the colvar and the extended
2505     degree of freedom due to thermal fluctuation.
2506     The spring force constant is calculated internally as $k_B T / \sigma^2$,
2507     where $\sigma$ is the value of \texttt{extendedFluctuation}.}
2508
2509 \item %
2510   \keydef
2511     {extendedTimeConstant}{%
2512     \texttt{colvar}}{%
2513     Oscillation period of the fictitious particle (fs)}{%
2514     positive decimal}{%
2515     \texttt{200}}{%
2516     Defines the inertial mass of the fictitious particle, by setting the
2517     oscillation period of the harmonic oscillator formed by the fictitious
2518     particle and the spring. The period
2519     should be much larger than the MD time step to ensure accurate integration
2520     of the extended particle's equation of motion.
2521     The fictitious mass is calculated internally as $k_B T (\tau/2 \pi \sigma)^2$,
2522     where $\tau$ is the period and $\sigma$ is the typical fluctuation (see above).}
2523
2524 \item %
2525   \keydef
2526     {extendedTemp}{%
2527     \texttt{colvar}}{%
2528     Temperature for the extended degree of freedom (K)}{%
2529     positive decimal}{%
2530     thermostat temperature}{%
2531     Temperature used for calculating the coupling force constant of the
2532     extended variable (see \texttt{extendedFluctuation}) and, if needed, as a
2533     target temperature for extended Langevin dynamics (see
2534     \texttt{extendedLangevinDamping}). This should normally be left at its
2535     default value.}
2536
2537 \item %
2538   \keydef
2539     {extendedLangevinDamping}{%
2540     \texttt{colvar}}{%
2541     Damping factor for extended Langevin dynamics
2542     (ps$^{-1}$)}{%
2543     positive decimal}{%
2544     \texttt{1.0}}{%
2545     If this is non-zero, the extended degree of freedom undergoes Langevin dynamics
2546     at temperature \texttt{extendedTemp}. The friction force is minus
2547     \texttt{extendedLangevinDamping} times the velocity. This is useful because
2548     the extended dynamics coordinate may heat up in the transient
2549     non-equilibrium regime of ABF. Use moderate damping values, to limit
2550     viscous friction (potentially slowing down diffusive sampling) and stochastic
2551     noise (increasing the variance of statistical measurements). In
2552     doubt, use the default value.}
2553 \end{itemize}
2554
2555
2556 \cvsubsec{Backward-compatibility}{sec:colvar_compatibility}
2557
2558 \begin{itemize}
2559
2560 \item %
2561   \keydef
2562     {subtractAppliedForce}{%
2563     \texttt{colvar}}{%
2564     Do not include biasing forces in the total force for this colvar}{%
2565     boolean}{%
2566     \texttt{off}}{%
2567     If the colvar supports total force calculation (see \ref{sec:cvc_sys_forces}), all forces applied to this colvar by biases will be removed from the total force.
2568     This keyword allows to recover some of the ``system force'' calculation available in the Colvars module     before version 2016-08-10.
2569     Please note that removal of \emph{all} other external forces (including biasing forces applied to a         different colvar) is \emph{no longer supported}, due to changes in the underlying simulation engines (primarily NAMD).
2570     This option may be useful when continuing a previous simulation where the removal of external/applied forces is essential.
2571     \emph{For all new simulations, the use of this option is not recommended.}
2572 }
2573
2574 \end{itemize}
2575
2576
2577 \cvsubsec{Statistical analysis}{sec:colvar_acf}
2578
2579 Run-time calculations of statistical properties that depend explicitly on time can be performed for individual collective variables.
2580 Currently, several types of time correlation functions, running averages and running standard deviations are implemented.
2581 For run-time computation of histograms, please see the histogram bias (\ref{sec:colvarbias_histogram}).
2582
2583 \begin{itemize}
2584
2585 \item %
2586   \labelkey{colvar|corrFunc}
2587   \keydef
2588     {corrFunc}{%
2589     \texttt{colvar}}{%
2590     Calculate a time correlation function?}{%
2591     boolean}{%
2592     \texttt{off}}{%
2593     Whether or not a time correlaction function should be calculated
2594     for this colvar.}
2595
2596 \item %
2597   \key
2598     {corrFuncWithColvar}{%
2599     \texttt{colvar}}{%
2600     Colvar name for the correlation function}{%
2601     string}{%
2602     By default, the auto-correlation function (ACF) of this colvar,
2603     $\xi_{i}$, is calculated.  When this option is specified, the
2604     correlation function is calculated instead with another colvar,
2605     $\xi_{j}$, which must be of the same type (scalar, vector, or
2606     quaternion) as $\xi_{i}$.}
2607
2608 \item%
2609   \keydef
2610     {corrFuncType}{%
2611     \texttt{colvar}}{%
2612     Type of the correlation function}{%
2613     \texttt{velocity}, \texttt{coordinate} or
2614     \texttt{coordinate\_p2}}{%
2615     \texttt{velocity}}{%
2616     With \texttt{coordinate} or \texttt{velocity}, the correlation
2617     function $C_{i,j}(t)$~= $\left\langle \Pi\left(\xi_{i}(t_{0}),
2618         \xi_{j}(t_{0}+t)\right) \right\rangle$ is calculated between
2619     the variables $\xi_{i}$ and $\xi_{j}$, or their velocities.
2620     $\Pi(\xi_{i}, \xi_{j})$ is the scalar product when calculated
2621     between scalar or vector values, whereas for quaternions it is the
2622     cosine between the two corresponding rotation axes.  With
2623     \texttt{coordinate\_p2}, the second order Legendre polynomial,
2624     $(3\cos(\theta)^{2}-1)/2$, is used instead of the cosine.}
2625
2626 \item %
2627   \keydef
2628     {corrFuncNormalize}{%
2629     \texttt{colvar}}{%
2630     Normalize the time correlation function?}{%
2631     boolean}{%
2632     \texttt{on}}{%
2633     If enabled, the value of the correlation function at $t$~= 0
2634     is normalized to 1; otherwise, it equals to $\left\langle
2635       O\left(\xi_{i}, \xi_{j}\right) \right\rangle$.}
2636
2637 \item %
2638   \keydef
2639     {corrFuncLength}{%
2640     \texttt{colvar}}{%
2641     Length of the time correlation function}{%
2642     positive integer}{%
2643     \texttt{1000}}{%
2644     Length (in number of points) of the time correlation function.}
2645
2646 \item %
2647   \keydef
2648     {corrFuncStride}{%
2649     \texttt{colvar}}{%
2650     Stride of the time correlation function}{%
2651     positive integer}{%
2652     \texttt{1}}{%
2653     Number of steps between two values of the time correlation function.}
2654
2655 \item %
2656   \keydef
2657     {corrFuncOffset}{%
2658     \texttt{colvar}}{%
2659     Offset of the time correlation function}{%
2660     positive integer}{%
2661     \texttt{0}}{%
2662     The starting time (in number of steps) of the time correlation
2663     function (default: $t$~= 0).  \textbf{Note:} \emph{the value at $t$~= 0 is always
2664     used for the normalization}.}
2665
2666 \item %
2667   \keydef
2668     {corrFuncOutputFile}{%
2669     \texttt{colvar}}{%
2670     Output file for the time correlation function}{%
2671     UNIX filename}{%
2672     \outputName\texttt{.$<$name$>$.corrfunc.dat}}{%
2673     The time correlation function is saved in this file.}
2674
2675 \item %
2676   \labelkey{colvar|runAve}
2677   \keydef
2678     {runAve}{%
2679     \texttt{colvar}}{%
2680     Calculate the running average and standard deviation}{%
2681     boolean}{%
2682     \texttt{off}}{%
2683     Whether or not the running average and standard deviation should
2684     be calculated for this colvar.}
2685
2686 \item %
2687   \keydef
2688     {runAveLength}{%
2689     \texttt{colvar}}{%
2690     Length of the running average window}{%
2691     positive integer}{%
2692     \texttt{1000}}{%
2693     Length (in number of points) of the running average window.}
2694
2695 \item %
2696   \keydef
2697     {runAveStride}{%
2698     \texttt{colvar}}{%
2699     Stride of the running average window values}{%
2700     positive integer}{%
2701     \texttt{1}}{%
2702     Number of steps between two values within the running average window.}
2703
2704 \item %
2705   \keydef
2706     {runAveOutputFile}{%
2707     \texttt{colvar}}{%
2708     Output file for the running average and standard deviation}{%
2709     UNIX filename}{%
2710     \outputName\texttt{.$<$name$>$.runave.traj}}{%
2711     The running average and standard deviation are saved in this file.}
2712
2713 \end{itemize}
2714
2715
2716 \cvsec{Selecting atoms}{sec:colvar_atom_groups}
2717
2718 To define collective variables, atoms are usually selected as groups.  Each group is defined using an identifier that is unique in the context of the specific colvar component (e.g.~for a distance component, the two groups are \texttt{group1} and \texttt{group2}).
2719 The identifier is followed by a brace-delimited block containing selection keywords and other parameters, including an optional \texttt{name}:
2720
2721 \begin{itemize}
2722 \item \key
2723   {name}{%
2724   atom group}{%
2725   Unique name for the atom group}{%
2726   string}{%
2727   This parameter defines a unique name for this atom group, which can be referred to
2728   in the definition of other atom groups (including in other colvars) by invoking
2729   \texttt{atomsOfGroup} as a selection keyword.}
2730 \end{itemize}
2731
2732
2733 \cvsubsec{Atom selection keywords}{sec:colvar_atom_groups_sel}
2734
2735 Selection keywords may be used individually or in combination with each other, and each can be repeated any number of times.
2736 Selection is incremental: each keyword adds the corresponding atoms to the selection, so that different sets of atoms can be combined.
2737 However, atoms included by multiple keywords are only counted once.
2738 Below is an example configuration for an atom group called ``\texttt{atoms}''.
2739 \textbf{Note: }\emph{this is an unusually varied combination of selection keywords, demonstrating how they can be combined together: most simulations only use one of them.}\\
2740
2741 {%
2742 % verbatim can't appear within commands
2743 \noindent\ttfamily atoms \{\\
2744 \\
2745 \-~~\# add atoms 1 and 3 to this group (note: the first atom in the system is 1)\\
2746 \-~~atomNumbers \{ \\
2747 \-~~~~1 3\\
2748 \-~~\}\\
2749 \\
2750 \-~~\# add atoms starting from 20 up to and including 50\\
2751 \-~~atomNumbersRange  20-50\\
2752 }
2753 \cvnamebasedonly{{%
2754 \noindent\ttfamily\\
2755 \-~~\# add all the atoms with occupancy 2 in the file atoms.pdb\\
2756 \-~~atomsFile             atoms.pdb\\
2757 \-~~atomsCol              O\\
2758 \-~~atomsColValue         2.0\\
2759 \\
2760 \-~~\# add all the C-alphas within residues 11 to 20 of segments "PR1" and "PR2"\\
2761 \-~~psfSegID              PR1 PR2\\
2762 \-~~atomNameResidueRange  CA 11-20\\
2763 \-~~atomNameResidueRange  CA 11-20\\
2764 }}
2765 {\noindent\ttfamily\\
2766 \-~~\# add index group (requires a .ndx file to be provided globally)\\
2767 \-~~indexGroup Water\\
2768 \}\\}
2769
2770
2771 The resulting selection includes atoms 1 and 3, those between 20 and 50,\cvnamebasedonly{ the $\mathrm{C}_{\alpha}$ atoms between residues 11 and 20 of the two segments \texttt{PR1} and \texttt{PR2},} and those in the index group called ``Water''.
2772 The indices of this group are read from the file provided by the global keyword \refkey{indexFile}{Colvars-global|indexFile}.
2773
2774 \cvvmdonly{In the current version, the Colvars module does not manipulate VMD atom selections directly: however, these can be converted to atom groups within the Colvars configuration string, using selection keywords such as \texttt{atomNumbers}.}
2775 The complete list of selection keywords available in \MDENGINE{} is:
2776
2777 \begin{itemize}
2778
2779 \item %
2780   \key
2781     {atomNumbers}{%
2782     atom group}{%
2783     List of atom numbers}{%
2784     space-separated list of positive integers}{%
2785     This option adds to the group all the atoms whose numbers are in
2786     the list.  \emph{The number of the first atom in the system is 1: to convert from a VMD selection, use ``atomselect get serial''.}
2787   }
2788
2789 \item %
2790   \key
2791     {indexGroup}{%
2792     atom group}{%
2793     Name of index group to be used (GROMACS format)}{%
2794     string}{%
2795     If the name of an index file has been provided by \texttt{indexFile}, this option allows to select one index group from that file: the atoms from that index group will be used to define the current group.}
2796
2797 \item %
2798   \key
2799     {atomsOfGroup}{%
2800     atom group}{%
2801     Name of group defined previously}{%
2802     string}{%
2803     Refers to a group defined previously using its user-defined \texttt{name}.
2804     This adds all atoms of that named group to the current group.}
2805
2806 \item %
2807   \key
2808     {atomNumbersRange}{%
2809     atom group}{%
2810     Atoms within a number range}{%
2811     $<$Starting number$>$-$<$Ending number$>$}{%
2812     This option includes in the group all atoms whose numbers are within the range specified.  \emph{The number of the first atom in the system is 1.}
2813   }
2814
2815 \cvnamebasedonly{
2816 \item %
2817   \key
2818     {atomNameResidueRange}{%
2819     atom group}{%
2820     Named atoms within a range of residue numbers}{%
2821     $<$Atom name$>$ $<$Starting residue$>$-$<$Ending residue$>$}{%
2822     This option adds to the group all the atoms with the provided
2823     name, within residues in the given range.}
2824
2825 \item %
2826   \key
2827     {psfSegID}{%
2828     atom group}{%
2829     PSF segment identifier}{%
2830     space-separated list of strings (max 4 characters)}{%
2831     This option sets the PSF segment identifier for
2832     \texttt{atomNameResidueRange}.  Multiple values may be provided,
2833     which correspond to multiple instances of
2834     \texttt{atomNameResidueRange}, in order of their occurrence.
2835     This option is only necessary if a PSF topology file is used.}
2836
2837 \item %
2838   \key
2839     {atomsFile}{%
2840     atom group}{%
2841     PDB file name for atom selection}{%
2842     UNIX filename}{%
2843     This option selects atoms from the PDB file provided and adds them
2844     to the group according to numerical flags in the column
2845     \texttt{atomsCol}.  \textbf{Note:} \emph{the sequence of atoms in the PDB file
2846     provided must match that in the system's topology}.}
2847
2848 \item %
2849   \key
2850     {atomsCol}{%
2851     atom group}{%
2852     PDB column to use for atom selection flags}{%
2853     \texttt{O}, \texttt{B}, \texttt{X}, \texttt{Y}, or \texttt{Z}}{%
2854     This option specifies which PDB column in \texttt{atomsFile} is used to determine which atoms are to be included in the group.
2855   }
2856
2857 \item %
2858   \key
2859     {atomsColValue}{%
2860     atom group}{%
2861     Atom selection flag in the PDB column}{%
2862     positive decimal}{%
2863     If defined, this value in \texttt{atomsCol} identifies atoms in \texttt{atomsFile} that are included in the group.
2864     If undefined, all atoms with a non-zero value in \texttt{atomsCol} are included.}
2865 }
2866
2867 \item %
2868   \key
2869     {dummyAtom}{%
2870     atom group}{%
2871     Dummy atom position (\AA{})}{%
2872     \texttt{(x, y, z)} triplet}{%
2873     Instead of selecting any atom, this option makes the group a virtual particle at a fixed position in space.  This is useful e.g.~to replace a group's center of geometry with a user-defined position.}
2874
2875 \end{itemize}
2876
2877 \cvsubsec{Moving frame of reference.}{sec:colvar_atom_groups_ref_frame}
2878
2879 The following options define an automatic calculation of an optimal translation (\texttt{centerReference}) or optimal rotation (\texttt{rotateReference}), that superimposes the positions of this group to a provided set of reference coordinates.
2880 This can allow, for example, to effectively remove from certain colvars the effects of molecular tumbling and of diffusion.
2881 Given the set of atomic positions $\mathbf{x}_{i}$, the colvar $\xi$ can be defined on a set of roto-translated positions $\mathbf{x}_{i}' = R(\mathbf{x}_{i} - \mathbf{x}^{\mathrm{C}}) + \mathbf{x}^{\mathrm{ref}}$.
2882 $\mathbf{x}^{\mathrm{C}}$ is the geometric center of the $\mathbf{x}_{i}$, $R$ is the optimal rotation matrix to the reference positions and $\mathbf{x}^{\mathrm{ref}}$ is the geometric center of the reference positions.
2883
2884 Components that are defined based on pairwise distances are naturally invariant under global roto-translations.
2885 Other components are instead affected by global rotations or translations: however, they can be made invariant if they are expressed in the frame of reference of a chosen group of atoms, using the \texttt{centerReference} and \texttt{rotateReference} options.
2886 Finally, a few components are defined by convention using a roto-translated frame (e.g. the minimal RMSD): for these components, \texttt{centerReference} and \texttt{rotateReference} are enabled by default.
2887 In typical applications, the default settings result in the expected behavior.
2888
2889 \paragraph*{Warning on rotating frames of reference and periodic boundary conditions.}
2890 \texttt{rotateReference} affects coordinates that depend on minimum-image distances in periodic boundary conditions (PBC).
2891 After rotation of the coordinates, the periodic cell vectors become irrelevant: the rotated system is effectively non-periodic.
2892 A safe way to handle this is to ensure that the relevant inter-group distance vectors remain smaller than the half-size of the periodic cell.
2893 If this is not desirable, one should avoid the rotating frame of reference, and apply orientational restraints to the reference group instead, in order to keep the orientation of the reference group consistent with the orientation of the periodic cell.
2894
2895 \paragraph*{Warning on rotating frames of reference and ABF.}
2896 Note that \texttt{centerReference} and \texttt{rotateReference} may affect the Jacobian derivative of colvar components in a way that is not taken into account by default.
2897 Be careful when using these options in ABF simulations or when using total force values.
2898
2899 \begin{itemize}
2900
2901 \item %
2902   \keydef
2903     {centerReference}{%
2904     atom group}{%
2905     Implicitly remove translations for this group}{%
2906     boolean}{%
2907     \texttt{off}}{%
2908     If this option is \texttt{on}, the center of geometry of the group will be aligned with that of the reference positions provided by \cvnamebasedonly{either} \texttt{refPositions} or \texttt{refPositionsFile}.
2909     Colvar components will only have access to the aligned positions.
2910 \textbf{Note}: unless otherwise specified, \texttt{rmsd} and \texttt{eigenvector} set this option to \texttt{on} \emph{by default}.
2911 }
2912
2913 \item %
2914   \keydef
2915     {rotateReference}{%
2916     atom group}{%
2917     Implicitly remove rotations for this group}{%
2918     boolean}{%
2919     \texttt{off}}{%
2920     If this option is \texttt{on}, the coordinates of this group will be optimally superimposed to the reference positions provided by \cvnamebasedonly{either} \texttt{refPositions} or \texttt{refPositionsFile}.
2921     The rotation will be performed around the center of geometry if \texttt{centerReference} is \texttt{on}, or around the origin otherwise.
2922     The algorithm used is the same employed by the \texttt{orientation} colvar component~\cite{Coutsias2004}.
2923     Forces applied to the atoms of this group will also be implicitly rotated back to the original frame.
2924     \textbf{Note}: unless otherwise specified, \texttt{rmsd} and \texttt{eigenvector} set this option to \texttt{on} \emph{by default}.
2925 }
2926
2927 \item %
2928   \labelkey{atom-group|refPositions}
2929   \key
2930     {refPositions}{%
2931     atom group}{%
2932     Reference positions for fitting (\AA)}{%
2933     space-separated list of \texttt{(x, y, z)} triplets}{%
2934     \label{key:colvars:atom_group:refPositions}
2935     This option provides a list of reference coordinates for \texttt{centerReference} and/or \texttt{rotateReference}, and is mutually exclusive with \texttt{refPositionsFile}.
2936     If only \texttt{centerReference} is \texttt{on}, the list may contain a single (x, y, z) triplet; if also \texttt{rotateReference} is \texttt{on}, the list should be as long as the atom group, and \emph{its order must match the order in which atoms were defined}.
2937 }
2938
2939 \item %
2940   \labelkey{atom-group|refPositionsFile}
2941   \key
2942     {refPositionsFile}{%
2943     atom group}{%
2944     File containing the reference positions for fitting}{%
2945     UNIX filename}{%
2946     \label{key:colvars:atom_group:refPositionsFile}
2947     This option provides a list of reference coordinates for \texttt{centerReference} and/or \texttt{rotateReference}, and is mutually exclusive with \texttt{refPositions}.
2948     The acceptable file format is XYZ, which is read in double precision\cvnamebasedonly{, or PDB; \emph{the latter is discouraged if the precision of the reference coordinates is a concern}}.
2949     Atomic positions are read differently depending on the following scenarios:
2950     \textbf{(i)} the file contains exactly as many records as the atoms in the group: all positions are read in sequence;
2951     \textbf{(ii)} (most common case) the file contains coordinates for the entire system: only the positions corresponding to the numeric indices of the atom group are read\cvnamebasedonly{;
2952       \textbf{(iii)} if the file is a PDB file and \texttt{refPositionsCol} is specified, positions are read according to the value of the column \texttt{refPositionsCol} (which may be the same as \texttt{atomsCol})}.
2953     In each case, atoms are read from the file \emph{in order of increasing number}.
2954 }
2955
2956
2957 \cvnamebasedonly{
2958 \item %
2959   \key
2960     {refPositionsCol}{%
2961     atom group}{%
2962     PDB column containing atom flags}{%
2963     \texttt{O}, \texttt{B}, \texttt{X}, \texttt{Y}, or \texttt{Z}}{%
2964     Like \texttt{atomsCol} for \texttt{atomsFile}, indicates which column to use to identify the atoms in \texttt{refPositionsFile} (if this is a PDB file).}
2965
2966 \item %
2967   \key
2968     {refPositionsColValue}{%
2969     atom group}{%
2970     Atom selection flag in the PDB column}{%
2971     positive decimal}{%
2972     Analogous to \texttt{atomsColValue}, but applied to \texttt{refPositionsCol}.}
2973 }
2974
2975 \item %
2976   \labelkey{atom-group|fittingGroup} %
2977   \keydef
2978     {fittingGroup}{%
2979     atom group}{%
2980     Use an alternate set of atoms to define the roto-translation}{%
2981     Block \texttt{fittingGroup \{ ... \}}}{%
2982     This group itself}{%
2983     If either \texttt{centerReference} or \texttt{rotateReference} is defined, this keyword defines an alternate atom group to calculate the optimal roto-translation.
2984     Use this option to define a continuous rotation if the structure of the group involved changes significantly (a typical symptom would be the message ``Warning: discontinuous rotation!'').
2985
2986 \cvnamebasedonly{
2987     The following example illustrates the syntax of \texttt{fittingGroup}: a group called ``\texttt{atoms}'' is defined, including 8 C$_{\alpha}$ atoms of a protein of 100 residues.
2988     An optimal roto-translation is calculated automatically by fitting the C$_{\alpha}$ trace of the rest of the protein onto the coordinates provided by a PDB file.}
2989
2990 {%
2991 \noindent\ttfamily
2992 \# Example: defining a group "atoms", with its coordinates expressed \\
2993 \# on a roto-translated frame of reference defined by a second group\\
2994 atoms \{\\
2995 \\
2996 \-~~psfSegID              PROT\\
2997 \-~~atomNameResidueRange  CA 41-48\\
2998 \\
2999 \-~~centerReference yes\\
3000 \-~~rotateReference yes\\
3001 \-~~fittingGroup \{\\
3002 \-~~~~\# define the frame by fitting the rest of the protein\\
3003 \-~~~~psfSegID              PROT PROT\\
3004 \-~~~~atomNameResidueRange  CA  1-40\\
3005 \-~~~~atomNameResidueRange  CA 49-100\\
3006 \-~~\} \\
3007 \-~~refPositionsFile all.pdb  \# can be the entire system\\
3008 \}\\}
3009 }
3010 \end{itemize}
3011
3012 The following two options have default values appropriate for the vast majority of applications, and are only provided to support rare, special cases.
3013 \begin{itemize}
3014
3015 \item %
3016   \keydef
3017     {enableFitGradients}{%
3018     atom group}{%
3019     Include the roto-translational contribution to colvar gradients}{%
3020     boolean}{%
3021     \texttt{on}}{%
3022     When either \texttt{centerReference} or \texttt{rotateReference} is on,
3023     the gradients of some colvars include terms proportional to
3024     $\partial{}R/\partial\mathbf{x}_{i}$ (rotational gradients) and
3025     $\partial\mathbf{x}^{\mathrm{C}}/\partial\mathbf{x}_{i}$ (translational gradients).
3026     By default, these terms are calculated and included in the total gradients;
3027     if this option is set to \texttt{off}, they are neglected.
3028     In the case of a minimum RMSD component, this flag is automatically disabled
3029     because the contributions of those derivatives to the gradients cancel out.
3030 }
3031
3032 \item %
3033   \keydef
3034     {enableForces}{%
3035     atom group}{%
3036     Apply forces from this colvar to this group}{%
3037     boolean}{%
3038     \texttt{on}}{%
3039     If this option is \texttt{off}, no forces are applied the atoms in the group.
3040     Other forces are not affected (i.e. those
3041     from the MD engine, from other colvars, and other external forces).
3042     For dummy atoms, this option is \texttt{off} by default.
3043  }
3044
3045 \end{itemize}
3046
3047
3048 \cvsubsec{Treatment of periodic boundary conditions.}{sec:colvar_atom_groups_wrapping}
3049
3050 \cvnamdonly{
3051  In simulations with periodic boundary conditions, NAMD maintains
3052   the coordinates of all the atoms within a molecule contiguous to
3053   each other (i.e.~there are no spurious ``jumps'' in the molecular
3054   bonds).  The Colvars module relies on this when calculating a group's
3055   center of geometry, but this condition may fail if the group spans
3056   different molecules.  In that case, writing the NAMD output and restart files
3057   using \texttt{wrapAll} or \texttt{wrapWater} could produce wrong results
3058   when a simulation run is continued from a previous one.  
3059   The user should then determine, according to which
3060   type of colvars are being calculated, whether \texttt{wrapAll} or
3061   \texttt{wrapWater} can be enabled.
3062
3063   In general, internal coordinate wrapping by NAMD does not affect the calculation of colvars if each atom group satisfies one or more of the following:
3064 }
3065 \cvlammpsonly{
3066 In simulations with periodic boundary conditions, many of the implemented colvar components rely on the fact that each position within a group of atoms is at the nearest periodic image from the center of geometry of the group itself.
3067 However, due to the internal wrapping of individual atomic positions done by LAMMPS, this assumption is broken if the group straddles one of the unit cell's boundaries.
3068 For this reason, within the Colvars module all coordinates are unwrapped by default to avoid discontinuities (see \texttt{unwrap} keyword in \ref{sec:colvars_mdengine_parameters}).
3069
3070 The user should determine whether maintaining the default value of \texttt{unwrap}, depending on the specifics of each system.
3071 In general, internal coordinate wrapping by LAMMPS does not affect the calculation of colvars if each atom group satisfies one or more of the following:
3072 }
3073 \cvvmdonly{
3074   When periodic boundary conditions are defined, the Colvars module requires that the coordinates of each molecular fragment are contiguous, without ``jumps'' when a fragment is partially wrapped near a periodic boundary.
3075   The Colvars module relies on this assumption when calculating a group's center of geometry, but the condition may fail if the group spans different molecules.
3076   In general, coordinate wrapping does not affect the calculation of colvars if each atom group satisfies one or more of the following:
3077 }
3078
3079 \begin{enumerate}
3080   \item[\emph{i)}] it is composed by only one atom;
3081   \item[\emph{ii)}] it is used by a colvar component which does not make use of its center of geometry, but only of pairwise distances (\texttt{distanceInv}, \texttt{coordNum}, \texttt{hBond}, \texttt{alpha}, \texttt{dihedralPC});
3082   \item[\emph{iii)}]  it is used by a colvar component that ignores the ill-defined Cartesian components of its center of mass (such as the $x$ and $y$ components of a membrane's center of mass modeled with \texttt{distanceZ})\cvnamdonly{;
3083   \item[\emph{iv)}] it has all of its atoms within the same molecular fragment%
3084 }.
3085 \end{enumerate}
3086 \cvvmdonly{If none of these conditions are met, wrapping may affect the calculation of collective variables: a possible solution is to use \texttt{pbc wrap} or \texttt{pbc unwrap} prior to processing a trajectory with the Colvars module.}
3087
3088
3089 \cvsubsec{Performance of a Colvars calculation based on group size.}{sec:colvar_atom_groups_scaling}
3090
3091 In simulations performed with message-passing programs (such as NAMD or LAMMPS), the calculation of energy and forces is distributed (i.e., parallelized) across multiple nodes, as well as over the processor cores of each node.
3092 When Colvars is enabled, certain atomic coordinates are collected on a single node, where the calculation of collective variables and of their biases is executed.
3093 This means that for simulations over large numbers of nodes, a Colvars calculation may produce a significant overhead, coming from the costs of transmitting atomic coordinates to one node and of processing them.
3094 \cvnamdonly{The latency-tolerant design and dynamic load balancing of NAMD may alleviate both factors, but a noticeable performance impact may be observed.}
3095
3096 Performance can be improved in multiple ways:
3097 \begin{itemize}
3098 \item The calculation of variables, components and biases can be distributed over the processor cores of the node where the Colvars module is executed.
3099   Currently, an equal weight is assigned to each colvar, or to each component of those colvars that include more than one component.
3100   The performance of simulations that use many colvars or components is improved automatically.
3101   For simulations that use a single large colvar, it may be advisable to partition it in multiple components, which will be then distributed across the available cores.
3102   \cvnamdonly{In NAMD, this feature is enabled in all binaries compiled using SMP builds of Charm++ with the CkLoop extension.}
3103   \cvlammpsonly{In LAMMPS, this feature is supported automatically when LAMMPS is compiled with OpenMP support.}
3104   If printed, the message ``SMP parallelism is available.'' indicates the availability of the option\cvvmdonly{ (will be supported in a future relase of VMD)}.
3105   If available, the option is turned on by default, but may be disabled using the keyword \refkey{smp}{Colvars-global|smp} if required for debugging.
3106
3107 \cvnamdonly{
3108   % Use the following command to identify them:
3109   % grep -B10 'provide(f_cvc_com_based' * |grep '\:\:'|grep '(std::string const &conf)'
3110 \item NAMD also offers a parallelized calculation of the centers of mass of groups of atoms.
3111   This option is on by default for all components that are simple functions of centers of mass, and is controlled by the keyword \refkey{scalable}{sec:cvc_common}.
3112   When supported, the message ``Will enable scalable calculation for group \ldots'' is printed for each group.
3113 }
3114
3115 \item As a general rule, the size of atom groups should be kept relatively small (up to a few thousands of atoms, depending on the size of the entire system in comparison).
3116 To gain an estimate of the computational cost of a large colvar, one can use a test calculation of the same colvar in VMD (hint: use the \texttt{time} Tcl command to measure the cost of running \texttt{cv update}).
3117 \end{itemize}
3118
3119
3120 \cvsec{Biasing and analysis methods}{sec:colvarbias}
3121
3122 All of the biasing and analysis methods implemented recognize the following options:
3123 \begin{itemize}
3124
3125 \item %
3126   \labelkey{colvarbias|name}
3127   \keydef
3128     {name}{%
3129     colvar bias}{%
3130     Identifier for the bias}{%
3131     string}{%
3132     \texttt{$<$type of bias$><$bias index$>$}}{%
3133     This string is used to identify the bias or analysis method in
3134     output messages and to name some output files.}
3135
3136 \item %
3137   \labelkey{colvarbias|colvars}
3138   \key
3139     {colvars}{%
3140     colvar bias}{%
3141     Collective variables involved}{%
3142     space-separated list of colvar names}{%
3143     This option selects by name all the colvars to which this bias or
3144     analysis will be applied.}
3145
3146 \item %
3147   \labelkey{colvarbias|outputEnergy}
3148   \keydef
3149     {outputEnergy}{%
3150     colvar bias}{%
3151     Write the current bias energy to the trajectory file}{%
3152     boolean}{%
3153     \texttt{off}}{%
3154     If this option is chosen and  \texttt{colvarsTrajFrequency} is not zero, the current value of the biasing energy will be written to the trajectory file during the simulation.
3155 }
3156
3157 \end{itemize}
3158
3159 In addition, restraint biases (\ref{sec:colvarbias_harmonic}, \ref{sec:colvarbias_harmonic_walls}, \ref{sec:colvarbias_linear}, ...) and metadynamics biases (\ref{sec:colvarbias_meta}) offer the following optional keywords, which allow the use of thermodynamic integration (TI) to compute potentials of mean force (PMFs).  In adaptive biasing force (ABF) biases (\ref{sec:colvarbias_abf}) the same keywords are not recognized because their functionality is always included.
3160
3161 \begin{itemize}
3162
3163 \item %
3164   \labelkey{colvarbias|writeTIPMF}
3165   \keydef
3166     {writeTIPMF}{%
3167     colvar bias}{%
3168     Write the PMF computed by thermodynamic integration}{%
3169     boolean}{%
3170     \texttt{off}}{%
3171     If the bias is applied to a variable that supports the calculation of total forces (see \refkey{outputTotalForce}{sec:colvar} and \ref{sec:cvc_sys_forces}), this option allows calculating the corresponding PMF by thermodynanic integration, and writing it to the file \outputName{}\texttt{.$<$name$>$.ti.pmf}, where \texttt{$<$name$>$} is the name of the bias.
3172     The total force includes the forces applied to the variable by all bias, except those from this bias itself.
3173     If any bias applies time-dependent forces besides the one using this option, an error is raised.
3174 }
3175
3176
3177 \item %
3178   \labelkey{colvarbias|writeTISamples}
3179   \keydef
3180     {writeTISamples}{%
3181     colvar bias}{%
3182     Write the free-energy gradient samples}{%
3183     boolean}{%
3184     \texttt{off}}{%
3185     This option allows to compute total forces for use with thermodynamic integration as done by the keyword \refkey{writeTIPMF}{sec:colvarbias}.
3186     The names of the files containing the variables' histogram and mean thermodyunamic forces are \outputName\texttt{.$<$name$>$.ti.count} and \outputName\texttt{.$<$name$>$.ti.grad}, respectively: these can be used by \texttt{abf\_integrate} or similar utility.
3187     This option on by default when \texttt{writeTIPMF} is on, but can be enabled separately if the bias is applied to more than one variable, making not possible the direct integration of the PMF at runtime.
3188     If any bias applies time-dependent forces besides the one using this option, an error is raised.
3189 }
3190
3191 \end{itemize}
3192
3193
3194 \cvsubsec{Adaptive Biasing Force}{sec:colvarbias_abf}
3195
3196 For a full description of the Adaptive Biasing Force method, see
3197 reference~\cite{Darve2008}. For details about this implementation,
3198 see references~\cite{Henin2004} and \cite{Henin2010}. \textbf{When
3199 publishing research that makes use of this functionality, please cite
3200 references~\cite{Darve2008} and \cite{Henin2010}.}
3201
3202 An alternate usage of this feature is the application of custom
3203 tabulated biasing potentials to one or more colvars. See
3204 \texttt{inputPrefix} and \texttt{updateBias} below.
3205
3206 Combining ABF with the extended Lagrangian feature (\ref{sec:colvar_extended})
3207 of the variables produces the extended-system ABF variant of the method
3208 (\ref{sec:eABF}).
3209
3210 ABF is based on the thermodynamic integration (TI) scheme for
3211 computing free energy profiles. The free energy as a function
3212 of a set of collective variables $\bm{\xi}=(\xi_{i})_{i\in[1,n]}$
3213 is defined from the canonical distribution of $\bm{\xi}$, ${\mathcal P}(\bm{\xi})$:
3214
3215 \begin{equation}
3216   \label{eq:free}
3217   A(\bm{\xi}) = -\frac{1}{\beta} \ln {\mathcal P}(\bm{\xi}) + A_0
3218 \end{equation}
3219
3220 In the TI formalism, the free energy is obtained from its gradient, 
3221 which is generally calculated in the form of the average of a force
3222 $\bm{F}_\xi$ exerted on $\bm{\xi}$, taken over an iso-$\bm{\xi}$ surface:
3223
3224 \begin{equation}
3225   \label{eq:gradient}
3226   \bm{\nabla}_\xi A(\bm{\xi}) = \left\langle -\bm{F}_\xi \right\rangle_{\bm{\xi}}
3227 \end{equation}
3228
3229 Several formulae that take the form of~(\ref{eq:gradient}) have been
3230 proposed.  This implementation relies partly on the classic
3231 formulation~\cite{Carter1989}, and partly on a more versatile scheme
3232 originating in a work by Ruiz-Montero et al.~\cite{Ruiz-Montero1997},
3233 generalized by den Otter~\cite{denOtter2000} and extended to multiple
3234 variables by Ciccotti et al.~\cite{Ciccotti2005}.  Consider a system
3235 subject to constraints of the form $\sigma_{k}(\vx) = 0$.  Let
3236 $(\bm{v}_{i})_{i\in[1,n]}$ be arbitrarily chosen vector fields
3237 ($\mathbb{R}^{3N}\rightarrow\mathbb{R}^{3N}$) verifying, for all $i$,
3238 $j$, and $k$:
3239
3240 \begin{eqnarray}
3241 \label{eq:ortho_gradient}
3242 \bm{v}_{i} \cdot \gradx \xi_{j}    & = & \delta_{ij}\\
3243 \label{eq:ortho_constraints}
3244 \bm{v}_{i} \cdot \gradx \sigma_{k} & = & 0
3245 \end{eqnarray}
3246
3247 then the following holds~\cite{Ciccotti2005}:
3248
3249 \begin{equation}
3250 \label{eq:gradient_vector}
3251 \frac{\partial A}{\partial \xi_{i}} = \left\langle \bm{v}_{i} \cdot \gradx V
3252 - k_B T \gradx \cdot \bm{v}_{i} \right\rangle_{\bm{\xi}}
3253 \end{equation}
3254
3255 where $V$ is the potential energy function.
3256 $\bm{v}_{i}$ can be interpreted as the direction along which the force
3257 acting on variable $\xi_{i}$ is measured, whereas the second term in the
3258 average corresponds to the geometric entropy contribution that appears
3259 as a Jacobian correction in the classic formalism~\cite{Carter1989}.
3260 Condition~(\ref{eq:ortho_gradient}) states that the direction along
3261 which the total force on $\xi_{i}$ is measured is orthogonal to the
3262 gradient of $\xi_{j}$, which means that the force measured on $\xi_{i}$
3263 does not act on $\xi_{j}$.
3264
3265 Equation~(\ref{eq:ortho_constraints}) implies that constraint forces
3266 are orthogonal to the directions along which the free energy gradient is
3267 measured, so that the measurement is effectively performed on unconstrained
3268 degrees of freedom.
3269 \cvnamdonly{In NAMD, constraints are typically applied to the lengths of
3270 bonds involving hydrogen atoms, for example in TIP3P water molecules (parameter \texttt{rigidBonds}\cvnamdugonly{, section~\ref{section:rigidBonds}}).}
3271
3272 In the framework of ABF,
3273 ${\bf F}_\xi$ is accumulated in bins of finite size $\delta \xi$,
3274 thereby providing an estimate of the free energy gradient
3275 according to equation~({\ref{eq:gradient}}).
3276 The biasing force applied along the collective variables
3277 to overcome free energy barriers is calculated as:
3278
3279 \begin{equation}
3280   \label{eq:abf}
3281   {\bf F}^{\rm ABF} = \alpha(N_\xi) \times \gradx \widetilde A(\bm{\xi})
3282 \end{equation}
3283
3284 where $\gradx \widetilde A$ denotes the current estimate of the
3285 free energy gradient at the current point $\bm{\xi}$ in the collective
3286 variable subspace, and $\alpha(N_\xi)$ is a scaling factor that is ramped
3287 from 0 to 1 as the local number of samples $N_\xi$ increases
3288 to prevent nonequilibrium effects in the early phase of the simulation,
3289 when the gradient estimate has a large variance.
3290 See the \texttt{fullSamples} parameter below for details.
3291
3292 As sampling of the phase space proceeds, the estimate
3293 $\gradx \widetilde A$ is progressively refined. The biasing
3294 force introduced in the equations of motion guarantees that in
3295 the bin centered around $\bm{\xi}$,
3296 the forces acting along the selected collective variables average
3297 to zero over time. Eventually, as the undelying free energy surface is canceled
3298 by the adaptive bias, evolution of the system along $\bm{\xi}$
3299 is governed mainly by diffusion.
3300 Although this implementation of ABF can in principle be used in 
3301 arbitrary dimension, a higher-dimension collective variable space is likely
3302 to result in sampling difficulties.
3303 Most commonly, the number of variables is one or two.
3304
3305
3306 \cvsubsubsec{ABF requirements on collective variables}{sec:colvarbias_abf_req}
3307
3308 The following conditions must be met for an ABF simulation to be possible and
3309 to produce an accurate estimate of the free energy profile.
3310 Note that these requirements do not apply when using the extended-system
3311 ABF method (\ref{sec:eABF}).
3312
3313 \begin{enumerate}
3314  \item \emph{Only linear combinations} of colvar components can be used in ABF calculations.
3315  \item \emph{Availability of total forces} is necessary. The following colvar components
3316 can be used in ABF calculations:
3317 \texttt{distance}, \texttt{distance\_xy}, \texttt{distance\_z}, \texttt{angle},
3318 \texttt{dihedral}, \texttt{gyration},  \texttt{rmsd} and \texttt{eigenvector}.
3319 Atom groups may not be replaced by dummy atoms, unless they are excluded
3320 from the force measurement by specifying \texttt{oneSiteTotalForce}, if available.
3321  \item \emph{Mutual orthogonality of colvars}. In a multidimensional ABF calculation,
3322 equation~(\ref{eq:ortho_gradient}) must be satisfied for any two colvars $\xi_{i}$ and $\xi_{j}$.
3323 Various cases fulfill this orthogonality condition:
3324 \begin{itemize}
3325  \item $\xi_{i}$ and $\xi_{j}$ are based on non-overlapping sets of atoms.
3326  \item atoms involved in the force measurement on $\xi_{i}$ do not participate in
3327 the definition of $\xi_{j}$. This can be obtained using the option \texttt{oneSiteTotalForce}
3328 of the \texttt{distance}, \texttt{angle}, and \texttt{dihedral} components
3329 (example: Ramachandran angles $\phi$, $\psi$).
3330  \item $\xi_{i}$ and $\xi_{j}$ are orthogonal by construction. Useful cases are the sum and
3331 difference of two components, or \texttt{distance\_z} and \texttt{distance\_xy} using the same axis.
3332 \end{itemize}
3333  \item \emph{Mutual orthogonality of components}: when several components are combined into a colvar,
3334 it is assumed that their vectors $\bm{v}_{i}$ (equation~(\ref{eq:gradient_vector}))
3335 are mutually orthogonal. The cases described for colvars in the previous paragraph apply.
3336 % (example: difference of distances).
3337  \item \emph{Orthogonality of colvars and constraints}: equation~\ref{eq:ortho_constraints} can
3338 be satisfied in two simple ways, if either no constrained atoms are involved in the force measurement
3339 (see point 3 above) or pairs of atoms joined by a constrained bond are part of an \textit{atom group}
3340 which only intervenes through its center (center of mass or geometric center) in the force measurement.
3341 In the latter case, the contributions of the two atoms to the left-hand side of equation~\ref{eq:ortho_constraints}
3342 cancel out. For example, all atoms of a rigid TIP3P water molecule can safely be included in an atom
3343 group used in a \texttt{distance} component.
3344 \end{enumerate}
3345
3346
3347 \cvsubsubsec{Parameters for ABF}{sec:colvarbias_abf_params}
3348
3349 ABF depends on parameters from collective variables to define the grid on which free
3350 energy gradients are computed. In the direction of each colvar, the grid ranges from
3351 \texttt{lowerBoundary} to \texttt{upperBoundary}, and the bin width (grid spacing)
3352 is set by the \refkey{width}{colvar|width} parameter.
3353 The following specific parameters can be set in the ABF configuration block:
3354
3355 \begin{itemize}
3356
3357 \item \dupkey{name}{\texttt{abf}}{sec:colvarbias}{biasing and analysis methods}
3358 \item \dupkey{colvars}{\texttt{abf}}{sec:colvarbias}{biasing and analysis methods}
3359 %\item \dupkey{outputEnergy}{sec:colvarbias}{biasing and analysis methods}
3360
3361 \item \keydef{fullSamples}{\texttt{abf}}{%
3362     Number of samples in a bin prior
3363     to application of the ABF}
3364   {positive integer}
3365   {200}
3366   {To avoid nonequilibrium effects due to large fluctuations of the force exerted along the
3367    colvars, it is recommended to apply a biasing force only after a the estimate has started
3368    converging. If \texttt{fullSamples} is non-zero, the applied biasing force is scaled by a factor
3369    $\alpha(N_\xi)$ between 0 and 1.
3370    If the number of samples $N_\xi$ in the current bin is higher than \texttt{fullSamples},
3371    the factor is one. If it is less than half of \texttt{fullSamples}, the factor is zero and
3372    no bias is applied. Between those two thresholds, the factor follows a linear ramp from
3373    0 to 1: $\alpha(N_\xi) =(2N_\xi/\mathrm{fullSamples})-1$}.
3374
3375 \item \keydef{maxForce}{\texttt{abf}}{%
3376     Maximum magnitude of the ABF force}
3377   {positive decimals (one per colvar)}
3378   {disabled}
3379   {This option enforces a cap on the magnitude of the biasing force effectively applied
3380    by this ABF bias on each colvar. This can be useful in the presence of singularities
3381    in the PMF such as hard walls, where the discretization of the average force becomes
3382    very inaccurate, causing the colvar's diffusion to get ``stuck'' at the singularity.
3383    To enable this cap, provide one non-negative value for each colvar. The unit of force
3384    is \cvnamdonly{kcal/mol}\cvvmdonly{kcal/mol}\cvlammpsonly{the unit of energy specified by \texttt{units}} divided by the colvar unit.}
3385
3386 \item \keydef{hideJacobian}{\texttt{abf}}{%
3387     Remove geometric entropy term from calculated
3388     free energy gradient?}
3389   {boolean}
3390   {\texttt{no}}
3391   {In a few special cases, most notably distance-based variables, an alternate definition of
3392     the potential of mean force is traditionally used, which excludes the Jacobian
3393     term describing the effect of geometric entropy on the distribution of the variable.
3394     This results, for example, in particle-particle potentials of mean force being flat
3395     at large separations.
3396     Setting this parameter to \texttt{yes} causes the output data to follow that convention,
3397     by removing this contribution from the output gradients while
3398     applying internally the corresponding correction to ensure uniform sampling.
3399     It is not allowed for colvars with multiple components.}
3400
3401 \item \keydef{outputFreq}{\texttt{abf}}{%
3402     Frequency (in timesteps) at which ABF data files are refreshed}
3403   {positive integer}
3404   {Colvars module restart frequency}
3405   {The files containing the free energy gradient estimate and sampling histogram
3406     (and the PMF in one-dimensional calculations) are written on disk at the given
3407     time interval.}
3408
3409 \item \keydef{historyFreq}{\texttt{abf}}{%
3410     Frequency (in timesteps) at which ABF history files are
3411   accumulated}
3412   {positive integer}
3413   {0}
3414   {If this number is non-zero, the free energy gradient estimate and sampling histogram
3415     (and the PMF in one-dimensional calculations) are appended to files on disk at
3416     the given time interval. History file names use the same prefix as output files, with
3417     ``\texttt{.hist}'' appended.}
3418
3419 \item \key{inputPrefix}{\texttt{abf}}{%
3420     Filename prefix for reading ABF data}
3421   {list of strings}
3422   {If this parameter is set, for each item in the list, ABF tries to read
3423     a gradient and a sampling files named \texttt{$<$inputPrefix$>$.grad}
3424     and \texttt{$<$inputPrefix$>$.count}. This is done at
3425     startup and sets the initial state of the ABF algorithm.
3426     The data from all provided files is combined appropriately.
3427     Also, the grid definition (min and max values, width) need not be the same
3428     that for the current run. This command is useful to piece together
3429     data from simulations in different regions of collective variable space,
3430     or change the colvar boundary values and widths. Note that it is not
3431     recommended to use it to switch to a smaller width, as that will leave
3432     some bins empty in the finer data grid.
3433     This option is NOT compatible with reading the data from a restart file\cvnamdonly{ (\texttt{colvarsInput} option of the NAMD config file)}\cvvmdonly{ (\texttt{cv load} command)}\cvlammpsonly{ (\texttt{input} keyword of the \texttt{fix ID group-ID colvars} command)}.}
3434
3435 \item \keydef{applyBias}{\texttt{abf}}{%
3436     Apply the ABF bias?}
3437   {boolean}
3438   {\texttt{yes}}
3439   { If this is set to no, the calculation proceeds normally but the adaptive
3440     biasing force is not applied. Data is still collected to compute
3441     the free energy gradient. This is mostly intended for testing purposes, and should
3442     not be used in routine simulations.
3443   }
3444
3445 \item \keydef{updateBias}{\texttt{abf}}{%
3446     Update the ABF bias?}
3447   {boolean}
3448   {\texttt{yes}}
3449   { If this is set to no, the initial biasing force (e.g. read from a restart file or
3450     through \texttt{inputPrefix}) is not updated during the simulation.
3451     As a result, a constant bias is applied. This can be used to apply a custom, tabulated
3452     biasing potential to any combination of colvars. To that effect, one should prepare
3453     a gradient file containing the gradient of the potential to be applied (negative
3454     of the bias force), and a count file containing only values greater than
3455     \texttt{fullSamples}. These files must match the grid parameters of the colvars.
3456   }
3457 \end{itemize}
3458
3459 \cvnamdonly{
3460 \cvsubsubsec{Multiple-replica ABF}{sec:colvarbias_abf_shared}
3461 \label{sec:mw-ABF}
3462
3463 \begin{itemize}
3464 \item \keydef{shared}{\texttt{abf}}{%
3465     Apply multiple-replica ABF, sharing force samples among the replicas?}
3466   {boolean}
3467   {\texttt{no}}
3468   { This is command requires that NAMD be compiled and executed with multiple-replica
3469     support.
3470     If \texttt{shared} is set to yes, the total force samples will be synchronized among all replicas
3471     at intervals defined by \texttt{sharedFreq}.
3472     This implements the multiple-walker ABF scheme described in \cite{Minoukadeh2010}; this
3473     implementation is documented in \cite{Comer2014c}.
3474     Thus, it is as if total force samples among all replicas are
3475     gathered in a single shared buffer, which why the algorithm is referred to as shared ABF.
3476     Shared ABF allows all replicas to benefit from the sampling done by other replicas and can lead to faster convergence of the biasing force.
3477   }
3478
3479 \item \keydef{sharedFreq}{\texttt{abf}}{%
3480     Frequency (in timesteps) at which force samples are synchronized among the replicas}
3481   {positive integer}
3482   {\texttt{outputFreq}}
3483   {
3484   In the current implementation of shared ABF, each replica maintains a separate
3485   buffer of total force samples that determine the biasing force.
3486   Every \texttt{sharedFreq} steps, the replicas communicate the samples that
3487   have been gathered since the last synchronization time, ensuring all replicas
3488   apply a similar biasing force.
3489   }
3490 \end{itemize}
3491 }
3492
3493 \cvsubsubsec{Output files}{sec:colvarbias_abf_output}
3494
3495 The ABF bias produces the following files, all in multicolumn text format:
3496 \begin{itemize}
3497 \item \outputName\texttt{.grad}: current estimate of the free energy gradient (grid),
3498   in multicolumn;
3499 \item \outputName\texttt{.count}: histogram of samples collected, on the same grid;
3500 \item \outputName\texttt{.pmf}: only for one-dimensional calculations, integrated
3501   free energy profile or PMF.
3502 \end{itemize}
3503
3504 If several ABF biases are defined concurrently, their name is inserted to produce
3505 unique filenames for output, as in \outputName\texttt{.abf1.grad}.
3506 This should not be done routinely and could lead to meaningless results:
3507 only do it if you know what you are doing!
3508
3509 If the colvar space has been partitioned into sections (\emph{windows}) in which independent
3510 ABF simulations have been run, the resulting data can be merged using the
3511 \texttt{inputPrefix} option described above (a run of 0 steps is enough).
3512
3513
3514 \cvsubsubsec{Post-processing: reconstructing a multidimensional free energy surface}{sec:colvarbias_abf_post}
3515
3516 If a one-dimensional calculation is performed, the estimated free energy
3517 gradient is automatically integrated and a potential of mean force is written
3518 under the file name \texttt{<outputName>.pmf}, in a plain text format that
3519 can be read by most data plotting and analysis programs (e.g. gnuplot).
3520
3521 In dimension 2 or greater, integrating the discretized gradient becomes non-trivial. The
3522 standalone utility \texttt{abf\_integrate} is provided to perform that task.
3523 \texttt{abf\_integrate} reads the gradient data and uses it to perform a Monte-Carlo (M-C)
3524 simulation in discretized collective variable space (specifically, on the same grid
3525 used by ABF to discretize the free energy gradient).
3526 By default, a history-dependent bias (similar in spirit to metadynamics) is used:
3527 at each M-C step, the bias at the current position is incremented by a preset amount
3528 (the \emph{hill height}).
3529 Upon convergence, this bias counteracts optimally the underlying gradient;
3530 it is negated to obtain the estimate of the free energy surface.
3531
3532 \texttt{abf\_integrate} is invoked using the command-line:\\
3533 {\small \noindent\ttfamily
3534 abf\_integrate <gradient\_file> [-n <nsteps>] [-t <temp>] [-m (0|1)] [-h <hill\_height>] [-f <factor>]
3535 }
3536
3537 The gradient file name is provided first, followed by other parameters in any order.
3538 They are described below, with their default value in square brackets:
3539 \begin{itemize}
3540 \setlength{\itemsep}{0pt}
3541 \item \texttt{-n}: number of M-C steps to be performed; by default, a minimal number of
3542 steps is chosen based on the size of the grid, and the integration runs until a convergence
3543 criterion is satisfied (based on the RMSD between the target gradient and the real PMF gradient)
3544 \item \texttt{-t}: temperature for M-C sampling (unrelated to the simulation temperature)
3545   [500~K]
3546 \item \texttt{-m}: use metadynamics-like biased sampling? (0 = false) [1]
3547 \item \texttt{-h}: increment for the history-dependent bias (``hill height'') [0.01~kcal/mol]
3548 \item \texttt{-f}: if non-zero, this factor is used to scale the increment stepwise in the 
3549   second half of the M-C sampling to refine the free energy estimate [0.5]
3550 \end{itemize}
3551
3552 Using the default values of all parameters should give reasonable results in most cases.
3553
3554 \bigskip
3555 \texttt{abf\_integrate} produces the following output files:
3556 \begin{itemize}
3557 \setlength{\itemsep}{0pt}
3558 \item \texttt{<gradient\_file>.pmf}: computed free energy surface
3559 \item \texttt{<gradient\_file>.histo}: histogram of M-C sampling (not
3560 usable in a straightforward way if the history-dependent bias has been applied)
3561 \item \texttt{<gradient\_file>.est}: estimated gradient of the calculated free energy surface
3562 (from finite differences)
3563 \item \texttt{<gradient\_file>.dev}: deviation between the user-provided numerical gradient
3564 and the actual gradient of the calculated free energy surface. The RMS norm of this vector
3565 field is used as a convergence criteria and displayed periodically during the integration.
3566 \end{itemize}
3567
3568 \textbf{Note:} Typically, the ``deviation'' vector field does not
3569 vanish as the integration converges. This happens because the
3570 numerical estimate of the gradient does not exactly derive from a
3571 potential, due to