--- trunk/tengDissertation/Methodology.tex	2006/06/11 02:53:15	2854
+++ trunk/tengDissertation/Methodology.tex	2006/06/29 23:56:11	2911
@@ -2,37 +2,39 @@ In order to mimic the experiments, which are usually p
 
 \section{\label{methodSection:rigidBodyIntegrators}Integrators for Rigid Body Motion in Molecular Dynamics}
 
-In order to mimic the experiments, which are usually performed under
+In order to mimic 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.
+as the canonical and isobaric-isothermal ensembles. 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 the
+issues.
 
-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
+Integration schemes for the rotational motion of the rigid molecules
+in the microcanonical ensemble have been extensively studied over
+the last two decades. Matubayasi developed a
+time-reversible integrator for rigid bodies in quaternion
+representation\cite{Matubayasi1999}. Although it is not symplectic, this integrator still
+demonstrates a better long-time energy conservation than Euler angle
+methods because of the time-reversible nature. Extending the
+Trotter-Suzuki factorization to general system with a flat phase
+space, Miller\cite{Miller2002} and his colleagues devised a novel
+symplectic, time-reversible and volume-preserving integrator in the
 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
+integrators in the 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.
+An alternative integration scheme utilizing the rotation matrix
+directly proposed by Dullweber, Leimkuhler and McLachlan (DLM) also
+preserved the same structural properties of the Hamiltonian
+propagator\cite{Dullweber1997}. In this section, the integration
+scheme of DLM method will be reviewed and extended to other
+ensembles.
 
 \subsection{\label{methodSection:DLM}Integrating the Equations of Motion: the
 DLM method}
@@ -47,7 +49,6 @@ for timesteps of length $h$.
 \item the error for a single time step is of order $\mathcal{O}\left(h^4\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, where $h= \delta t$:
 
@@ -62,10 +63,9 @@ velocity-Verlet style 2-part algorithm, where $h= \del
 {\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}
+\mathsf{Q}(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}
@@ -74,13 +74,13 @@ rotates both the rotation matrix ($\mathsf{A}$) and th
 / 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
+rotates both the rotation matrix $\mathsf{Q}$ 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, \\
+\mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
 {\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
 j}(0).
 \end{array}
@@ -100,25 +100,23 @@ All other rotations follow in a straightforward manner
 \end{array}
 \right).
 \end{equation}
-All other rotations follow in a straightforward manner.
+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
 
-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}}, \\
+    \times (\frac{\partial V}{\partial {\bf u}})_{u(t+h)}, \\
 %
-{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
+{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(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
+${\bf u}$ is automatically updated when the rotation matrix
+$\mathsf{Q}$ 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.
 
@@ -132,14 +130,24 @@ be advanced to the same time value.
 \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}.
+in Fig.~\ref{methodFig:timestep} where 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.
 
 \begin{figure}
 \centering
@@ -154,66 +162,39 @@ In Fig.~\ref{methodFig:timestep}, the resulting energy
 \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
+\dot{\mathsf{Q}} & = & \mathsf{Q} \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}
+rot}\left(\mathsf{Q}^{T} \cdot \frac{\partial V}{\partial
+\mathsf{Q}} \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
+where $\tau_T$ is the time constant for relaxation of the
+temperature to the target value, and the instantaneous temperature
+$T$ is given by
 \begin{equation}
-T = \frac{2 K}{f k_B}
+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}
+where $N_{\mathrm{orient}}$ is the number of molecules with
+orientational degrees of freedom. The integration of the equations of motion
+is carried out in a velocity-Verlet style 2 part algorithm:
 
-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\} ,\\
@@ -229,7 +210,7 @@ T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\b
     + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
     \chi(t) \right) ,\\
 %
-\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}
+\mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}
     \left(h * {\bf j}(t + h / 2)
     \overleftrightarrow{\mathsf{I}}^{-1} \right) ,\\
 %
@@ -237,21 +218,18 @@ T(t) &\leftarrow \left\{{\bf v}(t)\right\}, \left\{{\b
     + \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
+\overleftrightarrow{\mathsf{I}}^{-1})$ is the same symplectic Strang
+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{Q}$ 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:}
@@ -273,42 +251,37 @@ T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
     \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}
+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.
+$H_{\mathrm{NVT}}$, so the conserved quantity should be checked
+periodically to verify the quality of the integration.
 
 \subsection{\label{methodSection:NPTi}Constant-pressure integration with
-isotropic box deformations (NPTi)}
+isotropic box (NPTi)}
 
