--- trunk/tengDissertation/Methodology.tex	2006/04/03 18:07:54	2685
+++ trunk/tengDissertation/Methodology.tex	2006/06/11 02:53:15	2854
@@ -2,12 +2,744 @@
 
 \section{\label{methodSection:rigidBodyIntegrators}Integrators for Rigid Body Motion in Molecular Dynamics}
 
-\section{\label{methodSection:conservedQuantities}Integrators to Conserve Properties in Special Ensembles}
+In order to mimic the experiments, which are usually performed under
+constant temperature and/or pressure, extended Hamiltonian system
+methods have been developed to generate statistical ensembles, such
+as canonical ensemble and isobaric-isothermal ensemble \textit{etc}.
+In addition to the standard ensemble, specific ensembles have been
+developed to account for the anisotropy between the lateral and
+normal directions of membranes. The $NPAT$ ensemble, in which the
+normal pressure and the lateral surface area of the membrane are
+kept constant, and the $NP\gamma T$ ensemble, in which the normal
+pressure and the lateral surface tension are kept constant were
+proposed to address this issue.
 
-\section{\label{methodSection:hydrodynamics}Hydrodynamics}
+Integration schemes for rotational motion of the rigid molecules in
+microcanonical ensemble have been extensively studied in the last
+two decades. Matubayasi developed a time-reversible integrator for
+rigid bodies in quaternion representation. Although it is not
+symplectic, this integrator still demonstrates a better long-time
+energy conservation than traditional methods because of the
+time-reversible nature. Extending Trotter-Suzuki to general system
+with a flat phase space, Miller and his colleagues devised an novel
+symplectic, time-reversible and volume-preserving integrator in
+quaternion representation, which was shown to be superior to the
+Matubayasi's time-reversible integrator. However, all of the
+integrators in quaternion representation suffer from the
+computational penalty of constructing a rotation matrix from
+quaternions to evolve coordinates and velocities at every time step.
+An alternative integration scheme utilizing rotation matrix directly
+proposed by Dullweber, Leimkuhler and McLachlan (DLM) also preserved
+the same structural properties of the Hamiltonian flow. In this
+section, the integration scheme of DLM method will be reviewed and
+extended to other ensembles.
 
-\section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies}
+\subsection{\label{methodSection:DLM}Integrating the Equations of Motion: the
+DLM method}
 
-\section{\label{methodSection:coarseGrained}Coarse-Grained Modeling}
+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 $\mathcal{O}\left(h^4\right)$
+for timesteps of length $h$.
+\end{enumerate}
 
