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 1088 by mmeineke, Fri Mar 5 22:16:34 2004 UTC vs.
Revision 1089 by mmeineke, Mon Mar 8 22:15:54 2004 UTC

# Line 18 | Line 18 | source simulation package {\sc oopse}. It is important
18   \section{\label{oopseSec:foreword}Foreword}
19  
20   In this chapter, I present and detail the capabilities of the open
21 < source simulation package {\sc oopse}. It is important to note, that a
22 < simulation package of this size and scope would not have been possible
21 > source simulation program {\sc oopse}. It is important to note that a
22 > simulation program 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 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
27 > complete description to the program'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 the long
# Line 70 | Line 70 | code was not originally designed to simulate. Examples
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 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
77 < develop their own code or license and extend one of the commercial
78 < packages. What we have elected to do, is develop a package of
79 < simulation code capable of implementing the types of models upon which
80 < our research is based.
73 > code was not originally designed to simulate. Examples of techniques
74 > and energetics not commonly implemented include; dipole-dipole
75 > interactions, rigid body dynamics, and metallic potentials. When faced
76 > with these obstacles, a researcher must either develop their own code
77 > or license and extend one of the commercial packages. What we have
78 > elected to do is develop a body of simulation code capable of
79 > implementing the types of models upon which our research is based.
80  
81   In developing {\sc oopse}, we have adhered to the precepts of Open
82   Source development, and are releasing our source code with a
# Line 173 | Line 172 | maintained for each rigid body. At a minimum, the rota
172   each rigid body. In order to move between the space fixed and body
173   fixed coordinate axes, parameters describing the orientation must be
174   maintained for each rigid body. At a minimum, the rotation matrix
175 < (\textbf{A}) can be described by the three Euler angles ($\phi,
176 < \theta,$ and $\psi$), where the elements of \textbf{A} are composed of
175 > ($\mathsf{A}$) can be described by the three Euler angles ($\phi,
176 > \theta,$ and $\psi$), where the elements of $\mathsf{A}$ are composed of
177   trigonometric operations involving $\phi, \theta,$ and
178   $\psi$.\cite{Goldstein01} In order to avoid numerical instabilities
179   inherent in using the Euler angles, the four parameter ``quaternion''
180 < scheme is often used. The elements of \textbf{A} can be expressed as
180 > scheme is often used. The elements of $\mathsf{A}$ can be expressed as
181   arithmetic operations involving the four quaternions ($q_0, q_1, q_2,$
182   and $q_3$).\cite{allen87:csl} Use of quaternions also leads to
183   performance enhancements, particularly for very small
# Line 194 | Line 193 | molecule{
193  
194   \begin{lstlisting}[float,caption={[Defining rigid bodies]A sample definition of a rigid body},label={sch:rigidBody}]
195   molecule{
196 <  name = "TIP3P_water";
196 >  name = "TIP3P";
197 >  nAtoms = 3;
198 >  atom[0]{
199 >    type = "O_TIP3P";
200 >    position( 0.0, 0.0, -0.06556 );
201 >  }
202 >  atom[1]{
203 >    type = "H_TIP3P";
204 >    position( 0.0, 0.75695, 0.52032 );
205 >  }
206 >  atom[2]{
207 >    type = "H_TIP3P";
208 >    position( 0.0, -0.75695, 0.52032 );
209 >  }
210 >
211    nRigidBodies = 1;
212 <  rigidBody[0]{
213 <    nAtoms = 3;
214 <    atom[0]{
202 <      type = "O_TIP3P";
203 <      position( 0.0, 0.0, -0.06556 );    
204 <    }                                    
205 <    atom[1]{
206 <      type = "H_TIP3P";
207 <      position( 0.0, 0.75695, 0.52032 );
208 <    }
209 <    atom[2]{
210 <      type = "H_TIP3P";
211 <      position( 0.0, -0.75695, 0.52032 );
212 <    }
213 <    position( 0.0, 0.0, 0.0 );
214 <    orientation( 0.0, 0.0, 1.0 );
212 >  rigidBody[0]{
213 >    nMembers = 3;
214 >    members(0, 1, 2);
215    }
216   }
217   \end{lstlisting}
# Line 299 | Line 299 | point dipole interaction sites. By placing a dipole at
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 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
302 > point dipole interaction sites. By placing a dipole at the head
303 > group's center of mass, our model mimics the charge separation found
304 > in common phospholipid head groups such as
305 > phosphatidylcholine.\cite{Cevc87} Additionally, a large Lennard-Jones
306 > site is located at the pseudoatom's center of mass. The model is
307 > illustrated by the red atom in Fig.~\ref{oopseFig:lipidModel}. The
308 > water model we use to complement the dipoles of the lipids is our
309 > reparameterization of the soft sticky dipole (SSD) model of Ichiye
310   \emph{et al.}\cite{liu96:new_model}
311  
312   \begin{figure}
313   \centering
314 < \includegraphics[width=\linewidth]{lipidModel.eps}
314 > \includegraphics[width=\linewidth]{twoChainFig.eps}
315   \caption[A representation of a lipid model in {\sc duff}]{A representation of the lipid model. $\phi$ is the torsion angle, $\theta$ %
316 < is the bend angle, $\mu$ is the dipole moment of the head group, and n
317 < is the chain length.}
316 > is the bend angle, and $\mu$ is the dipole moment of the head group.}
317   \label{oopseFig:lipidModel}
318   \end{figure}
319  
# Line 338 | Line 337 | illustrated in Scheme \ref{sch:DUFF}.
337   integrating the equations of motion. A simulation using {\sc duff} is
338   illustrated in Scheme \ref{sch:DUFF}.
339  
340 < \begin{lstlisting}[float,caption={[Invocation of {\sc duff}]Sample \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}]
340 > \begin{lstlisting}[float,caption={[Invocation of {\sc duff}]A portion of a \texttt{.bass} file showing a simulation utilizing {\sc duff}},label={sch:DUFF}]
341  
342   #include "water.mdl"
343   #include "lipid.mdl"
# Line 580 | Line 579 | reference~\cite{Gezelter04}.
579   density corrected SSD models can be found in
580   reference~\cite{Gezelter04}.
581  
582 < \begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]An example file showing a simulation including {\sc ssd} water.},label={sch:ssd}]
582 > \begin{lstlisting}[float,caption={[A simulation of {\sc ssd} water]A portion of a \texttt{.bass} file showing a simulation including {\sc ssd} water.},label={sch:ssd}]
583  
584   #include "water.mdl"
585  
# Line 644 | Line 643 | These fits, are included in {\sc oopse}.
643   surrounding atom $i$ for both the density $\rho$ and pairwise $\phi$
644   interactions. Foiles \emph{et al}.~fit {\sc eam} potentials for the fcc
645   metals Cu, Ag, Au, Ni, Pd, Pt and alloys of these metals.\cite{FBD86}
646 < These fits, are included in {\sc oopse}.
646 > These fits are included in {\sc oopse}.
647  
648   \subsection{\label{oopseSec:pbc}Periodic Boundary Conditions}
649  
# Line 664 | Line 663 | Where $\mathbf{h}_j$ is the column vector of the $j$th
663   \begin{equation}
664   \mathsf{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
666 > Where $\mathbf{h}_{\alpha}$ is the column vector of the $\alpha$ 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 fluctuations when constraining
669   the pressure.
# Line 679 | Line 678 | cast each element to lie on the in the range $[-0.5,0.
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]$:
681 > cast each element to lie in the range $[-0.5,0.5]$:
682   \begin{equation}
683   s_{i}^{\prime}=s_{i}-\roundme(s_{i})
684   \end{equation}
685   Where $s_i$ is the $i$th element of $\mathbf{s}$, and
686 < $\roundme(s_i)$is given by
686 > $\roundme(s_i)$ is given by
687   \begin{equation}
688   \roundme(x) =
689          \begin{cases}
# Line 811 | Line 810 | output files.
810   entities are written out using quanternions, to save space in the
811   output files.
812  
813 < \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]
813 > \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]
814  
815   nAtoms
816   time; Hxx Hyx Hzx; Hxy Hyy Hzy; Hxz Hyz Hzz;
# Line 867 | Line 866 | statistics file is denoted with the \texttt{.stat} fil
866  
867   \section{\label{oopseSec:mechanics}Mechanics}
868  
870 –
871 – \section{\label{sec:mechanics}Mechanics}
872 –
869   \subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the
870   DLM method}
871  
# Line 879 | Line 875 | microcanonical (NVE) ensemble.\cite{}
875   (DLM).\cite{Dullweber1997} When there are no directional atoms or
876   rigid bodies present in the simulation, this integrator becomes the
877   standard velocity-Verlet integrator which is known to sample the
878 < microcanonical (NVE) ensemble.\cite{}
878 > microcanonical (NVE) ensemble.\cite{Frenkel1996}
879  
880   Previous integration methods for orientational motion have problems
881   that are avoided in the DLM method.  Direct propagation of the Euler
# Line 914 | Line 910 | where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian po
910   {\bf j}_i \right) +
911   V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right)
912   \end{equation}
913 < where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector
914 < and velocity of the center of mass of particle $i$, and ${\bf j}_i$
915 < and $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular
916 < momentum and moment of inertia tensor, respectively.  $\mathsf{A}_i$
913 > Where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector
914 > and velocity of the center of mass of particle $i$, and ${\bf j}_i$,
915 > $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular
916 > momentum and moment of inertia tensor respectively, and the
917 > superscript $T$ denotes the transpose of the vector.  $\mathsf{A}_i$
918   is the $3 \times 3$ rotation matrix describing the instantaneous
919   orientation of the particle.  $V$ is the potential energy function
920   which may depend on both the positions $\left\{{\bf r}\right\}$ and
# Line 998 | Line 995 | Monte Carlo applications, and
995   {\it symplectic}),
996   \item the integrator is time-{\it reversible}, making it suitable for Hybrid
997   Monte Carlo applications, and
998 < \item the error for a single time step is of order $O\left(h^3\right)$
998 > \item the error for a single time step is of order $\mathcal{O}\left(h^4\right)$
999   for timesteps of length $h$.
1000   \end{enumerate}
1001  
1002   The integration of the equations of motion is carried out in a
1003 < velocity-Verlet style 2-part algorithm:
1003 > velocity-Verlet style 2-part algorithm, where $h= \delta t$:
1004  
1005   {\tt moveA:}
1006 < \begin{eqnarray}
1007 < {\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
1008 < v}(t) + \frac{\delta t}{2} \left( {\bf f}(t) / m \right) \\
1009 < {\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t  {\bf
1010 < v}\left(t + \delta t / 2 \right) \\
1011 < {\bf j}\left(t + \delta t / 2 \right)  & \leftarrow & {\bf
1012 < j}(t) + \frac{\delta t}{2} {\bf \tau}^b(t)  \\
1013 < \mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left( \delta t
1014 < {\bf j}(t + \delta t / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1}
1015 < \right)
1016 < \end{eqnarray}
1006 > \begin{align*}
1007 > {\bf v}\left(t + h / 2\right)  &\leftarrow  {\bf v}(t)
1008 >        + \frac{h}{2} \left( {\bf f}(t) / m \right) \\
1009 > %
1010 > {\bf r}(t + h) &\leftarrow {\bf r}(t)
1011 >        + h  {\bf v}\left(t + h / 2 \right) \\
1012 > %
1013 > {\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
1014 >        + \frac{h}{2} {\bf \tau}^b(t)  \\
1015 > %
1016 > \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
1017 >        (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right)
1018 > \end{align*}
1019  
1020 < In this context, the $\mathrm{rot}$ function is the reversible product
1020 > In this context, the $\mathrm{rotate}$ function is the reversible product
1021   of the three body-fixed rotations,
1022   \begin{equation}
1023 < \mathrm{rot}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
1023 > \mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
1024   \mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y /
1025   2) \cdot \mathsf{G}_x(a_x /2)
1026   \end{equation}
# Line 1058 | Line 1057 | torques are calculated at the new positions and orient
1057   torques are calculated at the new positions and orientations
1058  
1059   {\tt doForces:}
1060 < \begin{eqnarray}
1061 < {\bf f}(t + \delta t) & \leftarrow & - \left(\frac{\partial V}{\partial {\bf
1062 < r}}\right)_{{\bf r}(t + \delta t)} \\
1063 < {\bf \tau}^{s}(t + \delta t) & \leftarrow & {\bf u}(t + \delta t)
1064 < \times \frac{\partial V}{\partial {\bf u}} \\
1065 < {\bf \tau}^{b}(t + \delta t) & \leftarrow & \mathsf{A}(t + \delta t)
1066 < \cdot {\bf \tau}^s(t + \delta t)
1067 < \end{eqnarray}
1060 > \begin{align*}
1061 > {\bf f}(t + h) &\leftarrow  
1062 >        - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)} \\
1063 > %
1064 > {\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h)
1065 >        \times \frac{\partial V}{\partial {\bf u}} \\
1066 > %
1067 > {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
1068 >        \cdot {\bf \tau}^s(t + h)
1069 > \end{align*}
1070  
1071   {\sc oopse} automatically updates ${\bf u}$ when the rotation matrix
1072   $\mathsf{A}$ is calculated in {\tt moveA}.  Once the forces and
# Line 1073 | Line 1074 | advanced to the same time value.
1074   advanced to the same time value.
1075  
1076   {\tt moveB:}
1077 < \begin{eqnarray}
1078 < {\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
1079 < v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
1080 < {\bf f}(t + \delta t) / m \right) \\
1081 < {\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
1082 < j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} {\bf
1083 < \tau}^b(t + \delta t)  
1083 < \end{eqnarray}
1077 > \begin{align*}
1078 > {\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2 \right)
1079 >        + \frac{h}{2} \left( {\bf f}(t + h) / m \right) \\
1080 > %
1081 > {\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t + h / 2 \right)
1082 >        + \frac{h}{2} {\bf \tau}^b(t + h)  
1083 > \end{align*}
1084  
1085   The matrix rotations used in the DLM method end up being more costly
1086   computationally than the simpler arithmetic quaternion
# Line 1088 | Line 1088 | figure \ref{timestep}.
1088   shows an average 7\% increase in computation time using the DLM method
1089   in place of quaternions. This cost is more than justified when
1090   comparing the energy conservation of the two methods as illustrated in
1091 < figure \ref{timestep}.
1091 > Fig.~\ref{timestep}.
1092  
1093   \begin{figure}
1094   \centering
# Line 1102 | Line 1102 | In figure \ref{timestep}, the resulting energy drift a
1102   \label{timestep}
1103   \end{figure}
1104  
1105 < In figure \ref{timestep}, the resulting energy drift at various time
1105 > In Fig.~\ref{timestep}, the resulting energy drift at various time
1106   steps for both the DLM and quaternion integration schemes is
1107   compared. All of the 1000 molecule water simulations started with the
1108   same configuration, and the only difference was the method for
# Line 1122 | Line 1122 | $\delta t$ & {\tt dt = 2.0;} & fs & none
1122   \begin{tabular}{llll}
1123   {\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf
1124   default value} \\  
1125 < $\delta t$ & {\tt dt = 2.0;} & fs & none
1125 > $h$ & {\tt dt = 2.0;} & fs & none
1126   \end{tabular}
1127   \end{center}
1128  
# Line 1136 | Line 1136 | NVE & microcanonical & {\tt ensemble = ``NVE''; } \\
1136   \begin{center}
1137   \begin{tabular}{lll}
1138   {\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\
1139 < NVE & microcanonical & {\tt ensemble = ``NVE''; } \\
1140 < NVT & canonical & {\tt ensemble = ``NVT''; } \\
1141 < NPTi & isobaric-isothermal (with isotropic volume changes) & {\tt
1142 < ensemble = ``NPTi'';} \\
1143 < NPTf & isobaric-isothermal (with changes to box shape) & {\tt
1144 < ensemble = ``NPTf'';} \\
1145 < NPTxyz & approximate isobaric-isothermal & {\tt ensemble =
1146 < ``NPTxyz'';} \\
1147 < &  (with separate barostats on each box dimension) &
1139 > NVE & microcanonical & {\tt ensemble = NVE; } \\
1140 > NVT & canonical & {\tt ensemble = NVT; } \\
1141 > NPTi & isobaric-isothermal & {\tt ensemble = NPTi;} \\
1142 >  &  (with isotropic volume changes) & \\
1143 > NPTf & isobaric-isothermal & {\tt ensemble = NPTf;} \\
1144 >  & (with changes to box shape) & \\
1145 > NPTxyz & approximate isobaric-isothermal & {\tt ensemble = NPTxyz;} \\
1146 > &  (with separate barostats on each box dimension) & \\
1147   \end{tabular}
1148   \end{center}
1149  
1150 < The relatively well-known Nos\'e-Hoover thermostat is implemented in
1151 < {\sc oopse}'s NVT integrator.  This method couples an extra degree of
1152 < freedom (the thermostat) to the kinetic energy of the system, and has
1153 < been shown to sample the canonical distribution in the system degrees
1154 < of freedom while conserving a quantity that is, to within a constant,
1155 < the Helmholtz free energy.
1150 > The relatively well-known Nos\'e-Hoover thermostat\cite{Hoover85} is
1151 > implemented in {\sc oopse}'s NVT integrator.  This method couples an
1152 > extra degree of freedom (the thermostat) to the kinetic energy of the
1153 > system, and has been shown to sample the canonical distribution in the
1154 > system degrees of freedom while conserving a quantity that is, to
1155 > within a constant, the Helmholtz free energy.\cite{melchionna93}
1156  
1157   NPT algorithms attempt to maintain constant pressure in the system by
1158   coupling the volume of the system to a barostat.  {\sc oopse} contains
1159   three different constant pressure algorithms.  The first two, NPTi and
1160   NPTf have been shown to conserve a quantity that is, to within a
1161 < constant, the Gibbs free energy.  The Melchionna modification to the
1162 < Hoover barostat is implemented in both NPTi and NPTf.  NPTi allows
1163 < only isotropic changes in the simulation box, while box {\it shape}
1164 < variations are allowed in NPTf.  The NPTxyz integrator has {\it not}
1165 < been shown to sample from the isobaric-isothermal ensemble.  It is
1166 < useful, however, in that it maintains orthogonality for the axes of
1167 < the simulation box while attempting to equalize pressure along the
1168 < three perpendicular directions in the box.
1161 > constant, the Gibbs free energy.\cite{melchionna93} The Melchionna
1162 > modification to the Hoover barostat is implemented in both NPTi and
1163 > NPTf.  NPTi allows only isotropic changes in the simulation box, while
1164 > box {\it shape} variations are allowed in NPTf.  The NPTxyz integrator
1165 > has {\it not} been shown to sample from the isobaric-isothermal
1166 > ensemble.  It is useful, however, in that it maintains orthogonality
1167 > for the axes of the simulation box while attempting to equalize
1168 > pressure along the three perpendicular directions in the box.
1169  
1170   Each of the extended system integrators requires additional keywords
1171   to set target values for the thermodynamic state variables that are
# Line 1174 | Line 1173 | variables.
1173   characteristic decay times for the dynamics of the extended
1174   variables.
1175  
1176 + \begin{center}
1177   \begin{tabular}{llll}
1178   {\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf
1179   default value} \\  
# Line 1182 | Line 1182 | $\tau_B$ & {\tt tauBarostat = 5e3;} & fs  & none \\
1182   $\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\
1183   $\tau_B$ & {\tt tauBarostat = 5e3;} & fs  & none \\
1184           & {\tt resetTime = 200;} & fs & none \\
1185 <         & {\tt useInitialExtendedSystemState = ``true'';} & logical &
1186 < false
1185 >         & {\tt useInitialExtendedSystemState = true;} & logical &
1186 > true
1187   \end{tabular}
1188 + \end{center}
1189  
1190   Two additional keywords can be used to either clear the extended
1191   system variables periodically ({\tt resetTime}), or to maintain the
# Line 1192 | Line 1193 | and their use in the integrators follows below.
1193   useInitialExtendedSystemState}).  More details on these variables
1194   and their use in the integrators follows below.
1195  
1196 < \subsubsection{\label{oopseSec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting}
1196 > \subsection{\label{oopseSec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting}
1197  
1198   The Nos\'e-Hoover equations of motion are given by\cite{Hoover85}
1199   \begin{eqnarray}
# Line 1239 | Line 1240 | part algorithm:
1240   part algorithm:
1241  
1242   {\tt moveA:}
1243 < \begin{eqnarray}
1244 < T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
1245 < {\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
1246 < v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
1247 < \chi(t)\right) \\
1248 < {\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf
1249 < v}\left(t + \delta t / 2 \right) \\
1250 < {\bf j}\left(t + \delta t / 2 \right)  & \leftarrow & {\bf
1251 < j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
1252 < \chi(t) \right) \\
1253 < \mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t *
1254 < {\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
1255 < \chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) +
1256 < \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
1257 < \right)
1258 < \end{eqnarray}
1243 > \begin{align*}
1244 > T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
1245 > %
1246 > {\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t)
1247 >        + \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
1248 >        \chi(t)\right) \\
1249 > %
1250 > {\bf r}(t + h) &\leftarrow {\bf r}(t)
1251 >        + h {\bf v}\left(t + h / 2 \right) \\
1252 > %
1253 > {\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
1254 >        + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
1255 >        \chi(t) \right) \\
1256 > %
1257 > \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}
1258 >        \left(h * {\bf j}(t + h / 2)
1259 >        \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
1260 > %
1261 > \chi\left(t + h / 2 \right) &\leftarrow \chi(t)
1262 >        + \frac{h}{2 \tau_T^2} \left( \frac{T(t)}
1263 >        {T_{\mathrm{target}}} - 1 \right)
1264 > \end{align*}
1265  
1266 < Here $\mathrm{rot}(\delta t * {\bf j}
1266 > Here $\mathrm{rotate}(h * {\bf j}
1267   \overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic Trotter
1268   factorization of the three rotation operations that was discussed in
1269   the section on the DLM integrator.  Note that this operation modifies
# Line 1272 | Line 1279 | advanced to the same time value.
1279   advanced to the same time value.
1280  
1281   {\tt moveB:}
1282 < \begin{eqnarray}
1283 < T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\},
1284 < \left\{{\bf j}(t + \delta t)\right\} \\
1285 < \chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t /
1286 < 2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta
1287 < t)}{T_{\mathrm{target}}} - 1 \right) \\
1288 < {\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
1289 < v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
1290 < \frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t)
1291 < \chi(t \delta t)\right) \\
1292 < {\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
1293 < j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf
1294 < \tau}^b(t + \delta t) - {\bf j}(t + \delta t)
1295 < \chi(t + \delta t) \right)
1296 < \end{eqnarray}
1282 > \begin{align*}
1283 > T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
1284 >        \left\{{\bf j}(t + h)\right\} \\
1285 > %
1286 > \chi\left(t + h \right) &\leftarrow \chi\left(t + h /
1287 >        2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)}
1288 >        {T_{\mathrm{target}}} - 1 \right) \\
1289 > %
1290 > {\bf v}\left(t + h \right)  &\leftarrow {\bf v}\left(t
1291 >        + h / 2 \right) + \frac{h}{2} \left(
1292 >        \frac{{\bf f}(t + h)}{m} - {\bf v}(t + h)
1293 >        \chi(t h)\right) \\
1294 > %
1295 > {\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t
1296 >        + h / 2 \right) + \frac{h}{2}
1297 >        \left( {\bf \tau}^b(t + h) - {\bf j}(t + h)
1298 >        \chi(t + h) \right)
1299 > \end{align*}
1300  
1301 < Since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ are required
1302 < to caclculate $T(t + \delta t)$ as well as $\chi(t + \delta t)$, they
1303 < indirectly depend on their own values at time $t + \delta t$.  {\tt
1304 < moveB} is therefore done in an iterative fashion until $\chi(t +
1305 < \delta t)$ becomes self-consistent.  The relative tolerance for the
1306 < self-consistency check defaults to a value of $\mbox{10}^{-6}$, but
1307 < {\sc oopse} will terminate the iteration after 4 loops even if the
1298 < consistency check has not been satisfied.
1301 > Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to caclculate
1302 > $T(t + h)$ as well as $\chi(t + h)$, they indirectly depend on their
1303 > own values at time $t + h$.  {\tt moveB} is therefore done in an
1304 > iterative fashion until $\chi(t + h)$ becomes self-consistent.  The
1305 > relative tolerance for the self-consistency check defaults to a value
1306 > of $\mbox{10}^{-6}$, but {\sc oopse} will terminate the iteration
1307 > after 4 loops even if the consistency check has not been satisfied.
1308  
1309   The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the
1310   extended system that is, to within a constant, identical to the
1311 < Helmholtz free energy,
1311 > Helmholtz free energy,\cite{melchionna93}
1312   \begin{equation}
1313   H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left(
1314   \frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime
1315   \right)
1316   \end{equation}
1317 < Poor choices of $\delta t$ or $\tau_T$ can result in non-conservation
1317 > Poor choices of $h$ or $\tau_T$ can result in non-conservation
1318   of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the
1319   last column of the {\tt .stat} file to allow checks on the quality of
1320   the integration.
# Line 1314 | Line 1323 | algorithms are given in section \ref{oopseSec:rattle}.
1323   {\tt moveB} portions of the algorithm.  Details on the constraint
1324   algorithms are given in section \ref{oopseSec:rattle}.
1325  
1326 < \subsubsection{\label{sec:NPTi}Constant-pressure integration with
1326 > \subsection{\label{sec:NPTi}Constant-pressure integration with
1327   isotropic box deformations (NPTi)}
1328  
1329   To carry out isobaric-isothermal ensemble calculations {\sc oopse}
# Line 1382 | Line 1391 | velocity-Verlet style 2 part algorithm:
1391   velocity-Verlet style 2 part algorithm:
1392  
1393   {\tt moveA:}
1394 < \begin{eqnarray}
1395 < T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
1396 < P(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(t)\right\} \\
1397 < {\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
1398 < v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
1399 < \left(\chi(t) + \eta(t) \right) \right) \\
1400 < {\bf j}\left(t + \delta t / 2 \right)  & \leftarrow & {\bf
1401 < j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
1402 < \chi(t) \right) \\
1403 < \mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t *
1404 < {\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
1405 < \chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) +
1406 < \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
1407 < \right) \\
1408 < \eta(t + \delta t / 2) & \leftarrow & \eta(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B
1409 < T(t) \tau_B^2} \left( P(t) - P_{\mathrm{target}} \right) \\
1410 < {\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf
1411 < v}\left(t + \delta t / 2 \right) + \eta(t + \delta t / 2)\left[ {\bf
1412 < r}(t + \delta t) - {\bf R}_0 \right] \right\} \\
1413 < \mathsf{H}(t + \delta t) & \leftarrow & e^{-\delta t \eta(t + \delta t
1414 < / 2)} \mathsf{H}(t)
1415 < \end{eqnarray}
1394 > \begin{align*}
1395 > T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
1396 > %
1397 > P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\} \\
1398 > %
1399 > {\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t)
1400 >        + \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
1401 >        \left(\chi(t) + \eta(t) \right) \right) \\
1402 > %
1403 > {\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
1404 >        + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
1405 >        \chi(t) \right) \\
1406 > %
1407 > \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
1408 >        {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1}
1409 >        \right) \\
1410 > %
1411 > \chi\left(t + h / 2 \right) &\leftarrow \chi(t) +
1412 >        \frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
1413 >        \right) \\
1414 > %
1415 > \eta(t + h / 2) &\leftarrow \eta(t) + \frac{h
1416 >        \mathcal{V}(t)}{2 N k_B T(t) \tau_B^2} \left( P(t)
1417 >        - P_{\mathrm{target}} \right) \\
1418 > %
1419 > {\bf r}(t + h) &\leftarrow {\bf r}(t) + h
1420 >        \left\{ {\bf v}\left(t + h / 2 \right)
1421 >        + \eta(t + h / 2)\left[ {\bf r}(t + h)
1422 >        - {\bf R}_0 \right] \right\} \\
1423 > %
1424 > \mathsf{H}(t + h) &\leftarrow e^{-h \eta(t + h / 2)}
1425 >        \mathsf{H}(t)
1426 > \end{align*}
1427  
1428   Most of these equations are identical to their counterparts in the NVT
1429 < integrator, but the propagation of positions to time $t + \delta t$
1429 > integrator, but the propagation of positions to time $t + h$
1430   depends on the positions at the same time.  {\sc oopse} carries out
1431   this step iteratively (with a limit of 5 passes through the iterative
1432   loop).  Also, the simulation box $\mathsf{H}$ is scaled uniformly for
1433   one full time step by an exponential factor that depends on the value
1434   of $\eta$ at time $t +
1435 < \delta t / 2$.  Reshaping the box uniformly also scales the volume of
1435 > h / 2$.  Reshaping the box uniformly also scales the volume of
1436   the box by
1437   \begin{equation}
1438 < \mathcal{V}(t + \delta t) \leftarrow e^{ - 3 \delta t \eta(t + \delta t /2)}
1438 > \mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /2)}
1439   \mathcal{V}(t)
1440   \end{equation}
1441  
# Line 1425 | Line 1445 | the same time value.
1445   the same time value.
1446  
1447   {\tt moveB:}
1448 < \begin{eqnarray}
1449 < T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\},
1450 < \left\{{\bf j}(t + \delta t)\right\} \\
1451 < P(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\},
1452 < \left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(t + \delta t)\right\} \\
1453 < \chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t /
1454 < 2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta
1455 < t)}{T_{\mathrm{target}}} - 1 \right) \\
1456 < \eta(t + \delta t) & \leftarrow & \eta(t + \delta t / 2) +
1457 < \frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2}
1458 < \left( P(t + \delta t) - P_{\mathrm{target}}
1459 < \right) \\
1460 < {\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
1461 < v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
1462 < \frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t)
1463 < (\chi(t + \delta t) + \eta(t + \delta t)) \right) \\
1464 < {\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
1465 < j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf
1466 < \tau}^b(t + \delta t) - {\bf j}(t + \delta t)
1467 < \chi(t + \delta t) \right)
1468 < \end{eqnarray}
1448 > \begin{align*}
1449 > T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
1450 >        \left\{{\bf j}(t + h)\right\} \\
1451 > %
1452 > P(t + h) &\leftarrow  \left\{{\bf r}(t + h)\right\},
1453 >        \left\{{\bf v}(t + h)\right\} \\
1454 > %
1455 > \chi\left(t + h \right) &\leftarrow \chi\left(t + h /
1456 >        2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)}
1457 >        {T_{\mathrm{target}}} - 1 \right) \\
1458 > %
1459 > \eta(t + h) &\leftarrow \eta(t + h / 2) +
1460 >        \frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h)
1461 >        \tau_B^2} \left( P(t + h) - P_{\mathrm{target}} \right) \\
1462 > %
1463 > {\bf v}\left(t + h \right)  &\leftarrow {\bf v}\left(t
1464 >        + h / 2 \right) + \frac{h}{2} \left(
1465 >        \frac{{\bf f}(t + h)}{m} - {\bf v}(t + h)
1466 >        (\chi(t + h) + \eta(t + h)) \right) \\
1467 > %
1468 > {\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t
1469 >        + h / 2 \right) + \frac{h}{2} \left( {\bf
1470 >        \tau}^b(t + h) - {\bf j}(t + h)
1471 >        \chi(t + h) \right)
1472 > \end{align*}
1473  
1474 < Once again, since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$
1475 < are required to caclculate $T(t + \delta t)$, $P(t + \delta t)$, $\chi(t +
1476 < \delta t)$, and $\eta(t + \delta t)$, they indirectly depend on their
1477 < own values at time $t + \delta t$.  {\tt moveB} is therefore done in
1478 < an iterative fashion until $\chi(t + \delta t)$ and $\eta(t + \delta
1479 < t)$ become self-consistent.  The relative tolerance for the
1480 < self-consistency check defaults to a value of $\mbox{10}^{-6}$, but
1457 < {\sc oopse} will terminate the iteration after 4 loops even if the
1474 > Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required
1475 > to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t +
1476 > h)$, they indirectly depend on their own values at time $t + h$.  {\tt
1477 > moveB} is therefore done in an iterative fashion until $\chi(t + h)$
1478 > and $\eta(t + h)$ become self-consistent.  The relative tolerance for
1479 > the self-consistency check defaults to a value of $\mbox{10}^{-6}$,
1480 > but {\sc oopse} will terminate the iteration after 4 loops even if the
1481   consistency check has not been satisfied.
1482  
1483   The Melchionna modification of the Nos\'e-Hoover-Andersen algorithm is
# Line 1481 | Line 1504 | algorithms are given in section \ref{oopseSec:rattle}.
1504   {\tt moveB} portions of the algorithm.  Details on the constraint
1505   algorithms are given in section \ref{oopseSec:rattle}.
1506  
1507 < \subsubsection{\label{sec:NPTf}Constant-pressure integration with a
1507 > \subsection{\label{sec:NPTf}Constant-pressure integration with a
1508   flexible box (NPTf)}
1509  
1510   There is a relatively simple generalization of the
# Line 1493 | Line 1516 | the box shape.  The equations of motion for this metho
1516   \begin{eqnarray}
1517   \dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right) \\
1518   \dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} +
1519 < \chi \mathsf{1}) {\bf v} \\
1519 > \chi \cdot \mathsf{1}) {\bf v} \\
1520   \dot{\mathsf{A}} & = & \mathsf{A} \cdot
1521   \mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\
1522   \dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1}
# Line 1501 | Line 1524 | V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\
1524   V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\
1525   \dot{\chi} & = & \frac{1}{\tau_{T}^2} \left(
1526   \frac{T}{T_{\mathrm{target}}} - 1 \right) \\
1527 < \dot{\overleftrightarrow{eta}} & = & \frac{1}{\tau_{B}^2 f k_B
1527 > \dot{\overleftrightarrow{\eta}} & = & \frac{1}{\tau_{B}^2 f k_B
1528   T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) \\
1529   \dot{\mathsf{H}} & = &  \overleftrightarrow{\eta} \cdot \mathsf{H}
1530   \label{eq:melchionna2}
# Line 1515 | Line 1538 | NPTi integration:
1538   NPTi integration:
1539  
1540   {\tt moveA:}
1541 < \begin{eqnarray}
1542 < T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
1543 < \overleftrightarrow{\mathsf{P}}(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(t)\right\} \\
1544 < {\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
1545 < v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} -
1546 < \left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot
1547 < {\bf v}(t) \right) \\
1548 < {\bf j}\left(t + \delta t / 2 \right)  & \leftarrow & {\bf
1549 < j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
1550 < \chi(t) \right) \\
1551 < \mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t *
1552 < {\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
1553 < \chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) +
1554 < \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
1555 < \right) \\
1556 < \overleftrightarrow{\eta}(t + \delta t / 2) & \leftarrow & \overleftrightarrow{\eta}(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B
1557 < T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t) - P_{\mathrm{target}}\mathsf{1} \right) \\
1558 < {\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf
1559 < v}\left(t + \delta t / 2 \right) + \overleftrightarrow{\eta}(t +
1560 < \delta t / 2) \cdot \left[ {\bf
1561 < r}(t + \delta t) - {\bf R}_0 \right] \right\} \\
1562 < \mathsf{H}(t + \delta t) & \leftarrow & \mathsf{H}(t) \cdot e^{-\delta t
1563 < \overleftrightarrow{\eta}(t + \delta t / 2)}
1564 < \end{eqnarray}
1541 > \begin{align*}
1542 > T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
1543 > %
1544 > \overleftrightarrow{\mathsf{P}}(t) &\leftarrow \left\{{\bf r}(t)\right\},
1545 >        \left\{{\bf v}(t)\right\} \\
1546 > %
1547 > {\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t)
1548 >        + \frac{h}{2} \left( \frac{{\bf f}(t)}{m} -
1549 >        \left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot
1550 >        {\bf v}(t) \right) \\
1551 > %
1552 > {\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
1553 >        + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
1554 >        \chi(t) \right) \\
1555 > %
1556 > \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
1557 >        {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1}
1558 >        \right) \\
1559 > %
1560 > \chi\left(t + h / 2 \right) &\leftarrow \chi(t) +
1561 >        \frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}}
1562 >        - 1 \right) \\
1563 > %
1564 > \overleftrightarrow{\eta}(t + h / 2) &\leftarrow
1565 >        \overleftrightarrow{\eta}(t) + \frac{h \mathcal{V}(t)}{2 N k_B
1566 >        T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t)
1567 >        - P_{\mathrm{target}}\mathsf{1} \right) \\
1568 > %
1569 > {\bf r}(t + h) &\leftarrow {\bf r}(t) + h \left\{ {\bf v}
1570 >        \left(t + h / 2 \right) + \overleftrightarrow{\eta}(t +
1571 >        h / 2) \cdot \left[ {\bf r}(t + h)
1572 >        - {\bf R}_0 \right] \right\} \\
1573 > %
1574 > \mathsf{H}(t + h) &\leftarrow \mathsf{H}(t) \cdot e^{-h
1575 >        \overleftrightarrow{\eta}(t + h / 2)}
1576 > \end{align*}
1577   {\sc oopse} uses a power series expansion truncated at second order
1578   for the exponential operation which scales the simulation box.
1579  
# Line 1546 | Line 1581 | NPTi integrator:
1581   NPTi integrator:
1582  
1583   {\tt moveB:}
1584 < \begin{eqnarray}
1585 < T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\},
1586 < \left\{{\bf j}(t + \delta t)\right\} \\
1587 < \overleftrightarrow{\mathsf{P}}(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\},
1588 < \left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(t + \delta t)\right\} \\
1589 < \chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t /
1590 < 2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta
1591 < t)}{T_{\mathrm{target}}} - 1 \right) \\
1592 < \overleftrightarrow{\eta}(t + \delta t) & \leftarrow & \overleftrightarrow{\eta}(t + \delta t / 2) +
1593 < \frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2}
1594 < \left( \overleftrightarrow{P}(t + \delta t) - P_{\mathrm{target}}\mathsf{1}
1595 < \right) \\
1596 < {\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
1597 < v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
1598 < \frac{{\bf f}(t + \delta t)}{m} -
1599 < (\chi(t + \delta t)\mathsf{1} + \overleftrightarrow{\eta}(t + \delta
1600 < t)) \right) \cdot {\bf v}(t + \delta t) \\
1601 < {\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
1602 < j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf
1603 < \tau}^b(t + \delta t) - {\bf j}(t + \delta t)
1604 < \chi(t + \delta t) \right)
1605 < \end{eqnarray}
1584 > \begin{align*}
1585 > T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
1586 >        \left\{{\bf j}(t + h)\right\} \\
1587 > %
1588 > \overleftrightarrow{\mathsf{P}}(t + h) &\leftarrow \left\{{\bf r}
1589 >        (t + h)\right\}, \left\{{\bf v}(t
1590 >        + h)\right\}, \left\{{\bf f}(t + h)\right\} \\
1591 > %
1592 > \chi\left(t + h \right) &\leftarrow \chi\left(t + h /
1593 >        2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+
1594 >        h)}{T_{\mathrm{target}}} - 1 \right) \\
1595 > %
1596 > \overleftrightarrow{\eta}(t + h) &\leftarrow
1597 >        \overleftrightarrow{\eta}(t + h / 2) +
1598 >        \frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h)
1599 >        \tau_B^2} \left( \overleftrightarrow{P}(t + h)
1600 >        - P_{\mathrm{target}}\mathsf{1} \right) \\
1601 > %
1602 > {\bf v}\left(t + h \right)  &\leftarrow {\bf v}\left(t
1603 >        + h / 2 \right) + \frac{h}{2} \left(
1604 >        \frac{{\bf f}(t + h)}{m} -
1605 >        (\chi(t + h)\mathsf{1} + \overleftrightarrow{\eta}(t
1606 >        + h)) \right) \cdot {\bf v}(t + h) \\
1607 > %
1608 > {\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t
1609 >        + h / 2 \right) + \frac{h}{2} \left( {\bf \tau}^b(t
1610 >        + h) - {\bf j}(t + h) \chi(t + h) \right)
1611 > \end{align*}
1612  
1613   The iterative schemes for both {\tt moveA} and {\tt moveB} are
1614   identical to those described for the NPTi integrator.
# Line 1585 | Line 1626 | electrostatic or Lennard-Jones cutoff radii.  It finds
1626   simulations.  Liquids have very small restoring forces in the
1627   off-diagonal directions, and the simulation box can very quickly form
1628   elongated and sheared geometries which become smaller than the
1629 < electrostatic or Lennard-Jones cutoff radii.  It finds most use in
1630 < simulating crystals or liquid crystals which assume non-orthorhombic
1631 < geometries.
1629 > electrostatic or Lennard-Jones cutoff radii.  The NPTf integrator
1630 > finds most use in simulating crystals or liquid crystals which assume
1631 > non-orthorhombic geometries.
1632  
1633 < \subsubsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)}
1633 > \subsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)}
1634  
1635   There is one additional extended system integrator which is somewhat
1636   simpler than the NPTf method described above.  In this case, the three
# Line 1911 | Line 1952 | Where $\beta$ are all of the unconstrained molecules i
1952   F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha} G_{\alpha}}
1953          {\sum_{\beta}\sum_i m_{\beta i}}
1954   \end{equation}
1955 < Where $\beta$ are all of the unconstrained molecules in the system.
1955 > Where $\beta$ are all of the unconstrained molecules in the
1956 > system. Similarly, the velocities of the unconstrained molecules must
1957 > also be scaled.
1958 > \begin{equation}
1959 > v_{\beta i} = v_{\beta i} + \sum_{\alpha}
1960 >        \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}
1961 > \end{equation}
1962  
1963   At the very beginning of the simulation, the molecules may not be at their
1964   constrained positions. To move a z-constrained molecule to its specified
# Line 1939 | Line 1986 | can be found in Table~\ref{oopseTb:gofrs}
1986   can be found in Table~\ref{oopseTb:gofrs}
1987  
1988   \begin{table}
1989 < \caption[The list of pair correlations in \texttt{staticProps}]{The different pair correlations in \texttt{staticProps} along with whether atom A or B must be directional.}
1989 > \caption[The list of pair correlations in \texttt{staticProps}]{THE DIFFERENT PAIR CORRELATIONS IN \texttt{staticProps}}
1990   \label{oopseTb:gofrs}
1991   \begin{center}
1992   \begin{tabular}{|l|c|c|}
# Line 1952 | Line 1999 | $\langle \cos \omega \rangle_{\text{AB}}(r)$ & Eq.~\re
1999   $\langle \cos \omega \rangle_{\text{AB}}(r)$ & Eq.~\ref{eq:cosOmegaOfR} &%
2000          both \\ \hline
2001   \end{tabular}
2002 + \begin{minipage}{\linewidth}
2003 + \centering
2004 + \vspace{2mm}
2005 + The third column specifies which atom, if any, need be a directional entity.
2006 + \end{minipage}
2007   \end{center}
2008   \end{table}
2009  
# Line 2120 | Line 2172 | $O(N/P)$ for atomic decomposition. Unfortunately all p
2172   Algorithmically simplest of the three methods is atomic decomposition
2173   where N particles in a simulation are split among P processors for the
2174   duration of the simulation. Computational cost scales as an optimal
2175 < $O(N/P)$ for atomic decomposition. Unfortunately all processors must
2176 < communicate positions and forces with all other processors at every
2177 < force evaluation, leading communication costs to scale as an
2178 < unfavorable $O(N)$, \emph{independent of the number of processors}. This
2179 < communication bottleneck led to the development of spatial and force
2180 < decomposition methods in which communication among processors scales
2181 < much more favorably. Spatial or domain decomposition divides the
2182 < physical spatial domain into 3D boxes in which each processor is
2183 < responsible for calculation of forces and positions of particles
2184 < located in its box. Particles are reassigned to different processors
2185 < as they move through simulation space. To calculate forces on a given
2186 < particle, a processor must know the positions of particles within some
2187 < cutoff radius located on nearby processors instead of the positions of
2188 < particles on all processors. Both communication between processors and
2189 < computation scale as $O(N/P)$ in the spatial method. However, spatial
2175 > $\mathcal{O}(N/P)$ for atomic decomposition. Unfortunately all
2176 > processors must communicate positions and forces with all other
2177 > processors at every force evaluation, leading communication costs to
2178 > scale as an unfavorable $\mathcal{O}(N)$, \emph{independent of the
2179 > number of processors}. This communication bottleneck led to the
2180 > development of spatial and force decomposition methods in which
2181 > communication among processors scales much more favorably. Spatial or
2182 > domain decomposition divides the physical spatial domain into 3D boxes
2183 > in which each processor is responsible for calculation of forces and
2184 > positions of particles located in its box. Particles are reassigned to
2185 > different processors as they move through simulation space. To
2186 > calculate forces on a given particle, a processor must know the
2187 > positions of particles within some cutoff radius located on nearby
2188 > processors instead of the positions of particles on all
2189 > processors. Both communication between processors and computation
2190 > scale as $\mathcal{O}(N/P)$ in the spatial method. However, spatial
2191   decomposition adds algorithmic complexity to the simulation code and
2192   is not very efficient for small N since the overall communication
2193 < scales as the surface to volume ratio $(N/P)^{2/3}$ in three
2194 < dimensions.
2193 > scales as the surface to volume ratio $\mathcal{O}(N/P)^{2/3}$ in
2194 > three dimensions.
2195  
2196   The parallelization method used in {\sc oopse} is the force
2197   decomposition method.  Force decomposition assigns particles to
# Line 2147 | Line 2200 | spatial method but still scales computationally as $O(
2200   and column processor groups. Forces are calculated on particles in a
2201   given row by particles located in that processors column
2202   assignment. Force decomposition is less complex to implement than the
2203 < spatial method but still scales computationally as $O(N/P)$ and scales
2204 < as $O(N/\sqrt{P})$ in communication cost. Plimpton has also found that
2205 < force decompositions scale more favorably than spatial decompositions
2206 < for systems up to 10,000 atoms and favorably compete with spatial
2207 < methods up to 100,000 atoms.\cite{plimpton95}
2203 > spatial method but still scales computationally as $\mathcal{O}(N/P)$
2204 > and scales as $\mathcal{O}(N/\sqrt{P})$ in communication
2205 > cost. Plimpton has also found that force decompositions scale more
2206 > favorably than spatial decompositions for systems up to 10,000 atoms
2207 > and favorably compete with spatial methods up to 100,000
2208 > atoms.\cite{plimpton95}
2209  
2210   \subsection{\label{oopseSec:memAlloc}Memory Issues in Trajectory Analysis}
2211  
# Line 2209 | Line 2263 | well.Documentation and source code for {\sc oopse} can
2263   These features are all brought together in a single open-source
2264   program. Allowing researchers to not only benefit from
2265   {\sc oopse}, but also contribute to {\sc oopse}'s development as
2266 < well.Documentation and source code for {\sc oopse} can be downloaded
2213 < from \texttt{http://www.openscience.org/oopse/}.
2266 > well.
2267  

Diff Legend

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