--- trunk/mattDisertation/oopse.tex	2004/03/05 17:21:20	1086
+++ trunk/mattDisertation/oopse.tex	2004/03/05 22:16:34	1087
@@ -867,106 +867,749 @@ statistics file is denoted with the \texttt{.stat} fil
 
 \section{\label{oopseSec:mechanics}Mechanics}
 
-\subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the Symplectic Step Integrator}
 
-Integration of the equations of motion was carried out using the
-symplectic splitting method proposed by Dullweber \emph{et
-al.}.\cite{Dullweber1997} The reason for the selection of this
-integrator, is the poor energy conservation of rigid body systems
-using quaternion dynamics. While quaternions work well for
-orientational motion in alternate ensembles, the microcanonical
-ensemble has a constant energy requirement that is quite sensitive to
-errors in the equations of motion. The original implementation of {\sc
-oopse} utilized quaternions for rotational motion propagation;
-however, a detailed investigation showed that they resulted in a
-steady drift in the total energy, something that has been observed by
-others.\cite{Laird97}
+\section{\label{sec:mechanics}Mechanics}
+
+\subsection{\label{oopseSec:integrate}Integrating the Equations of Motion: the
+DLM method}
+
+The default method for integrating the equations of motion in {\sc
+oopse} is a velocity-Verlet version of the symplectic splitting method
+proposed by Dullweber, Leimkuhler and McLachlan
+(DLM).\cite{Dullweber1997} When there are no directional atoms or
+rigid bodies present in the simulation, this integrator becomes the
+standard velocity-Verlet integrator which is known to sample the
+microcanonical (NVE) ensemble.\cite{}
+
+Previous integration methods for orientational motion have problems
+that are avoided in the DLM method.  Direct propagation of the Euler
+angles has a known $1/\sin\theta$ divergence in the equations of
+motion for $\phi$ and $\psi$,\cite{allen87:csl} leading to
+numerical instabilities any time one of the directional atoms or rigid
+bodies has an orientation near $\theta=0$ or $\theta=\pi$.  More
+modern quaternion-based integration methods have relatively poor
+energy conservation.  While quaternions work well for orientational
+motion in other ensembles, the microcanonical ensemble has a
+constant energy requirement that is quite sensitive to errors in the
+equations of motion.  An earlier implementation of {\sc oopse}
+utilized quaternions for propagation of rotational motion; however, a
+detailed investigation showed that they resulted in a steady drift in
+the total energy, something that has been observed by
+Laird {\it et al.}\cite{Laird97}      
 
 The key difference in the integration method proposed by Dullweber
-\emph{et al}.~({\sc dlm}) is that the entire rotation matrix is propagated from
-one time step to the next. In the past, this would not have been a
-feasible option, since the rotation matrix for a single body is nine
-elements long as opposed to three or four elements for Euler angles
-and quaternions respectively. System memory has become much less of an
-issue in recent times, and the {\sc dlm} method has used memory in
-exchange for substantial benefits in energy conservation.
+\emph{et al.} is that the entire $3 \times 3$ rotation matrix is
+propagated from one time step to the next. In the past, this would not
+have been feasible, since the rotation matrix for a single body has
+nine elements compared with the more memory-efficient methods (using
+three Euler angles or 4 quaternions).  Computer memory has become much
+less costly in recent years, and this can be translated into
+substantial benefits in energy conservation.
 
