Update default behavior and documentation for Drude
[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 by
869   introducing an auxiliary particle attached to each polarizable atom via a
870   zero-length harmonic spring.
871 The advantage with the Drude model is that it preserves the simple
872   particle-particle Coulomb electrostatic interaction employed in
873   nonpolarizable force fields, thus its implementation in \NAMD\ is more
874   straightforward than alternative models for polarization.
875 \NAMD\ performs the integration of Drude oscillators by employing a novel dual
876   Langevin thermostat to ``freeze'' the Drude oscillators while maintaining the
877   normal ``warm'' degrees of freedom at the desired
878   temperature~\cite{JIAN2011}.
879 Use of the Langevin thermostat enables better parallel scalability than the
880   earlier reported implementation which made use of a dual Nos\'e-Hoover
881   thermostat acting on, and within, each
882   nucleus-Drude pair~\cite{Lamoureux-2003a}.
883 Performance results show that the \NAMD\ implementation of the Drude model
884   maintains good parallel scalability with an increase in computational cost by
885   not more than twice that of using a nonpolarizable force
886   field~\cite{JIAN2011}.
887
888 Excessive ``hyperpolarization'' of Drude oscillators can be prevented by two
889   different schemes.
890 The default ``hard wall'' option reflects elongated springs back towards the
891   nucleus using a simple collision model.
892 Alternatively, the Drude oscillators can be supplemented by a flat-bottom
893   quartic restraining potential (usually with a large force constant).
894
895 The Drude polarizable force field requires some extensions to the CHARMM force
896   field.
897 An anisotropic spring term is added to account for out-of-plane forces from a
898   polarized atom and its covalently bonded neighbor with two more covalently
899   bonded neighbors (similar in structure to an improper bond).
900 The screened Coulomb correction of Thole is calculated between pairs of Drude
901   oscillators that would otherwise be excluded from nonbonded interaction and
902   optionally between non-excluded, nonbonded pairs of Drude oscillators
903   that are within a prescribed cutoff distance~\cite{Thole81,van1998molecular}.
904 Also included in the Drude force field are the use of off-centered massless
905   interaction sites, so called ``lone pairs'' (LPs), to avoid the limitations
906   of centrosymmetric-based Coulomb interactions~\cite{Harder2006}.
907 The coordinate of each LP site is constructed based on three ``host'' atoms.
908 The calculated forces on the massless LP must be transferred to the host atoms,
909   preserving total force and torque.
910 After an integration step of velocities and positions, the position of the LP
911   is updated based on the three host atoms, along with additional geometry
912   parameters that give displacement and in-plane and out-of-plane angles.
913 See our research web page (\url{http://www.ks.uiuc.edu/Research/Drude/}) for
914   additional details and parallel performance results.
915
916 \subsubsection{Required input files}
917
918 No additional files are required by \NAMD\ to use the Drude polarizable force
919   field.
920 However, it is presently beyond the capability of the \PSFGEN{} tool to
921   generate the PSF file needed to perform a simulation using the Drude model.
922 For now, CHARMM is needed to generate correct input files.
923
924 The CHARMM force field parameter files specific to the Drude model are
925   required.
926 The PDB file must also include the Drude particles (mass between 0.05 and 1.0)
927   and the LPs (mass 0).
928 The Drude particles always immediately follow their parent atom.
929 The PSF file augments the ``atom'' section with additional columns that include
930   the ``Thole'' and ``alpha'' parameters for the screened Coulomb interactions
931   of Thole.
932 The PSF file also requires additional sections that list the LPs, including
933   their host atoms and geometry parameters, and list the anisotropic
934   interaction terms, including their parameters.
935 A Drude-compatible PSF file is denoted by the keyword ``DRUDE'' given along the
936   top line.
937
938 \subsubsection{Standard output}
939
940 The \NAMD\ logging to standard output is extended to provide additional
941   temperature data on the cold and warm degrees of freedom.
942 Four additional quantities are listed on the {\tt ETITLE} and {\tt ENERGY}
943   lines:
944 \begin{description}
945 \item[{\tt DRUDECOM}] gives the temperature for the warm center-of-mass
946   degrees of freedom,
947 \item[{\tt DRUDEBOND}] gives the temperature for the cold Drude oscillator
948   degrees of freedom,
949 \item[{\tt DRCOMAVG}] gives the average temperature (averaged since the
950   previously reported temperature) for the warm center-of-mass degrees of
951   freedom,
952 \item[{\tt DRBONDAVG}] gives the average temperature (averaged since the
953   previously reported temperature) for the cold Drude oscillator degrees of
954   freedom.
955 \end{description}
956 The energies resulting from the Drude oscillators and the anisotropic
957   interactions are summed into the {\tt BOND} energy.
958 The energies resulting from the LPs and the screened Coulomb interactions of
959   Thole are summed into the {\tt ELECT} energy.
960
961 \subsubsection{Drude force field parameters}
962
963 The Drude model should be used with the Langevin thermostat enabled
964   ({\tt Langevin=on}).
965 Doing so permits the use of normal sized time steps (e.g., 1 fs).
966 The Drude model is also compatible with constant pressure simulation using the
967   Langevin piston.
968 Long-range electrostatics may be calculated using PME.
969 The nonbonded exclusions should generally be set to use either the 1-3
970   exclusion policy ({\tt exclude=1-3}) or the scaled 1-4 exclusion policy
971   ({\tt exclude=scaled1-4}).
972
973 The Drude water model (SWM4-NDP) is a 5-site model with four charge sites and
974   a negatively charged Drude particle~\cite{Lamoureux-2006a}, with the
975   particles ordered in the input files as oxygen, Drude particle, LP, hydrogen,
976   hydrogen.
977 The atoms in the water molecules should be constrained
978   ({\tt rigidBonds=water}), with use of the SETTLE algorithm recommended
979   ({\tt useSettle=on}).
980 Explicitly setting the water model ({\tt waterModel=swm4}) is optional.
981
982 \begin{itemize}
983
984 \item\NAMDCONFWDEF{drude}
985 {Perform integration of Drude oscillators?}
986 {{\tt on} or {\tt off}}
987 {{\tt off}}
988 {
989 The integration uses a dual Langevin thermostat to freeze the Drude
990   oscillators while maintaining the warm degrees of freedom at the desired
991   temperature.
992 Must also enable the Langevin thermostat.
993 If {\tt drude} is set to {\tt on}, then {\tt drudeTemp} must also be defined.
994 }
995
996 \item\NAMDCONF{drudeTemp}
997 {temperature for freezing the Drude oscillators (K)}
998 {non-negative decimal}
999 {
1000 For stability, the Drude oscillators must be kept at a very cold termpature.
1001 Using a Langevin thermostat, it is possible to set this temperature to 0 K.
1002 }
1003
1004 \item\NAMDCONF{drudeDamping}
1005 {damping coefficient for Drude oscillators (1/ps)}
1006 {positive decimal}
1007 {
1008 The Langevin coupling coefficient to be applied to the Drude oscillators.
1009 If not given, {\tt drudeDamping} is set to the value of {\tt langevinDamping},
1010   but values of as much as an order of magnitude greater may be appropriate.
1011 }
1012
1013 \item\NAMDCONFWDEF{drudeNbTholeCut}
1014 {nonbonded Thole interaction radius (\AA)}
1015 {positive decimal}
1016 {5.0}
1017 {
1018 If {\tt drudeNbTholeCut} is non-zero, the screened Coulomb correction of Thole
1019   is also calculated for non-excluded, nonbonded pairs of Drude oscillators
1020   that are within this radius of interaction.
1021 }
1022
1023 \item\NAMDCONFWDEF{drudeHardWall}
1024 {use collisions to correct hyperpolarization?}
1025 {on or off}
1026 {on}
1027 {
1028 Excessively elongated Drude oscillator bonds are avoided by reflective
1029   collisions induced at a fixed cutoff, {\tt drudeBondLen}.
1030 A large number of such events is usually indicative of unstable/unphysical
1031   dynamics and a simulation will stop if double the cutoff is exceeded.
1032 }
1033
1034 \item\NAMDCONFWDEF{drudeBondLen}
1035 {hyperpolarization cutoff (\AA)}
1036 {positive decimal}
1037 {0.25}
1038 {
1039 If using {\tt drudeHardWall on}, this is the distance at which collisions
1040   occur.
1041 Otherwise, this is the distance at which an additional quartic restraining
1042   potential is applied to each Drude oscillator.
1043 In this latter case, a value of 0.2~\AA{} (slightly smaller than default) is
1044   recommended.
1045 }
1046
1047 \item\NAMDCONFWDEF{drudeBondConst}
1048 {Drude oscillator restraining force constant}
1049 {positive decimal}
1050 {40000.0}
1051 {
1052 If {\tt drudeHardWall off} and {\tt drudeBondConst} is non-zero, an additional
1053   quartic restraining potential is applied to a Drude oscillator if its length
1054   exceeds {\tt drudeBondLen}.
1055 }
1056
1057 \end{itemize}
1058
1059
1060 \subsection{MARTINI Residue-Based Coarse-Grain Forcefield}
1061 %\label{sec:martini}
1062
1063 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.  
1064
1065 When preparing MARTINI simulations it is important to include only those
1066 dihedrals specified by the forcefield.  Using the ``auto dihedrals'' or 
1067 ``regenerate dihedrals'' feature of psfgen will create dihedrals for all
1068 possible sets of four bonded atoms.  This is incorrect for MARTINI and
1069 will result in energy jumps because the dihedral potential function is
1070 degenerate for the angles of 180 degrees allowed by cosine-based angles.
1071
1072 When using MARTINI the following configuration parameters should be set as indicated:
1073
1074 \begin{verbatim}
1075 cosAngles on
1076 martiniSwitching on
1077 dielectric 15.0
1078 PME off
1079 \end{verbatim}
1080
1081 \begin{itemize}
1082
1083 \item
1084 \NAMDCONFWDEF{cosAngles}{enable the MARTINI cosine-based angle potential function}{{\tt on} or {\tt off}}{{\tt off}}
1085 {Specifies whether or not the MARTINI forcefield is being used, specifically cosine-based angle potential function.
1086 The cosine-based potential will only be used for angles in CHARMM parameter files that specify the cos keyword.} 
1087
1088 \item
1089 \NAMDCONFWDEF{martiniSwitching}{enable the MARTINI Lennard-Jones switching function?}{{\tt on} or {\tt off}}{{\tt off}}
1090 {Specifies whether or not the MARTINI forcefield is being used, specifically the Lennard-Jones switching function.} 
1091
1092 \item
1093 \NAMDCONF{martiniDielAllow}{Allow dielectrics != 15.0 for use with MARTINI}{{\tt on} or {\tt off}}{{\tt off}}
1094 {Allows user to specify a {\tt dielectric} not equal to 15.0, ie a non-standard dielectric for MARTINI.}
1095
1096 \end{itemize}
1097
1098
1099 \subsection{Constraints and Restraints}
1100 %\label{section:config_add}
1101
1102 \subsubsection{Bond constraint parameters}
1103 \label{section:rigidBonds}
1104 \begin{itemize}
1105 \item
1106 \NAMDCONFWDEF{rigidBonds}{controls if and how ShakeH is used}{{\tt none},
1107 {\tt water}, {\tt all}}{{\tt none}} 
1108 {When {\tt water} is selected, the hydrogen-oxygen and hydrogen-hydrogen
1109 distances in waters are constrained to the nominal length or angle given
1110 in the parameter file, making the molecules completely rigid.
1111 When {\tt rigidBonds} is {\tt all}, waters are made rigid as described above
1112 and the bond between each hydrogen and the (one) atom to which
1113 it is bonded is similarly constrained.
1114 For the default case {\tt none}, no lengths are constrained.
1115 }
1116
1117 \item
1118 \NAMDCONFWDEF{rigidTolerance}{allowable bond-length error for ShakeH (\AA)}
1119 {positive decimal}{1.0e-8}
1120 {
1121 The ShakeH algorithm is assumed to have converged when all constrained
1122 bonds differ from the nominal bond length by less than this amount.
1123 }
1124
1125 \item
1126 \NAMDCONFWDEF{rigidIterations}{maximum ShakeH iterations}{positive integer}{100}
1127 {
1128 The maximum number of iterations ShakeH will perform before giving up
1129 on constraining the bond lengths.  If the bond lengths do not
1130 converge, a warning message is printed, and the atoms are left at the
1131 final value achieved by ShakeH.  
1132 Although the default value is 100, 
1133 convergence is usually reached after fewer than 10 iterations.
1134 }
1135
1136 \item
1137 \NAMDCONFWDEF{rigidDieOnError}{maximum ShakeH iterations}{{\tt on} or {\tt off}}
1138 {{\tt on}}
1139 {
1140 Exit and report an error if rigidTolerance is not achieved after rigidIterations.
1141 }
1142
1143 \item
1144 \NAMDCONFWDEF{useSettle}{Use SETTLE for waters.}{{\tt on} or {\tt off}}
1145 {{\tt on}}
1146 {
1147 If rigidBonds are enabled then use the non-iterative SETTLE algorithm to
1148 keep waters rigid rather than the slower SHAKE algorithm.
1149 }
1150 \end{itemize}
1151
1152 \subsubsection{Position restraint parameters}
1153
1154 The following describes the parameters for the position restraints feature of \NAMD.
1155 For historical reasons the term ``constraints'' has been carried over from X-PLOR.
1156 This feature allows a restraining potential to each atom of an arbitrary set during the simulation.
1157
1158 \begin{itemize}
1159
1160 \item
1161 \NAMDCONFWDEF{constraints}{are position restraints active?}{{\tt on} or {\tt off}}{{\tt off}}
1162 {Specifies whether or not position restraints are active.
1163 If it is set to {\tt off}, then no position restraints are computed.
1164 If it is set to {\tt on}, the potential $k\times(\mathbf{x} - \mathbf{x}_0)^p$ is applied to each atom.
1165 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}.
1166 }
1167
1168 \item
1169 \NAMDCONFWDEF{consexp}{exponent for position restraint energy function}{positive, even integer}{2}
1170 {Exponent to be use in the position restraint energy function.
1171 This value must be a positive integer, and only even values really make 
1172 sense.  This parameter is used only if {\tt constraints} is set to 
1173 {\tt on}.}
1174
1175 \item
1176 \NAMDCONF{consref}{PDB file containing restraint reference positions}{UNIX file name}
1177 {PDB file to use for reference positions for position restraints.
1178 Each atom that has a positive force constant will be restrained about the position specified in this file.}
1179
1180 \item
1181 \NAMDCONF{conskfile}{PDB file containing force constant values}{UNIX filename}
1182 {PDB file to use for force constants for 
1183 position restraints.}
1184
1185 \item
1186 \NAMDCONF{conskcol}{column of PDB file containing force constant}{{\tt X}, {\tt Y}, {\tt Z}, {\tt O}, or {\tt B}}
1187 {Column of the PDB file to use for the position restraint force constant.
1188 This parameter may specify any of the floating point fields of the PDB file, 
1189 either X, Y, Z, occupancy, or beta-coupling (temperature-coupling).
1190 Regardless of which column is used, a value of 0 indicates that the atom 
1191 qshould not be restrained.
1192 Otherwise, the value specified is used as the force constant for 
1193 that atom's restraining potential.}
1194
1195 \item
1196 \NAMDCONFWDEF{constraintScaling}{scaling factor for position restraint energy function}{positive}{1.0}
1197 {The position restraint energy function is multiplied by this parameter,
1198 making it possible to gradually turn off restraints during equilibration.
1199 This parameter is used only if {\tt constraints} is set to 
1200 {\tt on}.}
1201
1202 \item
1203 \NAMDCONFWDEF{selectConstraints}{Restrain only selected Cartesian components of the coordinates?}{{\tt on} or {\tt off}}{{\tt off}}
1204 {This option is useful to restrain the positions of atoms to a plane or a line in space. If active,
1205  this option will ensure that only selected Cartesian components of the coordinates are restrained.
1206  E.g.: Restraining the positions of atoms to their current z values with no restraints
1207  in x and y will allow the atoms to move in the x-y plane while retaining their original z-coordinate.
1208  Restraining the x and y values will lead to free motion only along the z coordinate.}
1209
1210 \item
1211 \NAMDCONFWDEF{selectConstrX}{Restrain X components of coordinates}{{\tt on} or {\tt off}}{{\tt off}}
1212 {Restrain the Cartesian x components of the positions.}
1213 \item
1214 \NAMDCONFWDEF{selectConstrY}{Restrain Y components of coordinates}{{\tt on} or {\tt off}}{{\tt off}}
1215 {Restrain the Cartesian y components of the positions.}
1216 \item
1217 \NAMDCONFWDEF{selectConstrZ}{Restrain Z components of coordinates}{{\tt on} or {\tt off}}{{\tt off}}
1218 {Restrain the Cartesian z components of the positions.}
1219
1220 \end{itemize}
1221
1222 \subsubsection{Fixed atoms parameters}
1223
1224 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.
1225
1226 \begin{itemize}
1227
1228 \item
1229 \NAMDCONFWDEF{fixedAtoms}{are there fixed atoms?}{{\tt on} or {\tt off}}{{\tt off}}
1230 {Specifies whether or not fixed atoms are present.} 
1231
1232 \item
1233 \NAMDCONFWDEF{fixedAtomsForces}{are forces between fixed atoms calculated?}{{\tt on} or {\tt off}}{{\tt off}}
1234 {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.
1235 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.
1236 The use of constant pressure with significant numbers of fixed atoms is not recommended.}
1237
1238 \item
1239 \NAMDCONFWDEF{fixedAtomsFile}{PDB file containing fixed atom parameters}
1240 {UNIX filename}{{\tt coordinates}}
1241 {PDB file to use for the fixed atom flags for each atom.  
1242 If this parameter is not specified, then 
1243 the PDB file specified by {\tt coordinates} is used.}
1244
1245 \item
1246 \NAMDCONFWDEF{fixedAtomsCol}{column of PDB containing fixed atom parameters}
1247 {{\tt X}, {\tt Y}, {\tt Z}, {\tt O}, or {\tt B}}{{\tt O}} 
1248 {Column of the PDB file to use for the containing fixed atom parameters for 
1249 each atom.  The coefficients can be read from any 
1250 floating point column of the PDB file.  
1251 A value of 0 indicates that the atom is not fixed.}
1252
1253 \end{itemize}
1254
1255 \subsubsection{Extra bond, angle, and dihedral restraints}
1256
1257 Additional bond, angle, and dihedral energy terms may be applied to system,
1258 allowing secondary or tertiary structure to be restrained, for example.
1259 Extra bonded terms are not considered part of the molecular structure
1260 and hence do not alter nonbonded exclusions.
1261 The energies from extra bonded terms are included with the normal
1262 bond, angle, and dihedral energies in \NAMD\ output.
1263
1264 All extra bonded terms are harmonic potentials of the form
1265 $U(x) = k (x-x_{ref})^2$ except dihedrals and impropers with
1266 a non-zero periodicity specified, which use
1267 $U(x) = k (1 + cos(n x - x_{ref}))$.
1268 The only difference between dihedrals and
1269 impropers is the output field that their potential energy is added to.
1270
1271 Due to a very old bug all NAMD releases prior to 2.13 have used the
1272 MARTINI cosine-based angle potential function for all extra angles.
1273 Since workflows may unknowingly depend on this undocumented behavior,
1274 cosine-based angles remain the default, but a warning is printed
1275 unless the desired behavior is specified via the new option
1276 extraBondsCosAngles (defaults to ``on'', set to ``off'' to use
1277 the normal harmonic angle potential function for all extra angles).
1278
1279 The extra bonded term implementation shares the parallel implementation
1280 of regular bonded terms in \NAMD, allowing large numbers of extra terms
1281 to be specified with minimal impact on parallel scalability.
1282 Extra bonded terms do not have to duplicate normal bonds/angles/dihedrals,
1283 but each extra bond/angle/dihedral should only involve nearby atoms.
1284 If the atoms involved are too far apart a bad global bond count will be
1285 reported in parallel runs.
1286
1287 Extra bonded terms are enabled via the following options:
1288
1289 \begin{itemize}
1290
1291 \item
1292 \NAMDCONFWDEF{extraBonds}{enable extra bonded terms?}{{\tt on} or {\tt off}}{{\tt off}}
1293 {Specifies whether or not extra bonded terms are present.} 
1294
1295 \item
1296 \NAMDCONFWDEF{extraBondsCosAngles}{are extra angles cosine-based?}{{\tt on} or {\tt off}}{{\tt on}}
1297 {Specifies whether or not all extra angles are cosine-based for consistency with previous versions.
1298 Set to {\tt off} to use normal harmonic angle potential for all extra angles.} 
1299
1300 \item
1301 \NAMDCONF{extraBondsFile}{file containing extra bonded terms}{file}
1302 {File containing extra bonded terms.  May be repeated for multiple files.} 
1303
1304 \end{itemize}
1305
1306 The extra bonds file(s) should contain lines of the following formats:
1307
1308 \begin{itemize}
1309 \item
1310 {\tt bond <atom> <atom> <k> <ref>}
1311 \item
1312 {\tt angle <atom> <atom> <atom> <k> <ref>}
1313 \item
1314 {\tt dihedral <atom> <atom> <atom> <atom> <k> <ref>}
1315 \item
1316 {\tt dihedral <atom> <atom> <atom> <atom> <k> <n> <ref>}
1317 \item
1318 {\tt improper <atom> <atom> <atom> <atom> <k> <ref>}
1319 \item
1320 {\tt improper <atom> <atom> <atom> <atom> <k> <n> <ref>}
1321 \item
1322 {\tt \# <comment ...>}
1323 \end{itemize}
1324
1325 In all cases {\tt <atom>} is a {\bf zero-based} atom index
1326 (the first atom has index 0),
1327 {\tt <ref>} is a reference distance in \AA\ (bond) or angle in degrees (others),
1328 and {\tt <k>} is a spring constant in the potential energy function
1329 $U(x) = k (x-x_{ref})^2$ or, for dihedrals and impropers with 
1330 periodicity {\tt <n>} specified and not 0, $U(x) = k (1 + cos(n x - x_{ref}))$.
1331 Note that $x_{ref}$ is only a minimum for the harmonic potential;
1332 the sinusoidal potential has minima at $(x_{ref} + 180)/n + i \times 360/n$.
1333
1334