-\section{\label{methodSection:moleculeScale}Molecular-Scale Modeling}
+The integration of the equations of motion is carried out in a
+velocity-Verlet style 2-part algorithm, where $h= \delta t$:
+
+{\tt moveA:}
+\begin{align*}
+{\bf v}\left(t + h / 2\right)  &\leftarrow  {\bf v}(t)
+    + \frac{h}{2} \left( {\bf f}(t) / m \right), \\
+%
+{\bf r}(t + h) &\leftarrow {\bf r}(t)
+    + h  {\bf v}\left(t + h / 2 \right), \\
+%
+{\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
+    + \frac{h}{2} {\bf \tau}^b(t), \\
+%
+\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
+    (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
+\end{align*}
+
+In this context, the $\mathrm{rotate}$ function is the reversible
+product of the three body-fixed rotations,
+\begin{equation}
+\mathrm{rotate}({\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{align*}
+{\bf f}(t + h) &\leftarrow
+    - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\
+%
+{\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h)
+    \times \frac{\partial V}{\partial {\bf u}}, \\
+%
+{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
+    \cdot {\bf \tau}^s(t + h).
+\end{align*}
+
+${\bf u}$ will be automatically updated 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{align*}
+{\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2
+\right)
+    + \frac{h}{2} \left( {\bf f}(t + h) / m \right), \\
+%
+{\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t + h / 2
+\right)
+    + \frac{h}{2} {\bf \tau}^b(t + h) .
+\end{align*}
+
+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 Fig.~\ref{methodFig:timestep}.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{timeStep.eps}
+\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{methodFig:timestep}
+\end{figure}
+
+In Fig.~\ref{methodFig: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 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 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.
+
+\subsection{\label{methodSection:NVT}Nos\'{e}-Hoover Thermostatting}
+
+The Nos\'e-Hoover equations of motion are given by\cite{Hoover1985}
+\begin{eqnarray}
+\dot{{\bf r}} & = & {\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}^2} \left(
+\frac{T}{T_{\mathrm{target}}} - 1 \right). \label{eq:nosehooverext}
+\end{equation}
+
+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{align*}
+T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\bf j}(t)\right\} ,\\
+%
+{\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t)
+    + \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
+    \chi(t)\right), \\
+%
+{\bf r}(t + h) &\leftarrow {\bf r}(t)
+    + h {\bf v}\left(t + h / 2 \right) ,\\
+%
+{\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
+    + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
+    \chi(t) \right) ,\\
+%
+\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}
+    \left(h * {\bf j}(t + h / 2)
+    \overleftrightarrow{\mathsf{I}}^{-1} \right) ,\\
+%
+\chi\left(t + h / 2 \right) &\leftarrow \chi(t)
+    + \frac{h}{2 \tau_T^2} \left( \frac{T(t)}
+    {T_{\mathrm{target}}} - 1 \right) .
+\end{align*}
+
+Here $\mathrm{rotate}(h * {\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{align*}
+T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
+    \left\{{\bf j}(t + h)\right\}, \\
+%
+\chi\left(t + h \right) &\leftarrow \chi\left(t + h /
+    2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)}
+    {T_{\mathrm{target}}} - 1 \right), \\
+%
+{\bf v}\left(t + h \right)  &\leftarrow {\bf v}\left(t
+    + h / 2 \right) + \frac{h}{2} \left(
+    \frac{{\bf f}(t + h)}{m} - {\bf v}(t + h)
+    \chi(t h)\right) ,\\
+%
+{\bf j}\left(t + h \right) &\leftarrow {\bf j}\left(t
+    + h / 2 \right) + \frac{h}{2}
+    \left( {\bf \tau}^b(t + h) - {\bf j}(t + h)
+    \chi(t + h) \right) .
+\end{align*}
+
+Since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required to
+caclculate $T(t + h)$ as well as $\chi(t + h)$, they indirectly
+depend on their own values at time $t + h$.  {\tt moveB} is
+therefore done in an iterative fashion until $\chi(t + h)$ becomes
+self-consistent.
+
+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,\cite{Melchionna1993}
+\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 $h$ 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.
+
+\subsection{\label{methodSection:NPTi}Constant-pressure integration with
+isotropic box deformations (NPTi)}
+
+Isobaric-isothermal ensemble integrator is implemented using the
+Melchionna modifications to the Nos\'e-Hoover-Andersen equations of
+motion,\cite{Melchionna1993}
+
+\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. 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{align*}
+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\} ,\\
+%
+{\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t)
+    + \frac{h}{2} \left( \frac{{\bf f}(t)}{m} - {\bf v}(t)
+    \left(\chi(t) + \eta(t) \right) \right), \\
+%
+{\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
+    + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
+    \chi(t) \right), \\
+%
+\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
+    {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1}
+    \right) ,\\
+%
+\chi\left(t + h / 2 \right) &\leftarrow \chi(t) +
+    \frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}} - 1
+    \right) ,\\
+%
+\eta(t + h / 2) &\leftarrow \eta(t) + \frac{h
+    \mathcal{V}(t)}{2 N k_B T(t) \tau_B^2} \left( P(t)
+    - P_{\mathrm{target}} \right), \\
+%
+{\bf r}(t + h) &\leftarrow {\bf r}(t) + h
+    \left\{ {\bf v}\left(t + h / 2 \right)
+    + \eta(t + h / 2)\left[ {\bf r}(t + h)
+    - {\bf R}_0 \right] \right\} ,\\
+%
+\mathsf{H}(t + h) &\leftarrow e^{-h \eta(t + h / 2)}
+    \mathsf{H}(t).
+\end{align*}
+
+Most of these equations are identical to their counterparts in the
+NVT integrator, but the propagation of positions to time $t + h$
+depends on the positions at the same time. 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 +
+h / 2$.  Reshaping the box uniformly also scales the volume of the
+box by
+\begin{equation}
+\mathcal{V}(t + h) \leftarrow e^{ - 3 h \eta(t + h /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{align*}
+T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
+    \left\{{\bf j}(t + h)\right\} ,\\
+%
+P(t + h) &\leftarrow  \left\{{\bf r}(t + h)\right\},
+    \left\{{\bf v}(t + h)\right\}, \\
+%
+\chi\left(t + h \right) &\leftarrow \chi\left(t + h /
+    2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+h)}
+    {T_{\mathrm{target}}} - 1 \right), \\
+%
+\eta(t + h) &\leftarrow \eta(t + h / 2) +
+    \frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h)
+    \tau_B^2} \left( P(t + h) - P_{\mathrm{target}} \right), \\
+%
+{\bf v}\left(t + h \right)  &\leftarrow {\bf v}\left(t
+    + h / 2 \right) + \frac{h}{2} \left(
+    \frac{{\bf f}(t + h)}{m} - {\bf v}(t + h)
+    (\chi(t + h) + \eta(t + h)) \right) ,\\
+%
+{\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t
+    + h / 2 \right) + \frac{h}{2} \left( {\bf
+    \tau}^b(t + h) - {\bf j}(t + h)
+    \chi(t + h) \right) .
+\end{align*}
+
+Once again, since ${\bf v}(t + h)$ and ${\bf j}(t + h)$ are required
+to caclculate $T(t + h)$, $P(t + h)$, $\chi(t + h)$, and $\eta(t +
+h)$, they indirectly depend on their own values at time $t + h$.
+{\tt moveB} is therefore done in an iterative fashion until $\chi(t
++ h)$ and $\eta(t + h)$ become self-consistent.
+
+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}
+
+\subsection{\label{methodSection: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 \cdot \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{align*}
+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\} ,\\
+%
+{\bf v}\left(t + h / 2\right)  &\leftarrow {\bf v}(t)
+    + \frac{h}{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 + h / 2 \right)  &\leftarrow {\bf j}(t)
+    + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
+    \chi(t) \right), \\
+%
+\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
+    {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1}
+    \right), \\
+%
+\chi\left(t + h / 2 \right) &\leftarrow \chi(t) +
+    \frac{h}{2 \tau_T^2} \left( \frac{T(t)}{T_{\mathrm{target}}}
+    - 1 \right), \\
+%
+\overleftrightarrow{\eta}(t + h / 2) &\leftarrow
+    \overleftrightarrow{\eta}(t) + \frac{h \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 + h) &\leftarrow {\bf r}(t) + h \left\{ {\bf v}
+    \left(t + h / 2 \right) + \overleftrightarrow{\eta}(t +
+    h / 2) \cdot \left[ {\bf r}(t + h)
+    - {\bf R}_0 \right] \right\}, \\
+%
+\mathsf{H}(t + h) &\leftarrow \mathsf{H}(t) \cdot e^{-h
+    \overleftrightarrow{\eta}(t + h / 2)} .
+\end{align*}
+Here, a power series expansion truncated at second order for the
+exponential operation is used to scale the simulation box.
+
+The {\tt moveB} portion of the algorithm is largely unchanged from
+the NPTi integrator:
+
+{\tt moveB:}
+\begin{align*}
+T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
+    \left\{{\bf j}(t + h)\right\}, \\
+%
+\overleftrightarrow{\mathsf{P}}(t + h) &\leftarrow \left\{{\bf r}
+    (t + h)\right\}, \left\{{\bf v}(t
+    + h)\right\}, \left\{{\bf f}(t + h)\right\} ,\\
+%
+\chi\left(t + h \right) &\leftarrow \chi\left(t + h /
+    2 \right) + \frac{h}{2 \tau_T^2} \left( \frac{T(t+
+    h)}{T_{\mathrm{target}}} - 1 \right), \\
+%
+\overleftrightarrow{\eta}(t + h) &\leftarrow
+    \overleftrightarrow{\eta}(t + h / 2) +
+    \frac{h \mathcal{V}(t + h)}{2 N k_B T(t + h)
+    \tau_B^2} \left( \overleftrightarrow{P}(t + h)
+    - P_{\mathrm{target}}\mathsf{1} \right) ,\\
+%
+{\bf v}\left(t + h \right)  &\leftarrow {\bf v}\left(t
+    + h / 2 \right) + \frac{h}{2} \left(
+    \frac{{\bf f}(t + h)}{m} -
+    (\chi(t + h)\mathsf{1} + \overleftrightarrow{\eta}(t
+    + h)) \right) \cdot {\bf v}(t + h), \\
+%
+{\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t
+    + h / 2 \right) + \frac{h}{2} \left( {\bf \tau}^b(t
+    + h) - {\bf j}(t + h) \chi(t + h) \right) .
+\end{align*}
+
+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{eqnarray*}
+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{eqnarray*}
+
+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.  The NPTf integrator
+finds most use in simulating crystals or liquid crystals which
+assume non-orthorhombic geometries.
+
+\subsubsection{\label{methodSection:NPAT}NPAT Ensemble}
+
+A comprehensive understanding of structure¨Cfunction relations of
+biological membrane system ultimately relies on structure and
+dynamics of lipid bilayer, which are strongly affected by the
+interfacial interaction between lipid molecules and surrounding
+media. One quantity to describe the interfacial interaction is so
+called the average surface area per lipid. Constat area and constant
+lateral pressure simulation can be achieved by extending the
+standard NPT ensemble with a different pressure control strategy
+
+\begin{equation}
+ \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll}
+                  \frac{{V(P_{\alpha \beta }  - P_{{\rm{target}}} )}}{{\tau_{\rm{B}}^{\rm{2}} fk_B T_{{\rm{target}}} }}
+                  & \mbox{if $ \alpha = \beta  = z$}\\
+                  0 & \mbox{otherwise}\\
+           \end{array}
+    \right.
+\end{equation}
+
+Note that the iterative schemes for NPAT are identical to those
+described for the NPTi integrator.
+
+\subsection{\label{methodSection:NPrT}NP$\gamma$T
+Ensemble}
+
+Theoretically, the surface tension $\gamma$ of a stress free
+membrane system should be zero since its surface free energy $G$ is
+minimum with respect to surface area $A$
+\[
+\gamma  = \frac{{\partial G}}{{\partial A}}.
+\]
+However, a surface tension of zero is not appropriate for relatively
+small patches of membrane. In order to eliminate the edge effect of
+the membrane simulation, a special ensemble, NP$\gamma$T, is
+proposed to maintain the lateral surface tension and normal
+pressure. The equation of motion for cell size control tensor,
+$\eta$, in $NP\gamma T$ is
+\begin{equation}
+ \dot {\overleftrightarrow{\eta}} _{\alpha \beta}=\left\{\begin{array}{ll}
+    - A_{xy} (\gamma _\alpha   - \gamma _{{\rm{target}}} ) & \mbox{$\alpha  = \beta  = x$ or $\alpha  = \beta  = y$}\\
+    \frac{{V(P_{\alpha \beta }  - P_{{\rm{target}}} )}}{{\tau _{\rm{B}}^{\rm{2}} fk_B T_{{\rm{target}}}}} & \mbox{$\alpha  = \beta  = z$} \\
+    0 & \mbox{$\alpha  \ne \beta$} \\
+       \end{array}
+    \right.
+\end{equation}
+where $ \gamma _{{\rm{target}}}$ is the external surface tension and
+the instantaneous surface tensor $\gamma _\alpha$ is given by
+\begin{equation}
+\gamma _\alpha   =  - h_z( \overleftrightarrow{P} _{\alpha \alpha }
+- P_{{\rm{target}}} )
+\label{methodEquation:instantaneousSurfaceTensor}
+\end{equation}
+
+There is one additional extended system integrator (NPTxyz), in
+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. It should be noted that the NPTxyz
+integrator is a special case of $NP\gamma T$ if the surface tension
+$\gamma$ is set to zero.
+
+\section{\label{methodSection:zcons}Z-Constraint Method}
+
+Based on the fluctuation-dissipation theorem, a force
+auto-correlation method was developed by Roux and Karplus to
+investigate the dynamics of ions inside ion channels\cite{Roux1991}.
+The time-dependent friction coefficient can be calculated from the
+deviation of the instantaneous force from its mean force.
+\begin{equation}
+\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T,
+\end{equation}
+where%
+\begin{equation}
+\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle.
+\end{equation}
+
+If the time-dependent friction decays rapidly, the static friction
+coefficient can be approximated by
+\begin{equation}
+\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta
+F(z,0)\rangle dt.
+\end{equation}
+Allowing diffusion constant to then be calculated through the
+Einstein relation:\cite{Marrink1994}
+\begin{equation}
+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}
+
+The Z-Constraint method, which fixes the z coordinates of the
+molecules with respect to the center of the mass of the system, has
+been a method suggested to obtain the forces required for the force
+auto-correlation calculation.\cite{Marrink1994} However, simply
+resetting the coordinate will move the center of the mass of the
+whole system. To avoid this problem, we reset the forces of
+z-constrained molecules as well as subtract the total constraint
+forces from the rest of the system after the force calculation at
+each time step instead of resetting the coordinate.
+
+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.
+Similarly, the velocities of the unconstrained molecules must also
+be scaled.
+\begin{equation}
+v_{\beta i} = v_{\beta i} + \sum_{\alpha}
+    \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}.
+\end{equation}
+
+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
+\begin{equation}
+U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2},%
+\end{equation}
+where $k_{\text{Harmonic}}$ is the harmonic force constant, $z(t)$
+is the current $z$ coordinate of the center of mass of the
+constrained molecule, and $z_{\text{cons}}$ is the constrained
+position. The harmonic force operating on the z-constrained molecule
+at time $t$ can be calculated by
+\begin{equation}
+F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}=
+    -k_{\text{Harmonic}}(z(t)-z_{\text{cons}}).
+\end{equation}