-The {\sc dlm} method allows for Verlet style integration of both
-linear and angular motion of rigid bodies. In the integration method,
-the orientational propagation involves a sequence of matrix
-evaluations to update the rotation matrix.\cite{Dullweber1997} These
-matrix rotations are more costly computationally than the simpler
-arithmetic quaternion propagation. With the same time step, a 1000 SSD
-particle simulation shows an average 7\% increase in computation time
-using the {\sc dlm} method in place of quaternions. This cost is more
-than justified when comparing the energy conservation of the two
-methods as illustrated in Fig.~\ref{timestep}.
+The basic equations of motion being integrated are derived from the
+Hamiltonian for conservative systems containing rigid bodies,
+\begin{equation}
+H = \sum_{i} \left( \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i +
+\frac{1}{2} {\bf j}_i^T \cdot \overleftrightarrow{\mathsf{I}}_i^{-1} \cdot
+{\bf j}_i \right) +
+V\left(\left\{{\bf r}\right\}, \left\{\mathsf{A}\right\}\right)
+\end{equation}
+where ${\bf r}_i$ and ${\bf v}_i$ are the cartesian position vector
+and velocity of the center of mass of particle $i$, and ${\bf j}_i$
+and $\overleftrightarrow{\mathsf{I}}_i$ are the body-fixed angular
+momentum and moment of inertia tensor, respectively.  $\mathsf{A}_i$
+is the $3 \times 3$ rotation matrix describing the instantaneous
+orientation of the particle.  $V$ is the potential energy function
+which may depend on both the positions $\left\{{\bf r}\right\}$ and
+orientations $\left\{\mathsf{A}\right\}$ of all particles.  The
+equations of motion for the particle centers of mass are derived from
+Hamilton's equations and are quite simple,
+\begin{eqnarray}
+\dot{{\bf r}} & = & {\bf v} \\
+\dot{{\bf v}} & = & \frac{{\bf f}}{m}
+\end{eqnarray}
+where ${\bf f}$ is the instantaneous force on the center of mass
+of the particle,
+\begin{equation}
+{\bf f} = - \frac{\partial}{\partial
+{\bf r}} V(\left\{{\bf r}(t)\right\}, \left\{\mathsf{A}(t)\right\}).
+\end{equation}
+
+The equations of motion for the orientational degrees of freedom are
+\begin{eqnarray}
+\dot{\mathsf{A}} & = & \mathsf{A} \cdot
+\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right) \\
+\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1}
+\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial
+V}{\partial \mathsf{A}} \right)
+\end{eqnarray}
+In these equations of motion, the $\mbox{skew}$ matrix of a vector
+${\bf v} = \left( v_1, v_2, v_3 \right)$ is defined:
+\begin{equation}
+\mbox{skew}\left( {\bf v} \right) := \left( 
+\begin{array}{ccc}
+0 & v_3 & - v_2 \\
+-v_3 & 0 & v_1 \\
+v_2 & -v_1 & 0 
+\end{array}
+\right)
+\end{equation}
+The $\mbox{rot}$ notation refers to the mapping of the $3 \times 3$
+rotation matrix to a vector of orientations by first computing the
+skew-symmetric part $\left(\mathsf{A} - \mathsf{A}^{T}\right)$ and
+then associating this with a length 3 vector by inverting the
+$\mbox{skew}$ function above:
+\begin{equation}
+\mbox{rot}\left(\mathsf{A}\right) := \mbox{ skew}^{-1}\left(\mathsf{A}
+- \mathsf{A}^{T} \right)
+\end{equation}
+Written this way, the $\mbox{rot}$ operation creates a set of
+conjugate angle coordinates to the body-fixed angular momenta
+represented by ${\bf j}$.  This equation of motion for angular momenta
+is equivalent to the more familiar body-fixed forms,
+\begin{eqnarray}
+\dot{j_{x}} & = & \tau^b_x(t)  +
+\left(\overleftrightarrow{\mathsf{I}}_{yy} - \overleftrightarrow{\mathsf{I}}_{zz} \right) j_y j_z \\
+\dot{j_{y}} & = & \tau^b_y(t) +
+\left(\overleftrightarrow{\mathsf{I}}_{zz} - \overleftrightarrow{\mathsf{I}}_{xx} \right) j_z j_x \\
+\dot{j_{z}} & = & \tau^b_z(t) +
+\left(\overleftrightarrow{\mathsf{I}}_{xx} - \overleftrightarrow{\mathsf{I}}_{yy} \right) j_x j_y 
+\end{eqnarray}
+which utilize the body-fixed torques, ${\bf \tau}^b$. Torques are
+most easily derived in the space-fixed frame, 
+\begin{equation}
+{\bf \tau}^b(t) = \mathsf{A}(t) \cdot {\bf \tau}^s(t)
+\end{equation}
+where the torques are either derived from the forces on the
+constituent atoms of the rigid body, or for directional atoms,
+directly from derivatives of the potential energy,
+\begin{equation}
+{\bf \tau}^s(t) = - \hat{\bf u}(t) \times \left( \frac{\partial}
+{\partial \hat{\bf u}} V\left(\left\{ {\bf r}(t) \right\}, \left\{
+\mathsf{A}(t) \right\}\right) \right).
+\end{equation}
+Here $\hat{\bf u}$ is a unit vector pointing along the principal axis
+of the particle in the space-fixed frame.
+
+The DLM method uses a Trotter factorization of the orientational
+propagator.  This has three effects:
+\begin{enumerate}
+\item the integrator is area-preserving in phase space (i.e. it is
+{\it symplectic}),
+\item the integrator is time-{\it reversible}, making it suitable for Hybrid
+Monte Carlo applications, and
+\item the error for a single time step is of order $O\left(h^3\right)$
+for timesteps of length $h$.
+\end{enumerate}
+
+The integration of the equations of motion is carried out in a
+velocity-Verlet style 2-part algorithm:
+
+{\tt moveA:}
+\begin{eqnarray}
+{\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
+v}(t) + \frac{\delta t}{2} \left( {\bf f}(t) / m \right) \\
+{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t  {\bf
+v}\left(t + \delta t / 2 \right) \\
+{\bf j}\left(t + \delta t / 2 \right)  & \leftarrow & {\bf
+j}(t) + \frac{\delta t}{2} {\bf \tau}^b(t)  \\
+\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left( \delta t
+{\bf j}(t + \delta t / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1}
+\right)
+\end{eqnarray}
+
+In this context, the $\mathrm{rot}$ function is the reversible product
+of the three body-fixed rotations,
+\begin{equation}
+\mathrm{rot}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
+\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y /
+2) \cdot \mathsf{G}_x(a_x /2)
+\end{equation}
+where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, rotates
+both the rotation matrix ($\mathsf{A}$) and the body-fixed angular
+momentum (${\bf j}$) by an angle $\theta$ around body-fixed axis
+$\alpha$,
+\begin{equation}
+\mathsf{G}_\alpha( \theta ) = \left\{
+\begin{array}{lcl}
+\mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T \\
+{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0)
+\end{array}
+\right.
+\end{equation}
+$\mathsf{R}_\alpha$ is a quadratic approximation to
+the single-axis rotation matrix.  For example, in the small-angle
+limit, the rotation matrix around the body-fixed x-axis can be
+approximated as
+\begin{equation}
+\mathsf{R}_x(\theta) \approx \left(
+\begin{array}{ccc}
+1 & 0 & 0 \\
+0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
+\theta^2 / 4} \\
+0 & \frac{\theta}{1+
+\theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}
+\end{array}
+\right).
+\end{equation}
+All other rotations follow in a straightforward manner.
 
