--- trunk/tengDissertation/Methodology.tex	2006/06/28 17:36:32	2904
+++ trunk/tengDissertation/Methodology.tex	2006/06/28 19:09:25	2905
@@ -48,7 +48,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$:
 
@@ -66,7 +65,6 @@ velocity-Verlet style 2-part algorithm, where $h= \del
 \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}
@@ -101,11 +99,10 @@ 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
@@ -117,7 +114,6 @@ torques are calculated at the new positions and orient
 {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
     \cdot {\bf \tau}^s(t + h).
 \end{align*}
-
 ${\bf u}$ is 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
@@ -133,14 +129,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
@@ -155,19 +161,6 @@ 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}
@@ -181,19 +174,17 @@ rot}\left(\mathsf{A}^{T} \cdot \frac{\partial V}{\part
 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}
+T = \frac{2 K}{f k_B}.
 \end{equation}
 Here, $f$ is the total number of degrees of freedom in the system,
 \begin{equation}
@@ -207,8 +198,7 @@ K = \sum_{i=1}^{N} \frac{1}{2} m_i {\bf v}_i^T \cdot {
 \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
+In Eq.~\ref{eq:nosehooverext}, $\tau_T$ is the time constant for
 relaxation of the temperature to the target value. The integration
 of the equations of motion is carried out in a velocity-Verlet style
 2 part algorithm:
@@ -236,7 +226,6 @@ 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
@@ -247,12 +236,10 @@ integrator.
 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.
+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.
 
-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\},
@@ -272,16 +259,13 @@ 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)
@@ -297,8 +281,7 @@ Nos\'e-Hoover-Andersen equations of motion,\cite{Melch
 
 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}
-
+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}, \\
@@ -315,7 +298,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
@@ -348,10 +330,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:
@@ -390,7 +371,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
@@ -402,7 +382,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
@@ -434,7 +413,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$.
@@ -535,11 +513,10 @@ 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.
+exponential operation is used to scale the simulation box. The {\tt
+moveB} portion of the algorithm is largely unchanged from the NPTi
+integrator:
 
-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\},
@@ -569,11 +546,9 @@ 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)
@@ -582,7 +557,6 @@ T_{\mathrm{target}}}{2}
 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
@@ -620,16 +594,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, 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
+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, 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,7 +616,6 @@ - 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
@@ -668,7 +638,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 +650,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 +658,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, we 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 +694,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