-Isobaric-isothermal ensemble integrator is implemented using the
-Melchionna modifications to the Nos\'e-Hoover-Andersen equations of
-motion,\cite{Melchionna1993}
-
+We can used an isobaric-isothermal ensemble integrator which 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
+\dot{\mathsf{Q}} & = & \mathsf{Q} \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}, \\
+rot}\left(\mathsf{Q}^{T} \cdot \frac{\partial
+V}{\partial \mathsf{Q}} \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
@@ -316,7 +289,6 @@ P_{\mathrm{target}} \right), \\
 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
@@ -349,10 +321,9 @@ P(t) = \frac{1}{3} \mathrm{Tr} \left(
 the Pressure tensor,
 \begin{equation}
 P(t) = \frac{1}{3} \mathrm{Tr} \left(
-\overleftrightarrow{\mathsf{P}}(t). \right)
+\overleftrightarrow{\mathsf{P}}(t) \right) .
 \end{equation}
-
-In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for
+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:
@@ -371,7 +342,7 @@ P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\b
     + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
     \chi(t) \right), \\
 %
-\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
+\mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left(h *
     {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1}
     \right) ,\\
 %
@@ -391,7 +362,6 @@ P(t) &\leftarrow \left\{{\bf r}(t)\right\}, \left\{{\b
 \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
@@ -403,7 +373,6 @@ box by
 \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
@@ -435,7 +404,6 @@ P(t + h) &\leftarrow  \left\{{\bf r}(t + h)\right\},
     \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$.
@@ -450,12 +418,9 @@ Poor choices of $\delta t$, $\tau_T$, or $\tau_B$ can 
 \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),
+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) +
@@ -476,12 +441,12 @@ method are
 \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
+\dot{\mathsf{Q}} & = & \mathsf{Q} \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} ,\\
+rot}\left(\mathsf{Q}^{T} \cdot \frac{\partial
+V}{\partial \mathsf{Q}} \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
@@ -514,7 +479,7 @@ r}(t)\right\},
     + \frac{h}{2} \left( {\bf \tau}^b(t) - {\bf j}(t)
     \chi(t) \right), \\
 %
-\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left(h *
+\mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left(h *
     {\bf j}(t + h / 2) \overleftrightarrow{\mathsf{I}}^{-1}
     \right), \\
 %
@@ -536,10 +501,9 @@ exponential operation is used to scale the simulation 
     \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:
+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*}
@@ -570,20 +534,17 @@ T(t + h) &\leftarrow \left\{{\bf v}(t + h)\right\},
     + 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:
+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}
+ & & + 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
@@ -592,16 +553,17 @@ assume non-orthorhombic geometries.
 finds most use in simulating crystals or liquid crystals which
 assume non-orthorhombic geometries.
 
-\subsubsection{\label{methodSection:NPAT}NPAT Ensemble}
+\subsection{\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
+A comprehensive understanding of relations between structures and
+functions in biological membrane system ultimately relies on
+structure and dynamics of lipid bilayers, which are strongly
+affected by the interfacial interaction between lipid molecules and
+surrounding media. One quantity used to describe the interfacial
+interaction is the average surface area per lipid.
+Constant area and constant lateral pressure simulations 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}
@@ -620,16 +582,13 @@ minimum with respect to surface area $A$
 
 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
+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 membrane simulations, a special
+ensemble, NP$\gamma$T, has been proposed to maintain the lateral
+surface tension and normal pressure. The equation of motion for the
+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$}\\
@@ -645,16 +604,15 @@ - P_{{\rm{target}}} )
 - 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.
+$\gamma$ is set to zero, and if $x$ and $y$ can move independently.
 
-\section{\label{methodSection:zcons}Z-Constraint Method}
+\section{\label{methodSection:zcons}The Z-Constraint Method}
 
 Based on the fluctuation-dissipation theorem, a force
 auto-correlation method was developed by Roux and Karplus to
@@ -668,7 +626,6 @@ 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}
@@ -681,7 +638,6 @@ D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B
 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
@@ -690,9 +646,8 @@ each time step instead of resetting the coordinate.
 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
+each time step instead of resetting the coordinate. After the force
+calculation, we define $G_\alpha$ as
 \begin{equation}
 G_{\alpha} = \sum_i F_{\alpha i}, \label{oopseEq:zc1}
 \end{equation}
@@ -727,7 +682,6 @@ v_{\beta i} = v_{\beta i} + \sum_{\alpha}
 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