+After the first part of the propagation, the forces and body-fixed
+torques are calculated at the new positions and orientations
+
+{\tt doForces:}
+\begin{eqnarray}
+{\bf f}(t + \delta t) & \leftarrow & - \left(\frac{\partial V}{\partial {\bf
+r}}\right)_{{\bf r}(t + \delta t)} \\
+{\bf \tau}^{s}(t + \delta t) & \leftarrow & {\bf u}(t + \delta t)
+\times \frac{\partial V}{\partial {\bf u}} \\
+{\bf \tau}^{b}(t + \delta t) & \leftarrow & \mathsf{A}(t + \delta t)
+\cdot {\bf \tau}^s(t + \delta t) 
+\end{eqnarray}
+
+{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix
+$\mathsf{A}$ is calculated in {\tt moveA}.  Once the forces and
+torques have been obtained at the new time step, the velocities can be
+advanced to the same time value.
+
+{\tt moveB:}
+\begin{eqnarray}
+{\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
+v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
+{\bf f}(t + \delta t) / m \right) \\
+{\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
+j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} {\bf
+\tau}^b(t + \delta t)  
+\end{eqnarray}
+
+The matrix rotations used in the DLM method end up being more costly
+computationally than the simpler arithmetic quaternion
+propagation. With the same time step, a 1000-molecule water simulation
+shows an average 7\% increase in computation time using the DLM method
+in place of quaternions. This cost is more than justified when
+comparing the energy conservation of the two methods as illustrated in
+figure \ref{timestep}.
+
 \begin{figure}
 \centering
 \includegraphics[width=\linewidth]{timeStep.eps}
-\caption[Energy conservation for quaternion versus {\sc dlm} dynamics]{Energy conservation using quaternion based integration versus 
-the {\sc dlm} method with
-increasing time step. For each time step, the dotted line is total
-energy using the {\sc dlm} integrator, and the solid line comes
-from the quaternion integrator. The larger time step plots are shifted
-up from the true energy baseline for clarity.}
+\caption[Energy conservation for quaternion versus DLM dynamics]{Energy conservation using quaternion based integration versus 
+the method proposed by Dullweber \emph{et al.} with increasing time
+step. For each time step, the dotted line is total energy using the
+DLM integrator, and the solid line comes from the quaternion
+integrator. The larger time step plots are shifted up from the true
+energy baseline for clarity.}
 \label{timestep}
 \end{figure}
 
-In Fig.~\ref{timestep}, the resulting energy drift at various time
-steps for both the {\sc dlm} and quaternion integration schemes
-is compared. All of the 1000 SSD particle simulations started with the
+In figure \ref{timestep}, the resulting energy drift at various time
+steps for both the DLM and quaternion integration schemes is
+compared. All of the 1000 molecule water simulations started with the
 same configuration, and the only difference was the method for
 handling rotational motion. At time steps of 0.1 and 0.5 fs, both
-methods for propagating particle rotation conserve energy fairly well,
+methods for propagating molecule rotation conserve energy fairly well,
 with the quaternion method showing a slight energy drift over time in
 the 0.5 fs time step simulation. At time steps of 1 and 2 fs, the
-energy conservation benefits of the {\sc dlm} method are clearly
+energy conservation benefits of the DLM method are clearly
 demonstrated. Thus, while maintaining the same degree of energy
 conservation, one can take considerably longer time steps, leading to
 an overall reduction in computation time.
 
-Energy drift in these SSD particle simulations was unnoticeable for
-time steps up to three femtoseconds. A slight energy drift on the
-order of 0.012 kcal/mol per nanosecond was observed at a time step of
-four femtoseconds, and as expected, this drift increases dramatically
-with increasing time step.
+There is only one specific keyword relevant to the default integrator,
+and that is the time step for integrating the equations of motion.
 
+\begin{center}
+\begin{tabular}{llll}
+{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf
+default value} \\  
+$\delta t$ & {\tt dt = 2.0;} & fs & none 
+\end{tabular}
+\end{center}
 
 \subsection{\label{sec:extended}Extended Systems for other Ensembles}
 
+{\sc oopse} implements a number of extended system integrators for
+sampling from other ensembles relevant to chemical physics.  The
+integrator can selected with the {\tt ensemble} keyword in the
+{\tt .bass} file:
 
-{\sc oopse} implements a 
+\begin{center}
+\begin{tabular}{lll}
+{\bf Integrator} & {\bf Ensemble} & {\bf {\tt .bass} line} \\
+NVE & microcanonical & {\tt ensemble = ``NVE''; } \\
+NVT & canonical & {\tt ensemble = ``NVT''; } \\
+NPTi & isobaric-isothermal (with isotropic volume changes) & {\tt
+ensemble = ``NPTi'';} \\
+NPTf & isobaric-isothermal (with changes to box shape) & {\tt
+ensemble = ``NPTf'';} \\
+NPTxyz & approximate isobaric-isothermal & {\tt ensemble =
+``NPTxyz'';} \\
+ &  (with separate barostats on each box dimension) &
+\end{tabular}
+\end{center}
 
+The relatively well-known Nos\'e-Hoover thermostat is implemented in
+{\sc oopse}'s NVT integrator.  This method couples an extra degree of
+freedom (the thermostat) to the kinetic energy of the system, and has
+been shown to sample the canonical distribution in the system degrees
+of freedom while conserving a quantity that is, to within a constant,
+the Helmholtz free energy.
 
-\subsection{\label{oopseSec:noseHooverThermo}Nose-Hoover Thermostatting}
+NPT algorithms attempt to maintain constant pressure in the system by
+coupling the volume of the system to a barostat.  {\sc oopse} contains
+three different constant pressure algorithms.  The first two, NPTi and
+NPTf have been shown to conserve a quantity that is, to within a
+constant, the Gibbs free energy.  The Melchionna modification to the
+Hoover barostat is implemented in both NPTi and NPTf.  NPTi allows
+only isotropic changes in the simulation box, while box {\it shape}
+variations are allowed in NPTf.  The NPTxyz integrator has {\it not}
+been shown to sample from the isobaric-isothermal ensemble.  It is
+useful, however, in that it maintains orthogonality for the axes of
+the simulation box while attempting to equalize pressure along the
+three perpendicular directions in the box.
 
-To mimic the effects of being in a constant temperature ({\sc nvt})
-ensemble, {\sc oopse} uses the Nose-Hoover extended system
-approach.\cite{Hoover85} In this method, the equations of motion for
-the particle positions and velocities are
+Each of the extended system integrators requires additional keywords
+to set target values for the thermodynamic state variables that are
+being held constant.  Keywords are also required to set the
+characteristic decay times for the dynamics of the extended
+variables.
+
+\begin{tabular}{llll}
+{\bf variable} & {\bf {\tt .bass} keyword} & {\bf units} & {\bf
+default value} \\  
+$T_{\mathrm{target}}$ & {\tt targetTemperature = 300;} &  K & none \\
+$P_{\mathrm{target}}$ & {\tt targetPressure = 1;} & atm & none \\
+$\tau_T$ & {\tt tauThermostat = 1e3;} & fs & none \\
+$\tau_B$ & {\tt tauBarostat = 5e3;} & fs  & none \\
+         & {\tt resetTime = 200;} & fs & none \\
+         & {\tt useInitialExtendedSystemState = ``true'';} & logical &
+false
+\end{tabular}
+
+Two additional keywords can be used to either clear the extended
+system variables periodically ({\tt resetTime}), or to maintain the
+state of the extended system variables between simulations ({\tt
+useInitialExtendedSystemState}).  More details on these variables
+and their use in the integrators follows below.
+
+\subsubsection{\label{oopseSec:noseHooverThermo}Nos\'{e}-Hoover Thermostatting}
+
+The Nos\'e-Hoover equations of motion are given by\cite{Hoover85}
 \begin{eqnarray}
 \dot{{\bf r}} & = & {\bf v} \\
-\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v}
+\dot{{\bf v}} & = & \frac{{\bf f}}{m} - \chi {\bf v} \\
+\dot{\mathsf{A}} & = & \mathsf{A} \cdot
+\mbox{ skew}\left(\overleftrightarrow{\mathsf{I}}^{-1} \cdot {\bf j}\right) \\
+\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{\mathsf{I}}^{-1}
+\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial
+V}{\partial \mathsf{A}} \right) - \chi {\bf j}
 \label{eq:nosehoovereom}
 \end{eqnarray}
 
 $\chi$ is an ``extra'' variable included in the extended system, and
 it is propagated using the first order equation of motion
 \begin{equation}
