--- trunk/tengDissertation/Methodology.tex	2006/06/06 14:56:36	2801
+++ trunk/tengDissertation/Methodology.tex	2006/06/11 02:53:15	2854
@@ -16,16 +16,16 @@ two decades. Matubayasi and Nakahara developed a time-
 
 Integration schemes for rotational motion of the rigid molecules in
 microcanonical ensemble have been extensively studied in the last
-two decades. Matubayasi and Nakahara 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
-time-reversible integrator of Matubayasi and Nakahara. However, all
-of the integrators in quaternion representation suffer from the
+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
@@ -117,7 +117,7 @@ torques are calculated at the new positions and orient
     \cdot {\bf \tau}^s(t + h).
 \end{align*}
 
-{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix
+${\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.
@@ -139,7 +139,7 @@ in Fig.~\ref{timestep}.
 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{timestep}.
+in Fig.~\ref{methodFig:timestep}.
 
 \begin{figure}
 \centering
@@ -278,10 +278,7 @@ self-consistent.  The relative tolerance for the self-
 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 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.
+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
@@ -299,9 +296,9 @@ To carry out isobaric-isothermal ensemble calculations
 \subsection{\label{methodSection: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{Melchionna1993}
+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), \\
@@ -356,11 +353,7 @@ relaxation of the pressure to the target value.  To se
 \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
+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:
 
@@ -401,12 +394,11 @@ depends on the positions at the same time.  {\sc oopse
 
 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.  {\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 + h / 2$.  Reshaping the
-box uniformly also scales the volume of the box by
+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)
@@ -448,10 +440,7 @@ + h)$ and $\eta(t + h)$ become self-consistent.  The r
 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 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.
++ 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,
@@ -473,10 +462,6 @@ Bond constraints are applied at the end of both the {\
 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}.
-
 \subsection{\label{methodSection:NPTf}Constant-pressure integration with a
 flexible box (NPTf)}
 
@@ -550,8 +535,8 @@ r}(t)\right\},
 \mathsf{H}(t + h) &\leftarrow \mathsf{H}(t) \cdot e^{-h
     \overleftrightarrow{\eta}(t + h / 2)} .
 \end{align*}
-{\sc oopse} uses a power series expansion truncated at second order
-for the exponential operation which scales the simulation box.
+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:
@@ -590,13 +575,14 @@ The NPTf integrator is known to conserve the following
 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}
+\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{equation}
+\end{eqnarray*}
 
 This integrator must be used with care, particularly in liquid
 simulations.  Liquids have very small restoring forces in the
@@ -606,8 +592,6 @@ assume non-orthorhombic geometries.
 finds most use in simulating crystals or liquid crystals which
 assume non-orthorhombic geometries.
 
-\subsection{\label{methodSection:otherSpecialEnsembles}Other Special Ensembles}
-
 \subsubsection{\label{methodSection:NPAT}NPAT Ensemble}
 
 A comprehensive understanding of structure�Cfunction relations of
@@ -631,7 +615,8 @@ described for the NPTi integrator.
 Note that the iterative schemes for NPAT are identical to those
 described for the NPTi integrator.
 
-\subsubsection{\label{methodSection:NPrT}NP$\gamma$T Ensemble}
+\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
@@ -669,20 +654,92 @@ $\gamma$ is set to zero.
 integrator is a special case of $NP\gamma T$ if the surface tension
 $\gamma$ is set to zero.
 
-%\section{\label{methodSection:constraintMethod}Constraint Method}
+\section{\label{methodSection:zcons}Z-Constraint Method}
 
-%\subsection{\label{methodSection:bondConstraint}Bond Constraint for Rigid Body}
-
-%\subsection{\label{methodSection:zcons}Z-constraint Method}
-
-\section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies}
-
-\subsection{\label{methodSection:temperature}Temperature Control}
-
-\subsection{\label{methodSection:pressureControl}Pressure Control}
+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}
 
-\section{\label{methodSection:hydrodynamics}Hydrodynamics}
+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}
 
-%\section{\label{methodSection:coarseGrained}Coarse-Grained Modeling}
+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.
 
-%\section{\label{methodSection:moleculeScale}Molecular-Scale Modeling}
+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}
