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

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.

\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^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$:

{\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*}

{\sc oopse} automatically updates ${\bf u}$ 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 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.

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

To carry out isobaric-isothermal ensemble calculations {\sc oopse}
implements 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.  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
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.  {\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
\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 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.

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}

Bond constraints are applied at the end of both the {\tt moveA} and
{\tt moveB} portions of the algorithm.

\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*}
{\sc oopse} uses a power series expansion truncated at second order
for the exponential operation which scales 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{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}
\mathrm{Tr}\left[\overleftrightarrow{\eta}(t)\right]^2 \tau_B^2.
\end{equation}

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:otherSpecialEnsembles}Other Special Ensembles}

\subsubsection{\label{methodSection:NPAT}\textbf{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.

\subsubsection{\label{methodSection:NPrT}\textbf{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}


\section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies}

%\subsection{\label{methodSection:temperature}Temperature Control}

%\subsection{\label{methodSection:pressureControl}Pressure Control}

%\section{\label{methodSection:hydrodynamics}Hydrodynamics}

%applications of langevin dynamics
As an excellent alternative to newtonian dynamics, Langevin
dynamics, which mimics a simple heat bath with stochastic and
dissipative forces, has been applied in a variety of studies. The
stochastic treatment of the solvent enables us to carry out
substantially longer time simulation. Implicit solvent Langevin
dynamics simulation of met-enkephalin not only outperforms explicit
solvent simulation on computation efficiency, but also agrees very
well with explicit solvent simulation on dynamics
properties\cite{Shen2002}. Recently, applying Langevin dynamics with
UNRES model, Liow and his coworkers suggest that protein folding
pathways can be possibly exploited within a reasonable amount of
time\cite{Liwo2005}. The stochastic nature of the Langevin dynamics
also enhances the sampling of the system and increases the
probability of crossing energy barrier\cite{Banerjee2004, Cui2003}.
Combining Langevin dynamics with Kramers's theory, Klimov and
Thirumalai identified the free-energy barrier by studying the
viscosity dependence of the protein folding rates\cite{Klimov1997}.
In order to account for solvent induced interactions missing from
implicit solvent model, Kaya incorporated desolvation free energy
barrier into implicit coarse-grained solvent model in protein
folding/unfolding study and discovered a higher free energy barrier
between the native and denatured states. Because of its stability
against noise, Langevin dynamics is very suitable for studying
remagnetization processes in various
systems\cite{Palacios1998,Berkov2002,Denisov2003}. For instance, the
oscillation power spectrum of nanoparticles from Langevin dynamics
simulation has the same peak frequencies for different wave
vectors,which recovers the property of magnetic excitations in small
finite structures\cite{Berkov2005a}. In an attempt to reduce the
computational cost of simulation, multiple time stepping (MTS)
methods have been introduced and have been of great interest to
macromolecule and protein community\cite{Tuckerman1992}. Relying on
the observation that forces between distant atoms generally
demonstrate slower fluctuations than forces between close atoms, MTS
method are generally implemented by evaluating the slowly
fluctuating forces less frequently than the fast ones.
Unfortunately, nonlinear instability resulting from increasing
timestep in MTS simulation have became a critical obstruction
preventing the long time simulation. Due to the coupling to the heat
bath, Langevin dynamics has been shown to be able to damp out the
resonance artifact more efficiently\cite{Sandu1999}.

It is very important to develop stable and efficient methods to
integrate the equations of motion of orientational degrees of
freedom. Euler angles are the nature choice to describe the
rotational degrees of freedom. However, due to its singularity, the
numerical integration of corresponding equations of motion is very
inefficient and inaccurate. Although an alternative integrator using
different sets of Euler angles can overcome this
difficulty\cite{Ryckaert1977, Andersen1983}, the computational
penalty and the lost of angular momentum conservation still remain.
In 1977, a singularity free representation utilizing quaternions was
developed by Evans\cite{Evans1977}. Unfortunately, this approach
suffer from the nonseparable Hamiltonian resulted from quaternion
representation, which prevents the symplectic algorithm to be
utilized. Another different approach is to apply holonomic
constraints to the atoms belonging to the rigid
body\cite{Barojas1973}. Each atom moves independently under the
normal forces deriving from potential energy and constraint forces
which are used to guarantee the rigidness. However, due to their
iterative nature, SHAKE and Rattle algorithm converge very slowly
when the number of constraint increases.

The break through in geometric literature suggests that, in order to
develop a long-term integration scheme, one should preserve the
geometric structure of the flow. 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. 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
is RSHAKE\cite{Kol1997}, in which a conjugate momentum to rotation
matrix is introduced to re-formulate the Hamiltonian's equation and
the Hamiltonian is evolved in a constraint manifold by iteratively
satisfying the orthogonality constraint. However, RSHAKE is
inefficient because of the iterative procedure. An extremely
efficient integration scheme in rotation matrix representation,
which also preserves the same structural properties of the
Hamiltonian flow as Miller's integrator, is proposed by Dullweber,
Leimkuhler and McLachlan (DLM)\cite{Dullweber1997}.

%review langevin/browninan dynamics for arbitrarily shaped rigid body
Combining Langevin or Brownian dynamics with rigid body dynamics,
one can study the slow processes in biomolecular systems. Modeling
the DNA as a chain of rigid spheres beads, which subject to harmonic
potentials as well as excluded volume potentials, Mielke and his
coworkers discover rapid superhelical stress generations from the
stochastic simulation of twin supercoiling DNA with response to
induced torques\cite{Mielke2004}. Membrane fusion is another key
biological process which controls a variety of physiological
functions, such as release of neurotransmitters \textit{etc}. A
typical fusion event happens on the time scale of millisecond, which
is impracticable to study using all atomistic model with newtonian
mechanics. With the help of coarse-grained rigid body model and
stochastic dynamics, the fusion pathways were exploited by many
researchers\cite{Noguchi2001,Noguchi2002,Shillcock2005}. Due to the
difficulty of numerical integration of anisotropy rotation, most of
the rigid body models are simply modeled by sphere, cylinder,
ellipsoid or other regular shapes in stochastic simulations. In an
effort to account for the diffusion anisotropy of the arbitrary
particles, Fernandes and de la Torre improved the original Brownian
dynamics simulation algorithm\cite{Ermak1978,Allison1991} by
incorporating a generalized $6\times6$ diffusion tensor and
introducing a simple rotation evolution scheme consisting of three
consecutive rotations\cite{Fernandes2002}. Unfortunately, unexpected
error and bias are introduced into the system due to the arbitrary
order of applying the noncommuting rotation
operators\cite{Beard2003}. Based on the observation the momentum
relaxation time is much less than the time step, one may ignore the
inertia in Brownian dynamics. However, assumption of the zero
average acceleration is not always true for cooperative motion which
is common in protein motion. An inertial Brownian dynamics (IBD) was
proposed to address this issue by adding an inertial correction
term\cite{Beard2003}. As a complement to IBD which has a lower bound
in time step because of the inertial relaxation time, long-time-step
inertial dynamics (LTID) can be used to investigate the inertial
behavior of the polymer segments in low friction
regime\cite{Beard2003}. LTID can also deal with the rotational
dynamics for nonskew bodies without translation-rotation coupling by
separating the translation and rotation motion and taking advantage
of the analytical solution of hydrodynamics properties. However,
typical nonskew bodies like cylinder and ellipsoid are inadequate to
represent most complex macromolecule assemblies. These intricate
molecules have been represented by a set of beads and their
hydrodynamics properties can be calculated using variant
hydrodynamic interaction tensors.

The goal of the present work is to develop a Langevin dynamics
algorithm for arbitrary rigid particles by integrating the accurate
estimation of friction tensor from hydrodynamics theory into the
sophisticated rigid body dynamics.

\subsection{\label{introSection:frictionTensor} Friction Tensor}
Theoretically, the friction kernel can be determined using velocity
autocorrelation function. However, this approach become impractical
when the system become more and more complicate. Instead, various
approaches based on hydrodynamics have been developed to calculate
the friction coefficients. The friction effect is isotropic in
Equation, $\zeta$ can be taken as a scalar. In general, friction
tensor $\Xi$ is a $6\times 6$ matrix given by
\[
\Xi  = \left( {\begin{array}{*{20}c}
   {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
   {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
\end{array}} \right).
\]
Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are translational friction
tensor and rotational resistance (friction) tensor respectively,
while ${\Xi^{tr} }$ is translation-rotation coupling tensor and $
{\Xi^{rt} }$ is rotation-translation coupling tensor. When a
particle moves in a fluid, it may experience friction force or
torque along the opposite direction of the velocity or angular
velocity,
\[
\left( \begin{array}{l}
 F_R  \\
 \tau _R  \\
 \end{array} \right) =  - \left( {\begin{array}{*{20}c}
   {\Xi ^{tt} } & {\Xi ^{rt} }  \\
   {\Xi ^{tr} } & {\Xi ^{rr} }  \\
\end{array}} \right)\left( \begin{array}{l}
 v \\
 w \\
 \end{array} \right)
\]
where $F_r$ is the friction force and $\tau _R$ is the friction
toque.

\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}}

For a spherical particle, the translational and rotational friction
constant can be calculated from Stoke's law,
\[
\Xi ^{tt}  = \left( {\begin{array}{*{20}c}
   {6\pi \eta R} & 0 & 0  \\
   0 & {6\pi \eta R} & 0  \\
   0 & 0 & {6\pi \eta R}  \\
\end{array}} \right)
\]
and
\[
\Xi ^{rr}  = \left( {\begin{array}{*{20}c}
   {8\pi \eta R^3 } & 0 & 0  \\
   0 & {8\pi \eta R^3 } & 0  \\
   0 & 0 & {8\pi \eta R^3 }  \\
\end{array}} \right)
\]
where $\eta$ is the viscosity of the solvent and $R$ is the
hydrodynamics radius.

Other non-spherical shape, such as cylinder and ellipsoid
\textit{etc}, are widely used as reference for developing new
hydrodynamics theory, because their properties can be calculated
exactly. In 1936, Perrin extended Stokes's law to general ellipsoid,
also called a triaxial ellipsoid, which is given in Cartesian
coordinates by\cite{Perrin1934, Perrin1936}
\[
\frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
}} = 1
\]
where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
due to the complexity of the elliptic integral, only the ellipsoid
with the restriction of two axes having to be equal, \textit{i.e.}
prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
exactly. Introducing an elliptic integral parameter $S$ for prolate,
\[
S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
} }}{b},
\]
and oblate,
\[
S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
}}{a}
\],
one can write down the translational and rotational resistance
tensors
\[
\begin{array}{l}
 \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}} \\
 \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S + 2a}} \\
 \end{array},
\]
and
\[
\begin{array}{l}
 \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}} \\
 \Xi _b^{rr}  = \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}} \\
 \end{array}.
\]

\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shape}}

Unlike spherical and other regular shaped molecules, there is not
analytical solution for friction tensor of any arbitrary shaped
rigid molecules. The ellipsoid of revolution model and general
triaxial ellipsoid model have been used to approximate the
hydrodynamic properties of rigid bodies. However, since the mapping
from all possible ellipsoidal space, $r$-space, to all possible
combination of rotational diffusion coefficients, $D$-space is not
unique\cite{Wegener1979} as well as the intrinsic coupling between
translational and rotational motion of rigid body, general ellipsoid
is not always suitable for modeling arbitrarily shaped rigid
molecule. A number of studies have been devoted to determine the
friction tensor for irregularly shaped rigid bodies using more
advanced method where the molecule of interest was modeled by
combinations of spheres(beads)\cite{Carrasco1999} and the
hydrodynamics properties of the molecule can be calculated using the
hydrodynamic interaction tensor. Let us consider a rigid assembly of
$N$ beads immersed in a continuous medium. Due to hydrodynamics
interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different
than its unperturbed velocity $v_i$,
\[
v'_i  = v_i  - \sum\limits_{j \ne i} {T_{ij} F_j }
\]
where $F_i$ is the frictional force, and $T_{ij}$ is the
hydrodynamic interaction tensor. The friction force of $i$th bead is
proportional to its ``net'' velocity
\begin{equation}
F_i  = \zeta _i v_i  - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
\label{introEquation:tensorExpression}
\end{equation}
This equation is the basis for deriving the hydrodynamic tensor. In
1930, Oseen and Burgers gave a simple solution to Equation
\ref{introEquation:tensorExpression}
\begin{equation}
T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor}
\end{equation}
Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$.
A second order expression for element of different size was
introduced by Rotne and Prager\cite{Rotne1969} and improved by
Garc\'{i}a de la Torre and Bloomfield\cite{Torre1977},
\begin{equation}
T_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I +
\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma
_i^2  + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} -
\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
\label{introEquation:RPTensorNonOverlapped}
\end{equation}
Both of the Equation \ref{introEquation:oseenTensor} and Equation
\ref{introEquation:RPTensorNonOverlapped} have an assumption $R_{ij}
\ge \sigma _i  + \sigma _j$. An alternative expression for
overlapping beads with the same radius, $\sigma$, is given by
\begin{equation}
T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
\frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I +
\frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
\label{introEquation:RPTensorOverlapped}
\end{equation}

To calculate the resistance tensor at an arbitrary origin $O$, we
construct a $3N \times 3N$ matrix consisting of $N \times N$
$B_{ij}$ blocks
\begin{equation}
B = \left( {\begin{array}{*{20}c}
   {B_{11} } &  \ldots  & {B_{1N} }  \\
    \vdots  &  \ddots  &  \vdots   \\
   {B_{N1} } &  \cdots  & {B_{NN} }  \\
\end{array}} \right),
\end{equation}
where $B_{ij}$ is given by
\[
B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
)T_{ij}
\]
where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
$B$, we obtain

\[
C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
   {C_{11} } &  \ldots  & {C_{1N} }  \\
    \vdots  &  \ddots  &  \vdots   \\
   {C_{N1} } &  \cdots  & {C_{NN} }  \\
\end{array}} \right)
\]
, which can be partitioned into $N \times N$ $3 \times 3$ block
$C_{ij}$. With the help of $C_{ij}$ and skew matrix $U_i$
\[
U_i  = \left( {\begin{array}{*{20}c}
   0 & { - z_i } & {y_i }  \\
   {z_i } & 0 & { - x_i }  \\
   { - y_i } & {x_i } & 0  \\
\end{array}} \right)
\]
where $x_i$, $y_i$, $z_i$ are the components of the vector joining
bead $i$ and origin $O$. Hence, the elements of resistance tensor at
arbitrary origin $O$ can be written as
\begin{equation}
\begin{array}{l}
 \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\
 \Xi _{}^{tr}  = \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
 \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j  \\
 \end{array}
\label{introEquation:ResistanceTensorArbitraryOrigin}
\end{equation}

The resistance tensor depends on the origin to which they refer. The
proper location for applying friction force is the center of
resistance (reaction), at which the trace of rotational resistance
tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
resistance is defined as an unique point of the rigid body at which
the translation-rotation coupling tensor are symmetric,
\begin{equation}
\Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
\label{introEquation:definitionCR}
\end{equation}
Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
we can easily find out that the translational resistance tensor is
origin independent, while the rotational resistance tensor and
translation-rotation coupling resistance tensor depend on the
origin. Given resistance tensor at an arbitrary origin $O$, and a
vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
obtain the resistance tensor at $P$ by
\begin{equation}
\begin{array}{l}
 \Xi _P^{tt}  = \Xi _O^{tt}  \\
 \Xi _P^{tr}  = \Xi _P^{rt}  = \Xi _O^{tr}  - U_{OP} \Xi _O^{tt}  \\
 \Xi _P^{rr}  = \Xi _O^{rr}  - U_{OP} \Xi _O^{tt} U_{OP}  + \Xi _O^{tr} U_{OP}  - U_{OP} \Xi _O^{{tr} ^{^T }}  \\
 \end{array}
 \label{introEquation:resistanceTensorTransformation}
\end{equation}
where
\[
U_{OP}  = \left( {\begin{array}{*{20}c}
   0 & { - z_{OP} } & {y_{OP} }  \\
   {z_i } & 0 & { - x_{OP} }  \\
   { - y_{OP} } & {x_{OP} } & 0  \\
\end{array}} \right)
\]
Using Equations \ref{introEquation:definitionCR} and
\ref{introEquation:resistanceTensorTransformation}, one can locate
the position of center of resistance,
\begin{eqnarray*}
 \left( \begin{array}{l}
 x_{OR}  \\
 y_{OR}  \\
 z_{OR}  \\
 \end{array} \right) & = &\left( {\begin{array}{*{20}c}
   {(\Xi _O^{rr} )_{yy}  + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} }  \\
   { - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz}  + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} }  \\
   { - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx}  + (\Xi _O^{rr} )_{yy} }  \\