-\dot{\chi} = \frac{1}{\tau_{T}} \left( \frac{T}{T_{target}} - 1 \right)
+\dot{\chi} = \frac{1}{\tau_{T}^2} \left( \frac{T}{T_{\mathrm{target}}} - 1 \right).
 \label{eq:nosehooverext}
 \end{equation}
-where $T_{target}$ is the target temperature for the simulation, and
-$\tau_{T}$ is a time constant for the thermostat.  
 
-To select the Nose-Hoover {\sc nvt} ensemble, the {\tt ensemble = NVT;} 
-command would be used in the simulation's {\sc bass} file.  There is
-some subtlety in choosing values for $\tau_{T}$, and it is usually set
-to values of a few ps.  Within a {\sc bass} file, $\tau_{T}$ could be
-set to 1 ps using the {\tt tauThermostat = 1000; } command.
+The instantaneous temperature $T$ is proportional to the total kinetic
+energy (both translational and orientational) and is given by
+\begin{equation}
+T = \frac{2 K}{f k_B}
+\end{equation}
+Here, $f$ is the total number of degrees of freedom in the system,
+\begin{equation}
+f = 3 N + 3 N_{\mathrm{orient}} - N_{\mathrm{constraints}}
+\end{equation}
+and $K$ is the total kinetic energy,
+\begin{equation}
+K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {\bf v}_i +
+\sum_{i=1}^{N_{\mathrm{orient}}}  \frac{1}{2} {\bf j}_i^T \cdot
+\overleftrightarrow{\mathsf{I}}_i^{-1} \cdot {\bf j}_i
+\end{equation}
 
