1 \section{Force Field Parameters}
2 \label{section:forcefield}
5 \subsection{Potential energy functions}
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.
32 \subsubsection{Bonded potential energy terms}
33 %\label{sec:bonded}
35 The bonded potential terms involve 2--, 3--, and 4--body interactions
36 of covalently bonded atoms.
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.
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.
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.
97 \subsubsection{Nonbonded potential energy terms}
98 %\label{sec:nonbonded}
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.
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.
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
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.
156 \subsection{Non-bonded interactions}
157 \label{section:electdesc}
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.
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.
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}
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}.
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.
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}
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}.
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}
248 \subsubsection{Non-bonded force field parameters}
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.}
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.}
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.}
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.}
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 }
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 %}
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 %}
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}.}
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 }
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.}
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.}
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.}
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.}
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 }
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 }
453 \end{itemize}
456 \subsubsection{PME parameters}
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.
462 \begin{itemize}
464 \item
465 \NAMDCONFWDEF{PME}{Use particle mesh Ewald for electrostatics?}{{\tt yes} or {\tt no}}{{\tt no}}
466 {Turns on particle mesh Ewald.}
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.}
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.}
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).}
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).}
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).}
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).}
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.}
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.}
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.}
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.}
521 \end{itemize}
524 \subsubsection{MSM parameters}
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.
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.
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.
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.
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}.
604 \begin{itemize}
606 \item
607 \NAMDCONFWDEF{MSM}{Use multilevel summation method for electrostatics?}{{\tt yes} or {\tt no}}{{\tt no}}
608 {Turns on multilevel summation method.}
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)$.}
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.
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 }
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 }
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 }
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.}
689 \item
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.}
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}.}
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}.}
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.}
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.}
719 \end{itemize}
722 \subsubsection{Full direct parameters}
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.
736 \begin{itemize}
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.}
743 \end{itemize}
746 \subsubsection{Tabulated nonbonded interaction parameters}
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.
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}
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}
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}
816 The following three parameters are required for tabulated energies.
818 \begin{itemize}
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.}
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.}
830 \item
831 \NAMDCONF{tableInterpType}{cubic or linear interpolation}{{\tt cubic} or {\tt linear}}
832 {Specifies the order for interpolating between energy table entries.}
834 \end{itemize}
837 \subsection{Water Models}
838 \label{section:water_models}
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.
847 \begin{itemize}
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. }
862 \end{itemize}
865 \subsection{Drude polarizable force field}
866 \label{section:drude_forcefield}
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}.
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).
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.
916 \subsubsection{Required input files}
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.
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.
938 \subsubsection{Standard output}
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.
961 \subsubsection{Drude force field parameters}
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}).
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.
982 \begin{itemize}
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 }
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 }
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 }
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 }
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 }
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 }
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 }
1057 \end{itemize}
1060 \subsection{MARTINI Residue-Based Coarse-Grain Forcefield}
1061 %\label{sec:martini}
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.
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.
1072 When using MARTINI the following configuration parameters should be set as indicated:
1074 \begin{verbatim}
1075 cosAngles on
1076 martiniSwitching on
1077 dielectric 15.0
1078 PME off
1079 \end{verbatim}
1081 \begin{itemize}
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.}
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.}
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.}
1096 \end{itemize}
1099 \subsection{Constraints and Restraints}
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 }
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 }
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 }
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 }
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}
1152 \subsubsection{Position restraint parameters}
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.
1158 \begin{itemize}
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 }
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}.}
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.}
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.}
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.}
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}.}
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.}
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.}
1220 \end{itemize}
1222 \subsubsection{Fixed atoms parameters}
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.
1226 \begin{itemize}
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.}
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.}
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.}
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.}
1253 \end{itemize}
1255 \subsubsection{Extra bond, angle, and dihedral restraints}
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.
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.
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).
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.
1287 Extra bonded terms are enabled via the following options:
1289 \begin{itemize}
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.}
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.}
1300 \item
1301 \NAMDCONF{extraBondsFile}{file containing extra bonded terms}{file}
1302 {File containing extra bonded terms.  May be repeated for multiple files.}
1304 \end{itemize}
1306 The extra bonds file(s) should contain lines of the following formats:
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}
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$.