e5dd11e9b77c9377c26143c88a0aa12c5266ca49
[namd.git] / ug / ug_forcefield.tex
1 \section{Force Field Parameters}
2 \label{section:forcefield}
3
4
5 \subsection{Potential energy functions}
6
7 Evaluating the force is the most computationally demanding
8 part of molecular dynamics.
9 The force is the negative gradient of a scalar potential energy function,
10 \begin{equation}
11 \vec{F}(\vec{r}) = -\nabla U(\vec{r}),
12 \end{equation}
13 and, for systems of biomolecules,
14 this potential function involves the summing,
15 \begin{equation}
16 U(\vec{r}) = \sum U_{\text{bonded}}(\vec{r})
17   + \sum U_{\text{nonbonded}}(\vec{r}),
18 \end{equation}
19 over a large number of bonded and nonbonded terms.
20 The bonded potential terms involve 2--, 3--, and 4--body interactions
21 of covalently bonded atoms,
22 with $O(N)$ terms in the summation.
23 The nonbonded potential terms involve interactions
24 between all pairs of atoms
25 (usually excluding pairs of atoms already involved in a bonded term),
26 with $O(N^2)$ terms in the summation,
27 although fast evaluation techniques are used to
28 compute good approximations to their contribution to the potential
29 with $O(N)$ or $O(N \log N)$ computational cost.
30
31
32 \subsubsection{Bonded potential energy terms}
33 %\label{sec:bonded}
34
35 The bonded potential terms involve 2--, 3--, and 4--body interactions
36 of covalently bonded atoms.
37
38 The 2--body spring bond potential
39 describes the harmonic vibrational motion
40 between an $(i,j)$--pair of covalently bonded atoms,
41 \begin{equation}
42 U_{\text{bond}} = k (r_{ij} - r_0)^2,
43 \end{equation}
44 where $r_{ij} = \|\vec{r}_j - \vec{r}_i\|$ gives the distance
45 between the atoms,
46 $r_0$ is the equilibrium distance,
47 and $k$ is the spring constant.
48
49 The 3--body angular bond potential
50 describes the angular vibrational motion
51 occurring between an $(i,j,k)$--triple of covalently bonded atoms,
52 \begin{equation}
53 U_{\text{angle}} = k_{\theta} (\theta - \theta_0)^2
54                  + k_{\text{ub}} (r_{ik} - r_{\text{ub}})^2,
55 \end{equation}
56 where, in the first term,
57 $\theta$ is the angle in radians between vectors
58 $\vec{r}_{ij} = \vec{r}_j - \vec{r}_i$
59 and $\vec{r}_{kj} = \vec{r}_j - \vec{r}_k$,
60 $\theta_0$ is the equilibrium angle,
61 and $k_{\theta}$ is the angle constant.
62 The second term is the Urey--Bradley term
63 used to describe a 
64 (noncovalent) spring between the outer $i$ and $k$ atoms,
65 active when constant $k_{\text{ub}} \neq 0$,
66 where, like the spring bond,
67 $r_{ik} = \|\vec{r}_k - \vec{r}_i\|$ gives the distance between
68 the pair of atoms and
69 $r_{\text{ub}}$ is the equilibrium distance.
70
71 The 4--body torsion angle (also known as dihedral angle) potential
72 describes the angular spring between the planes formed
73 by the first three and last three atoms of
74 a consecutively bonded $(i,j,k,l)$--quadruple of atoms,
75 \begin{equation}
76 U_{\text{tors}} =
77   \begin{cases}
78   k (1 + \cos(n \psi + \phi)) & \text{if $n > 0$,} \\
79   k (\psi - \phi)^2           & \text{if $n = 0$,}
80   \end{cases}
81 \end{equation}
82 where $\psi$ is the angle in radians between
83 the $(i,j,k)$--plane and the $(j,k,l)$--plane.
84 The integer constant $n$ is nonnegative and indicates the periodicity.
85 For $n > 0$, $\phi$ is the phase shift angle
86 and $k$ is the multiplicative constant.
87 For $n = 0$, $\phi$ acts as an equilibrium angle
88 and the units of $k$ change to $\text{potential}/\text{rad}^2$.
89 A given $(i,j,k,l)$--quadruple of atoms might contribute
90 multiple terms to the potential,
91 each with its own parameterization.
92 The use of multiple terms for a torsion angle allows for
93 complex angular variation of the potential,
94 effectively a truncated Fourier series.
95
96
97 \subsubsection{Nonbonded potential energy terms}
98 %\label{sec:nonbonded}
99
100 The nonbonded potential terms involve interactions
101 between all $(i,j)$--pairs of atoms,
102 usually excluding pairs of atoms already involved in a bonded term.
103 Even using a fast evaluation methods
104 the cost of computing the nonbonded potentials dominates the work
105 required for each time step of an MD simulation.
106
107 The Lennard--Jones potential
108 accounts for the weak dipole attraction between distant atoms and
109 the hard core repulsion as atoms become close,
110 \begin{equation}
111 U_{\text{LJ}} = (-E_{\text{min}}) \left[ 
112     \left( \frac{R_{\text{min}}}{r_{ij}} \right)^{12} -
113     2 \left( \frac{R_{\text{min}}}{r_{ij}} \right)^{6} \right],
114 \end{equation}
115 where $r_{ij} = \|\vec{r}_j - \vec{r}_i\|$ gives the distance
116 between the pair of atoms.
117 The parameter $E_{\text{min}} = U_{\text{LJ}}(R_{\text{min}})$ is 
118 the minimum of the potential term
119 ($E_{\text{min}} < 0$, which means that $-E_{\text{min}}$ is the well-depth).
120 The Lennard--Jones potential approaches 0 rapidly as $r_{ij}$
121 increases, so it is usually truncated (smoothly shifted) to 0
122 past a cutoff radius, requiring $O(N)$ computational cost.
123
124 The electrostatic potential
125 is repulsive for atomic charges with the same sign
126 and attractive for atomic charges with opposite signs,
127 \begin{equation}
128 U_{\text{elec}} = \epsilon_{14} \frac{C q_i q_j}{\epsilon_0 r_{ij}},
129 \end{equation}
130 where $r_{ij} = \|\vec{r}_j - \vec{r}_i\|$ gives the distance
131 between the pair of atoms,
132 and $q_i$ and $q_j$ are the charges on the respective atoms.
133 Coulomb's constant $C$ and the dielectric constant $\epsilon_0$
134 are fixed for all electrostatic interactions.
135 The parameter $\epsilon_{14}$ is a unitless scaling factor 
136 whose value is 1,
137 except for a modified 1--4 interaction,
138 where the pair of atoms is separated by a sequence
139 of three covalent bonds (so that the atoms might
140 also be involved in a torsion angle interaction),
141 in which case $\epsilon_{14} = \varepsilon$,
142 for a fixed constant $0 \leq \varepsilon \leq 1$.
143 Although the electrostatic potential may be computed with
144 a cutoff like the Lennard--Jones potential,
145 the $1/r$ potential approaches 0 much more
146 slowly than the $1/r^6$ potential,
147 so neglecting the long range electrostatic terms
148 can degrade qualitative results,
149 especially for highly charged systems.
150 There are other fast evaluation methods that approximate
151 the contribution to the long range electrostatic terms
152 that require $O(N)$ or $O(N \log N)$ computational cost,
153 depending on the method.
154
155
156 \subsection{Non-bonded interactions}
157 \label{section:electdesc}
158
159 \NAMD\ has a number of options that control the way that non-bonded
160 interactions are calculated.  These options are interrelated and
161 can be quite confusing, so this section attempts to explain the
162 behavior of the non-bonded interactions and how to use these
163 parameters.
164
165 \subsubsection{Van der Waals interactions}
166 The simplest non-bonded 
167 interaction is the van der Waals interaction.  In 
168 \NAMD, van der Waals interactions are always truncated at the 
169 cutoff distance, specified by {\tt cutoff}.  
170 The main option that effects van der Waals interactions
171 is the {\tt switching} parameter.  With this option set to {\tt on},
172 a smooth switching function will be used to truncate the
173 van der Waals potential energy smoothly at the cutoff distance.  
174 A graph of the van der Waals 
175 potential with this switching function is shown in Figure 
176 \ref{fig:switching}.  If {\tt switching} is set to {\tt off}, the 
177 van der Waals energy is just abruptly truncated at the cutoff 
178 distance, so that energy may not be conserved.  
179
180 \begin{figure}[htb]
181   \center{\includegraphics{figures/switching}}
182   \caption[Graph of van der Waals potential with and without switching]
183   {\small Graph of van der Waals potential with and without the
184   application of the switching function.  With the switching function
185   active, the potential is smoothly reduced to 0 at the cutoff distance.
186   Without the switching function, there is a discontinuity where the
187   potential is truncated.}
188   \label{fig:switching}
189 \end{figure}
190
191 The switching function used is based on the X-PLOR switching
192 function.  The parameter {\tt switchdist} specifies the distance
193 at which the switching function should start taking effect to
194 bring the van der Waals potential to 0 smoothly at the cutoff distance.  
195 Thus, the value of {\tt switchdist} must always be less than that 
196 of {\tt cutoff}.
197
198
199 \subsubsection{Electrostatic interactions}
200 The handling of electrostatics is slightly
201 more complicated due to the incorporation of multiple timestepping for full
202 electrostatic interactions.  There are two cases to consider, one where
203 full electrostatics is employed and the other where electrostatics
204 are truncated at a given distance.
205 \prettypar
206 First let us consider the latter case, where electrostatics are truncated at
207 the cutoff distance.  Using this scheme, all electrostatic interactions
208 beyond a specified distance are ignored, or assumed to be zero.  If
209 {\tt switching} is set to {\tt on}, rather than having a discontinuity
210 in the potential
211 at the cutoff distance, a shifting function is applied to the electrostatic
212 potential as shown in Figure \ref{fig:shifting}.  As this figure shows, the
213 shifting function shifts the entire potential curve so that the curve
214 intersects the x-axis at the cutoff distance.  This shifting function
215 is based on the
216 shifting function used by X-PLOR.
217
218 \begin{figure}[htb]
219   \center{\includegraphics{figures/shifting}}
220   \caption[Graph of electrostatic potential with and without shifting function]
221   {\small Graph showing an electrostatic potential with and without the
222   application of the shifting function.}
223   \label{fig:shifting}
224 \end{figure}
225
226 Next, consider the case where full electrostatics are calculated.  In this
227 case, the electrostatic interactions are not truncated at any distance.  In
228 this scheme, the {\tt cutoff} parameter has a slightly different meaning
229 for the electrostatic interactions --- it represents
230 the {\it local interaction distance\/}, or distance within which electrostatic
231 pairs will be directly calculated every timestep.  Outside of this distance,
232 interactions will be calculated only periodically.  These forces
233 will be applied using a multiple timestep integration scheme as described in
234 Section \ref{section:mts}.
235
236 \begin{figure}[htb]
237   \center{\includegraphics{figures/fmaOn}}
238   \caption[Graph of electrostatic split between short and long range forces]
239   {\small Graph showing an electrostatic potential 
240   when full electrostatics are used within \NAMD, 
241   with one curve portion calculated directly 
242   and the other calculated using PME.}
243   \label{fig:fmaOn}
244 \end{figure}
245
246
247
248 \subsubsection{Non-bonded force field parameters}
249
250 \begin{itemize}
251 %\item
252 %\NAMDCONF{eleccutoff}%
253 %{local interaction distance for electrostatic calculations (\AA)}%
254 %{positive decimal}%
255 %{Specifies the distance at which
256 %electrostatic interactions are truncated.
257 %If {\tt eleccutoff} is defined, it supersedes {\tt cutoff}.
258 %If {\tt eleccutoff} is not defined, then \verb }cutoff} {\em must}
259 %be defined.
260 %See Section \ref{section:electdesc} for a further discussion
261 %of this configuration value.}
262
263 %\item
264 %\NAMDCONF{vdwcutoff}%
265 %{local interaction distance for van der Waals calculations (\AA)}%
266 %{positive decimal}%
267 %{This value specifies the distance at which
268 %van der Waals interactions are truncated.
269 %If {\tt vdwcutoff} is defined, it supersedes {\tt cutoff}.
270 %If {\tt vdwcutoff} is not defined, then \verb }cutoff} {\em must}
271 %be defined.
272 %See Section \ref{section:electdesc} for a further discussion
273 %of this configuration value.}
274
275 \item
276 \NAMDCONF{cutoff}
277 {local interaction distance common to both electrostatic 
278 and van der Waals calculations (\AA)}
279 {positive decimal}
280 {%This value can substitute for either {\tt eleccutoff}
281 %or {\tt vdwcutoff} if either of those is undefined.
282 See Section \ref{section:electdesc} for more information.}
283
284 \item
285 \NAMDCONFWDEF{switching}{use switching function?}{{\tt on} or {\tt off}}
286 {{\tt on}}
287 {If {\tt switching} is
288 specified to be {\tt off}, then a truncated cutoff is performed.
289 If {\tt switching} is turned {\tt on}, then smoothing functions
290 are applied to both the electrostatics and van der Waals forces.
291 For a complete description of the non-bonded force parameters see
292 Section \ref{section:electdesc}.  If {\tt switching} is set to
293 {\tt on}, then {\tt switchdist} must also be defined.}
294
295 \item
296 \NAMDCONFWDEF{vdwForceSwitching}{use force switching for VDW?}{{\tt on} or {\tt off}}
297 {{\tt off}}
298 {If both {\tt switching} and {\tt vdwForceSwitching} are set to {\tt on},
299 then CHARMM force switching is used for van der Waals forces.
300 }
301
302 %\item
303 %\NAMDCONF{elecswitchdist}{distance at which to activate switching function for %electrostatic calculations (\AA)}{positive decimal $\leq$ {\tt eleccutoff}}
304 %{Distance at which the switching function
305 %used to smooth the truncation of
306 %electrostatic forces should begin to take effect.  
307 %This parameter only has meaning if {\tt switching} is 
308 %set to {\tt on}.  
309 %The value of {\tt elecswitchdist} must be less than
310 %or equal to the value of {\tt eleccutoff}, since the switching function
311 %is only applied on the range from {\tt elecswitchdist} to {\tt eleccutoff}.
312 %If {\tt elecswitchdist} is defined, it supersedes {\tt switchdist}.
313 %If {\tt elecswitchdist} is not defined and {\tt switching} is
314 %{\tt on}, then \verb }switchdist} {\em must} be defined.
315 %For a complete description of the non-bonded force parameters, see
316 %Section \ref{section:electdesc}.
317 %}
318
319 %\item
320 %\NAMDCONF{vdwswitchdist}%
321 %{distance at which to activate switching function 
322 %for van der Waals calculations (\AA)}%
323 %{positive decimal $\leq$ {\tt vdwcutoff}}%
324 %{Distance at which the switching function
325 %used to smooth the truncation of
326 %van der Waals forces should begin to take effect.  
327 %This parameter only has meaning if {\tt switching} is 
328 %set to {\tt on}.  
329 %The value of {\tt vdwswitchdist} must be less than
330 %or equal to the value of {\tt vdwcutoff}, since the switching function
331 %is only applied on the range from {\tt vdwswitchdist} to {\tt vdwcutoff}.
332 %If {\tt vdwswitchdist} is defined, it supersedes {\tt switchdist}.
333 %If {\tt vdwswitchdist} is not defined and {\tt switching} is
334 %{\tt on}, then \verb }switchdist} {\em must} be defined.
335 %For a complete description of the non-bonded force parameters, see
336 %Section \ref{section:electdesc}.
337 %}
338
339 \item
340 \NAMDCONF{switchdist}
341 {distance at which to activate switching/splitting function 
342 for electrostatic and van der Waals calculations (\AA)}
343 {positive decimal $\leq$ {\tt cutoff}}
344 {Distance at which the switching function
345 should begin to take effect.  
346 This parameter only has meaning if {\tt switching} is 
347 set to {\tt on}.  
348 The value of {\tt switchdist} must be less than
349 or equal to the value of {\tt cutoff}, since the switching function
350 is only applied on the range from {\tt switchdist} to {\tt cutoff}.  
351 For a complete description of the non-bonded force parameters see
352 Section \ref{section:electdesc}.}
353
354 \item
355 \NAMDCONF{exclude}
356 {non-bonded exclusion policy to use}
357 {{\tt none}, {\tt 1-2}, {\tt 1-3}, {\tt 1-4}, or {\tt scaled1-4}}
358 {\label{param:exclude}
359 %% This parameter is {\it required\/} for every simulation.
360 This parameter specifies which pairs of bonded atoms should
361 be excluded from non-bonded
362 interactions.  With the value of {\tt none}, no bonded pairs of atoms 
363 will be excluded.  With the value of {\tt 1-2}, all atom pairs that
364 are directly connected via a linear bond will be excluded.  With the
365 value of {\tt 1-3}, all {\tt 1-2} pairs will be excluded along with
366 all pairs of atoms that are bonded to a common
367 third atom (i.e., if atom A is bonded to atom B and atom B is bonded
368 to atom C, then the atom pair A-C would be excluded).
369 With the value of {\tt 1-4}, all {\tt 1-3} pairs will be excluded along
370 with all pairs connected by a set of two bonds (i.e., if atom A is bonded
371 to atom B, and atom B is bonded to atom C, and atom C is bonded to
372 atom D, then the atom pair A-D would be excluded).  With the value
373 of {\tt scaled1-4}, all {\tt 1-3} pairs are excluded and all pairs
374 that match the {\tt 1-4} criteria are modified.  The electrostatic
375 interactions for such pairs are modified by the constant factor
376 defined by {\tt 1-4scaling}.  
377 The van der Waals interactions are modified
378 by using the special 1-4 parameters defined in the parameter files.
379 The value of {\tt scaled1-4} is necessary to enable the modified
380 1-4 VDW parameters present in the CHARMM parameter files.
381 }
382
383 \item
384 \NAMDCONFWDEF{1-4scaling}{scaling factor for 1-4 electrostatic interactions}
385 {0 $\leq$ decimal $\leq$ 1}{1.0}
386 {Scaling factor for 1-4 electrostatic interactions.
387 This factor is only used when the
388 {\tt exclude} parameter is set to {\tt scaled1-4}.  In this case, this
389 factor is used to modify the electrostatic interactions between 1-4 atom
390 pairs.  If the {\tt exclude} parameter is set to anything but 
391 {\tt scaled1-4}, this parameter has no effect regardless of its value.}
392
393 \item
394 \NAMDCONFWDEF{dielectric}{dielectric constant for system}
395 {decimal $\geq$ 1.0}{1.0}
396 {Dielectric constant for the system.  A value of 1.0 implies no modification
397 of the electrostatic interactions.  Any larger value will lessen the
398 electrostatic forces acting in the system.}
399
400 \item
401 \NAMDCONFWDEF{nonbondedScaling}{scaling factor for nonbonded forces}
402 {decimal $\geq$ 0.0}{1.0}
403 {Scaling factor for electrostatic and van der Waals forces.
404 A value of 1.0 implies no modification of the interactions.
405 Any smaller value will lessen the
406 nonbonded forces acting in the system.}
407
408 \item
409 \NAMDCONFWDEF{vdwGeometricSigma}{use geometric mean to combine L-J sigmas}
410 {{\tt yes} or {\tt no}}{{\tt no}}
411 {Use geometric mean, as required by \index{OPLS} OPLS, rather than
412 traditional arithmetic mean when combining Lennard-Jones sigma parameters
413 for different atom types.}
414
415 \item
416 \NAMDCONFWDEF{limitdist}
417 {maximum distance between pairs for limiting interaction strength(\AA)}
418 {non-negative decimal}
419 {{\tt 0.}}
420 {
421 The electrostatic and van der Waals potential functions diverge
422 as the distance between two atoms approaches zero.
423 The potential for atoms closer than {\tt limitdist} is instead
424 treated as $a r^2 + c$ with parameters chosen to match the
425 force and potential at {\tt limitdist}.
426 This option should primarily be useful for alchemical free energy
427 perturbation calculations, since it makes the process of creating
428 and destroying atoms far less drastic energetically.
429 The larger the value of {\tt limitdist} the more the maximum force
430 between atoms will be reduced.
431 In order to not alter the other interactions in the simulation,
432 {\tt limitdist} should be less than the closest approach
433 of any non-bonded pair of atoms; 1.3\,\AA\ appears to satisfy this
434 for typical simulations but the user is encouraged to experiment.
435 There should be no performance impact from enabling this feature.
436 }
437
438 \item
439 \NAMDCONFWDEF{LJcorrection}
440 {Apply long-range corrections to the system energy and virial to
441 account for neglected vdW forces?}{{\tt yes} or {\tt no}}{{\tt no}}
442 {Apply an analytical correction to the reported vdW energy and virial
443 that is equal to the amount lost due to switching and cutoff of the LJ
444 potential. The correction will use the average of vdW parameters for
445 all particles in the system and assume a constant, homogeneous
446 distribution of particles beyond the switching distance. See 
447 \cite{Shirts2007} for details (the equations used in the NAMD
448 implementation are slightly different due to the use of a different
449 switching function). Periodic boundary conditions are required to make
450 use of tail corrections.
451 }
452
453 \end{itemize}
454
455
456 \subsubsection{PME parameters}
457
458 PME stands for Particle Mesh Ewald and is an efficient
459 full electrostatics method for use with periodic boundary conditions.
460 None of the parameters should affect energy conservation, although they may affect the accuracy of the results and momentum conservation.
461
462 \begin{itemize}
463
464 \item
465 \NAMDCONFWDEF{PME}{Use particle mesh Ewald for electrostatics?}{{\tt yes} or {\tt no}}{{\tt no}}
466 {Turns on particle mesh Ewald.}
467
468 \item
469 \NAMDCONFWDEF{PMETolerance}{PME direct space tolerance}{positive decimal}{$10^{-6}$}
470 {Affects the value of the Ewald coefficient and the overall accuracy of the results.}
471
472 \item
473 \NAMDCONFWDEF{PMEInterpOrder}{PME interpolation order}{positive integer}{4 (cubic)}
474 {Charges are interpolated onto the grid and forces are interpolated off using this many points, equal to the order of the interpolation function plus one.}
475
476 \item
477 \NAMDCONF{PMEGridSpacing}{maximum space between grid points}{positive real}
478 {The grid spacing partially determines the accuracy and efficiency of PME.
479 If any of the grid sizes below are not set, then PMEGridSpacing must be set
480 (recommended value is 1.0 \AA ) and will be used to calculate them.
481 If a grid size is set, then the grid spacing must be
482 at least PMEGridSpacing (if set, or a very large default of 1.5).}
483
484 \item
485 \NAMDCONF{PMEGridSizeX}{number of grid points in x dimension}{positive integer}
486 {The grid size partially determines the accuracy and efficiency of PME.
487 For speed, {\tt PMEGridSizeX} should have only small integer factors (2, 3 and 5).}
488
489 \item
490 \NAMDCONF{PMEGridSizeY}{number of grid points in y dimension}{positive integer}
491 {The grid size partially determines the accuracy and efficiency of PME.
492 For speed, {\tt PMEGridSizeY} should have only small integer factors (2, 3 and 5).}
493
494 \item
495 \NAMDCONF{PMEGridSizeZ}{number of grid points in z dimension}{positive integer}
496 {The grid size partially determines the accuracy and efficiency of PME.
497 For speed, {\tt PMEGridSizeZ} should have only small integer factors (2, 3 and 5).}
498
499 \item
500 \NAMDCONFWDEF{PMEProcessors}{processors for FFT and reciprocal sum}{positive integer}{larger of x and y grid sizes up to all available processors}
501 {For best performance on some systems and machines, it may be necessary to
502 restrict the amount of parallelism used.  Experiment with this parameter if
503 your parallel performance is poor when PME is used.}
504
505 \item
506 \NAMDCONFWDEF{FFTWEstimate}{Use estimates to optimize FFT?}{{\tt yes} or {\tt no}}{{\tt no}}
507 {Do not optimize FFT based on measurements, but on FFTW rules of thumb.
508 This reduces startup time, but may affect performance.}
509
510 \item
511 \NAMDCONFWDEF{FFTWUseWisdom}{Use FFTW wisdom archive file?}{{\tt yes} or {\tt no}}{{\tt yes}}
512 {Try to reduce startup time when possible by reading FFTW ``wisdom'' from a file, and saving wisdom generated by performance measurements to the same file for future use.
513 This will reduce startup time when running the same size PME grid on the same number of processors as a previous run using the same file.}
514
515 \item
516 \NAMDCONFWDEF{FFTWWisdomFile}{name of file for FFTW wisdom archive}{file name}{FFTW\_NAMD\_{\em version}\_{\em platform}.txt}
517 {File where FFTW wisdom is read and saved.
518 If you only run on one platform this may be useful to reduce startup times for all runs.
519 The default is likely sufficient, as it is version and platform specific.}
520
521 \end{itemize}
522
523
524 \subsubsection{MSM parameters}
525
526 The multilevel summation method (MSM)~\cite{HARD2015}
527 is an alternative to PME for calculating full electrostatic interactions.
528 The use of the FFT in PME has two drawbacks:
529 (1) it generally requires the use of periodic boundary conditions, 
530 in which the simulation describes an infinite three-dimensional lattice,
531 with each lattice cell containing a copy of the simulated system, and
532 (2) calculation of the FFT becomes a considerable performance bottleneck
533 to the parallel scalability of MD simulations, due to the many-to-many
534 communication pattern employed.
535 MSM avoids the use of the FFT in its calculation,
536 instead employing the nested interpolation in real space
537 of softened pair potentials, 
538 which permits in addition to periodic boundary conditions 
539 the use of
540 semi-periodic boundaries, in which there is periodicity along 
541 just one or two basis vectors, 
542 or non-periodic boundaries, in which the simulation is 
543 performed in a vacuum. 
544 Also, better parallel scaling has been observed with MSM 
545 when scaling a sufficiently large system to a large number of processors. 
546 See the MSM research web page (\url{http://www.ks.uiuc.edu/Research/msm/}) 
547 for more information. 
548
549 In order to use the MSM, 
550 one need only specify ``MSM on'' in the configuration file. 
551 For production use, 
552 we presently recommend using the default 
553 ``MSMQuality 0''
554 ($C^1$ cubic interpolation with $C^2$ Taylor splitting), 
555 which has been validated to correctly reproduce
556 the PME results~\cite{HARD2015}. 
557 At this time, we discourage use of the higher order interpolation schemes 
558 (Hermite, quintic, etc.), 
559 as they are still under development. 
560 With cubic interpolation, MSM now gets roughly half the performance of PME. 
561 Comparable performance and better scaling for MSM 
562 have been observed with the optimizations described
563 in Ref.~\cite{HARD2015}, which will be available shortly. 
564
565 For now, \NAMD's implementation of the MSM
566 does not calculate the long-range electrostatic 
567 contribution to the virial, so use with a barostat for 
568 constant pressure simulation is inappropriate. 
569 (Note that the experiments in Ref.~\cite{HARD2015} 
570 involving constant pressure simulation with MSM 
571 made use of a custom version that is incompatible with 
572 some other \NAMD\ features, so is not yet available.) 
573 The performance of PME is generally still better for smaller systems 
574 with smaller processor counts. 
575 MSM is the only efficient method in \NAMD\ for calculating 
576 full electrostatics for simulations with semi-periodic or
577 non-periodic boundaries. 
578
579 The periodicity is defined through setting the cell basis vectors 
580 appropriately, as discussed in Sec.~\ref{section:dynamics}. 
581 The cutoff distance, discussed earlier in this section, 
582 also determines the splitting distance between the 
583 MSM short-range part, calculated exactly, and long-range part, 
584 interpolated from the grid hierarchy; 
585 this splitting distance is the primary control for 
586 accuracy for a given interpolation and splitting, 
587 although most simulations will likely want to keep the 
588 cutoff set to the CHARMM-prescribed value of 12~\AA. 
589
590 The configuration options specific to MSM are listed below. 
591 A simulation employing non-periodic boundaries in one or more 
592 dimensions might have atoms that attempt to drift beyond
593 the predetermined extent of the grid. 
594 In the case that an atom does drift beyond the grid, 
595 the simulation will be halted prematurely with an error message. 
596 Several options listed below deal with defining the extent of the 
597 grid along non-periodic dimensions beyond what can be automatically 
598 determined by the initial coordinates. 
599 It is also recommended for non-periodic simulation to 
600 configure boundary restraints to contain the atoms, for instance,
601 through Tcl boundary forces in Sec.~\ref{section:tclBC}. 
602
603
604 \begin{itemize}
605
606 \item
607 \NAMDCONFWDEF{MSM}{Use multilevel summation method for electrostatics?}{{\tt yes} or {\tt no}}{{\tt no}}
608 {Turns on multilevel summation method.}
609
610 \item
611 \NAMDCONFWDEF{MSMGridSpacing}{spacing between finest level grid points (\AA)}{positive real}{2.5}
612 {The grid spacing determines in part the accuracy and efficiency of MSM. 
613 An error versus cost analysis shows that the best tradeoff is setting 
614 the grid spacing to a value close to the inter-particle spacing. 
615 The default value works well in practice for atomic scale simulation.
616 This value will be exact along non-periodic dimensions. 
617 For periodic dimensions, the grid spacing must evenly divide the 
618 basis vector length; the actual spacing for a desired grid spacing $h$ 
619 is guaranteed to be within the interval
620 $\bigl[ \frac{4}{5} h, \frac{6}{5} h \bigr)$.}
621
622 \item
623 \NAMDCONFWDEF{MSMQuality}{select the approximation quality}{$0,1,2,3,4$}{0}
624 {This parameter offers a simplified way to select higher order
625 interpolation and splitting for MSM.  The available choices are: 
626 \begin{itemize}
627 \item 0 sets $C^1$ cubic ($p=3$) interpolation with $C^2$ Taylor splitting,
628 \item 1 sets $C^1$ Hermite ($p=4$) interpolation with $C^3$ Taylor splitting,
629 \item 2 sets $C^1$ quintic ($p=5$) interpolation with $C^3$ Taylor splitting,
630 \item 3 sets $C^1$ septic ($p=7$) interpolation with $C^4$ Taylor splitting,
631 \item 4 sets $C^1$ nonic ($p=9$) interpolation with $C^5$ Taylor splitting.
632 \end{itemize}
633 \emph{We presently recommend using the default selection, 
634 which has been validated to correctly reproduce
635 the PME results~\cite{HARD2015}, 
636 and discourage use of the higher order interpolation schemes, 
637 as they are still under development.}
638 With cubic interpolation, MSM now gets roughly half the performance of PME. 
639 Comparable performance and better scaling for MSM 
640 have been observed with the optimizations described
641 in Ref.~\cite{HARD2015}, which will be available shortly. 
642
643 There is generally a tradeoff between quality and performance. 
644 Empirical results show that the $C^1$ interpolation schemes offer a little
645 better accuracy than the alternative
646 interpolation schemes that have greater continuity. 
647 Also, better accuracy has been observed by using 
648 a splitting function with $C^{\lceil (p+1) / 2 \rceil}$ continuity 
649 where $p$ is the order of the interpolant. 
650 }
651
652 \item
653 \NAMDCONFWDEF{MSMApprox}{select the interpolant}{$0,1,\dots,7$}{0}
654 {Select the interpolation scheme:
655 \begin{itemize}
656 \item 0 sets $C^1$ cubic ($p=3$) interpolation,
657 \item 1 sets $C^1$ quintic ($p=5$) interpolation,
658 \item 2 sets $C^2$ quintic ($p=5$) interpolation,
659 \item 3 sets $C^1$ septic ($p=7$) interpolation,
660 \item 4 sets $C^3$ septic ($p=7$) interpolation,
661 \item 5 sets $C^1$ nonic ($p=9$) interpolation,
662 \item 6 sets $C^4$ nonic ($p=9$) interpolation,
663 \item 7 sets $C^1$ Hermite ($p=4$) interpolation.
664 \end{itemize}
665 }
666
667 \item
668 \NAMDCONFWDEF{MSMSplit}{select the splitting}{$0,1,\dots,6$}{0}
669 {Select the splitting function: 
670 \begin{itemize}
671 \item 0 sets $C^2$ Taylor splitting,
672 \item 1 sets $C^3$ Taylor splitting,
673 \item 2 sets $C^4$ Taylor splitting,
674 \item 3 sets $C^5$ Taylor splitting,
675 \item 4 sets $C^6$ Taylor splitting,
676 \item 5 sets $C^7$ Taylor splitting,
677 \item 6 sets $C^8$ Taylor splitting.
678 \end{itemize}
679 }
680
681 \item
682 \NAMDCONFWDEF{MSMLevels}{maximum number of levels}{non-negative integer}{0}
683 {Set the maximum number of levels to use in the grid hierarchy. 
684 Although setting slightly lower than the default might (or might not) 
685 improve performance and/or accuracy for non-periodic simulation, 
686 it is generally best to leave this at the default value "0" which will
687 then automatically adjust the levels to the size of the given system.}
688
689 \item
690 \NAMDCONFWDEF{MSMPadding}{grid padding (\AA)}{non-negative real}{2.5}
691 {The grid padding applies only to non-periodic dimensions, for which 
692 the extent of the grid is automatically determined by the 
693 maximum and minimum of the initial coordinates plus the padding value.}
694
695 \item
696 \NAMDCONF{MSMxmin, MSMymin, MSMzmin}{minimum x-, y-, z-coordinate (\AA)}{real}
697 {Set independently the minimum x-, y-, or z-coordinates of 
698 the simulation.  This parameter is applicable only to non-periodic dimensions. 
699 It is useful in conjunction with setting a boundary restraining force 
700 with Tcl boundary forces in Sec.~\ref{section:tclBC}.}
701
702 \item
703 \NAMDCONF{MSMxmax, MSMymax, MSMzmax}{maximum x-, y-, z-coordinate (\AA)}{real}
704 {Set independently the maximum x-, y-, or z-coordinates of 
705 the simulation.  This parameter is applicable only to non-periodic dimensions. 
706 It is useful in conjunction with setting a boundary restraining force 
707 with Tcl boundary forces in Sec.~\ref{section:tclBC}.}
708
709 \item
710 \NAMDCONFWDEF{MSMBlockSizeX, MSMBlockSizeY, MSMBlockSizeZ}{block size for grid decomposition}{positive integer}{8}
711 {Tune parallel performance by adjusting the block size used for parallel 
712 domain decomposition of the grid.  Recommended to keep the default.}
713
714 \item
715 \NAMDCONFWDEF{MSMSerial}{Use serial long-range solver?}{{\tt yes} or {\tt no}}{{\tt no}}
716 {Enable instead the slow serial long-range solver. 
717 Intended to be used only for testing and diagnostic purposes.}
718
719 \end{itemize}
720
721
722 \subsubsection{Full direct parameters}
723
724 The direct computation of electrostatics 
725 is not intended to be used during 
726 real calculations, but rather as a testing or 
727 comparison measure.  Because of the ${\mathcal O}(N^2)$ 
728 computational complexity for performing 
729 direct calculations, this is {\it much} 
730 slower than using PME or MSM to compute full 
731 electrostatics for large systems.
732 In the case of periodic boundary conditions,
733 the nearest image convention is used rather than a
734 full Ewald sum.
735
736 \begin{itemize}
737
738 \item
739 \NAMDCONFWDEF{FullDirect}{calculate full electrostatics directly?}{{\tt yes} or {\tt no}}{{\tt no}}
740 {Specifies whether or not direct computation of 
741 full electrostatics should be performed.}
742
743 \end{itemize}
744
745
746 \subsubsection{Tabulated nonbonded interaction parameters}
747
748 In order to support coarse grained models and semiconductor force fields,
749 the tabulated energies feature replaces the normal van der Waals potential
750 for specified pairs of atom types with one interpolated from user-supplied
751 energy tables.  The electrostatic potential is not altered.
752
753 Pairs of atom types to which the modified interactions apply are specified
754 in a CHARMM parameter file by an {\tt NBTABLE} section consisting of lines
755 with two atom types and a corresponding interaction type name.
756 For example, tabulated interactions for SI-O, O-O, and SI-SI pairs would
757 be specified in a parameter file as:
758 \begin{verbatim}
759 NBTABLE
760 SI O SIO
761 O O OO
762 SI SI SISI
763 \end{verbatim}
764
765 Each interaction type must correspond to an entry in the energy table file.
766 The table file consists of a header formatted as:
767 \begin{verbatim}
768 # multiple comment lines
769 <number_of_tables> <table_spacing (A)> <maximum_distance (A)>
770 \end{verbatim}
771 followed by {\em number\_of\_tables} energy tables formatted as:
772 \begin{verbatim}
773 TYPE <interaction type name>
774 0 <energy (kcal/mol)> <force (kcal/mol/A)>
775 <table_spacing> <energy (kcal/mol)> <force (kcal/mol/A)>
776 <2*table_spacing> <energy (kcal/mol)> <force (kcal/mol/A)>
777 <3*table_spacing> <energy (kcal/mol)> <force (kcal/mol/A)>
778 ...
779 <maximum_distance - 3*table_spacing> <energy (kcal/mol)> <force (kcal/mol/A)>
780 <maximum_distance - 2*table_spacing> <energy (kcal/mol)> <force (kcal/mol/A)>
781 <maximum_distance - table_spacing> <energy (kcal/mol)> <force (kcal/mol/A)>
782 \end{verbatim}
783
784 The table entry at {\em maximum\_distance} will match the energy of the
785 previous entry but have a force of zero.  The maximum distance must be at
786 least equal to the nonbonded cutoff distance and entries beyond the cutoff
787 distance will be ignored.  For the above example with a cutoff of 12 \AA\
788 the table file could look like:
789 \begin{verbatim}
790 # parameters for silicon dioxide
791 3 0.01 14.0
792 TYPE SIO
793 0 5.092449e+26 3.055469e+31
794 0.01 5.092449e+14 3.055469e+17
795 0.02 7.956951e+12 2.387085e+15
796 0.03 6.985526e+11 1.397105e+14
797 ...
798 13.98 0.000000e+00 -0.000000e+00
799 13.99 0.000000e+00 -0.000000e+00
800 TYPE OO
801 0 1.832907e+27 1.099744e+32
802 0.01 1.832907e+15 1.099744e+18
803 0.02 2.863917e+13 8.591751e+15
804 0.03 2.514276e+12 5.028551e+14
805 ...
806 13.98 0.000000e+00 -0.000000e+00
807 13.99 0.000000e+00 -0.000000e+00
808 TYPE SISI
809 0 0.000000e+00 -0.000000e+00
810 0.01 0.000000e+00 -0.000000e+00
811 ...
812 13.98 0.000000e+00 -0.000000e+00
813 13.99 0.000000e+00 -0.000000e+00
814 \end{verbatim}
815
816 The following three parameters are required for tabulated energies.
817
818 \begin{itemize}
819
820 \item
821 \NAMDCONFWDEF{tabulatedEnergies}{use tabulated energies}{{\tt yes} or {\tt no}}{{\tt no}}
822 {Specifies whether or not tabulated energies will be used for
823 van der Waals interactions between specified pairs of atom types.}
824
825 \item
826 \NAMDCONF{tabulatedEnergiesFile}{file containing energy table}{file name}
827 {Provides one energy table for each interaction type in parameter file.
828 See format above.}
829
830 \item
831 \NAMDCONF{tableInterpType}{cubic or linear interpolation}{{\tt cubic} or {\tt linear}}
832 {Specifies the order for interpolating between energy table entries.}
833
834 \end{itemize}
835
836
837 \subsection{Water Models}
838 \label{section:water_models}
839
840 \NAMD~currently supports the 3-site TIP3P water model, 
841 the 4-site TIP4P water model,
842 and the 5-site SWM4-NDP water model 
843 (from the Drude force field)~\cite{Lamoureux-2006a}.  
844 TIP3P is the current default water model.  
845 Usage of alternative water models is described below. 
846
847 \begin{itemize}
848
849   \item
850     \NAMDCONFWDEF{waterModel}{using which water model?}{%
851 {\tt tip3}, {\tt tip4}, {\tt swm4}}{ {\tt tip3}}
852     {Specifies the water model to be used.
853 When using the TIP3P water model, the ordering of atoms within each
854 TIP3P water molecule must be oxygen, hydrogen, hydrogen.  
855 When using the TIP4P water model, the ordering of atoms within each
856 TIP4P water molecule must be oxygen, hydrogen, hydrogen, lone pair.
857 When using the SWM4-NDP water model, 
858 the ordering of atoms within each SWM4-NDP water molecule
859 must be oxygen, Drude particle, lone pair, hydrogen, hydrogen.  
860 Alternative orderings will fail. } 
861
862 \end{itemize}
863
864
865 \subsection{Drude polarizable force field}
866 \label{section:drude_forcefield}
867
868 The Drude oscillator model represents induced electronic polarization 
869 by introducing an auxiliary particle attached to each polarizable 
870 atom via a harmonic spring.  The advantage with the Drude model is 
871 that it preserves the simple particle-particle Coulomb electrostatic 
872 interaction employed in nonpolarizable force fields, 
873 thus its implementation in \NAMD\ is more straightforward than 
874 alternative models for polarization.  
875 \NAMD\ performs the integration of Drude oscillators 
876 by employing a novel dual Langevin thermostat to freeze the Drude 
877 oscillators while maintaining the warm degrees of freedom 
878 at the desired temperature~\cite{JIAN2011}.  
879 Use of the Langevin thermostat enables better parallel scalability 
880 than the earlier reported implementation which made use of 
881 a dual Nos\'e-Hoover thermostat acting on, and within,
882 each nucleus-Drude pair~\cite{Lamoureux-2003a}.  
883 Performance results 
884 show that the \NAMD\ implementation of the Drude model maintains good 
885 parallel scalability, 
886 with an increase in computational cost by not more than twice 
887 that of using a nonpolarizable force field~\cite{JIAN2011}.  
888
889 The Drude polarizable force field requires some extensions to the
890 CHARMM force field.  
891 The Drude oscillators differ from typical spring bonds only in that they 
892 have an equilibrium length of zero.  
893 The Drude oscillators are optionally supplemented by a maximal 
894 bond length parameter, beyond which a quartic restraining potential 
895 is also applied.  
896 The force field is also extended by an anisotropic spring term 
897 that accounts for out-of-plane forces from a polarized atom and its 
898 covalently bonded neighbor with two more covalently bonded 
899 neighbors (similar in structure to an improper bond).  
900 The screened Coulomb correction of Thole is calculated between 
901 pairs of Drude oscillators that would otherwise be excluded from 
902 nonbonded interaction and 
903 optionally between non-excluded, nonbonded pairs of Drude oscillators 
904 that are within a prescribed cutoff distance~\cite{Thole81,van1998molecular}.  
905
906 Also included in the Drude force field 
907 are the use of off-centered massless interaction sites, 
908 so called ``lone pairs'' (LPs),
909 to avoid the limitations 
910 of centrosymmetric-based Coulomb interactions~\cite{Harder2006}.  
911 The coordinate of each LP site is constructed based on three ``guide'' atoms.
912 The calculated forces on the massless LP must be transferred 
913 to the guide atoms, 
914 preserving total force and torque.  
915 After an integration step of velocities and positions, 
916 the position of the LP is updated based on the three guide atoms, 
917 along with additional geometry parameters that give displacement 
918 and in-plane and out-of-plane angles.  
919 See our research web page
920 (\url{http://www.ks.uiuc.edu/Research/Drude/})
921 for additional details and parallel performance results.  
922
923
924 \subsubsection{Required input files}
925
926 No additional files are required by \NAMD\ to use the 
927 Drude polarizable force field.  
928 However, it is presently beyond the capability of 
929 the Psfgen tool to generate the PSF file needed to perform 
930 a simulation using the Drude model.  
931 For now, CHARMM is needed to generate correct input files.  
932
933 The CHARMM force field parameter files
934 specific to the Drude model are required.  
935 The PDB file must also include 
936 the Drude particles (mass between 0.1 and 1.0) 
937 and the LPs (mass 0).  
938 The Drude particles always immediately follow their parent atom.  
939 The PSF file augments the ``atom'' section with 
940 additional columns that include the 
941 ``Thole'' and ``alpha'' parameters for the screened Coulomb 
942 interactions of Thole.  
943 The PSF file also requires additional sections that 
944 list the LPs, including their guide atoms and geometry parameters, 
945 and list the anisotropic interaction terms, including their parameters.  
946 A Drude-compatible PSF file is denoted by the 
947 keyword ``DRUDE'' given along the top line.  
948
949
950 \subsubsection{Standard output}
951
952 The \NAMD\ logging to standard output is extended to provide additional 
953 temperature data on the cold and warm degrees of freedom.  
954 Four additional quantities are listed 
955 on the {\tt ETITLE} and {\tt ENERGY} lines:
956 \begin{description}
957 \item[{\tt DRUDECOM}] gives the temperature for the warm center-of-mass 
958 degrees of freedom,
959 \item[{\tt DRUDEBOND}] gives the temperature for the cold Drude oscillator 
960 degrees of freedom,
961 \item[{\tt DRCOMAVG}] gives the average temperature 
962 (averaged since the previously reported temperature) 
963 for the warm center-of-mass degrees of freedom, 
964 \item[{\tt DRBONDAVG}] gives the average temperature 
965 (averaged since the previously reported temperature) 
966 for the cold Drude oscillator degrees of freedom.
967 \end{description}
968 The energies resulting from the Drude oscillators and the 
969 anisotropic interactions are summed into the {\tt BOND} energy.  
970 The energies resulting from the LPs and the screened Coulomb 
971 interactions of Thole are summed into the {\tt ELECT} energy.  
972
973
974 \subsubsection{Drude force field parameters}
975
976 The Drude model should be used with the Langevin thermostat enabled
977 ({\tt Langevin=on}).  
978 Doing so permits the use of normal sized time steps (e.g., 1 fs).  
979 The Drude model is also compatible with constant pressure simulation 
980 using the Langevin piston.  
981 Long-range electrostatics may be calculated using PME.  
982 The nonbonded exclusions should generally be set to use either the
983 1-3 exclusion policy ({\tt exclude=1-3})
984 or the scaled 1-4 exclusion policy ({\tt exclude=scaled1-4}).  
985
986 The Drude water model (SWM4-NDP) is a 5-site model
987 with four charge sites and 
988 a negatively charged Drude particle~\cite{Lamoureux-2006a}, 
989 with the particles ordered in the input files as 
990 oxygen, Drude particle, LP, hydrogen, hydrogen.  
991 The atoms in the water molecules should be 
992 constrained ({\tt rigidBonds=water}),
993 with use of the SETTLE algorithm recommended ({\tt useSettle=on}).  
994 Explicitly setting the water model ({\tt waterModel=swm4}) is optional.  
995
996 \begin{itemize}
997
998 \item\NAMDCONFWDEF{drude}
999 {Perform integration of Drude oscillators?}
1000 {{\tt on} or {\tt off}}
1001 {{\tt off}}
1002 {The integration uses a dual Langevin thermostat to freeze the Drude 
1003 oscillators while maintaining the warm degrees of freedom at the
1004 desired temperature.  Must also enable the Langevin thermostat.
1005 If {\tt drude} is set to {\tt on}, then {\tt drudeTemp} must also be defined.}
1006
1007 \item\NAMDCONF{drudeTemp}
1008 {temperature for freezing the Drude oscillators (K)}
1009 {non-negative decimal}
1010 {For stability, the Drude oscillators must be kept at a very cold termpature.
1011 Using a Langevin thermostat, it is possible to set this temperature to 0 K.}
1012
1013 \item\NAMDCONF{drudeDamping}
1014 {damping coefficient for Drude oscillators (1/ps)}
1015 {positive decimal}
1016 {The Langevin coupling coefficient to be applied to the Drude oscillators.  
1017 If not given, {\tt drudeDamping} is set to the value of {\tt langevinDamping}.}
1018
1019 \item\NAMDCONF{drudeBondLen}
1020 {Drude oscillator bond length, beyond which to apply restraint (\AA)}
1021 {positive decimal}
1022 {An additional quartic restraining potential is applied to a Drude 
1023 oscillator if its length exceeds {\tt drudeBondLen}.  
1024 The recommended value is 0.2 \AA, fitted from QM calculations.}
1025
1026 \item\NAMDCONF{drudeBondConst}
1027 {Drude oscillator restraining force constant}
1028 {positive decimal}
1029 {If {\tt drudeBondConst} is defined,
1030 an additional quartic restraining potential is applied to a Drude 
1031 oscillator if its length exceeds {\tt drudeBondLen}.
1032 The recommended value is 40000, fitted from QM calculations.}
1033
1034 \item\NAMDCONF{drudeNbTholeCut}
1035 {nonbonded Thole interaction radius (\AA)}
1036 {positive decimal}
1037 {If {\tt drudeNbTholeCut} is defined, 
1038 the screened Coulomb correction of Thole is also calculated
1039 for non-excluded, nonbonded pairs of Drude oscillators that are
1040 within this radius of interaction.
1041 The recommended value is 5.0 \AA.}
1042
1043 \end{itemize}
1044
1045
1046 \subsection{MARTINI Residue-Based Coarse-Grain Forcefield}
1047 %\label{sec:martini}
1048
1049 The MARTINI forcefield for residue-based coarse-grain models allows simulation of several tens of atoms as only several large coarse-grained particles~\cite{MARR2004, MARR2007, MONT2008}.  In the MARTINI model, each protein residue is represented by a backbone bead and usually one or more sidechain beads.  
1050
1051 When preparing MARTINI simulations it is important to include only those
1052 dihedrals specified by the forcefield.  Using the ``auto dihedrals'' or 
1053 ``regenerate dihedrals'' feature of psfgen will create dihedrals for all
1054 possible sets of four bonded atoms.  This is incorrect for MARTINI and
1055 will result in energy jumps because the dihedral potential function is
1056 degenerate for the angles of 180 degrees allowed by cosine-based angles.
1057
1058 When using MARTINI the following configuration parameters should be set as indicated:
1059
1060 \begin{verbatim}
1061 cosAngles on
1062 martiniSwitching on
1063 dielectric 15.0
1064 PME off
1065 \end{verbatim}
1066
1067 \begin{itemize}
1068
1069 \item
1070 \NAMDCONFWDEF{cosAngles}{enable the MARTINI cosine-based angle potential function}{{\tt on} or {\tt off}}{{\tt off}}
1071 {Specifies whether or not the MARTINI forcefield is being used, specifically cosine-based angle potential function.
1072 The cosine-based potential will only be used for angles in CHARMM parameter files that specify the cos keyword.} 
1073
1074 \item
1075 \NAMDCONFWDEF{martiniSwitching}{enable the MARTINI Lennard-Jones switching function?}{{\tt on} or {\tt off}}{{\tt off}}
1076 {Specifies whether or not the MARTINI forcefield is being used, specifically the Lennard-Jones switching function.} 
1077
1078 \item
1079 \NAMDCONF{martiniDielAllow}{Allow dielectrics != 15.0 for use with MARTINI}{{\tt on} or {\tt off}}{{\tt off}}
1080 {Allows user to specify a {\tt dielectric} not equal to 15.0, ie a non-standard dielectric for MARTINI.}
1081
1082 \end{itemize}
1083
1084
1085 \subsection{Constraints and Restraints}
1086 %\label{section:config_add}
1087
1088 \subsubsection{Bond constraint parameters}
1089 \label{section:rigidBonds}
1090 \begin{itemize}
1091 \item
1092 \NAMDCONFWDEF{rigidBonds}{controls if and how ShakeH is used}{{\tt none},
1093 {\tt water}, {\tt all}}{{\tt none}} 
1094 {When {\tt water} is selected, the hydrogen-oxygen and hydrogen-hydrogen
1095 distances in waters are constrained to the nominal length or angle given
1096 in the parameter file, making the molecules completely rigid.
1097 When {\tt rigidBonds} is {\tt all}, waters are made rigid as described above
1098 and the bond between each hydrogen and the (one) atom to which
1099 it is bonded is similarly constrained.
1100 For the default case {\tt none}, no lengths are constrained.
1101 }
1102
1103 \item
1104 \NAMDCONFWDEF{rigidTolerance}{allowable bond-length error for ShakeH (\AA)}
1105 {positive decimal}{1.0e-8}
1106 {
1107 The ShakeH algorithm is assumed to have converged when all constrained
1108 bonds differ from the nominal bond length by less than this amount.
1109 }
1110
1111 \item
1112 \NAMDCONFWDEF{rigidIterations}{maximum ShakeH iterations}{positive integer}{100}
1113 {
1114 The maximum number of iterations ShakeH will perform before giving up
1115 on constraining the bond lengths.  If the bond lengths do not
1116 converge, a warning message is printed, and the atoms are left at the
1117 final value achieved by ShakeH.  
1118 Although the default value is 100, 
1119 convergence is usually reached after fewer than 10 iterations.
1120 }
1121
1122 \item
1123 \NAMDCONFWDEF{rigidDieOnError}{maximum ShakeH iterations}{{\tt on} or {\tt off}}
1124 {{\tt on}}
1125 {
1126 Exit and report an error if rigidTolerance is not achieved after rigidIterations.
1127 }
1128
1129 \item
1130 \NAMDCONFWDEF{useSettle}{Use SETTLE for waters.}{{\tt on} or {\tt off}}
1131 {{\tt on}}
1132 {
1133 If rigidBonds are enabled then use the non-iterative SETTLE algorithm to
1134 keep waters rigid rather than the slower SHAKE algorithm.
1135 }
1136 \end{itemize}
1137
1138 \subsubsection{Position restraint parameters}
1139
1140 The following describes the parameters for the position restraints feature of \NAMD.
1141 For historical reasons the term ``constraints'' has been carried over from X-PLOR.
1142 This feature allows a restraining potential to each atom of an arbitrary set during the simulation.
1143
1144 \begin{itemize}
1145
1146 \item
1147 \NAMDCONFWDEF{constraints}{are position restraints active?}{{\tt on} or {\tt off}}{{\tt off}}
1148 {Specifies whether or not position restraints are active.
1149 If it is set to {\tt off}, then no position restraints are computed.
1150 If it is set to {\tt on}, the potential $k\times(\mathbf{x} - \mathbf{x}_0)^p$ is applied to each atom.
1151 Per-atom values for $k$ can be defined by either {\tt conskfile} or {\tt conskcol}, for $\mathbf{x}_0$ by {\tt consref}, and for $p$ by {\tt consexp}.
1152 }
1153
1154 \item
1155 \NAMDCONFWDEF{consexp}{exponent for position restraint energy function}{positive, even integer}{2}
1156 {Exponent to be use in the position restraint energy function.
1157 This value must be a positive integer, and only even values really make 
1158 sense.  This parameter is used only if {\tt constraints} is set to 
1159 {\tt on}.}
1160
1161 \item
1162 \NAMDCONF{consref}{PDB file containing restraint reference positions}{UNIX file name}
1163 {PDB file to use for reference positions for position restraints.
1164 Each atom that has a positive force constant will be restrained about the position specified in this file.}
1165
1166 \item
1167 \NAMDCONF{conskfile}{PDB file containing force constant values}{UNIX filename}
1168 {PDB file to use for force constants for 
1169 position restraints.}
1170
1171 \item
1172 \NAMDCONF{conskcol}{column of PDB file containing force constant}{{\tt X}, {\tt Y}, {\tt Z}, {\tt O}, or {\tt B}}
1173 {Column of the PDB file to use for the position restraint force constant.
1174 This parameter may specify any of the floating point fields of the PDB file, 
1175 either X, Y, Z, occupancy, or beta-coupling (temperature-coupling).
1176 Regardless of which column is used, a value of 0 indicates that the atom 
1177 qshould not be restrained.
1178 Otherwise, the value specified is used as the force constant for 
1179 that atom's restraining potential.}
1180
1181 \item
1182 \NAMDCONFWDEF{constraintScaling}{scaling factor for position restraint energy function}{positive}{1.0}
1183 {The position restraint energy function is multiplied by this parameter,
1184 making it possible to gradually turn off restraints during equilibration.
1185 This parameter is used only if {\tt constraints} is set to 
1186 {\tt on}.}
1187
1188 \item
1189 \NAMDCONFWDEF{selectConstraints}{Restrain only selected Cartesian components of the coordinates?}{{\tt on} or {\tt off}}{{\tt off}}
1190 {This option is useful to restrain the positions of atoms to a plane or a line in space. If active,
1191  this option will ensure that only selected Cartesian components of the coordinates are restrained.
1192  E.g.: Restraining the positions of atoms to their current z values with no restraints
1193  in x and y will allow the atoms to move in the x-y plane while retaining their original z-coordinate.
1194  Restraining the x and y values will lead to free motion only along the z coordinate.}
1195
1196 \item
1197 \NAMDCONFWDEF{selectConstrX}{Restrain X components of coordinates}{{\tt on} or {\tt off}}{{\tt off}}
1198 {Restrain the Cartesian x components of the positions.}
1199 \item
1200 \NAMDCONFWDEF{selectConstrY}{Restrain Y components of coordinates}{{\tt on} or {\tt off}}{{\tt off}}
1201 {Restrain the Cartesian y components of the positions.}
1202 \item
1203 \NAMDCONFWDEF{selectConstrZ}{Restrain Z components of coordinates}{{\tt on} or {\tt off}}{{\tt off}}
1204 {Restrain the Cartesian z components of the positions.}
1205
1206 \end{itemize}
1207
1208 \subsubsection{Fixed atoms parameters}
1209
1210 Atoms may be held fixed during a simulation.  \NAMD\ avoids calculating most interactions in which all affected atoms are fixed unless {\tt fixedAtomsForces} is specified.
1211
1212 \begin{itemize}
1213
1214 \item
1215 \NAMDCONFWDEF{fixedAtoms}{are there fixed atoms?}{{\tt on} or {\tt off}}{{\tt off}}
1216 {Specifies whether or not fixed atoms are present.} 
1217
1218 \item
1219 \NAMDCONFWDEF{fixedAtomsForces}{are forces between fixed atoms calculated?}{{\tt on} or {\tt off}}{{\tt off}}
1220 {Specifies whether or not forces between fixed atoms are calculated.  This option is required to turn fixed atoms off in the middle of a simulation.
1221 These forces will affect the pressure calculation, and you should leave this option off when using constant pressure if the coordinates of the fixed atoms have not been minimized.
1222 The use of constant pressure with significant numbers of fixed atoms is not recommended.}
1223
1224 \item
1225 \NAMDCONFWDEF{fixedAtomsFile}{PDB file containing fixed atom parameters}
1226 {UNIX filename}{{\tt coordinates}}
1227 {PDB file to use for the fixed atom flags for each atom.  
1228 If this parameter is not specified, then 
1229 the PDB file specified by {\tt coordinates} is used.}
1230
1231 \item
1232 \NAMDCONFWDEF{fixedAtomsCol}{column of PDB containing fixed atom parameters}
1233 {{\tt X}, {\tt Y}, {\tt Z}, {\tt O}, or {\tt B}}{{\tt O}} 
1234 {Column of the PDB file to use for the containing fixed atom parameters for 
1235 each atom.  The coefficients can be read from any 
1236 floating point column of the PDB file.  
1237 A value of 0 indicates that the atom is not fixed.}
1238
1239 \end{itemize}
1240
1241 \subsubsection{Extra bond, angle, and dihedral restraints}
1242
1243 Additional bond, angle, and dihedral energy terms may be applied to system,
1244 allowing secondary or tertiary structure to be restrained, for example.
1245 Extra bonded terms are not considered part of the molecular structure
1246 and hence do not alter nonbonded exclusions.
1247 The energies from extra bonded terms are included with the normal
1248 bond, angle, and dihedral energies in \NAMD\ output.
1249
1250 All extra bonded terms are harmonic potentials of the form
1251 $U(x) = k (x-x_{ref})^2$ except dihedrals and impropers with
1252 a non-zero periodicity specified, which use
1253 $U(x) = k (1 + cos(n x - x_{ref}))$.
1254 The only difference between dihedrals and
1255 impropers is the output field that their potential energy is added to.
1256
1257 Due to a very old bug all NAMD releases prior to 2.13 have used the
1258 MARTINI cosine-based angle potential function for all extra angles.
1259 Since workflows may unknowingly depend on this undocumented behavior,
1260 cosine-based angles remain the default, but a warning is printed
1261 unless the desired behavior is specified via the new option
1262 extraBondsCosAngles (defaults to ``on'', set to ``off'' to use
1263 the normal harmonic angle potential function for all extra angles).
1264
1265 The extra bonded term implementation shares the parallel implementation
1266 of regular bonded terms in \NAMD, allowing large numbers of extra terms
1267 to be specified with minimal impact on parallel scalability.
1268 Extra bonded terms do not have to duplicate normal bonds/angles/dihedrals,
1269 but each extra bond/angle/dihedral should only involve nearby atoms.
1270 If the atoms involved are too far apart a bad global bond count will be
1271 reported in parallel runs.
1272
1273 Extra bonded terms are enabled via the following options:
1274
1275 \begin{itemize}
1276
1277 \item
1278 \NAMDCONFWDEF{extraBonds}{enable extra bonded terms?}{{\tt on} or {\tt off}}{{\tt off}}
1279 {Specifies whether or not extra bonded terms are present.} 
1280
1281 \item
1282 \NAMDCONFWDEF{extraBondsCosAngles}{are extra angles cosine-based?}{{\tt on} or {\tt off}}{{\tt on}}
1283 {Specifies whether or not all extra angles are cosine-based for consistency with previous versions.
1284 Set to {\tt off} to use normal harmonic angle potential for all extra angles.} 
1285
1286 \item
1287 \NAMDCONF{extraBondsFile}{file containing extra bonded terms}{file}
1288 {File containing extra bonded terms.  May be repeated for multiple files.} 
1289
1290 \end{itemize}
1291
1292 The extra bonds file(s) should contain lines of the following formats:
1293
1294 \begin{itemize}
1295 \item
1296 {\tt bond <atom> <atom> <k> <ref>}
1297 \item
1298 {\tt angle <atom> <atom> <atom> <k> <ref>}
1299 \item
1300 {\tt dihedral <atom> <atom> <atom> <atom> <k> <ref>}
1301 \item
1302 {\tt dihedral <atom> <atom> <atom> <atom> <k> <n> <ref>}
1303 \item
1304 {\tt improper <atom> <atom> <atom> <atom> <k> <ref>}
1305 \item
1306 {\tt improper <atom> <atom> <atom> <atom> <k> <n> <ref>}
1307 \item
1308 {\tt \# <comment ...>}
1309 \end{itemize}
1310
1311 In all cases {\tt <atom>} is a {\bf zero-based} atom index
1312 (the first atom has index 0),
1313 {\tt <ref>} is a reference distance in \AA\ (bond) or angle in degrees (others),
1314 and {\tt <k>} is a spring constant in the potential energy function
1315 $U(x) = k (x-x_{ref})^2$ or, for dihedrals and impropers with 
1316 periodicity {\tt <n>} specified and not 0, $U(x) = k (1 + cos(n x - x_{ref}))$.
1317 Note that $x_{ref}$ is only a minimum for the harmonic potential;
1318 the sinusoidal potential has minima at $(x_{ref} + 180)/n + i \times 360/n$.
1319
1320