+In eq.(\ref{eq:nosehooverext}), $\tau_T$ is the time constant for
+relaxation of the temperature to the target value.  To set values for
+$\tau_T$ or $T_{\mathrm{target}}$ in a simulation, one would use the
+{\tt tauThermostat} and {\tt targetTemperature} keywords in the {\tt
+.bass} file.  The units for {\tt tauThermostat} are fs, and the units
+for the {\tt targetTemperature} are degrees K.   The integration of
+the equations of motion is carried out in a velocity-Verlet style 2
+part algorithm:
+
+{\tt moveA:}
+\begin{eqnarray}
+T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
+{\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
+v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
+\chi(t)\right) \\
+{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t {\bf
+v}\left(t + \delta t / 2 \right) \\
+{\bf j}\left(t + \delta t / 2 \right)  & \leftarrow & {\bf
+j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
+\chi(t) \right) \\
+\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t *
+{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
+\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) +
+\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
+\right) 
+\end{eqnarray}
+
+Here $\mathrm{rot}(\delta t * {\bf j}
+\overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic Trotter
+factorization of the three rotation operations that was discussed in
+the section on the DLM integrator.  Note that this operation modifies
+both the rotation matrix $\mathsf{A}$ and the angular momentum ${\bf
+j}$.  {\tt moveA} propagates velocities by a half time step, and
+positional degrees of freedom by a full time step.  The new positions
+(and orientations) are then used to calculate a new set of forces and
+torques in exactly the same way they are calculated in the {\tt
+doForces} portion of the DLM integrator.
+
+Once the forces and torques have been obtained at the new time step,
+the temperature, velocities, and the extended system variable can be
+advanced to the same time value.
+
+{\tt moveB:}
+\begin{eqnarray}
+T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\},
+\left\{{\bf j}(t + \delta t)\right\} \\
+\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t /
+2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta
+t)}{T_{\mathrm{target}}} - 1 \right) \\
+{\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
+v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
+\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t)
+\chi(t \delta t)\right) \\
+{\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
+j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf
+\tau}^b(t + \delta t) - {\bf j}(t + \delta t)
+\chi(t + \delta t) \right) 
+\end{eqnarray}
+
+Since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$ are required
+to caclculate $T(t + \delta t)$ as well as $\chi(t + \delta t)$, they
+indirectly depend on their own values at time $t + \delta t$.  {\tt
+moveB} is therefore done in an iterative fashion until $\chi(t +
+\delta t)$ becomes self-consistent.  The relative tolerance for the
+self-consistency check defaults to a value of $\mbox{10}^{-6}$, but
+{\sc oopse} will terminate the iteration after 4 loops even if the
+consistency check has not been satisfied.
+
+The Nos\'e-Hoover algorithm is known to conserve a Hamiltonian for the
+extended system that is, to within a constant, identical to the
+Helmholtz free energy,
+\begin{equation}
+H_{\mathrm{NVT}} = V + K + f k_B T_{\mathrm{target}} \left(
+\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime
+\right)
+\end{equation}
+Poor choices of $\delta t$ or $\tau_T$ can result in non-conservation
+of $H_{\mathrm{NVT}}$, so the conserved quantity is maintained in the
+last column of the {\tt .stat} file to allow checks on the quality of
+the integration.
+
+Bond constraints are applied at the end of both the {\tt moveA} and
+{\tt moveB} portions of the algorithm.  Details on the constraint
+algorithms are given in section \ref{oopseSec:rattle}.
+
+\subsubsection{\label{sec:NPTi}Constant-pressure integration with 
+isotropic box deformations (NPTi)}
+
+To carry out isobaric-isothermal ensemble calculations {\sc oopse}
+implements the Melchionna modifications to the Nos\'e-Hoover-Andersen
+equations of motion,\cite{melchionna93}
+
+\begin{eqnarray}
+\dot{{\bf r}} & = & {\bf v} + \eta \left( {\bf r} - {\bf R}_0 \right) \\
+\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\eta + \chi) {\bf v} \\
+\dot{\mathsf{A}} & = & \mathsf{A} \cdot
+\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\
+\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1}
+\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial
+V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\
+\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left(
+\frac{T}{T_{\mathrm{target}}} - 1 \right) \\
+\dot{\eta} & = & \frac{1}{\tau_{B}^2 f k_B T_{\mathrm{target}}} V \left( P -
+P_{\mathrm{target}} \right) \\
+\dot{\mathcal{V}} & = & 3 \mathcal{V} \eta 
+\label{eq:melchionna1}
+\end{eqnarray}
+
+$\chi$ and $\eta$ are the ``extra'' degrees of freedom in the extended
+system.  $\chi$ is a thermostat, and it has the same function as it
+does in the Nos\'e-Hoover NVT integrator.  $\eta$ is a barostat which
+controls changes to the volume of the simulation box.  ${\bf R}_0$ is
+the location of the center of mass for the entire system, and
+$\mathcal{V}$ is the volume of the simulation box.  At any time, the
+volume can be calculated from the determinant of the matrix which
+describes the box shape:
+\begin{equation}
+\mathcal{V} = \det(\mathsf{H})
+\end{equation}
+
+The NPTi integrator requires an instantaneous pressure. This quantity
+is calculated via the pressure tensor,
+\begin{equation}
+\overleftrightarrow{\mathsf{P}}(t) = \frac{1}{\mathcal{V}(t)} \left(
+\sum_{i=1}^{N} m_i {\bf v}_i(t) \otimes {\bf v}_i(t) \right) +
+\overleftrightarrow{\mathsf{W}}(t)
+\end{equation}
+The kinetic contribution to the pressure tensor utilizes the {\it
+outer} product of the velocities denoted by the $\otimes$ symbol.  The
+stress tensor is calculated from another outer product of the
+inter-atomic separation vectors (${\bf r}_{ij} = {\bf r}_j - {\bf
+r}_i$) with the forces between the same two atoms,
+\begin{equation}
+\overleftrightarrow{\mathsf{W}}(t) = \sum_{i} \sum_{j>i} {\bf r}_{ij}(t)
+\otimes {\bf f}_{ij}(t)
+\end{equation}
+The instantaneous pressure is then simply obtained from the trace of
+the Pressure tensor,
+\begin{equation}
+P(t) = \frac{1}{3} \mathrm{Tr} \left( \overleftrightarrow{\mathsf{P}}(t)
+\right)
+\end{equation}
+
+In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for
+relaxation of the pressure to the target value.  To set values for
+$\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the
+{\tt tauBarostat} and {\tt targetPressure} keywords in the {\tt .bass}
+file.  The units for {\tt tauBarostat} are fs, and the units for the
+{\tt targetPressure} are atmospheres.  Like in the NVT integrator, the
+integration of the equations of motion is carried out in a
+velocity-Verlet style 2 part algorithm:
+
+{\tt moveA:}
+\begin{eqnarray}
+T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
+P(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(t)\right\} \\
+{\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
+v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
+\left(\chi(t) + \eta(t) \right) \right) \\
+{\bf j}\left(t + \delta t / 2 \right)  & \leftarrow & {\bf
+j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
+\chi(t) \right) \\
+\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t *
+{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
+\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) +
+\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
+\right) \\
+\eta(t + \delta t / 2) & \leftarrow & \eta(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B
+T(t) \tau_B^2} \left( P(t) - P_{\mathrm{target}} \right) \\ 
+{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf
+v}\left(t + \delta t / 2 \right) + \eta(t + \delta t / 2)\left[ {\bf
+r}(t + \delta t) - {\bf R}_0 \right] \right\} \\
+\mathsf{H}(t + \delta t) & \leftarrow & e^{-\delta t \eta(t + \delta t
+/ 2)} \mathsf{H}(t)
+\end{eqnarray}
+
+Most of these equations are identical to their counterparts in the NVT
+integrator, but the propagation of positions to time $t + \delta t$
+depends on the positions at the same time.  {\sc oopse} carries out
+this step iteratively (with a limit of 5 passes through the iterative
+loop).  Also, the simulation box $\mathsf{H}$ is scaled uniformly for
+one full time step by an exponential factor that depends on the value
+of $\eta$ at time $t +
+\delta t / 2$.  Reshaping the box uniformly also scales the volume of
+the box by
+\begin{equation}
+\mathcal{V}(t + \delta t) \leftarrow e^{ - 3 \delta t \eta(t + \delta t /2)}
+\mathcal{V}(t)
+\end{equation}
+
+The {\tt doForces} step for the NPTi integrator is exactly the same as
+in both the DLM and NVT integrators.  Once the forces and torques have
+been obtained at the new time step, the velocities can be advanced to
+the same time value.
+
+{\tt moveB:}
+\begin{eqnarray}
+T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\},
+\left\{{\bf j}(t + \delta t)\right\} \\
+P(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\},
+\left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(t + \delta t)\right\} \\
+\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t /
+2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta
+t)}{T_{\mathrm{target}}} - 1 \right) \\
+\eta(t + \delta t) & \leftarrow & \eta(t + \delta t / 2) +
+\frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2}
+\left( P(t + \delta t) - P_{\mathrm{target}} 
+\right) \\ 
+{\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
+v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
+\frac{{\bf f}(t + \delta t)}{m} - {\bf v}(t + \delta t)
+(\chi(t + \delta t) + \eta(t + \delta t)) \right) \\
+{\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
+j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf
+\tau}^b(t + \delta t) - {\bf j}(t + \delta t)
+\chi(t + \delta t) \right) 
+\end{eqnarray}
+
+Once again, since ${\bf v}(t + \delta t)$ and ${\bf j}(t + \delta t)$
+are required to caclculate $T(t + \delta t)$, $P(t + \delta t)$, $\chi(t +
+\delta t)$, and $\eta(t + \delta t)$, they indirectly depend on their
+own values at time $t + \delta t$.  {\tt moveB} is therefore done in
+an iterative fashion until $\chi(t + \delta t)$ and $\eta(t + \delta
+t)$ become self-consistent.  The relative tolerance for the
+self-consistency check defaults to a value of $\mbox{10}^{-6}$, but
+{\sc oopse} will terminate the iteration after 4 loops even if the
+consistency check has not been satisfied.
+
+The Melchionna modification of the Nos\'e-Hoover-Andersen algorithm is
+known to conserve a Hamiltonian for the extended system that is, to
+within a constant, identical to the Gibbs free energy,
+\begin{equation}
+H_{\mathrm{NPTi}} = V + K + f k_B T_{\mathrm{target}} \left(
+\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime
+\right) + P_{\mathrm{target}} \mathcal{V}(t).
+\end{equation}
+Poor choices of $\delta t$, $\tau_T$, or $\tau_B$ can result in
+non-conservation of $H_{\mathrm{NPTi}}$, so the conserved quantity is
+maintained in the last column of the {\tt .stat} file to allow checks
+on the quality of the integration.  It is also known that this
+algorithm samples the equilibrium distribution for the enthalpy
+(including contributions for the thermostat and barostat), 
+\begin{equation}
+H_{\mathrm{NPTi}} = V + K + \frac{f k_B T_{\mathrm{target}}}{2} \left(
+\chi^2 \tau_T^2 + \eta^2 \tau_B^2 \right) +  P_{\mathrm{target}}
+\mathcal{V}(t). 
+\end{equation}
+
+Bond constraints are applied at the end of both the {\tt moveA} and
+{\tt moveB} portions of the algorithm.  Details on the constraint
+algorithms are given in section \ref{oopseSec:rattle}.
+
+\subsubsection{\label{sec:NPTf}Constant-pressure integration with a
+flexible box (NPTf)} 
+
+There is a relatively simple generalization of the
+Nos\'e-Hoover-Andersen method to include changes in the simulation box
+{\it shape} as well as in the volume of the box.  This method utilizes
+the full $3 \times 3$ pressure tensor and introduces a tensor of
+extended variables ($\overleftrightarrow{\eta}$) to control changes to
+the box shape.  The equations of motion for this method are
+\begin{eqnarray}
+\dot{{\bf r}} & = & {\bf v} + \overleftrightarrow{\eta} \cdot \left( {\bf r} - {\bf R}_0 \right) \\
+\dot{{\bf v}} & = & \frac{{\bf f}}{m} - (\overleftrightarrow{\eta} +
+\chi \mathsf{1}) {\bf v} \\
+\dot{\mathsf{A}} & = & \mathsf{A} \cdot
+\mbox{ skew}\left(\overleftrightarrow{I}^{-1} \cdot {\bf j}\right) \\
+\dot{{\bf j}} & = & {\bf j} \times \left( \overleftrightarrow{I}^{-1}
+\cdot {\bf j} \right) - \mbox{ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial
+V}{\partial \mathsf{A}} \right) - \chi {\bf j} \\
+\dot{\chi} & = & \frac{1}{\tau_{T}^2} \left(
+\frac{T}{T_{\mathrm{target}}} - 1 \right) \\
+\dot{\overleftrightarrow{eta}} & = & \frac{1}{\tau_{B}^2 f k_B
+T_{\mathrm{target}}} V \left( \overleftrightarrow{\mathsf{P}} - P_{\mathrm{target}}\mathsf{1} \right) \\
+\dot{\mathsf{H}} & = &  \overleftrightarrow{\eta} \cdot \mathsf{H} 
+\label{eq:melchionna2}
+\end{eqnarray}
+
+Here, $\mathsf{1}$ is the unit matrix and $\overleftrightarrow{\mathsf{P}}$
+is the pressure tensor.  Again, the volume, $\mathcal{V} = \det
+\mathsf{H}$. 
+
+The propagation of the equations of motion is nearly identical to the
+NPTi integration:
+
+{\tt moveA:}
+\begin{eqnarray}
+T(t) & \leftarrow & \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} \\
+\overleftrightarrow{\mathsf{P}}(t) & \leftarrow & \left\{{\bf r}(t)\right\}, \left\{{\bf v}(t)\right\}, \left\{{\bf f}(t)\right\} \\
+{\bf v}\left(t + \delta t / 2\right)  & \leftarrow & {\bf
+v}(t) + \frac{\delta t}{2} \left( \frac{{\bf f}(t)}{m} - 
+\left(\chi(t)\mathsf{1} + \overleftrightarrow{\eta}(t) \right) \cdot
+{\bf v}(t) \right) \\ 
+{\bf j}\left(t + \delta t / 2 \right)  & \leftarrow & {\bf
+j}(t) + \frac{\delta t}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
+\chi(t) \right) \\
+\mathsf{A}(t + \delta t) & \leftarrow & \mathrm{rot}\left(\delta t *
+{\bf j}(t + \delta t / 2) \overleftrightarrow{\mathsf{I}}^{-1} \right) \\
+\chi\left(t + \delta t / 2 \right) & \leftarrow & \chi(t) +
+\frac{\delta t}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
+\right) \\
+\overleftrightarrow{\eta}(t + \delta t / 2) & \leftarrow & \overleftrightarrow{\eta}(t) + \frac{\delta t \mathcal{V}(t)}{2 N k_B
+T(t) \tau_B^2} \left( \overleftrightarrow{\mathsf{P}}(t) - P_{\mathrm{target}}\mathsf{1} \right) \\ 
+{\bf r}(t + \delta t) & \leftarrow & {\bf r}(t) + \delta t \left\{ {\bf
+v}\left(t + \delta t / 2 \right) + \overleftrightarrow{\eta}(t +
+\delta t / 2) \cdot \left[ {\bf
+r}(t + \delta t) - {\bf R}_0 \right] \right\} \\
+\mathsf{H}(t + \delta t) & \leftarrow & \mathsf{H}(t) \cdot e^{-\delta t
+\overleftrightarrow{\eta}(t + \delta t / 2)} 
+\end{eqnarray}
+{\sc oopse} uses a power series expansion truncated at second order
+for the exponential operation which scales the simulation box.
+
+The {\tt moveB} portion of the algorithm is largely unchanged from the
+NPTi integrator:
+
+{\tt moveB:}
+\begin{eqnarray}
+T(t + \delta t) & \leftarrow & \left\{{\bf v}(t + \delta t)\right\},
+\left\{{\bf j}(t + \delta t)\right\} \\
+\overleftrightarrow{\mathsf{P}}(t + \delta t) & \leftarrow & \left\{{\bf r}(t + \delta t)\right\},
+\left\{{\bf v}(t + \delta t)\right\}, \left\{{\bf f}(t + \delta t)\right\} \\
+\chi\left(t + \delta t \right) & \leftarrow & \chi\left(t + \delta t /
+2 \right) + \frac{\delta t}{2 \tau_T^2} \left( \frac{T(t+\delta
+t)}{T_{\mathrm{target}}} - 1 \right) \\
+\overleftrightarrow{\eta}(t + \delta t) & \leftarrow & \overleftrightarrow{\eta}(t + \delta t / 2) +
+\frac{\delta t \mathcal{V}(t + \delta t)}{2 N k_B T(t + \delta t) \tau_B^2}
+\left( \overleftrightarrow{P}(t + \delta t) - P_{\mathrm{target}}\mathsf{1} 
+\right) \\ 
+{\bf v}\left(t + \delta t \right)  & \leftarrow & {\bf
+v}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left(
+\frac{{\bf f}(t + \delta t)}{m} - 
+(\chi(t + \delta t)\mathsf{1} + \overleftrightarrow{\eta}(t + \delta
+t)) \right) \cdot {\bf v}(t + \delta t) \\
+{\bf j}\left(t + \delta t \right)  & \leftarrow & {\bf
+j}\left(t + \delta t / 2 \right) + \frac{\delta t}{2} \left( {\bf
+\tau}^b(t + \delta t) - {\bf j}(t + \delta t)
+\chi(t + \delta t) \right) 
+\end{eqnarray}
+
+The iterative schemes for both {\tt moveA} and {\tt moveB} are
+identical to those described for the NPTi integrator.
+
+The NPTf integrator is known to conserve the following Hamiltonian:
+\begin{equation}
+H_{\mathrm{NPTf}} = V + K + f k_B T_{\mathrm{target}} \left(
+\frac{\tau_{T}^2 \chi^2(t)}{2} + \int_{0}^{t} \chi(t^\prime) dt^\prime
+\right) + P_{\mathrm{target}} \mathcal{V}(t) + \frac{f k_B
+T_{\mathrm{target}}}{2}
+\mathrm{Tr}\left[\overleftrightarrow{\eta}(t)\right]^2 \tau_B^2.
+\end{equation}
+
+This integrator must be used with care, particularly in liquid
+simulations.  Liquids have very small restoring forces in the
+off-diagonal directions, and the simulation box can very quickly form
+elongated and sheared geometries which become smaller than the
+electrostatic or Lennard-Jones cutoff radii.  It finds most use in
+simulating crystals or liquid crystals which assume non-orthorhombic
+geometries.
+
+\subsubsection{\label{nptxyz}Constant pressure in 3 axes (NPTxyz)}
+
+There is one additional extended system integrator which is somewhat
+simpler than the NPTf method described above.  In this case, the three
+axes have independent barostats which each attempt to preserve the
+target pressure along the box walls perpendicular to that particular
+axis.  The lengths of the box axes are allowed to fluctuate
+independently, but the angle between the box axes does not change.
+The equations of motion are identical to those described above, but
+only the {\it diagonal} elements of $\overleftrightarrow{\eta}$ are
+computed.  The off-diagonal elements are set to zero (even when the
+pressure tensor has non-zero off-diagonal elements).
+
+It should be noted that the NPTxyz integrator is {\it not} known to
+preserve any Hamiltonian of interest to the chemical physics
+community.  The integrator is extremely useful, however, in generating
+initial conditions for other integration methods.  It {\it is} suitable
+for use with liquid simulations, or in cases where there is
+orientational anisotropy in the system (i.e. in lipid bilayer
+simulations).
+
 \subsection{\label{oopseSec:rattle}The {\sc rattle} Method for Bond 
 	Constraints}
 
@@ -1222,11 +1865,12 @@ coefficient can be approximated by
 If the time-dependent friction decays rapidly, the static friction
 coefficient can be approximated by
 \begin{equation}
-\xi^{static}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt
+\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta F(z,0)\rangle dt
 \end{equation}
-Therefore, the diffusion constant can then be estimated by
+Allowing diffusion constant to then be calculated through the
+Einstein relation:\cite{Marrink94}
 \begin{equation}
-D(z)=\frac{k_{B}T}{\xi^{static}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
+D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
 }\langle\delta F(z,t)\delta F(z,0)\rangle dt}%
 \end{equation}
 
@@ -1236,16 +1880,39 @@ resetting the coordinate, we reset the forces of z-con
 auto-correlation calculation.\cite{Marrink94} However, simply resetting the
 coordinate will move the center of the mass of the whole system. To
 avoid this problem, a new method was used in {\sc oopse}. Instead of
-resetting the coordinate, we reset the forces of z-constraint
+resetting the coordinate, we reset the forces of z-constrained
 molecules as well as subtract the total constraint forces from the
-rest of the system after force calculation at each time step.
-\begin{align}
-F_{\alpha i}&=0\\
-V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{i}M_{_{\alpha i}}}\\
-F_{\alpha i}&=F_{\alpha i}-\frac{M_{_{\alpha i}}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}\sum\limits_{\beta}F_{\beta}\\
-V_{\alpha i}&=V_{\alpha i}-\frac{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}V_{\alpha i}}{\sum\limits_{\alpha}\sum\limits_{i}M_{_{\alpha i}}}
-\end{align}
+rest of the system after the force calculation at each time step.
 
