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

Comparing trunk/oopsePaper/oopsePaper.tex (file contents):
Revision 1121 by mmeineke, Mon Apr 19 21:00:24 2004 UTC vs.
Revision 1155 by mmeineke, Tue May 11 19:36:28 2004 UTC

# Line 2 | Line 2
2   \usepackage{amsmath}
3   \usepackage{amssymb}
4   \usepackage{endfloat}
5 %\usepackage{berkeley}
5   \usepackage{listings}
6 + \usepackage{palatino}
7   \usepackage{graphicx}
8   \usepackage[ref]{overcite}
9   \usepackage{setspace}
# Line 18 | Line 18
18   \renewcommand\citemid{\ } % no comma in optional reference note
19  
20   \begin{document}
21 < \lstset{language=C,float,frame=tblr,frameround=tttt}
21 > \lstset{language=C,frame=TB,basicstyle=\small,basicstyle=\ttfamily, %
22 >        xleftmargin=0.5in, xrightmargin=0.5in,captionpos=b, %
23 >        abovecaptionskip=0.5cm, belowcaptionskip=0.5cm}
24   \renewcommand{\lstlistingname}{Scheme}
25   \title{{\sc oopse}: An Open Source Object-Oriented Parallel Simulation
26   Engine for Molecular Dynamics}
27  
28 < \author{Matthew A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher J. Fennell and J. Daniel Gezelter\\
28 > \author{Matthew A. Meineke, Charles F. Vardeman II, Teng Lin,\\
29 > Christopher J. Fennell and J. Daniel Gezelter\\
30   Department of Chemistry and Biochemistry\\
31   University of Notre Dame\\
32   Notre Dame, Indiana 46556}
# Line 33 | Line 36 | package ({\sc oopse}) that can perform molecular dynam
36  
37   \begin{abstract}
38   We detail the capabilities of a new open-source parallel simulation
39 < package ({\sc oopse}) that can perform molecular dynamics simulations
40 < on atom types that are missing from other popular packages.  In
41 < particular, {\sc oopse} is capable of performing orientational
39 < dynamics on dipolar systems, and it can handle simulations of metallic
39 > progrm for MD ({\sc oopse}) that can work with  atom types that are missing from other popular packages.  In
40 > particular, {\sc oopse} is capable of performing efficient orientational
41 > dynamics on dipolar or rigid body systems, and it can handle simulations of metallic
42   systems using the embedded atom method ({\sc eam}).
43   \end{abstract}
44  
# Line 75 | Line 77 | We have structured this chapter to first discuss the e
77   those of the community of scientists who contribute back to the
78   project.
79  
80 < We have structured this chapter to first discuss the empirical energy
80 > We have structured this paper to first discuss the empirical energy
81   functions that {\sc oopse } implements in
82   Sec.~\ref{oopseSec:empiricalEnergy}. Following that is a discussion of
83   the various input and output files associated with the package
84   (Sec.~\ref{oopseSec:IOfiles}). Sec.~\ref{oopseSec:mechanics}
85   elucidates the various Molecular Dynamics algorithms {\sc oopse}
86   implements in the integration of the Newtonian equations of
87 < motion. Basic analysis of the trajectories obtained from the
86 < simulation is discussed in Sec.~\ref{oopseSec:props}. Program design
87 > motion.  Program design
88   considerations are presented in Sec.~\ref{oopseSec:design}. And
89   lastly, Sec.~\ref{oopseSec:conclusion} concludes the chapter.
90  
91 < \section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions}
91 > \section{\label{oopseSec:IOfiles}Concepts \& Files}
92  
93 + \subsection{{\sc bass} and Model Files}
94 +
95 + Every {\sc oopse} simulation begins with a Bizarre Atom Simulation
96 + Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed
97 + by {\sc oopse} at runtime. The {\sc bass} file allows for the user to
98 + completely describe the system they wish to simulate, as well as tailor
99 + {\sc oopse}'s behavior during the simulation. {\sc bass} files are
100 + denoted with the extension
101 + \texttt{.bass}, an example file is shown in
102 + Scheme~\ref{sch:bassExample}.
103 +
104 + \begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}]
105 +
106 + molecule{
107 +  name = "Ar";
108 +  nAtoms = 1;
109 +  atom[0]{
110 +    type="Ar";
111 +    position( 0.0, 0.0, 0.0 );
112 +  }
113 + }
114 +
115 + nComponents = 1;
116 + component{
117 +  type = "Ar";
118 +  nMol = 108;
119 + }
120 +
121 + initialConfig = "./argon.init";
122 +
123 + forceField = "LJ";
124 + ensemble = "NVE"; // specify the simulation ensemble
125 + dt = 1.0;         // the time step for integration
126 + runTime = 1e3;    // the total simulation run time
127 + sampleTime = 100; // trajectory file frequency
128 + statusTime = 50;  // statistics file frequency
129 +
130 + \end{lstlisting}
131 +
132 + Within the \texttt{.bass} file it is necessary to provide a complete
133 + description of the molecule before it is actually placed in the
134 + simulation. The {\sc bass} syntax was originally developed with this
135 + goal in mind, and allows for the specification of all the atoms in a
136 + molecular prototype, as well as any bonds, bends, or torsions. These
137 + descriptions can become lengthy for complex molecules, and it would be
138 + inconvenient to duplicate the simulation at the beginning of each {\sc
139 + bass} script. Addressing this issue {\sc bass} allows for the
140 + inclusion of model files at the top of a \texttt{.bass} file. These
141 + model files, denoted with the \texttt{.mdl} extension, allow the user
142 + to describe a molecular prototype once, then simply include it into
143 + each simulation containing that molecule. Returning to the example in
144 + Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would
145 + be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would
146 + become Scheme~\ref{sch:bassExPrime}.
147 +
148 + \begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}]
149 +
150 + molecule{
151 +  name = "Ar";
152 +  nAtoms = 1;
153 +  atom[0]{
154 +    type="Ar";
155 +    position( 0.0, 0.0, 0.0 );
156 +  }
157 + }
158 +
159 + \end{lstlisting}
160 +
161 + \begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}]
162 +
163 + #include "argon.mdl"
164 +
165 + nComponents = 1;
166 + component{
167 +  type = "Ar";
168 +  nMol = 108;
169 + }
170 +
171 + initialConfig = "./argon.init";
172 +
173 + forceField = "LJ";
174 + ensemble = "NVE";
175 + dt = 1.0;
176 + runTime = 1e3;
177 + sampleTime = 100;
178 + statusTime = 50;
179 +
180 + \end{lstlisting}
181 +
182   \subsection{\label{oopseSec:atomsMolecules}Atoms, Molecules and Rigid Bodies}
183  
184   The basic unit of an {\sc oopse} simulation is the atom. The
# Line 198 | Line 288 | molecule{
288    rigidBody[0]{
289      nMembers = 3;
290      members(0, 1, 2);
291 +  }
292 + }
293 + \end{lstlisting}
294 +
295 + \subsection{\label{sec:miscConcepts}Putting a Script Together}
296 +
297 + The actual creation of a {\sc bass} script requires several key components. The first  part of the script needs to be the declaration of all of the molecule prototypes used in the simulation. This is typically done through the inclusion of {\tt .mdl} files. Only the molecules actually present in the simulation need to be declared, however {\sc bass} allows for the declaration of more molecules than are needed. This gives the user the ability to build up a library of commonly used molecules into a single {\tt .mdl} file.
298 +
299 + Once all prototypes are declared, the ordering of the rest of the script is less stringent. Typically, the next to follow the molecular prototypes are the component statements. These statements specify which molecules are present within the simulation. The number of components must first be declared before the first component block statement (an example is seen in Sch.~\ref{sch:bassExPrime}).
300 +
301 + \subsection{\label{oopseSec:coordFiles}Coordinate Files}
302 +
303 + The standard format for storage of a systems coordinates is a modified
304 + xyz-file syntax, the exact details of which can be seen in
305 + Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information
306 + is stored in the \texttt{.bass} and \texttt{.mdl} files, the
307 + coordinate files are simply the complete set of coordinates for each
308 + atom at a given simulation time. One important note, although the
309 + simulation propagates the complete rotation matrix, directional
310 + entities are written out using quanternions, to save space in the
311 + output files.
312 +
313 + \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 $\mathsf{H}$ column vectors. It is important to note, that for extended system ensembles, additional information pertinent to the integrators may be stored on this line as well. The next lines are the atomic coordinates for all atoms in the system. First is the name followed by position, velocity, quanternions, and lastly, body fixed angular momentum.},label=sch:dumpFormat]
314 +
315 + nAtoms
316 + time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz;
317 + Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz
318 + Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz
319 + etc...
320 +
321 + \end{lstlisting}
322 +
323 +
324 + There are three major files used by {\sc oopse} written in the
325 + coordinate format, they are as follows: the initialization file
326 + (\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and
327 + the final coordinates of the simulation (\texttt{.eor}). The initialization file is
328 + necessary for {\sc oopse} to start the simulation with the proper
329 + coordinates, and is generated before the simulation run. The
330 + trajectory file is created at the beginning of the simulation, and is
331 + used to store snapshots of the simulation at regular intervals. The
332 + first frame is a duplication of the
333 + \texttt{.init} file, and each subsequent frame is appended to the file
334 + at an interval specified in the \texttt{.bass} file with the
335 + \texttt{sampleTime} flag. The final coordinate file is the end of run file. The
336 + \texttt{.eor} file stores the final configuration of the system for a
337 + given simulation. The file is updated at the same time as the
338 + \texttt{.dump} file, however, it only contains the most recent
339 + frame. In this way, an \texttt{.eor} file may be used as the
340 + initialization file to a second simulation in order to continue a
341 + simulation or recover one from a processor that has crashed during the
342 + course of the run.
343 +
344 + \subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates}
345 +
346 + As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization
347 + file is needed to provide the starting coordinates for a
348 + simulation.  Several helper programs are provided with {\sc oopse} to illustrate possible build routes. However, as each simulation is different, system creation is left to the end user. The {\tt .init} file must list the atoms in the correct order or {\sc oopse} will give an atom mismatch error.
349 +
350 + The correct ordering of the atoms relies on the ordering of atoms and molecules within the model and {\sc bass} scripts. {\sc oopse} expects the order to comply with the following guidelines:
351 + \begin{enumerate}
352 + \item All of the molecules of the first declared component are given before proceeding to the molecules of the second component, and so on for all declared components.
353 + \item The ordering of the atoms for each molecule follows the order declared in the molecule's declaration within the model file.
354 + \end{enumerate}
355 + An example is given in Scheme~\ref{sch:initEx1} resulting in the {\tt .init} file shown in Scheme~\ref{sch:initEx2}.
356 +
357 + \begin{lstlisting}[float,caption={This scheme illustrates the declaration of the $\text{I}_2$ molecule and the HCl molecule. The two molecules are then included into a simulation.}, label=sch:initEx1]
358 +
359 + molecule{
360 +  name = "I2";
361 +  nAtoms = 2;
362 +  atom[0]{
363 +    type = "I";
364 +  }
365 +  atom[1]{
366 +    type = "I";
367 +  }
368 +  nBonds = 1;
369 +  bond[0]{
370 +    members( 0, 1);
371 +  }
372 + }
373 +
374 + molecule{
375 +  name = "HCl"
376 +  nAtoms = 2;
377 +  atom[0]{
378 +    type = "H";
379 +  }
380 +  atom[1]{
381 +    type = "Cl";
382 +  }
383 +  nBonds = 1;
384 +  bond[0]{
385 +    members( 0, 1);
386    }
387   }
388 +
389 + nComponents = 2;
390 + component{
391 +  type = "HCl";
392 +  nMol = 4;
393 + }
394 + component{
395 +  type = "I2";
396 +  nMol = 1;
397 + }
398 +
399 + initialConfig = "mixture.init";
400 +
401   \end{lstlisting}
402  
403 + \begin{lstlisting}[float,caption={This is the contents of the {\tt mixture.init} file matching the declarations in Scheme~\ref{sch:initEx1}. Note that even though $\text{I}_2$ is declared before HCl, the {\tt .init} file follows the order in which the components were included.},label=sch:initEx2]
404 +
405 + 10
406 + 0.0;  10.0  0.0  0.0;  0.0  10.0  0.0;  0.0  0.0  10.0;
407 + H  ...
408 + Cl ...
409 + H  ...
410 + Cl ...
411 + H  ...
412 + Cl ...
413 + H  ...
414 + Cl ...
415 + I  ...
416 + I  ...
417 +
418 + \end{lstlisting}
419 +
420 +
421 + \subsection{The Statistics File}
422 +
423 + The last output file generated by {\sc oopse} is the statistics
424 + file. This file records such statistical quantities as the
425 + instantaneous temperature, volume, pressure, etc. It is written out
426 + with the frequency specified in the \texttt{.bass} file with the
427 + \texttt{statusTime} keyword. The file allows the user to observe the
428 + system variables as a function of simulation time while the simulation
429 + is in progress. One useful function the statistics file serves is to
430 + monitor the conserved quantity of a given simulation ensemble, this
431 + allows the user to observe the stability of the integrator. The
432 + statistics file is denoted with the \texttt{.stat} file extension.
433 +
434 +
435 + \section{\label{oopseSec:empiricalEnergy}The Empirical Energy Functions}
436 +
437 + \
438   \subsection{\label{sec:LJPot}The Lennard Jones Force Field}
439  
440   The most basic force field implemented in {\sc oopse} is the
# Line 691 | Line 924 | the inter-atomic forces.
924   the inter-atomic forces.
925  
926  
694 \section{\label{oopseSec:IOfiles}Input and Output Files}
927  
696 \subsection{{\sc bass} and Model Files}
697
698 Every {\sc oopse} simulation begins with a Bizarre Atom Simulation
699 Syntax ({\sc bass}) file. {\sc bass} is a script syntax that is parsed
700 by {\sc oopse} at runtime. The {\sc bass} file allows for the user to
701 completely describe the system they wish to simulate, as well as tailor
702 {\sc oopse}'s behavior during the simulation. {\sc bass} files are
703 denoted with the extension
704 \texttt{.bass}, an example file is shown in
705 Scheme~\ref{sch:bassExample}.
706
707 \begin{lstlisting}[float,caption={[An example of a complete {\sc bass} file] An example showing a complete {\sc bass} file.},label={sch:bassExample}]
708
709 molecule{
710  name = "Ar";
711  nAtoms = 1;
712  atom[0]{
713    type="Ar";
714    position( 0.0, 0.0, 0.0 );
715  }
716 }
717
718 nComponents = 1;
719 component{
720  type = "Ar";
721  nMol = 108;
722 }
723
724 initialConfig = "./argon.init";
725
726 forceField = "LJ";
727 ensemble = "NVE"; // specify the simulation ensemble
728 dt = 1.0;         // the time step for integration
729 runTime = 1e3;    // the total simulation run time
730 sampleTime = 100; // trajectory file frequency
731 statusTime = 50;  // statistics file frequency
732
733 \end{lstlisting}
734
735 Within the \texttt{.bass} file it is necessary to provide a complete
736 description of the molecule before it is actually placed in the
737 simulation. The {\sc bass} syntax was originally developed with this
738 goal in mind, and allows for the specification of all the atoms in a
739 molecular prototype, as well as any bonds, bends, or torsions. These
740 descriptions can become lengthy for complex molecules, and it would be
741 inconvenient to duplicate the simulation at the beginning of each {\sc
742 bass} script. Addressing this issue {\sc bass} allows for the
743 inclusion of model files at the top of a \texttt{.bass} file. These
744 model files, denoted with the \texttt{.mdl} extension, allow the user
745 to describe a molecular prototype once, then simply include it into
746 each simulation containing that molecule. Returning to the example in
747 Scheme~\ref{sch:bassExample}, the \texttt{.mdl} file's contents would
748 be Scheme~\ref{sch:mdlExample}, and the new \texttt{.bass} file would
749 become Scheme~\ref{sch:bassExPrime}.
750
751 \begin{lstlisting}[float,caption={An example \texttt{.mdl} file.},label={sch:mdlExample}]
752
753 molecule{
754  name = "Ar";
755  nAtoms = 1;
756  atom[0]{
757    type="Ar";
758    position( 0.0, 0.0, 0.0 );
759  }
760 }
761
762 \end{lstlisting}
763
764 \begin{lstlisting}[float,caption={Revised {\sc bass} example.},label={sch:bassExPrime}]
765
766 #include "argon.mdl"
767
768 nComponents = 1;
769 component{
770  type = "Ar";
771  nMol = 108;
772 }
773
774 initialConfig = "./argon.init";
775
776 forceField = "LJ";
777 ensemble = "NVE";
778 dt = 1.0;
779 runTime = 1e3;
780 sampleTime = 100;
781 statusTime = 50;
782
783 \end{lstlisting}
784
785 \subsection{\label{oopseSec:coordFiles}Coordinate Files}
786
787 The standard format for storage of a systems coordinates is a modified
788 xyz-file syntax, the exact details of which can be seen in
789 Scheme~\ref{sch:dumpFormat}. As all bonding and molecular information
790 is stored in the \texttt{.bass} and \texttt{.mdl} files, the
791 coordinate files are simply the complete set of coordinates for each
792 atom at a given simulation time. One important note, although the
793 simulation propagates the complete rotation matrix, directional
794 entities are written out using quanternions, to save space in the
795 output files.
796
797 \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 $\mathsf{H}$ column vectors. It is important to note, that for extended system ensembles, additional information pertinent to the integrators may be stored on this line as well. 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 velocities.},label=sch:dumpFormat]
798
799 nAtoms
800 time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz;
801 Name1 x y z vx vy vz q0 q1 q2 q3 jx jy jz
802 Name2 x y z vx vy vz q0 q1 q2 q3 jx jy jz
803 etc...
804
805 \end{lstlisting}
806
807
808 There are three major files used by {\sc oopse} written in the
809 coordinate format, they are as follows: the initialization file
810 (\texttt{.init}), the simulation trajectory file (\texttt{.dump}), and
811 the final coordinates of the simulation. The initialization file is
812 necessary for {\sc oopse} to start the simulation with the proper
813 coordinates, and is generated before the simulation run. The
814 trajectory file is created at the beginning of the simulation, and is
815 used to store snapshots of the simulation at regular intervals. The
816 first frame is a duplication of the
817 \texttt{.init} file, and each subsequent frame is appended to the file
818 at an interval specified in the \texttt{.bass} file with the
819 \texttt{sampleTime} flag. The final coordinate file is the end of run file. The
820 \texttt{.eor} file stores the final configuration of the system for a
821 given simulation. The file is updated at the same time as the
822 \texttt{.dump} file, however, it only contains the most recent
823 frame. In this way, an \texttt{.eor} file may be used as the
824 initialization file to a second simulation in order to continue a
825 simulation or recover one from a processor that has crashed during the
826 course of the run.
827
828 \subsection{\label{oopseSec:initCoords}Generation of Initial Coordinates}
829
830 As was stated in Sec.~\ref{oopseSec:coordFiles}, an initialization
831 file is needed to provide the starting coordinates for a
832 simulation. The {\sc oopse} package provides several system building
833 programs to aid in the creation of the \texttt{.init}
834 file. The programs use {\sc bass}, and will recognize
835 arguments and parameters in the \texttt{.bass} file that would
836 otherwise be ignored by the simulation.
837
838 \subsection{The Statistics File}
839
840 The last output file generated by {\sc oopse} is the statistics
841 file. This file records such statistical quantities as the
842 instantaneous temperature, volume, pressure, etc. It is written out
843 with the frequency specified in the \texttt{.bass} file with the
844 \texttt{statusTime} keyword. The file allows the user to observe the
845 system variables as a function of simulation time while the simulation
846 is in progress. One useful function the statistics file serves is to
847 monitor the conserved quantity of a given simulation ensemble, this
848 allows the user to observe the stability of the integrator. The
849 statistics file is denoted with the \texttt{.stat} file extension.
850
928   \section{\label{oopseSec:mechanics}Mechanics}
929  
930   \subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the
# Line 1635 | Line 1712 | simulations).
1712   orientational anisotropy in the system (i.e. in lipid bilayer
1713   simulations).
1714  
1715 < \subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond
1715 > \subsection{\label{sec:constraints}Constraint Methods}
1716 >
1717 > \subsubsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond
1718          Constraints}
1719  
1720   In order to satisfy the constraints of fixed bond lengths within {\sc
1721   oopse}, we have implemented the {\sc rattle} algorithm of
1722   Andersen.\cite{andersen83} The algorithm is a velocity verlet
1723   formulation of the {\sc shake} method\cite{ryckaert77} of iteratively
1724 < solving the Lagrange multipliers of constraint. The system of Lagrange
1646 < multipliers allows one to reformulate the equations of motion with
1647 < explicit constraint forces.\cite{fowles99:lagrange}
1724 > solving the Lagrange multipliers of constraint.
1725  
1726 < Consider a system described by coordinates $q_1$ and $q_2$ subject to an
1650 < equation of constraint:
1651 < \begin{equation}
1652 < \sigma(q_1, q_2,t) = 0
1653 < \label{oopseEq:lm1}
1654 < \end{equation}
1655 < The Lagrange formulation of the equations of motion can be written:
1656 < \begin{equation}
1657 < \delta\int_{t_1}^{t_2}L\, dt =
1658 <        \int_{t_1}^{t_2} \sum_i \biggl [ \frac{\partial L}{\partial q_i}
1659 <        - \frac{d}{dt}\biggl(\frac{\partial L}{\partial \dot{q}_i}
1660 <        \biggr ) \biggr] \delta q_i \, dt = 0.
1661 < \label{oopseEq:lm2}
1662 < \end{equation}
1663 < Here, $\delta q_i$ is not independent for each $q$, as $q_1$ and $q_2$
1664 < are linked by $\sigma$. However, $\sigma$ is fixed at any given
1665 < instant of time, giving:
1666 < \begin{align}
1667 < \delta\sigma &= \biggl( \frac{\partial\sigma}{\partial q_1} \delta q_1 %
1668 <        + \frac{\partial\sigma}{\partial q_2} \delta q_2 \biggr) = 0 ,\\
1669 < %
1670 < \frac{\partial\sigma}{\partial q_1} \delta q_1 &= %
1671 <        - \frac{\partial\sigma}{\partial q_2} \delta q_2, \\
1672 < %
1673 < \delta q_2 &= - \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / %
1674 <        \frac{\partial\sigma}{\partial q_2} \biggr) \delta q_1.
1675 < \end{align}
1676 < Substituted back into Eq.~\ref{oopseEq:lm2},
1677 < \begin{equation}
1678 < \int_{t_1}^{t_2}\biggl [ \biggl(\frac{\partial L}{\partial q_1}
1679 <        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1680 <        \biggr)
1681 <        - \biggl( \frac{\partial L}{\partial q_1}
1682 <        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1683 <        \biggr) \biggl(\frac{\partial\sigma}{\partial q_1} \bigg / %
1684 <        \frac{\partial\sigma}{\partial q_2} \biggr)\biggr] \delta q_1 \, dt = 0.
1685 < \label{oopseEq:lm3}
1686 < \end{equation}
1687 < Leading to,
1688 < \begin{equation}
1689 < \frac{\biggl(\frac{\partial L}{\partial q_1}
1690 <        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1691 <        \biggr)}{\frac{\partial\sigma}{\partial q_1}} =
1692 < \frac{\biggl(\frac{\partial L}{\partial q_2}
1693 <        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_2}
1694 <        \biggr)}{\frac{\partial\sigma}{\partial q_2}}.
1695 < \label{oopseEq:lm4}
1696 < \end{equation}
1697 < This relation can only be statisfied, if both are equal to a single
1698 < function $-\lambda(t)$,
1699 < \begin{align}
1700 < \frac{\biggl(\frac{\partial L}{\partial q_1}
1701 <        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1702 <        \biggr)}{\frac{\partial\sigma}{\partial q_1}} &= -\lambda(t), \\
1703 < %
1704 < \frac{\partial L}{\partial q_1}
1705 <        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1} &=
1706 <         -\lambda(t)\,\frac{\partial\sigma}{\partial q_1} ,\\
1707 < %
1708 < \frac{\partial L}{\partial q_1}
1709 <        - \frac{d}{dt}\,\frac{\partial L}{\partial \dot{q}_1}
1710 <         + \mathcal{G}_i &= 0,
1711 < \end{align}
1712 < where $\mathcal{G}_i$, the force of constraint on $i$, is:
1713 < \begin{equation}
1714 < \mathcal{G}_i = \lambda(t)\,\frac{\partial\sigma}{\partial q_1}.
1715 < \label{oopseEq:lm5}
1716 < \end{equation}
1717 <
1718 < In a simulation, this would involve the solution of a set of $(m + n)$
1719 < number of equations. Where $m$ is the number of constraints, and $n$
1720 < is the number of constrained coordinates. In practice, this is not
1721 < done, as the matrix inversion necessary to solve the system of
1722 < equations would be very time consuming to solve. Additionally, the
1723 < numerical error in the solution of the set of $\lambda$'s would be
1724 < compounded by the error inherent in propagating by the Velocity Verlet
1725 < algorithm ($\Delta t^4$). The Verlet propagation error is negligible
1726 < in an unconstrained system, as one is interested in the statistics of
1727 < the run, and not that the run be numerically exact to the ``true''
1728 < integration. This relates back to the ergodic hypothesis that a time
1729 < integral of a valid trajectory will still give the correct ensemble
1730 < average. However, in the case of constraints, if the equations of
1731 < motion leave the ``true'' trajectory, they are departing from the
1732 < constrained surface. The method that is used, is to iteratively solve
1733 < for $\lambda(t)$ at each time step.
1734 <
1735 < In {\sc rattle} the equations of motion are modified subject to the
1736 < following two constraints:
1737 < \begin{align}
1738 < \sigma_{ij}[\mathbf{r}(t)] \equiv
1739 <        [ \mathbf{r}_i(t) - \mathbf{r}_j(t)]^2  - d_{ij}^2 &= 0 %
1740 <        \label{oopseEq:c1}, \\
1741 < %
1742 < [\mathbf{\dot{r}}_i(t) - \mathbf{\dot{r}}_j(t)] \cdot
1743 <        [\mathbf{r}_i(t) - \mathbf{r}_j(t)] &= 0 .\label{oopseEq:c2}
1744 < \end{align}
1745 < Eq.~\ref{oopseEq:c1} is the set of bond constraints, where $d_{ij}$ is
1746 < the constrained distance between atom $i$ and
1747 < $j$. Eq.~\ref{oopseEq:c2} constrains the velocities of $i$ and $j$ to
1748 < be perpendicular to the bond vector, so that the bond can neither grow
1749 < nor shrink. The constrained dynamics equations become:
1750 < \begin{equation}
1751 < m_i \mathbf{\ddot{r}}_i = \mathbf{F}_i + \mathbf{\mathcal{G}}_i,
1752 < \label{oopseEq:r1}
1753 < \end{equation}
1754 < where,$\mathbf{\mathcal{G}}_i$ are the forces of constraint on $i$,
1755 < and are defined:
1756 < \begin{equation}
1757 < \mathbf{\mathcal{G}}_i = - \sum_j \lambda_{ij}(t)\,\nabla \sigma_{ij}.
1758 < \label{oopseEq:r2}
1759 < \end{equation}
1760 <
1761 < In Velocity Verlet, if $\Delta t = h$, the propagation can be written:
1762 < \begin{align}
1763 < \mathbf{r}_i(t+h) &=
1764 <        \mathbf{r}_i(t) + h\mathbf{\dot{r}}(t) +
1765 <        \frac{h^2}{2m_i}\,\Bigl[ \mathbf{F}_i(t) +
1766 <        \mathbf{\mathcal{G}}_{Ri}(t) \Bigr] \label{oopseEq:vv1}, \\
1767 < %
1768 < \mathbf{\dot{r}}_i(t+h) &=
1769 <        \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i}
1770 <        \Bigl[ \mathbf{F}_i(t) + \mathbf{\mathcal{G}}_{Ri}(t) +
1771 <        \mathbf{F}_i(t+h) + \mathbf{\mathcal{G}}_{Vi}(t+h) \Bigr] ,%
1772 <        \label{oopseEq:vv2}
1773 < \end{align}
1774 < where:
1775 < \begin{align}
1776 < \mathbf{\mathcal{G}}_{Ri}(t) &=
1777 <        -2 \sum_j \lambda_{Rij}(t) \mathbf{r}_{ij}(t) ,\\
1778 < %
1779 < \mathbf{\mathcal{G}}_{Vi}(t+h) &=
1780 <        -2 \sum_j \lambda_{Vij}(t+h) \mathbf{r}(t+h).
1781 < \end{align}
1782 < Next, define:
1783 < \begin{align}
1784 < g_{ij} &= h \lambda_{Rij}(t) ,\\
1785 < k_{ij} &= h \lambda_{Vij}(t+h), \\
1786 < \mathbf{q}_i &= \mathbf{\dot{r}}_i(t) + \frac{h}{2m_i} \mathbf{F}_i(t)
1787 <        - \frac{1}{m_i}\sum_j g_{ij}\mathbf{r}_{ij}(t).
1788 < \end{align}
1789 < Using these definitions, Eq.~\ref{oopseEq:vv1} and \ref{oopseEq:vv2}
1790 < can be rewritten as,
1791 < \begin{align}
1792 < \mathbf{r}_i(t+h) &= \mathbf{r}_i(t) + h \mathbf{q}_i ,\\
1793 < %
1794 < \mathbf{\dot{r}}(t+h) &= \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h)
1795 <        -\frac{1}{m_i}\sum_j k_{ij} \mathbf{r}_{ij}(t+h).
1796 < \end{align}
1797 <
1798 < To integrate the equations of motion, the {\sc rattle} algorithm first
1799 < solves for $\mathbf{r}(t+h)$. Let,
1800 < \begin{equation}
1801 < \mathbf{q}_i = \mathbf{\dot{r}}(t) + \frac{h}{2m_i}\mathbf{F}_i(t).
1802 < \end{equation}
1803 < Here $\mathbf{q}_i$ corresponds to an initial unconstrained move. Next
1804 < pick a constraint $j$, and let,
1805 < \begin{equation}
1806 < \mathbf{s} = \mathbf{r}_i(t) + h\mathbf{q}_i(t)
1807 <        - \mathbf{r}_j(t) + h\mathbf{q}_j(t).
1808 < \label{oopseEq:ra1}
1809 < \end{equation}
1810 < If
1811 < \begin{equation}
1812 < \Big| |\mathbf{s}|^2 - d_{ij}^2 \Big| > \text{tolerance},
1813 < \end{equation}
1814 < then the constraint is unsatisfied, and corrections are made to the
1815 < positions. First we define a test corrected configuration as,
1816 < \begin{align}
1817 < \mathbf{r}_i^T(t+h) = \mathbf{r}_i(t) + h\biggl[\mathbf{q}_i -
1818 <        g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_i} \biggr] ,\\
1819 < %
1820 < \mathbf{r}_j^T(t+h) = \mathbf{r}_j(t) + h\biggl[\mathbf{q}_j +
1821 <        g_{ij}\,\frac{\mathbf{r}_{ij}(t)}{m_j} \biggr].
1822 < \end{align}
1823 < And we chose $g_{ij}$ such that, $|\mathbf{r}_i^T - \mathbf{r}_j^T|^2
1824 < = d_{ij}^2$. Solving the quadratic for $g_{ij}$ we obtain the
1825 < approximation,
1826 < \begin{equation}
1827 < g_{ij} = \frac{(s^2 - d^2)}{2h[\mathbf{s}\cdot\mathbf{r}_{ij}(t)]
1828 <        (\frac{1}{m_i} + \frac{1}{m_j})}.
1829 < \end{equation}
1830 < Although not an exact solution for $g_{ij}$, as this is an iterative
1831 < scheme overall, the eventual solution will converge. With a trial
1832 < $g_{ij}$, the new $\mathbf{q}$'s become,
1833 < \begin{align}
1834 < \mathbf{q}_i &= \mathbf{q}^{\text{old}}_i - g_{ij}\,
1835 <        \frac{\mathbf{r}_{ij}(t)}{m_i} ,\\
1836 < %
1837 < \mathbf{q}_j &= \mathbf{q}^{\text{old}}_j + g_{ij}\,
1838 <        \frac{\mathbf{r}_{ij}(t)}{m_j} .
1839 < \end{align}
1840 < The whole algorithm is then repeated from Eq.~\ref{oopseEq:ra1} until
1841 < all constraints are satisfied.
1842 <
1843 < The second step of {\sc rattle}, is to then update the velocities. The
1844 < step starts with,
1845 < \begin{equation}
1846 < \mathbf{\dot{r}}_i(t+h) = \mathbf{q}_i + \frac{h}{2m_i}\mathbf{F}_i(t+h).
1847 < \end{equation}
1848 < Next we pick a constraint $j$, and calculate the dot product $\ell$.
1849 < \begin{equation}
1850 < \ell = \mathbf{r}_{ij}(t+h) \cdot \mathbf{\dot{r}}_{ij}(t+h).
1851 < \label{oopseEq:rv1}
1852 < \end{equation}
1853 < Here if constraint Eq.~\ref{oopseEq:c2} holds, $\ell$ should be
1854 < zero. Therefore if $\ell$ is greater than some tolerance, then
1855 < corrections are made to the $i$ and $j$ velocities.
1856 < \begin{align}
1857 < \mathbf{\dot{r}}_i^T &= \mathbf{\dot{r}}_i(t+h) - k_{ij}
1858 <        \frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_i}, \\
1859 < %
1860 < \mathbf{\dot{r}}_j^T &= \mathbf{\dot{r}}_j(t+h) + k_{ij}
1861 <        \frac{\mathbf{\dot{r}}_{ij}(t+h)}{m_j}.
1862 < \end{align}
1863 < Like in the previous step, we select a value for $k_{ij}$ such that
1864 < $\ell$ is zero.
1865 < \begin{equation}
1866 < k_{ij} = \frac{\ell}{d^2_{ij}(\frac{1}{m_i} + \frac{1}{m_j})}.
1867 < \end{equation}
1868 < The test velocities, $\mathbf{\dot{r}}^T_i$ and
1869 < $\mathbf{\dot{r}}^T_j$, then replace their respective velocities, and
1870 < the algorithm is iterated from Eq.~\ref{oopseEq:rv1} until all
1871 < constraints are satisfied.
1726 > \subsubsection{\label{oopseSec:zcons}Z-Constraint Method}
1727  
1873
1874 \subsection{\label{oopseSec:zcons}Z-Constraint Method}
1875
1728   Based on the fluctuation-dissipation theorem, a force auto-correlation
1729   method was developed by Roux and Karplus to investigate the dynamics
1730   of ions inside ion channels.\cite{Roux91} The time-dependent friction
# Line 1959 | Line 1811 | F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\part
1811          -k_{\text{Harmonic}}(z(t)-z_{\text{cons}}).
1812   \end{equation}
1813  
1962 \section{\label{oopseSec:props}Trajectory Analysis}
1963
1964 \subsection{\label{oopseSec:staticProps}Static Property Analysis}
1965
1966 The static properties of the trajectories are analyzed with the
1967 program \texttt{staticProps}. The code is capable of calculating a
1968 number of pair correlations between species A and B. Some of which
1969 only apply to directional entities. The summary of pair correlations
1970 can be found in Table~\ref{oopseTb:gofrs}
1971
1972 \begin{table}
1973 \caption{THE DIFFERENT PAIR CORRELATIONS IN \texttt{staticProps}}
1974 \label{oopseTb:gofrs}
1975 \begin{center}
1976 \begin{tabular}{|l|c|c|}
1977 \hline
1978 Name      & Equation & Directional Atom \\ \hline
1979 $g_{\text{AB}}(r)$              & Eq.~\ref{eq:gofr}         & neither \\ \hline
1980 $g_{\text{AB}}(r, \cos \theta)$ & Eq.~\ref{eq:gofrCosTheta} & A \\ \hline
1981 $g_{\text{AB}}(r, \cos \omega)$ & Eq.~\ref{eq:gofrCosOmega} & both \\ \hline
1982 $g_{\text{AB}}(x, y, z)$        & Eq.~\ref{eq:gofrXYZ}      & neither \\ \hline
1983 $\langle \cos \omega \rangle_{\text{AB}}(r)$ & Eq.~\ref{eq:cosOmegaOfR} &%
1984        both \\ \hline
1985 \end{tabular}
1986 \begin{minipage}{\linewidth}
1987 \centering
1988 \vspace{2mm}
1989 The third column specifies which atom, if any, need be a directional entity.
1990 \end{minipage}
1991 \end{center}
1992 \end{table}
1993
1994 The first pair correlation, $g_{\text{AB}}(r)$, is defined as follows:
1995 \begin{equation}
1996 g_{\text{AB}}(r) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle %%
1997        \sum_{i \in \text{A}} \sum_{j \in \text{B}} %%
1998        \delta( r - |\mathbf{r}_{ij}|) \rangle, \label{eq:gofr}
1999 \end{equation}
2000 where $\mathbf{r}_{ij}$ is the vector
2001 \begin{equation*}
2002 \mathbf{r}_{ij} = \mathbf{r}_j - \mathbf{r}_i, \notag
2003 \end{equation*}
2004 and $\frac{V}{N_{\text{A}}N_{\text{B}}}$ normalizes the average over
2005 the expected pair density at a given $r$.
2006
2007 The next two pair correlations, $g_{\text{AB}}(r, \cos \theta)$ and
2008 $g_{\text{AB}}(r, \cos \omega)$, are similar in that they are both two
2009 dimensional histograms. Both use $r$ for the primary axis then a
2010 $\cos$ for the secondary axis ($\cos \theta$ for
2011 Eq.~\ref{eq:gofrCosTheta} and $\cos \omega$ for
2012 Eq.~\ref{eq:gofrCosOmega}). This allows for the investigator to
2013 correlate alignment on directional entities. $g_{\text{AB}}(r, \cos
2014 \theta)$ is defined as follows:
2015 \begin{equation}
2016 g_{\text{AB}}(r, \cos \theta) = \frac{V}{N_{\text{A}}N_{\text{B}}}\langle  
2017 \sum_{i \in \text{A}} \sum_{j \in \text{B}}  
2018 \delta( \cos \theta - \cos \theta_{ij})
2019 \delta( r - |\mathbf{r}_{ij}|) \rangle.
2020 \label{eq:gofrCosTheta}
2021 \end{equation}
2022 Here
2023 \begin{equation*}
2024 \cos \theta_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{r}}_{ij},
2025 \end{equation*}
2026 where $\mathbf{\hat{i}}$ is the unit directional vector of species $i$
2027 and $\mathbf{\hat{r}}_{ij}$ is the unit vector associated with vector
2028 $\mathbf{r}_{ij}$.
2029
2030 The second two dimensional histogram is of the form:
2031 \begin{equation}
2032 g_{\text{AB}}(r, \cos \omega) =
2033        \frac{V}{N_{\text{A}}N_{\text{B}}}\langle
2034        \sum_{i \in \text{A}} \sum_{j \in \text{B}}
2035        \delta( \cos \omega - \cos \omega_{ij})
2036        \delta( r - |\mathbf{r}_{ij}|) \rangle. \label{eq:gofrCosOmega}
2037 \end{equation}
2038 Here
2039 \begin{equation*}
2040 \cos \omega_{ij} = \mathbf{\hat{i}} \cdot \mathbf{\hat{j}}.
2041 \end{equation*}
2042 Again, $\mathbf{\hat{i}}$ and $\mathbf{\hat{j}}$ are the unit
2043 directional vectors of species $i$ and $j$.
2044
2045 The static analysis code is also cable of calculating a three
2046 dimensional pair correlation of the form:
2047 \begin{equation}\label{eq:gofrXYZ}
2048 g_{\text{AB}}(x, y, z) =
2049        \frac{V}{N_{\text{A}}N_{\text{B}}}\langle
2050        \sum_{i \in \text{A}} \sum_{j \in \text{B}}
2051        \delta( x - x_{ij})
2052        \delta( y - y_{ij})
2053        \delta( z - z_{ij}) \rangle,
2054 \end{equation}
2055 where $x_{ij}$, $y_{ij}$, and $z_{ij}$ are the $x$, $y$, and $z$
2056 components respectively of vector $\mathbf{r}_{ij}$.
2057
2058 The final pair correlation is similar to
2059 Eq.~\ref{eq:gofrCosOmega}. $\langle \cos \omega
2060 \rangle_{\text{AB}}(r)$ is calculated in the following way:
2061 \begin{equation}\label{eq:cosOmegaOfR}
2062 \langle \cos \omega \rangle_{\text{AB}}(r)  =
2063        \langle \sum_{i \in \text{A}} \sum_{j \in \text{B}}
2064        (\cos \omega_{ij}) \delta( r - |\mathbf{r}_{ij}|) \rangle.
2065 \end{equation}
2066 Here $\cos \omega_{ij}$ is defined in the same way as in
2067 Eq.~\ref{eq:gofrCosOmega}. This equation is a single dimensional pair
2068 correlation that gives the average correlation of two directional
2069 entities as a function of their distance from each other.
2070
2071 \subsection{\label{dynamicProps}Dynamic Property Analysis}
2072
2073 The dynamic properties of a trajectory are calculated with the program
2074 \texttt{dynamicProps}. The program calculates the following properties:
2075 \begin{gather}
2076 \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle, \label{eq:rms}\\
2077 \langle \mathbf{v}(t) \cdot \mathbf{v}(0) \rangle, \label{eq:velCorr} \\
2078 \langle \mathbf{j}(t) \cdot \mathbf{j}(0) \rangle. \label{eq:angularVelCorr}
2079 \end{gather}
2080
2081 Eq.~\ref{eq:rms} is the root mean square displacement function. Which
2082 allows one to observe the average displacement of an atom as a
2083 function of time. The quantity is useful when calculating diffusion
2084 coefficients because of the Einstein Relation, which is valid at long
2085 times.\cite{allen87:csl}
2086 \begin{equation}
2087 2tD = \langle | \mathbf{r}(t) - \mathbf{r}(0) |^2 \rangle.
2088 \label{oopseEq:einstein}
2089 \end{equation}
2090
2091 Eq.~\ref{eq:velCorr} and \ref{eq:angularVelCorr} are the translational
2092 velocity and angular velocity correlation functions respectively. The
2093 latter is only applicable to directional species in the
2094 simulation. The velocity autocorrelation functions are useful when
2095 determining vibrational information about the system of interest.
2096
1814   \section{\label{oopseSec:design}Program Design}
1815  
1816   \subsection{\label{sec:architecture} {\sc oopse} Architecture}
# Line 2191 | Line 1908 | atoms.\cite{plimpton95}
1908   and favorably compete with spatial methods up to 100,000
1909   atoms.\cite{plimpton95}
1910  
2194 \subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis}
2195
2196 For large simulations, the trajectory files can sometimes reach sizes
2197 in excess of several gigabytes. In order to effectively analyze that
2198 amount of data, two memory management schemes have been devised for
2199 \texttt{staticProps} and for \texttt{dynamicProps}. The first scheme,
2200 developed for \texttt{staticProps}, is the simplest. As each frame's
2201 statistics are calculated independent of each other, memory is
2202 allocated for each frame, then freed once correlation calculations are
2203 complete for the snapshot. To prevent multiple passes through a
2204 potentially large file, \texttt{staticProps} is capable of calculating
2205 all requested correlations per frame with only a single pair loop in
2206 each frame and a single read of the file.
2207
2208 The second, more advanced memory scheme, is used by
2209 \texttt{dynamicProps}. Here, the program must have multiple frames in
2210 memory to calculate time dependent correlations. In order to prevent a
2211 situation where the program runs out of memory due to large
2212 trajectories, the user is able to specify that the trajectory be read
2213 in blocks. The number of frames in each block is specified by the
2214 user, and upon reading a block of the trajectory,
2215 \texttt{dynamicProps} will calculate all of the time correlation frame
2216 pairs within the block. After in-block correlations are complete, a
2217 second block of the trajectory is read, and the cross correlations are
2218 calculated between the two blocks. This second block is then freed and
2219 then incremented and the process repeated until the end of the
2220 trajectory. Once the end is reached, the first block is freed then
2221 incremented, and the again the internal time correlations are
2222 calculated. The algorithm with the second block is then repeated with
2223 the new origin block, until all frame pairs have been correlated in
2224 time. This process is illustrated in
2225 Fig.~\ref{oopseFig:dynamicPropsMemory}.
2226
2227 %\begin{figure}
2228 %\centering
2229 %\includegraphics[width=\linewidth]{dynamicPropsMem.eps}
2230 %\caption[A representation of the block correlations in \texttt{dynamicProps}]{This diagram illustrates the memory management used by \texttt{dynamicProps}, which follows the scheme: $\sum^{N_{\text{memory blocks}}}_{i=1}[ \operatorname{self}(i) + \sum^{N_{\text{memory blocks}}}_{j>i} \operatorname{cross}(i,j)]$. The shaded region represents the self correlation of the memory block, and the open blocks are read one at a time and the cross correlations between blocks are calculated.}
2231 %\label{oopseFig:dynamicPropsMemory}
2232 %\end{figure}
2233
1911   \section{\label{oopseSec:conclusion}Conclusion}
1912  
1913   We have presented the design and implementation of our open source
# Line 2252 | Line 1929 | The authors would like to thank espresso for fueling t
1929  
1930   \newpage
1931   \section{Acknowledgments}
1932 < The authors would like to thank espresso for fueling this work, and
2256 < would also like to send a special acknowledgement to single malt
2257 < scotch for its wonderful calming effects and its ability to make the
2258 < troubles of the world float away.
1932 > The authors would like to thank the Notre Dame BoB computer cluster where much of this project was tested. Additionally, the authors would like to acknowledge their funding from {\LARGE FIX ME}.
1933  
1934   \bibliographystyle{achemso}
1935   \bibliography{oopsePaper}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines