trunk/tengDissertation/Methodology.tex

\chapter{\label{chapt:methodology}METHODOLOGY}

\section{\label{methodSection:rigidBodyIntegrators}Integrators for Rigid Body Motion in Molecular Dynamics}

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 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 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 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 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}

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^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, 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{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}
\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{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{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}
\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}})_{u(t+h)}, \\
%
{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h)
    \cdot {\bf \tau}^s(t + h).
\end{align*}
${\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.

{\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} where the resulting energy drifts at
various time steps for both the DLM and quaternion integration
schemes are 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
\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}

\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{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{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}
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}.
\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}
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:

{\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{Q}(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 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:}
\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 should be checked
periodically to verify the quality of the integration.

\subsection{\label{methodSection:NPTi}Constant-pressure integration with
isotropic box (NPTi)}

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{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{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
\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{Q}(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}
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{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{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
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{Q}(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.

\subsection{\label{methodSection:NPAT}NPAT Ensemble}

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}
                  \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$,
\begin{equation}
\gamma  = \frac{{\partial G}}{{\partial A}}=0.
\end{equation}
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$}\\
    \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, and if $x$ and $y$ can move independently.

\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
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, we 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}
Revision:	2941
Committed:	Mon Jul 17 20:01:05 2006 UTC (17 years, 11 months ago) by tim
Content type:	application/x-tex
File size:	30077 byte(s)
Log Message:	references corrections