\end{array}} \right)^{ - 1}  \\
  & & \left( \begin{array}{l}
 (\Xi _O^{tr} )_{yz}  - (\Xi _O^{tr} )_{zy}  \\
 (\Xi _O^{tr} )_{zx}  - (\Xi _O^{tr} )_{xz}  \\
 (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
 \end{array} \right) \\
\end{eqnarray*}

where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
joining center of resistance $R$ and origin $O$.

\subsection{Langevin dynamics for rigid particles of arbitrary shape}

Consider a Langevin equation of motions in generalized coordinates
\begin{equation}
M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)}  + F_{r,i} (t)
\label{LDGeneralizedForm}
\end{equation}
where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
and moment of inertial) matrix and $V_i$ is a generalized velocity,
$V_i = V_i(v_i,\omega _i)$. The right side of Eq.
(\ref{LDGeneralizedForm}) consists of three generalized forces in
lab-fixed frame, systematic force $F_{s,i}$, dissipative force
$F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
system in Newtownian mechanics typically refers to lab-fixed frame,
it is also convenient to handle the rotation of rigid body in
body-fixed frame. Thus the friction and random forces are calculated
in body-fixed frame and converted back to lab-fixed frame by:
\[
\begin{array}{l}
 F_{f,i}^l (t) = A^T F_{f,i}^b (t), \\
 F_{r,i}^l (t) = A^T F_{r,i}^b (t) \\
 \end{array}.
\]
Here, the body-fixed friction force $F_{r,i}^b$ is proportional to
the body-fixed velocity at center of resistance $v_{R,i}^b$ and
angular velocity $\omega _i$,
\begin{equation}
F_{r,i}^b (t) = \left( \begin{array}{l}
 f_{r,i}^b (t) \\
 \tau _{r,i}^b (t) \\
 \end{array} \right) =  - \left( {\begin{array}{*{20}c}
   {\Xi _{R,t} } & {\Xi _{R,c}^T }  \\
   {\Xi _{R,c} } & {\Xi _{R,r} }  \\
\end{array}} \right)\left( \begin{array}{l}
 v_{R,i}^b (t) \\
 \omega _i (t) \\
 \end{array} \right),
\end{equation}
while the random force $F_{r,i}^l$ is a Gaussian stochastic variable
with zero mean and variance
\begin{equation}
\left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
\left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
2k_B T\Xi _R \delta (t - t').
\end{equation}
The equation of motion for $v_i$ can be written as
\begin{equation}
m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
f_{r,i}^l (t)
\end{equation}
Since the frictional force is applied at the center of resistance
which generally does not coincide with the center of mass, an extra
torque is exerted at the center of mass. Thus, the net body-fixed
frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is
given by
\begin{equation}
\tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b
\end{equation}
where $r_{MR}$ is the vector from the center of mass to the center
of the resistance. Instead of integrating angular velocity in
lab-fixed frame, we consider the equation of motion of angular
momentum in body-fixed frame
\begin{equation}
\dot \pi _i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b
(t) + \tau _{r,i}^b(t)
\end{equation}

Embedding the friction terms into force and torque, one can
integrate the langevin equations of motion for rigid body of
arbitrary shape in a velocity-Verlet style 2-part algorithm, where
$h= \delta t$:

{\tt part one:}
\begin{align*}
 v_i (t + h/2) &\leftarrow v_i (t) + \frac{{hf_{t,i}^l (t)}}{{2m_i }} \\
 \pi _i (t + h/2) &\leftarrow \pi _i (t) + \frac{{h\tau _{t,i}^b (t)}}{2} \\
 r_i (t + h) &\leftarrow r_i (t) + hv_i (t + h/2) \\
 A_i (t + h) &\leftarrow rotate\left( {h\pi _i (t + h/2)I^{ - 1} } \right) \\
\end{align*}
In this context, the $\mathrm{rotate}$ function is the reversible
product of five consecutive 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 ($\pi$) 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 friction and random
forces are generated at the center of resistance in body-fixed frame
and converted back into lab-fixed frame
\[
f_{t,i}^l (t + h) =  - \left( {\frac{{\partial V}}{{\partial r_i }}}
\right)_{r_i (t + h)}  + A_i^T (t + h)[F_{f,i}^b (t + h) + F_{r,i}^b
(t + h)],
\]
while the system torque in lab-fixed frame is transformed into
body-fixed frame,
\[
\tau _{t,i}^b (t + h) = A\tau _{s,i}^l (t) + \tau _{n,i}^b (t) +
\tau _{r,i}^b (t).
\]
Once the forces and torques have been obtained at the new time step,
the velocities can be advanced to the same time value.

{\tt part two:}
\begin{align*}
 v_i (t) &\leftarrow v_i (t + h/2) + \frac{{hf_{t,i}^l (t + h)}}{{2m_i }} \\
 \pi _i (t) &\leftarrow \pi _i (t + h/2) + \frac{{h\tau _{t,i}^b (t + h)}}{2} \\
\end{align*}

\subsection{Results and discussion}
Revision:	2850
Committed:	Sun Jun 11 01:55:48 2006 UTC (18 years ago) by tim
Content type:	application/x-tex
File size:	56405 byte(s)
Log Message:	move friction tensor section to chapter2