ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/mattDisertation/oopse.tex
(Generate patch)

Comparing trunk/mattDisertation/oopse.tex (file contents):
Revision 1054 by mmeineke, Mon Feb 16 20:35:58 2004 UTC vs.
Revision 1071 by mmeineke, Thu Feb 26 21:55:43 2004 UTC

# Line 22 | Line 22 | F.~Vardeman II, Teng Lin, Christopher J.~Fennell and J
22   simulation package of this size and scope would not have been possible
23   without the collaborative efforts of my colleagues: Charles
24   F.~Vardeman II, Teng Lin, Christopher J.~Fennell and J.~Daniel
25 < Gezelter. Although my contributions to {\sc oopse} are significant,
26 < consideration of my work apart from the others, would not give a
25 > Gezelter. Although my contributions to {\sc oopse} are major,
26 > consideration of my work apart from the others would not give a
27   complete description to the package's capabilities. As such, all
28   contributions to {\sc oopse} to date are presented in this chapter.
29  
30 < Charles Vardeman is responsible for the parallelization of {\sc oopse}
31 < (Sec.~\ref{oopseSec:parallelization}) as well as the inclusion of the
32 < embedded-atom potential (Sec.~\ref{oopseSec:eam}). Teng Lin's
33 < contributions include refinement of the periodic boundary conditions
30 > Charles Vardeman is responsible for the parallelization of the long
31 > range forces in {\sc oopse} (Sec.~\ref{oopseSec:parallelization}) as
32 > well as the inclusion of the embedded-atom potential for transition
33 > metals (Sec.~\ref{oopseSec:eam}). Teng Lin's contributions include
34 > refinement of the periodic boundary conditions
35   (Sec.~\ref{oopseSec:pbc}), the z-constraint method
36   (Sec.~\ref{oopseSec:zcons}), refinement of the property analysis
37   programs (Sec.~\ref{oopseSec:props}), and development in the extended
38 < state integrators (Sec.~\ref{oopseSec:noseHooverThermo}). Christopher
38 > system integrators (Sec.~\ref{oopseSec:noseHooverThermo}). Christopher
39   Fennell worked on the symplectic integrator
40   (Sec.~\ref{oopseSec:integrate}) and the refinement of the {\sc ssd}
41   water model (Sec.~\ref{oopseSec:SSD}). Daniel Gezelter lent his
# Line 48 | Line 49 | the property analysis (Sec.~\ref{oopseSec:props}) and
49   of the Lennard-Jones (Sec.~\ref{sec:LJPot}) and {\sc duff}
50   (Sec.~\ref{oopseSec:DUFF}) force fields, and initial implementation of
51   the property analysis (Sec.~\ref{oopseSec:props}) and system
52 < initialization (Sec.~\ref{oopseSec:initCoords}) utility programs.
52 > initialization (Sec.~\ref{oopseSec:initCoords}) utility programs. {\sc
53 > oopse}, like many other Molecular Dynamics programs, is a work in
54 > progress, and will continue to be so for many graduate student
55 > lifetimes.
56  
57   \section{\label{sec:intro}Introduction}
58  
# Line 66 | Line 70 | researchers try to develop techniques or energetic mod
70  
71   Despite their utility, problems with these packages arise when
72   researchers try to develop techniques or energetic models that the
73 < code was not originally designed to do. Examples of uncommonly
73 > code was not originally designed to simulate. Examples of uncommonly
74   implemented techniques and energetics include; dipole-dipole
75   interactions, rigid body dynamics, and metallic embedded
76   potentials. When faced with these obstacles, a researcher must either
# Line 75 | Line 79 | our research is based.
79   simulation code capable of implementing the types of models upon which
80   our research is based.
81  
82 < Having written {\sc oopse} we are implementing the concept of Open
83 < Source development, and releasing our source code into the public
84 < domain. It is our intent that by doing so, other researchers might
85 < benefit from our work, and add their own contributions to the
86 < package. The license under which {\sc oopse} is distributed allows any
87 < researcher to download and modify the source code for their own
88 < use. In this way further development of {\sc oopse} is not limited to
89 < only the models of interest to ourselves, but also those of the
90 < community of scientists who contribute back to the project.
82 > In developing {\sc oopse}, we have adhered to the precepts of Open
83 > Source development, and are releasing our source code with a
84 > permissive license. It is our intent that by doing so, other
85 > researchers might benefit from our work, and add their own
86 > contributions to the package. The license under which {\sc oopse} is
87 > distributed allows any researcher to download and modify the source
88 > code for their own use. In this way further development of {\sc oopse}
89 > is not limited to only the models of interest to ourselves, but also
90 > those of the community of scientists who contribute back to the
91 > project.
92  
93   We have structured this chapter to first discuss the empirical energy
94   functions that {\sc oopse } implements in
95   Sec.~\ref{oopseSec:empiricalEnergy}. Following that is a discussion of
96   the various input and output files associated with the package
97 < (Sec.~\ref{oopseSec:IOfiles}). In Sec.~\ref{oopseSec:mechanics}
97 > (Sec.~\ref{oopseSec:IOfiles}). Sec.~\ref{oopseSec:mechanics}
98   elucidates the various Molecular Dynamics algorithms {\sc oopse}
99   implements in the integration of the Newtonian equations of
100   motion. Basic analysis of the trajectories obtained from the
101   simulation is discussed in Sec.~\ref{oopseSec:props}. Program design
102 < considerations as well as the software distribution license is
103 < presented in Sec.~\ref{oopseSec:design}. And lastly,
99 < Sec.~\ref{oopseSec:conclusion} concludes the chapter.
102 > considerations are presented in Sec.~\ref{oopseSec:design}. And
103 > lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter.
104  
105   \section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions}
106  
# Line 109 | Line 113 | dipoles). Charges, permanent dipoles, and Lennard-Jone
113   methyl and carbonyl groups. The atoms are also capable of having
114   directional components associated with them (\emph{e.g.}~permanent
115   dipoles). Charges, permanent dipoles, and Lennard-Jones parameters for
116 < a given atom type are set in the force field parameter files
116 > a given atom type are set in the force field parameter files.
117  
118   \begin{lstlisting}[float,caption={[Specifier for molecules and atoms] A sample specification of an Ar molecule},label=sch:AtmMole]
119   molecule{
# Line 129 | Line 133 | internal interactions (\emph{i.e.}~bonds, bends, and t
133   identities of all the atoms and rigid bodies associated with
134   themselves, and are responsible for the evaluation of their own
135   internal interactions (\emph{i.e.}~bonds, bends, and torsions). Scheme
136 < \ref{sch:AtmMole} shows how one creates a molecule in a
136 > \ref{sch:AtmMole} shows how one creates a molecule in a ``model'' or
137   \texttt{.mdl} file. The position of the atoms given in the
138   declaration are relative to the origin of the molecule, and is used
139   when creating a system containing the molecule.
# Line 139 | Line 143 | potential and move collectively.\cite{Goldstein01} The
143   to handle rigid body dynamics. Rigid bodies are non-spherical
144   particles or collections of particles that have a constant internal
145   potential and move collectively.\cite{Goldstein01} They are not
146 < included in most simulation packages because of the requirement to
147 < propagate the orientational degrees of freedom. Until recently,
148 < integrators which propagate orientational motion have been lacking.
146 > included in most simulation packages because of the algorithmic
147 > complexity involved in propagating orientational degrees of
148 > freedom. Until recently, integrators which propagate orientational
149 > motion have been much worse than those available for translational
150 > motion.
151  
152   Moving a rigid body involves determination of both the force and
153   torque applied by the surroundings, which directly affect the
# Line 150 | Line 156 | Accumulation of the total torque on the rigid body is
156   first be calculated for all the internal particles. The total force on
157   the rigid body is simply the sum of these external forces.
158   Accumulation of the total torque on the rigid body is more complex
159 < than the force in that it is the torque applied on the center of mass
160 < that dictates rotational motion. The torque on rigid body {\it i} is
159 > than the force because the torque is applied to the center of mass of
160 > the rigid body. The torque on rigid body $i$ is
161   \begin{equation}
162   \boldsymbol{\tau}_i=
163 <        \sum_{a}(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia}
164 <        + \boldsymbol{\tau}_{ia},
163 >        \sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia}
164 >        + \boldsymbol{\tau}_{ia}\biggr]
165   \label{eq:torqueAccumulate}
166   \end{equation}
167   where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and
# Line 164 | Line 170 | The summation of the total torque is done in the body
170   position of, and torque on the component particles of the rigid body.
171  
172   The summation of the total torque is done in the body fixed axis of
173 < the rigid body. In order to move between the space fixed and body
173 > each rigid body. In order to move between the space fixed and body
174   fixed coordinate axes, parameters describing the orientation must be
175   maintained for each rigid body. At a minimum, the rotation matrix
176   (\textbf{A}) can be described by the three Euler angles ($\phi,
# Line 179 | Line 185 | systems.\cite{Evans77}
185   systems.\cite{Evans77}
186  
187   {\sc oopse} utilizes a relatively new scheme that propagates the
188 < entire nine parameter rotation matrix internally. Further discussion
188 > entire nine parameter rotation matrix. Further discussion
189   on this choice can be found in Sec.~\ref{oopseSec:integrate}. An example
190   definition of a rigid body can be seen in Scheme
191   \ref{sch:rigidBody}. The positions in the atom definitions are the
# Line 233 | Line 239 | Lennard-Jones force field.
239  
240   \begin{lstlisting}[float,caption={[Invocation of the Lennard-Jones force field] A sample system using the Lennard-Jones force field.},label={sch:LJFF}]
241  
236 – /*
237 – * The Ar molecule is specified
238 – * external to the.bass file
239 – */
240 –
242   #include "argon.mdl"
243  
244   nComponents = 1;
# Line 246 | Line 247 | component{
247    nMol = 108;
248   }
249  
249 – /*
250 – * The initial configuration is generated
251 – * before the simulation is invoked.
252 – */
253 –
250   initialConfig = "./argon.init";
251  
252   forceField = "LJ";
# Line 259 | Line 255 | keep the pair evaluations to a manageable number, {\sc
255   Because this potential is calculated between all pairs, the force
256   evaluation can become computationally expensive for large systems. To
257   keep the pair evaluations to a manageable number, {\sc oopse} employs
258 < a cut-off radius.\cite{allen87:csl} The cutoff radius is set to be
258 > a cut-off radius.\cite{allen87:csl} The cutoff radius can either be
259 > specified in the \texttt{.bass} file, or left as its default value of
260   $2.5\sigma_{ii}$, where $\sigma_{ii}$ is the largest Lennard-Jones
261   length parameter present in the simulation. Truncating the calculation
262   at $r_{\text{cut}}$ introduces a discontinuity into the potential
263 < energy. To offset this discontinuity, the energy value at
264 < $r_{\text{cut}}$ is subtracted from the potential. This causes the
265 < potential to go to zero smoothly at the cut-off radius.
263 > energy and the force. To offset this discontinuity in the potential,
264 > the energy value at $r_{\text{cut}}$ is subtracted from the
265 > potential. This causes the potential to go to zero smoothly at the
266 > cut-off radius, and preserves conservation of energy in integrating
267 > the equations of motion.
268  
269   Interactions between dissimilar particles requires the generation of
270   cross term parameters for $\sigma$ and $\epsilon$. These are
# Line 280 | Line 279 | and
279   \epsilon_{ij} = \sqrt{\epsilon_{ii} \epsilon_{jj}}
280   \label{eq:epsilonMix}
281   \end{equation}
283 –
284 –
282  
283   \subsection{\label{oopseSec:DUFF}Dipolar Unified-Atom Force Field}
284  
# Line 296 | Line 293 | interaction sites. This simplification cuts the length
293   charges. Charge-neutral distributions were replaced with dipoles,
294   while most atoms and groups of atoms were reduced to Lennard-Jones
295   interaction sites. This simplification cuts the length scale of long
296 < range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, allowing us
297 < to avoid the computationally expensive Ewald sum. Instead, we can use
298 < neighbor-lists, reaction field, and cutoff radii for the dipolar
299 < interactions.
296 > range interactions from $\frac{1}{r}$ to $\frac{1}{r^3}$, and allows
297 > us to avoid the computationally expensive Ewald sum. Instead, we can
298 > use neighbor-lists and cutoff radii for the dipolar interactions, or
299 > include a reaction field to mimic larger range interactions.
300  
301   As an example, lipid head-groups in {\sc duff} are represented as
302 < point dipole interaction sites. By placing a dipole of 20.6~Debye at
303 < the head group center of mass, our model mimics the head group of
304 < phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones
305 < site is located at the pseudoatom's center of mass. The model is
306 < illustrated by the dark grey atom in Fig.~\ref{oopseFig:lipidModel}. The
307 < water model we use to complement the dipoles of the lipids is our
308 < reparameterization of the soft sticky dipole (SSD) model of Ichiye
302 > point dipole interaction sites. By placing a dipole at the head group
303 > center of mass, our model mimics the charge separation found in common
304 > phospholipids such as phosphatidylcholine.\cite{Cevc87} Additionally,
305 > a large Lennard-Jones site is located at the pseudoatom's center of
306 > mass. The model is illustrated by the red atom in
307 > Fig.~\ref{oopseFig:lipidModel}. The water model we use to complement
308 > the dipoles of the lipids is our reparameterization of the soft sticky
309 > dipole (SSD) model of Ichiye
310   \emph{et al.}\cite{liu96:new_model}
311  
312   \begin{figure}
# Line 328 | Line 326 | it generalizes the types of atoms in an alkyl chain to
326   equilibria using Gibbs ensemble Monte Carlo simulation
327   techniques.\cite{Siepmann1998} One of the advantages of TraPPE is that
328   it generalizes the types of atoms in an alkyl chain to keep the number
329 < of pseudoatoms to a minimum; the parameters for an atom such as
329 > of pseudoatoms to a minimum; the parameters for a unified atom such as
330   $\text{CH}_2$ do not change depending on what species are bonded to
331   it.
332  
333   TraPPE also constrains all bonds to be of fixed length. Typically,
334   bond vibrations are the fastest motions in a molecular dynamic
335   simulation. Small time steps between force evaluations must be used to
336 < ensure adequate sampling of the bond potential to ensure conservation
337 < of energy. By constraining the bond lengths, larger time steps may be
338 < used when integrating the equations of motion. A simulation using {\sc
339 < duff} is illustrated in Scheme \ref{sch:DUFF}.
336 > ensure adequate energy conservation in the bond degrees of freedom. By
337 > constraining the bond lengths, larger time steps may be used when
338 > integrating the equations of motion. A simulation using {\sc duff} is
339 > illustrated in Scheme \ref{sch:DUFF}.
340  
341   \begin{lstlisting}[float,caption={[Invocation of {\sc duff}]Sample \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}]
342  
# Line 407 | Line 405 | V_{\text{torsion}}(\phi) = c_1[1 + \cos \phi]
405          + c_3[1 + \cos(3\phi)]
406   \label{eq:origTorsionPot}
407   \end{equation}
408 < Here $\phi$ is the angle defined by four bonded neighbors $i$,
411 < $j$, $k$, and $l$ (again, see Fig.~\ref{oopseFig:lipidModel}). For
412 < computational efficiency, the torsion potential has been recast after
413 < the method of {\sc charmm},\cite{Brooks83} in which the angle series is
414 < converted to a power series of the form:
408 > Where:
409   \begin{equation}
410 + \cos\phi = (\hat{\mathbf{r}}_{ij} \times \hat{\mathbf{r}}_{jk}) \cdot
411 +        (\hat{\mathbf{r}}_{jk} \times \hat{\mathbf{r}}_{kl})
412 + \label{eq:torsPhi}
413 + \end{equation}
414 + Here, $\hat{\mathbf{r}}_{\alpha\beta}$ are the set of unit bond
415 + vectors between atoms $i$, $j$, $k$, and $l$. For computational
416 + efficiency, the torsion potential has been recast after the method of
417 + {\sc charmm},\cite{Brooks83} in which the angle series is converted to
418 + a power series of the form:
419 + \begin{equation}
420   V_{\text{torsion}}(\phi) =  
421          k_3 \cos^3 \phi + k_2 \cos^2 \phi + k_1 \cos \phi + k_0
422   \label{eq:torsionPot}
# Line 452 | Line 456 | V_{\text{dipole}}(\mathbf{r}_{ij},\boldsymbol{\Omega}_
456          \boldsymbol{\Omega}_{j}) = \frac{|\mu_i||\mu_j|}{4\pi\epsilon_{0}r_{ij}^{3}} \biggl[
457          \boldsymbol{\hat{u}}_{i} \cdot \boldsymbol{\hat{u}}_{j}
458          -
459 <        \frac{3(\boldsymbol{\hat{u}}_i \cdot \mathbf{r}_{ij}) %
460 <                (\boldsymbol{\hat{u}}_j \cdot \mathbf{r}_{ij}) }
457 <                {r^{2}_{ij}} \biggr]
459 >        3(\boldsymbol{\hat{u}}_i \cdot \hat{\mathbf{r}}_{ij}) %
460 >                (\boldsymbol{\hat{u}}_j \cdot \hat{\mathbf{r}}_{ij}) \biggr]
461   \label{eq:dipolePot}
462   \end{equation}
463   Here $\mathbf{r}_{ij}$ is the vector starting at atom $i$ pointing
# Line 466 | Line 469 | unit vector pointing along $\mathbf{r}_{ij}$
469   unit vector pointing along $\mathbf{r}_{ij}$
470   ($\boldsymbol{\hat{r}}_{ij}=\mathbf{r}_{ij}/|\mathbf{r}_{ij}|$).
471  
472 + To improve computational efficiency of the dipole-dipole interactions,
473 + {\sc oopse} employs an electrostatic cutoff radius. This parameter can
474 + be set in the \texttt{.bass} file, and controls the length scale over
475 + which dipole interactions are felt. To compensate for the
476 + discontinuity in the potential and the forces at the cutoff radius, we
477 + have implemented a switching function to smoothly scale the
478 + dipole-dipole interaction at the cutoff.
479 + \begin{equation}
480 + S(r_{ij}) =
481 +        \begin{cases}
482 +        1 & \text{if $r_{ij} \le r_t$},\\
483 +        \frac{(r_{\text{cut}} + 2r_{ij} - 3r_t)(r_{\text{cut}} - r_{ij})^2}
484 +        {(r_{\text{cut}} - r_t)^2}
485 +        & \text{if $r_t < r_{ij} \le r_{\text{cut}}$}, \\
486 +        0 & \text{if $r_{ij} > r_{\text{cut}}$.}
487 +        \end{cases}
488 + \label{eq:dipoleSwitching}
489 + \end{equation}
490 + Here $S(r_{ij})$ scales the potential at a given $r_{ij}$, and $r_t$
491 + is the taper radius some given thickness less than the electrostatic
492 + cutoff. The switching thickness can be set in the \texttt{.bass} file.
493  
494 < \subsubsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF}
494 > \subsection{\label{oopseSec:SSD}The {\sc duff} Water Models: SSD/E and SSD/RF}
495  
496   In the interest of computational efficiency, the default solvent used
497   by {\sc oopse} is the extended Soft Sticky Dipole (SSD/E) water
# Line 527 | Line 551 | articles.\cite{liu96:new_model,liu96:monte_carlo,chand
551   can be found in the original SSD
552   articles.\cite{liu96:new_model,liu96:monte_carlo,chandra99:ssd_md,Ichiye03}
553  
554 < Since SSD is a single-point {\it dipolar} model, the force
554 > Since SSD/E is a single-point {\it dipolar} model, the force
555   calculations are simplified significantly relative to the standard
556   {\it charged} multi-point models. In the original Monte Carlo
557   simulations using this model, Ichiye {\it et al.} reported that using
558   SSD decreased computer time by a factor of 6-7 compared to other
559   models.\cite{liu96:new_model} What is most impressive is that these savings
560   did not come at the expense of accurate depiction of the liquid state
561 < properties.  Indeed, SSD maintains reasonable agreement with the Soper
561 > properties.  Indeed, SSD/E maintains reasonable agreement with the Head-Gordon
562   diffraction data for the structural features of liquid
563 < water.\cite{Soper86,liu96:new_model} Additionally, the dynamical properties
564 < exhibited by SSD agree with experiment better than those of more
563 > water.\cite{hura00,liu96:new_model} Additionally, the dynamical properties
564 > exhibited by SSD/E agree with experiment better than those of more
565   computationally expensive models (like TIP3P and
566   SPC/E).\cite{chandra99:ssd_md} The combination of speed and accurate depiction
567 < of solvent properties makes SSD a very attractive model for the
567 > of solvent properties makes SSD/E a very attractive model for the
568   simulation of large scale biochemical simulations.
569  
570   Recent constant pressure simulations revealed issues in the original
# Line 551 | Line 575 | model. Solvent parameters can be easily modified in an
575   of a reaction field long-range interaction correction is desired, it
576   is recommended that the parameters be modified to those of the SSD/RF
577   model. Solvent parameters can be easily modified in an accompanying
578 < {\sc BASS} file as illustrated in the scheme below. A table of the
578 > \texttt{.bass} file as illustrated in the scheme below. A table of the
579   parameter values and the drawbacks and benefits of the different
580   density corrected SSD models can be found in
581   reference~\cite{Gezelter04}.
# Line 571 | Line 595 | forceField = "DUFF";
595   forceField = "DUFF";
596  
597   /*
574 – * The reactionField flag toggles reaction
575 – * field corrections.
576 – */
577 –
578 – reactionField = false; // defaults to false
579 – dielectric = 80.0; // dielectric for reaction field
580 –
581 – /*
598   * The following two flags set the cutoff
599   * radius for the electrostatic forces
600   * as well as the skin thickness of the switching
# Line 597 | Line 613 | describe bonding transition metal
613   capacity to simulate metallic systems, including some that have
614   parallel computational abilities\cite{plimpton93}. Potentials that
615   describe bonding transition metal
616 < systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have a
616 > systems\cite{Finnis84,Ercolessi88,Chen90,Qi99,Ercolessi02} have an
617   attractive interaction which models  ``Embedding''
618   a positively charged metal ion in the electron density due to the
619   free valance ``sea'' of electrons created by the surrounding atoms in
620 < the system. A mostly repulsive pairwise part of the potential
620 > the system. A mostly-repulsive pairwise part of the potential
621   describes the interaction of the positively charged metal core ions
622   with one another. A particular potential description called the
623   Embedded Atom Method\cite{Daw84,FBD86,johnson89,Lu97}({\sc eam}) that has
624   particularly wide adoption has been selected for inclusion in {\sc oopse}. A
625 < good review of {\sc eam} and other metallic potential formulations was done
625 > good review of {\sc eam} and other metallic potential formulations was written
626   by Voter.\cite{voter}
627  
628   The {\sc eam} potential has the form:
# Line 614 | Line 630 | V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i}
630   V & = & \sum_{i} F_{i}\left[\rho_{i}\right] + \sum_{i} \sum_{j \neq i}
631   \phi_{ij}({\bf r}_{ij})  \\
632   \rho_{i}  & = & \sum_{j \neq i} f_{j}({\bf r}_{ij})
633 < \end{eqnarray}S
618 <
633 > \end{eqnarray}
634   where $F_{i} $ is the embedding function that equates the energy required to embed a
635   positively-charged core ion $i$ into a linear superposition of
636   spherically averaged atomic electron densities given by
# Line 634 | Line 649 | interactions. Foiles et al. fit EAM potentials for fcc
649  
650   \newcommand{\roundme}{\operatorname{round}}
651  
652 < \textit{Periodic boundary conditions} are widely used to simulate truly
653 < macroscopic systems with a relatively small number of particles. The
654 < simulation box is replicated throughout space to form an infinite lattice.
655 < During the simulation, when a particle moves in the primary cell, its image in
656 < other boxes move in exactly the same direction with exactly the same
657 < orientation.Thus, as a particle leaves the primary cell, one of its images
658 < will enter through the opposite face.If the simulation box is large enough to
659 < avoid \textquotedblleft feeling\textquotedblright\ the symmetries of the
660 < periodic lattice, surface effects can be ignored. Cubic, orthorhombic and
661 < parallelepiped are the available periodic cells In OOPSE. We use a matrix to
662 < describe the property of the simulation box. Therefore, both the size and
648 < shape of the simulation box can be changed during the simulation. The
649 < transformation from box space vector $\mathbf{s}$ to its corresponding real
650 < space vector $\mathbf{r}$ is defined by
652 > \textit{Periodic boundary conditions} are widely used to simulate bulk properties with a relatively small number of particles. The
653 > simulation box is replicated throughout space to form an infinite
654 > lattice.  During the simulation, when a particle moves in the primary
655 > cell, its image in other cells move in exactly the same direction with
656 > exactly the same orientation. Thus, as a particle leaves the primary
657 > cell, one of its images will enter through the opposite face. If the
658 > simulation box is large enough to avoid ``feeling'' the symmetries of
659 > the periodic lattice, surface effects can be ignored. The available
660 > periodic cells in OOPSE are cubic, orthorhombic and parallelepiped. We
661 > use a $3 \times 3$ matrix, $\mathbf{H}$, to describe the shape and
662 > size of the simulation box. $\mathbf{H}$ is defined:
663   \begin{equation}
664 < \mathbf{r}=\underline{\mathbf{H}}\cdot\mathbf{s}%
664 > \mathbf{H} = ( \mathbf{h}_x, \mathbf{h}_y, \mathbf{h}_z )
665   \end{equation}
666 + Where $\mathbf{h}_j$ is the column vector of the $j$th axis of the
667 + box.  During the course of the simulation both the size and shape of
668 + the box can be changed to allow volume fluctations when constraining
669 + the pressure.
670  
671 <
672 < where $H=(h_{x},h_{y},h_{z})$ is a transformation matrix made up of the three
673 < box axis vectors. $h_{x},h_{y}$ and $h_{z}$ represent the three sides of the
674 < simulation box respectively.
675 <
676 < To find the minimum image of a vector $\mathbf{r}$, we convert the real vector
677 < to its corresponding vector in box space first, \bigskip%
671 > A real space vector, $\mathbf{r}$ can be transformed in to a box space
672 > vector, $\mathbf{s}$, and back through the following transformations:
673 > \begin{align}
674 > \mathbf{s} &= \mathbf{H}^{-1} \mathbf{r} \\
675 > \mathbf{r} &= \mathbf{H} \mathbf{s}
676 > \end{align}
677 > The vector $\mathbf{s}$ is now a vector expressed as the number of box
678 > lengths in the $\mathbf{h}_x$, $\mathbf{h}_y$, and $\mathbf{h}_z$
679 > directions. To find the minimum image of a vector $\mathbf{r}$, we
680 > first convert it to its corresponding vector in box space, and then,
681 > cast each element to lie on the in the range $[-0.5,0.5]$:
682   \begin{equation}
663 – \mathbf{s}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{r}%
664 – \end{equation}
665 – And then, each element of $\mathbf{s}$ is wrapped to lie between -0.5 to 0.5,
666 – \begin{equation}
683   s_{i}^{\prime}=s_{i}-\roundme(s_{i})
684   \end{equation}
685 < where
686 <
671 < %
672 <
685 > Where $s_i$ is the $i$th element of $\mathbf{s}$, and
686 > $\roundme(s_i)$is given by
687   \begin{equation}
688 < \roundme(x)=\left\{
689 < \begin{array}{cc}%
690 < \lfloor{x+0.5}\rfloor & \text{if \ }x\geqslant 0 \\
691 < \lceil{x-0.5}\rceil & \text{otherwise}%
692 < \end{array}
679 < \right.
688 > \roundme(x) =
689 >        \begin{cases}
690 >        \lfloor x+0.5 \rfloor & \text{if $x \ge 0$} \\
691 >        \lceil x-0.5 \rceil & \text{if $x < 0$ }
692 >        \end{cases}
693   \end{equation}
694 + Here $\lfloor x \rfloor$ is the floor operator, and gives the largest
695 + integer value that is not greater than $x$, and $\lceil x \rceil$ is
696 + the ceiling operator, and gives the smallest integer that is not less
697 + than $x$.  For example, $\roundme(3.6)=4$, $\roundme(3.1)=3$,
698 + $\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$.
699  
682 –
683 – For example, $\roundme(3.6)=4$,$\roundme(3.1)=3$, $\roundme(-3.6)=-4$, $\roundme(-3.1)=-3$.
684 –
700   Finally, we obtain the minimum image coordinates $\mathbf{r}^{\prime}$ by
701 < transforming back to real space,%
687 <
701 > transforming back to real space,
702   \begin{equation}
703 < \mathbf{r}^{\prime}=\underline{\mathbf{H}}^{-1}\cdot\mathbf{s}^{\prime}%
703 > \mathbf{r}^{\prime}=\mathbf{H}^{-1}\mathbf{s}^{\prime}%
704   \end{equation}
705 + In this way, particles are allowed to diffuse freely in $\mathbf{r}$,
706 + but their minimum images, $\mathbf{r}^{\prime}$ are used to compute
707 + the interatomic forces.
708  
709  
710   \section{\label{oopseSec:IOfiles}Input and Output Files}
711  
712   \subsection{{\sc bass} and Model Files}
713  
714 < Every {\sc oopse} simulation begins with a {\sc bass} file. {\sc bass}
715 < (\underline{B}izarre \underline{A}tom \underline{S}imulation
716 < \underline{S}yntax) is a script syntax that is parsed by {\sc oopse} at
717 < runtime. The {\sc bass} file allows for the user to completely describe the
718 < system they are to simulate, as well as tailor {\sc oopse}'s behavior during
719 < the simulation. {\sc bass} files are denoted with the extension
714 > Every {\sc oopse} simulation begins with a Bizarre Atom Simulation
715 > Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed
716 > by {\sc oopse} at runtime. The {\sc bass} file allows for the user to
717 > completely describe the system they wish to simulate, as well as tailor
718 > {\sc oopse}'s behavior during the simulation. {\sc bass} files are
719 > denoted with the extension
720   \texttt{.bass}, an example file is shown in
721 < Fig.~\ref{fig:bassExample}.
721 > Scheme~\ref{sch:bassExample}.
722  
723 < \begin{figure}
707 < \centering
708 < \framebox[\linewidth]{\rule{0cm}{0.75\linewidth}I'm a {\sc bass} file!}
709 < \caption{Here is an example \texttt{.bass} file}
710 < \label{fig:bassExample}
711 < \end{figure}
723 > \begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}]
724  
725 + molecule{
726 +  name = "Ar";
727 +  nAtoms = 1;
728 +  atom[0]{
729 +    type="Ar";
730 +    position( 0.0, 0.0, 0.0 );
731 +  }
732 + }
733 +
734 + nComponents = 1;
735 + component{
736 +  type = "Ar";
737 +  nMol = 108;
738 + }
739 +
740 + initialConfig = "./argon.init";
741 +
742 + forceField = "LJ";
743 + ensemble = "NVE"; // specify the simulation enesemble
744 + dt = 1.0;         // the time step for integration
745 + runTime = 1e3;    // the total simulation run time
746 + sampleTime = 100; // trajectory file frequency
747 + statusTime = 50;  // statistics file frequency
748 +
749 + \end{lstlisting}
750 +
751   Within the \texttt{.bass} file it is necessary to provide a complete
752   description of the molecule before it is actually placed in the
753 < simulation. The {\sc bass} syntax was originally developed with this goal in
754 < mind, and allows for the specification of all the atoms in a molecular
755 < prototype, as well as any bonds, bends, or torsions. These
753 > simulation. The {\sc bass} syntax was originally developed with this
754 > goal in mind, and allows for the specification of all the atoms in a
755 > molecular prototype, as well as any bonds, bends, or torsions. These
756   descriptions can become lengthy for complex molecules, and it would be
757 < inconvenient to duplicate the simulation at the beginning of each {\sc bass}
758 < script. Addressing this issue {\sc bass} allows for the inclusion of model
759 < files at the top of a \texttt{.bass} file. These model files, denoted
760 < with the \texttt{.mdl} extension, allow the user to describe a
761 < molecular prototype once, then simply include it into each simulation
762 < containing that molecule.
757 > inconvenient to duplicate the simulation at the beginning of each {\sc
758 > bass} script. Addressing this issue {\sc bass} allows for the
759 > inclusion of model files at the top of a \texttt{.bass} file. These
760 > model files, denoted with the \texttt{.mdl} extension, allow the user
761 > to describe a molecular prototype once, then simply include it into
762 > each simulation containing that molecule. Returning to the example in
763 > Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would
764 > be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would
765 > become Scheme~\ref{sch:bassExPrime}.
766 >
767 > \begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}]
768 >
769 > molecule{
770 >  name = "Ar";
771 >  nAtoms = 1;
772 >  atom[0]{
773 >    type="Ar";
774 >    position( 0.0, 0.0, 0.0 );
775 >  }
776 > }
777  
778 + \end{lstlisting}
779 +
780 + \begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}]
781 +
782 + #include "argon.mdl"
783 +
784 + molecule{
785 +  name = "Ar";
786 +  nAtoms = 1;
787 +  atom[0]{
788 +    type="Ar";
789 +    position( 0.0, 0.0, 0.0 );
790 +  }
791 + }
792 +
793 + nComponents = 1;
794 + component{
795 +  type = "Ar";
796 +  nMol = 108;
797 + }
798 +
799 + initialConfig = "./argon.init";
800 +
801 + forceField = "LJ";
802 + ensemble = "NVE";
803 + dt = 1.0;
804 + runTime = 1e3;
805 + sampleTime = 100;
806 + statusTime = 50;
807 +
808 + \end{lstlisting}
809 +
810   \subsection{\label{oopseSec:coordFiles}Coordinate Files}
811  
812   The standard format for storage of a systems coordinates is a modified
813   xyz-file syntax, the exact details of which can be seen in
814 < App.~\ref{appCoordFormat}. As all bonding and molecular information is
815 < stored in the \texttt{.bass} and \texttt{.mdl} files, the coordinate
816 < files are simply the complete set of coordinates for each atom at a
817 < given simulation time.
814 > Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information
815 > is stored in the \texttt{.bass} and \texttt{.mdl} files, the
816 > coordinate files are simply the complete set of coordinates for each
817 > atom at a given simulation time. One important note, although the
818 > simulation propagates the complete rotation matrix, directional
819 > entities are written out using quanternions, to save space in the
820 > output files.
821  
822 < There are three major files used by {\sc oopse} written in the coordinate
823 < format, they are as follows: the initialization file, the simulation
824 < trajectory file, and the final coordinates of the simulation. The
825 < initialization file is necessary for {\sc oopse} to start the simulation
826 < with the proper coordinates. It is typically denoted with the
827 < extension \texttt{.init}. The trajectory file is created at the
828 < beginning of the simulation, and is used to store snapshots of the
829 < simulation at regular intervals. The first frame is a duplication of
830 < the \texttt{.init} file, and each subsequent frame is appended to the
831 < file at an interval specified in the \texttt{.bass} file. The
832 < trajectory file is given the extension \texttt{.dump}. The final
833 < coordinate file is the end of run or \texttt{.eor} file. The
822 > \begin{lstlisting}[float,caption={[The format of the coordinate files]Shows the format of the coordinate files. The fist line is the number of atoms. The second line begins with the time stamp followed by the three $\mathbf{H}$ column vectors. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly angular momentum.},label=sch:dumpFormat]
823 >
824 > nAtoms
825 > time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz;
826 > Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz
827 > Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz
828 > etc...
829 >
830 > \end{lstlisting}
831 >
832 >
833 > There are three major files used by {\sc oopse} written in the
834 > coordinate format, they are as follows: the initialization file
835 > (\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and
836 > the final coordinates of the simulation. The initialization file is
837 > necessary for {\sc oopse} to start the simulation with the proper
838 > coordinates, and is generated before the simulation run. The
839 > trajectory file is created at the beginning of the simulation, and is
840 > used to store snapshots of the simulation at regular intervals. The
841 > first frame is a duplication of the
842 > \texttt{.init} file, and each subsequent frame is appended to the file
843 > at an interval specified in the \texttt{.bass} file with the
844 > \texttt{sampleTime} flag. The final coordinate file is the end of run file. The
845   \texttt{.eor} file stores the final configuration of the system for a
846   given simulation. The file is updated at the same time as the
847 < \texttt{.dump} file. However, it only contains the most recent
847 > \texttt{.dump} file, however, it only contains the most recent
848   frame. In this way, an \texttt{.eor} file may be used as the
849 < initialization file to a second simulation in order to continue or
850 < recover the previous simulation.
849 > initialization file to a second simulation in order to continue a
850 > simulation or recover one from a processor that has crashed during the
851 > course of the run.
852  
853   \subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates}
854  
855 < As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization file
856 < is needed to provide the starting coordinates for a simulation. The
857 < {\sc oopse} package provides a program called \texttt{sysBuilder} to aid in
858 < the creation of the \texttt{.init} file. \texttt{sysBuilder} is {\sc bass}
859 < aware, and will recognize arguments and parameters in the
860 < \texttt{.bass} file that would otherwise be ignored by the
861 < simulation. The program itself is under continual development, and is
763 < offered here as a helper tool only.
855 > As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization
856 > file is needed to provide the starting coordinates for a
857 > simulation. The {\sc oopse} package provides a program called
858 > \texttt{sysBuilder} to aid in the creation of the \texttt{.init}
859 > file. \texttt{sysBuilder} uses {\sc bass}, and will recognize
860 > arguments and parameters in the \texttt{.bass} file that would
861 > otherwise be ignored by the simulation.
862  
863   \subsection{The Statistics File}
864  
865 < The last output file generated by {\sc oopse} is the statistics file. This
866 < file records such statistical quantities as the instantaneous
867 < temperature, volume, pressure, etc. It is written out with the
868 < frequency specified in the \texttt{.bass} file. The file allows the
869 < user to observe the system variables as a function of simulation time
870 < while the simulation is in progress. One useful function the
871 < statistics file serves is to monitor the conserved quantity of a given
872 < simulation ensemble, this allows the user to observe the stability of
873 < the integrator. The statistics file is denoted with the \texttt{.stat}
874 < file extension.
865 > The last output file generated by {\sc oopse} is the statistics
866 > file. This file records such statistical quantities as the
867 > instantaneous temperature, volume, pressure, etc. It is written out
868 > with the frequency specified in the \texttt{.bass} file with the
869 > \texttt{statusTime} keyword. The file allows the user to observe the
870 > system variables as a function of simulation time while the simulation
871 > is in progress. One useful function the statistics file serves is to
872 > monitor the conserved quantity of a given simulation ensemble, this
873 > allows the user to observe the stability of the integrator. The
874 > statistics file is denoted with the \texttt{.stat} file extension.
875  
876   \section{\label{oopseSec:mechanics}Mechanics}
877  
# Line 781 | Line 879 | symplectic splitting method proposed by Dullweber \emp
879  
880   Integration of the equations of motion was carried out using the
881   symplectic splitting method proposed by Dullweber \emph{et
882 < al.}.\cite{Dullweber1997} The reason for this integrator selection
883 < deals with poor energy conservation of rigid body systems using
884 < quaternions. While quaternions work well for orientational motion in
885 < alternate ensembles, the microcanonical ensemble has a constant energy
886 < requirement that is quite sensitive to errors in the equations of
887 < motion. The original implementation of this code utilized quaternions
888 < for rotational motion propagation; however, a detailed investigation
889 < showed that they resulted in a steady drift in the total energy,
890 < something that has been observed by others.\cite{Laird97}
882 > al.}.\cite{Dullweber1997} The reason for the selection of this
883 > integrator, is the poor energy conservation of rigid body systems
884 > using quaternion dynamics. While quaternions work well for
885 > orientational motion in alternate ensembles, the microcanonical
886 > ensemble has a constant energy requirement that is quite sensitive to
887 > errors in the equations of motion. The original implementation of {\sc
888 > oopse} utilized quaternions for rotational motion propagation;
889 > however, a detailed investigation showed that they resulted in a
890 > steady drift in the total energy, something that has been observed by
891 > others.\cite{Laird97}
892  
893   The key difference in the integration method proposed by Dullweber
894 < \emph{et al.} is that the entire rotation matrix is propagated from
895 < one time step to the next. In the past, this would not have been as
896 < feasible a option, being that the rotation matrix for a single body is
897 < nine elements long as opposed to 3 or 4 elements for Euler angles and
898 < quaternions respectively. System memory has become much less of an
899 < issue in recent times, and this has resulted in substantial benefits
900 < in energy conservation. There is still the issue of 5 or 6 additional
802 < elements for describing the orientation of each particle, which will
803 < increase dump files substantially. Simply translating the rotation
804 < matrix into its component Euler angles or quaternions for storage
805 < purposes relieves this burden.
894 > \emph{et al}.~({\sc dlm}) is that the entire rotation matrix is propagated from
895 > one time step to the next. In the past, this would not have been a
896 > feasible option, since the rotation matrix for a single body is nine
897 > elements long as opposed to three or four elements for Euler angles
898 > and quaternions respectively. System memory has become much less of an
899 > issue in recent times, and the {\sc dlm} method has used memory in
900 > exchange for substantial benefits in energy conservation.
901  
902 < The symplectic splitting method allows for Verlet style integration of
903 < both linear and angular motion of rigid bodies. In the integration
904 < method, the orientational propagation involves a sequence of matrix
902 > The {\sc dlm} method allows for Verlet style integration of both
903 > linear and angular motion of rigid bodies. In the integration method,
904 > the orientational propagation involves a sequence of matrix
905   evaluations to update the rotation matrix.\cite{Dullweber1997} These
906 < matrix rotations end up being more costly computationally than the
907 < simpler arithmetic quaternion propagation. With the same time step, a
908 < 1000 SSD particle simulation shows an average 7\% increase in
909 < computation time using the symplectic step method in place of
910 < quaternions. This cost is more than justified when comparing the
911 < energy conservation of the two methods as illustrated in figure
817 < \ref{timestep}.
906 > matrix rotations are more costly computationally than the simpler
907 > arithmetic quaternion propagation. With the same time step, a 1000 SSD
908 > particle simulation shows an average 7\% increase in computation time
909 > using the {\sc dlm} method in place of quaternions. This cost is more
910 > than justified when comparing the energy conservation of the two
911 > methods as illustrated in Fig.~\ref{timestep}.
912  
913   \begin{figure}
914   \centering
915   \includegraphics[width=\linewidth]{timeStep.eps}
916 < \caption{Energy conservation using quaternion based integration versus
917 < the symplectic step method proposed by Dullweber \emph{et al.} with
916 > \caption[Energy conservation for quaternion versus {\sc dlm} dynamics]{Energy conservation using quaternion based integration versus
917 > the {\sc dlm} method with
918   increasing time step. For each time step, the dotted line is total
919 < energy using the symplectic step integrator, and the solid line comes
919 > energy using the {\sc dlm} integrator, and the solid line comes
920   from the quaternion integrator. The larger time step plots are shifted
921   up from the true energy baseline for clarity.}
922   \label{timestep}
923   \end{figure}
924  
925 < In figure \ref{timestep}, the resulting energy drift at various time
926 < steps for both the symplectic step and quaternion integration schemes
925 > In Fig.~\ref{timestep}, the resulting energy drift at various time
926 > steps for both the {\sc dlm} and quaternion integration schemes
927   is compared. All of the 1000 SSD particle simulations started with the
928   same configuration, and the only difference was the method for
929   handling rotational motion. At time steps of 0.1 and 0.5 fs, both
930   methods for propagating particle rotation conserve energy fairly well,
931   with the quaternion method showing a slight energy drift over time in
932   the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
933 < energy conservation benefits of the symplectic step method are clearly
933 > energy conservation benefits of the {\sc dlm} method are clearly
934   demonstrated. Thus, while maintaining the same degree of energy
935   conservation, one can take considerably longer time steps, leading to
936   an overall reduction in computation time.
# Line 845 | Line 939 | four femtoseconds, and as expected, this drift increas
939   time steps up to three femtoseconds. A slight energy drift on the
940   order of 0.012 kcal/mol per nanosecond was observed at a time step of
941   four femtoseconds, and as expected, this drift increases dramatically
942 < with increasing time step. To insure accuracy in the constant energy
849 < simulations, time steps were set at 2 fs and kept at this value for
850 < constant pressure simulations as well.
942 > with increasing time step.
943  
944  
945   \subsection{\label{sec:extended}Extended Systems for other Ensembles}
# Line 856 | Line 948 | constant pressure simulations as well.
948   {\sc oopse} implements a
949  
950  
951 < \subsubsection{\label{oopseSec:noseHooverThermo}Nose-Hoover Thermostatting}
951 > \subsection{\label{oopseSec:noseHooverThermo}Nose-Hoover Thermostatting}
952  
953   To mimic the effects of being in a constant temperature ({\sc nvt})
954   ensemble, {\sc oopse} uses the Nose-Hoover extended system
# Line 883 | Line 975 | set to 1 ps using the {\tt tauThermostat = 1000; } com
975   to values of a few ps.  Within a {\sc bass} file, $\tau_{T}$ could be
976   set to 1 ps using the {\tt tauThermostat = 1000; } command.
977  
978 + \subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond
979 +        Constraints}
980  
981 + In order to satisfy the constraints of fixed bond lengths within {\sc
982 + oopse}, we have implemented the {\sc rattle} algorithm of
983 + Andersen.\cite{andersen83} The algorithm is a velocity verlet
984 + formulation of the {\sc shake} method\cite{ryckaert77} of iteratively
985 + solving the Lagrange multipliers of constraint. The system of lagrange
986 + multipliers allows one to reformulate the equations of motion with
987 + explicit constraint forces on the equations of
988 + motion.\cite{fowles99:lagrange}
989 +
990 + Consider a system described by qoordinates $q_1$ and $q_2$ subject to an
991 + equation of constraint:
992 + \begin{equation}
993 + \sigma(q_1, q_2,t) = 0
994 + \label{oopseEq:lm1}
995 + \end{equation}
996 + The Lagrange formulation of the equations of motion can be written:
997 + \begin{equation}
998 + \delta\int_{t_1}^{t_2}L\, dt =
999 +        \int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i}
1000 +        - \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i}
1001 +        \biggr ) \biggr] \delta q_i \, dt = 0
1002 + \label{oopseEq:lm2}
1003 + \end{equation}
1004 + Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$
1005 + are linked by $\sigma$. However, $\sigma$ is fixed at any given
1006 + instant of time, giving:
1007 + \begin{align}
1008 + \delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 %
1009 +        + \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 \\
1010 + %
1011 + \frac{\partial\sigma}{\partial q_1} \delta q_1 &= %
1012 +        - \frac{\partial\sigma}{\partial q_2} \delta q_2 \\
1013 + %
1014 + \delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / %
1015 +        \frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1
1016 + \end{align}
1017 + Substituted back into Eq.~\ref{oopseEq:lm2},
1018 + \begin{equation}
1019 + \int_{t_1}^{t_2}\biggl [ \biggl(\frac{\partial L}{\partial q_1}
1020 +        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1021 +        \biggr)
1022 +        - \biggl( \frac{\partial L}{\partial q_1}
1023 +        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1024 +        \biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / %
1025 +        \frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0
1026 + \label{oopseEq:lm3}
1027 + \end{equation}
1028 + Leading to,
1029 + \begin{equation}
1030 + \frac{\biggl(\frac{\partial L}{\partial q_1}
1031 +        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1032 +        \biggr)}{\frac{\partial\sigma}{\partial q_1}} =
1033 + \frac{\biggl(\frac{\partial L}{\partial q_2}
1034 +        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2}
1035 +        \biggr)}{\frac{\partial\sigma}{\partial q_2}}
1036 + \label{oopseEq:lm4}
1037 + \end{equation}
1038 + This relation can only be statisfied, if both are equal to a single
1039 + function $-\lambda(t)$,
1040 + \begin{align}
1041 + \frac{\biggl(\frac{\partial L}{\partial q_1}
1042 +        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1043 +        \biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t) \\
1044 + %
1045 + \frac{\partial L}{\partial q_1}
1046 +        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &=
1047 +         -\lambda(t)\,\frac{\partial\sigma}{\partial q_1} \\
1048 + %
1049 + \frac{\partial L}{\partial q_1}
1050 +        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1051 +         + \mathcal{G}_i &= 0
1052 + \end{align}
1053 + Where $\mathcal{G}_i$, the force of constraint on $i$, is:
1054 + \begin{equation}
1055 + \mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1}
1056 + \label{oopseEq:lm5}
1057 + \end{equation}
1058 +
1059 + In a simulation, this would involve the solution of a set of $(m + n)$
1060 + number of equations. Where $m$ is the number of constraints, and $n$
1061 + is the number of constrained coordinates. In practice, this is not
1062 + done, as the matrix inversion neccassary to solve the system of
1063 + equations would be very time consuming to solve. Additionally, the
1064 + numerical error in the solution of the set of $\lambda$'s would be
1065 + compounded by the error inherent in propagating by the Velocity Verlet
1066 + algorithm ($\Delta t^4$). The verlet propagation error is negligible
1067 + in an unconstrained system, as one is interested in the statisitics of
1068 + the run, and not that the run be numerically exact to the ``true''
1069 + integration. This relates back to the ergodic hypothesis that a time
1070 + integral of a valid trajectory will still give the correct enesemble
1071 + average. However, in the case of constraints, if the equations of
1072 + motion leave the ``true'' trajectory, they are departing from the
1073 + constrained surface. The method that is used, is to iteratively solve
1074 + for $\lambda(t)$ at each time step.
1075 +
1076 + In {\sc rattle} the equations of motion are modified subject to the
1077 + following two constraints:
1078 + \begin{align}
1079 + \sigma_{ij}[\mathbf{r}(t)] \equiv
1080 +        [ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2  - d_{ij}^2 &= 0 %
1081 +        \label{oopseEq:c1} \\
1082 + %
1083 + [\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot
1084 +        [\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 \label{oopseEq:c2}
1085 + \end{align}
1086 + Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is
1087 + the constrained distance between atom $i$ and
1088 + $j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to
1089 + be perpindicular to the bond vector, so that the bond can neither grow
1090 + nor shrink. The constrained dynamics equations become:
1091 + \begin{equation}
1092 + m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i
1093 + \label{oopseEq:r1}
1094 + \end{equation}
1095 + Where,
1096 + \begin{equation}
1097 + \mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij}
1098 + \label{oopseEq:r2}
1099 + \end{equation}
1100 +
1101 + In Velocity Verlet, if $\Delta t = h$, the propagation can be written:
1102 + \begin{align}
1103 + \mathbf{r}_i(t+h) &=
1104 +        \mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) +
1105 +        \frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) +
1106 +        \mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1} \\
1107 + %
1108 + \mathbf{\dot{r}}_i(t+h) &=
1109 +        \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i}
1110 +        \Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) +
1111 +        \mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] %
1112 +        \label{oopseEq:vv2}
1113 + \end{align}
1114 +
1115 +
1116 +
1117   \subsection{\label{oopseSec:zcons}Z-Constraint Method}
1118  
1119 < Based on fluctuation-dissipation theorem,\bigskip\ force auto-correlation
1119 > Based on fluctuation-dissipation theorem, a force auto-correlation
1120   method was developed to investigate the dynamics of ions inside the ion
1121   channels.\cite{Roux91} Time-dependent friction coefficient can be calculated
1122   from the deviation of the instantaneous force from its mean force.

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines