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