+After the force calculation, define $G_\alpha$ as
+\begin{equation}
+G_{\alpha} = \sum_i F_{\alpha i}
+\label{oopseEq:zc1}
+\end{equation}
+Where $F_{\alpha i}$ is the force in the z direction of atom $i$ in
+z-constrained molecule $\alpha$. The forces of the z constrained
+molecule are then set to:
+\begin{equation}
+F_{\alpha i} = F_{\alpha i} - 
+	\frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}}
+\end{equation}
+Here, $m_{\alpha i}$ is the mass of atom $i$ in the z-constrained
+molecule. Having rescaled the forces, the velocities must also be
+rescaled to subtract out any center of mass velocity in the z
+direction.
+\begin{equation}
+v_{\alpha i} = v_{\alpha i} -
+	\frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}
+\end{equation}
+Where $v_{\alpha i}$ is the velocity of atom $i$ in the z direction.
+Lastly, all of the accumulated z constrained forces must be subtracted
+from the system to keep the system center of mass from drifting.
+\begin{equation}
+F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha} G_{\alpha}}
+	{\sum_{\beta}\sum_i m_{\beta i}}
+\end{equation}
+Where $\beta$ are all of the unconstrained molecules in the system.
+
 At the very beginning of the simulation, the molecules may not be at their
 constrained positions. To move a z-constrained molecule to its specified
 position, a simple harmonic potential is used