trunk/langevin/langevin.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
%\documentclass[aps,prb,preprint]{revtex4}
\documentclass[11pt]{article}
\usepackage{endfloat}
\usepackage{amsmath,bm}
\usepackage{amssymb}
\usepackage{times}
\usepackage{mathptmx}
\usepackage{setspace}
\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{bibentry}
\usepackage{mathrsfs}
\usepackage[ref]{overcite}
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
9.0in \textwidth 6.5in \brokenpenalty=10000
\renewcommand{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note

\begin{document}

\title{An algorithm for performing Langevin dynamics on rigid bodies of arbitrary shape }

\author{Teng Lin, Xiuquan Sun and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail:
gezelter@nd.edu} \\
Department of Chemistry and Biochemistry\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle \doublespacing

\begin{abstract}

\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

%applications of langevin dynamics
As 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
simulations. Implicit solvent Langevin dynamics simulations of
met-enkephalin not only outperform explicit solvent simulations for
computational efficiency, but also agrees very well with explicit
solvent simulations for dynamical properties.\cite{Shen2002}
Recently, applying Langevin dynamics with the UNRES model, Liow and
his coworkers suggest that protein folding pathways can be possibly
explored 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 barriers.\cite{Banerjee2004, Cui2003} Combining Langevin
dynamics with Kramers's theory, Klimov and Thirumalai identified
free-energy barriers 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 studies
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}

%review rigid body dynamics
Rigid bodies are frequently involved in the modeling of different
areas, from engineering, physics, to chemistry. For example,
missiles and vehicle are usually modeled by rigid bodies.  The
movement of the objects in 3D gaming engine or other physics
simulator is governed by the rigid body dynamics. In molecular
simulation, rigid body is used to simplify the model in
protein-protein docking study{\cite{Gray2003}}.

It is very important to develop stable and efficient methods to
integrate the equations of motion for orientational degrees of
freedom. Euler angles are the natural choice to describe the
rotational degrees of freedom. However, due to $\frac {1}{sin
\theta}$ singularities, the numerical integration of corresponding
equations of these motion is very inefficient and inaccurate.
Although an alternative integrator using multiple sets of Euler
angles can overcome this difficulty\cite{Barojas1973}, the
computational penalty and the loss of angular momentum conservation
still remain. A singularity-free representation utilizing
quaternions was developed by Evans in 1977.\cite{Evans1977}
Unfortunately, this approach used a nonseparable Hamiltonian
resulting from the quaternion representation, which prevented the
symplectic algorithm from being utilized. Another different approach
is to apply holonomic constraints to the atoms belonging to the
rigid body. 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,
the SHAKE and Rattle algorithms also converge very slowly when the
number of constraints increases.\cite{Ryckaert1977, Andersen1983}

A break-through in geometric literature suggests that, in order to
develop a long-term integration scheme, one should preserve the
symplectic structure of the propagator. By introducing a conjugate
momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's
equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was
proposed to evolve the Hamiltonian system in a constraint manifold
by iteratively satisfying the orthogonality constraint $Q^T Q = 1$.
An alternative method using the quaternion representation was
developed by Omelyan.\cite{Omelyan1998} However, both of these
methods are iterative and inefficient. In this section, we descibe a
symplectic Lie-Poisson integrator for rigid bodies developed by
Dullweber and his coworkers\cite{Dullweber1997} in depth.

%review langevin/browninan dynamics for arbitrarily shaped rigid body
Combining Langevin or Brownian dynamics with rigid body dynamics,
one can study slow processes in biomolecular systems. Modeling DNA
as a chain of rigid beads, which are subject to harmonic potentials
as well as excluded volume potentials, Mielke and his coworkers
discovered 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 a millisecond, which is impractical to
study using atomistic models with newtonian mechanics. With the help
of coarse-grained rigid body model and stochastic dynamics, the
fusion pathways were explored by many
researchers.\cite{Noguchi2001,Noguchi2002,Shillcock2005} Due to the
difficulty of numerical integration of anisotropic rotation, most of
the rigid body models are simply modeled using spheres, cylinders,
ellipsoids or other regular shapes in stochastic simulations. In an
effort to account for the diffusion anisotropy of 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
errors and biases 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, the assumption of 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{Beard2000} 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{Beard2000} 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 cylinders and ellipsoids are inadequate
to represent most complex macromolecule assemblies. These intricate
molecules have been represented by a set of beads and their
hydrodynamic properties can be calculated using variants on the
standard hydrodynamic interaction tensors.

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

\section{Computational Methods{\label{methodSec}}}

\subsection{\label{introSection:frictionTensor}Friction Tensor}
Theoretically, the friction kernel can be determined using the
velocity autocorrelation function. However, this approach becomes
impractical when the system becomes more and more complicated.
Instead, various approaches based on hydrodynamics have been
developed to calculate the friction coefficients. In general, the
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 $3 \times 3$
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
torque.

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

For a spherical particle with slip boundary conditions, 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
hydrodynamic radius.

Other non-spherical shapes, such as cylinders and ellipsoids, are
widely used as references for developing new hydrodynamics theory,
because their properties can be calculated exactly. In 1936, Perrin
extended Stokes's law to general ellipsoids, 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 being 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
ellipsoids :
\[
S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
} }}{b},
\]
and oblate ellipsoids:
\[
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{eqnarray*}
 \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{eqnarray*}
and
\begin{eqnarray*}
 \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{eqnarray*}

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

Unlike spherical and other simply shaped molecules, there is no
analytical solution for the friction tensor for arbitrarily 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 spaces, $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 bodies, general
ellipsoids are not always suitable for modeling arbitrarily shaped
rigid molecules. A number of studies have been devoted to
determining the friction tensor for irregularly shaped rigid bodies
using more advanced methods where the molecule of interest was
modeled by a combinations of spheres\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 hydrodynamic
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
Eq.~\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 Eq.~\ref{introEquation:oseenTensor} and
Eq.~\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 the Kronecker delta function. Inverting the
$B$ matrix, 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 the 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$, the elements of resistance tensor at
arbitrary origin $O$ can be written as
\begin{eqnarray}
 \Xi _{}^{tt}  & = & \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\
 \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. \notag \\
\label{introEquation:ResistanceTensorArbitraryOrigin}
\end{eqnarray}
The resistance tensor depends on the origin to which they refer. The
proper location for applying the friction force is the center of
resistance (or center of reaction), at which the trace of rotational
resistance tensor, $ \Xi ^{rr}$ reaches a minimum value.
Mathematically, the center of resistance is defined as an unique
point of the rigid body at which the translation-rotation coupling
tensors are symmetric,
\begin{equation}
\Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
\label{introEquation:definitionCR}
\end{equation}
From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
we can easily derive that the translational resistance tensor is
origin independent, while the rotational resistance tensor and
translation-rotation coupling resistance tensor depend on the
origin. Given the 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 Eq.~\ref{introEquation:definitionCR} and
Eq.~\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\label{LDRB}}

Consider the Langevin equations of motion 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) = Q^T F_{f,i}^b (t), \\
 F_{r,i}^l (t) = Q^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'). \label{randomForce}
\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 the angular velocity in
lab-fixed frame, we consider the equation of angular momentum in
body-fixed frame
\begin{equation}
\dot j_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 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}}, \\
%
{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h)
    \cdot {\bf \tau}^s(t + h).
\end{align*}
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*}

\section{Results and Discussion}

The Langevin algorithm described in previous section has been
implemented in {\sc oopse}\cite{Meineke2005} and applied to studies
of the static and dynamic properties in several systems.

\subsection{Temperature Control}

As shown in Eq.~\ref{randomForce}, random collisions associated with
the solvent's thermal motions is controlled by the external
temperature. The capability to maintain the temperature of the whole
system was usually used to measure the stability and efficiency of
the algorithm. In order to verify the stability of this new
algorithm, a series of simulations are performed on system
consisiting of 256 SSD water molecules with different viscosities.
The initial configuration for the simulations is taken from a 1ns
NVT simulation with a cubic box of 19.7166~\AA. All simulation are
carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns
with reference temperature at 300~K. The average temperature as a
function of $\eta$ is shown in Table \ref{langevin:viscosity} where
the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 -
1$ poise. The better temperature control at higher viscosity can be
explained by the finite size effect and relative slow relaxation
rate at lower viscosity regime.
\begin{table}
\caption{AVERAGE TEMPERATURES FROM LANGEVIN DYNAMICS SIMULATIONS OF
SSD WATER MOLECULES WITH REFERENCE TEMPERATURE AT 300~K.}
\label{langevin:viscosity}
\begin{center}
\begin{tabular}{lll}
  \hline
  $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\
  \hline
  1    & 300.47 & 10.99 \\
  0.1  & 301.19 & 11.136 \\
  0.01 & 303.04 & 11.796 \\
  \hline
\end{tabular}
\end{center}
\end{table}

Another set of calculations were performed to study the efficiency of
temperature control using different temperature coupling schemes.
The starting configuration is cooled to 173~K and evolved using NVE,
NVT, and Langevin dynamic with time step of 2 fs.
Fig.~\ref{langevin:temperature} shows the heating curve obtained as
the systems reach equilibrium. The orange curve in
Fig.~\ref{langevin:temperature} represents the simulation using
Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps
which gives reasonable tight coupling, while the blue one from
Langevin dynamics with viscosity of 0.1 poise demonstrates a faster
scaling to the desire temperature. When $ \eta = 0$, Langevin dynamics becomes normal
NVE (see orange curve in Fig.~\ref{langevin:temperature}) which
loses the temperature control ability.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{temperature.pdf}
\caption[Plot of Temperature Fluctuation Versus Time]{Plot of
temperature fluctuation versus time.} \label{langevin:temperature}
\end{figure}

\subsection{Langevin Dynamics simulation {\it vs} NVE simulations}

To validate our langevin dynamics simulation. We performed several NVE
simulations with explicit solvents for different shaped
molecules. There are one solute molecule and 1929 solvent molecules in
NVE simulation. The parameters are shown in table
\ref{tab:parameters}. The force field between spheres is standard
Lennard-Jones, and ellipsoids interact with other ellipsoids and
spheres with generalized Gay-Berne potential. All simulations are
carried out at 300 K and 1 Atm. The time step is 25 ns, and a
switching function was applied to all potentials to smoothly turn off
the interactions between a range of $22$ and $25$ \AA.  The switching
function was the standard (cubic) function,
\begin{equation}
s(r) = 
        \begin{cases}
        1 & \text{if $r \le r_{\text{sw}}$},\\
        \frac{(r_{\text{cut}} + 2r - 3r_{\text{sw}})(r_{\text{cut}} - r)^2}
        {(r_{\text{cut}} - r_{\text{sw}})^3} 
        & \text{if $r_{\text{sw}} < r \le r_{\text{cut}}$}, \\
        0 & \text{if $r > r_{\text{cut}}$.}
        \end{cases}
\label{eq:switchingFunc}
\end{equation}
We have computed translational diffusion constants for lipid molecules
from the mean-square displacement,
\begin{equation}
D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle,
\end{equation}
of the solute molecules. Translational diffusion constants for the
different shaped molecules are shown in table
\ref{tab:translation}.  We have also computed orientational correlation
times for different shaped molecules from fits of the second-order
Legendre polynomial correlation function,
\begin{equation}
C_{\ell}(t)  =  \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
\mu}_{i}(0) \right) 
\end{equation}
the results are shown in table \ref{tab:rotation}. We used einstein
format of the pressure correlation function,
\begin{equation}
C_{\ell}(t)  =  \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
\mu}_{i}(0) \right) 
\end{equation}
to estimate the viscosity of the systems from NVE simulations. The
viscosity can also be calculated by Green-Kubo pressure correlaton
function,
\begin{equation}
C_{\ell}(t)  =  \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
\mu}_{i}(0) \right) 
\end{equation}
However, this method converges slowly, and the statistics are not good
enough to give us a very accurate value. The langevin dynamics
simulations for different shaped molecules are performed at the same
conditions as the NVE simulations with viscosity estimated from NVE
simulations. To get better statistics, 1024 non-interacting solute
molecules are put into one simulation box for each langevin
simulation, this is equal to 1024 simulations for single solute
systems. The diffusion constants and rotation relaxation times for
different shaped molecules are shown in table \ref{tab:translation}
and \ref{tab:rotation} to compare to the results calculated from NVE
simulations. The theoretical values for sphere is calculated from the
Stokes-Einstein law, the theoretical values for ellipsoid is
calculated from Perrin's fomula, the theoretical values for dumbbell
is calculated from StinXX and Davis theory. The exact method is
applied to the langevin dynamics simulations for sphere and ellipsoid,
the bead model is applied to the simulation for dumbbell molecule, and
the rough shell model is applied to ellipsoid, dumbbell, banana and
lipid molecules. The results from all the langevin dynamics
simulations, including exact, bead model and rough shell, match the
theoretical values perfectly for all different shaped molecules. This
indicates that our simulation package for langevin dynamics is working
well. The approxiate methods ( bead model and rough shell model) are
accurate enough for the current simulations. The goal of the langevin
dynamics theory is to replace the explicit solvents by the friction
forces. We compared the dynamic properties of different shaped
molecules in langevin dynamics simulations with that in NVE
simulations. The results are reasonable close. Overall, the
translational diffusion constants calculated from langevin dynamics
simulations are very close to the values from the NVE simulation. For
sphere and lipid molecules, the diffusion constants are a little bit
off from the NVE simulation results. One possible reason is that the
calculation of the viscosity is very difficult to be accurate. Another
possible reason is that although we save very frequently during the
NVE simulations and run pretty long time simulations, there is only
one solute molecule in the system which makes the calculation for the
diffusion constant difficult. The sphere molecule behaves as a free
rotor in the solvent, so there is no rotation relaxation time
calculated from NVE simulations. The rotation relaxation time is not
very close to the NVE simulations results. The banana and lipid
molecules match the NVE simulations results pretty well. The mismatch
between langevin dynamics and NVE simulation for ellipsoid is possibly
caused by the slip boundary condition. For dumbbell, the mismatch is
caused by the size of the solvent molecule is pretty large compared to
dumbbell molecule in NVE simulations.

According to our simulations, the langevin dynamics is a reliable
theory to apply to replace the explicit solvents, especially for the
translation properties. For large molecules, the rotation properties
are also mimiced reasonablly well.

\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{}
\begin{tabular}{llccccccc}
\hline 
  & & Sphere & Ellipsoid & Dumbbell(2 spheres) & Banana(3 ellpsoids) &
Lipid(head) & lipid(tail) & Solvent \\
\hline 
$d$ (\AA) & & 6.5 & 4.6  & 6.5 &  4.2 & 6.5 & 4.6 & 4.7 \\
$l$ (\AA) & & $= d$ & 13.8 & $=d$ & 11.2 & $=d$ & 13.8 & 4.7 \\
$\epsilon^s$ (kcal/mol) & & 0.8 & 0.8 & 0.8 & 0.8 & 0.185 & 0.8 & 0.8 \\
$\epsilon_r$ (well-depth aspect ratio)& & 1 & 0.2 & 1 & 0.2 & 1 & 0.2 & 1 \\
$m$ (amu) & & 190 & 200 & 190 & 240 & 196 & 760 & 72.06 \\
%$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$) & & & & \\
%\multicolumn{2}c{$I_{xx}$} & 1125 & 45000 & N/A \\
%\multicolumn{2}c{$I_{yy}$} & 1125 & 45000 & N/A \\
%\multicolumn{2}c{$I_{zz}$} &  0 &    9000 & N/A \\
%$\mu$ (Debye) & & varied & 0 & 0 \\
\end{tabular}
\label{tab:parameters}
\end{center}
\end{minipage}
\end{table*}

\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{}
\begin{tabular}{lccccc}
\hline
 & & & & &Translation \\
\hline 
 & NVE &  & Theoretical & Langevin & \\
\hline 
 & $\eta$ & D & D & method & D \\
\hline
sphere & 3.480159e-03 & 1.643135e-04 & 1.942779e-04 & exact & 1.982283e-04 \\
ellipsoid & 2.551262e-03 & 2.437492e-04 & 2.335756e-04 & exact & 2.374905e-04 \\
 & 2.551262e-03  & 2.437492e-04 & 2.335756e-04 & rough shell & 2.284088e-04 \\
dumbell & 2.41276e-03  & 2.129432e-04 & 2.090239e-04 & bead model & 2.148098e-04 \\
 & 2.41276e-03 & 2.129432e-04 & 2.090239e-04 & rough shell & 2.013219e-04 \\
banana & 2.9846e-03 & 1.527819e-04 &  & rough shell & 1.54807e-04 \\
lipid & 3.488661e-03 & 0.9562979e-04 &  & rough shell & 1.320987e-04 \\
\end{tabular}
\label{tab:translation}
\end{center}
\end{minipage}
\end{table*}

\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{}
\begin{tabular}{lccccc}
\hline
 & & & & &Rotation \\
\hline 
 & NVE &  & Theoretical & Langevin & \\
\hline 
 & $\eta$ & $\tau_0$ & $\tau_0$ & method & $\tau_0$ \\
\hline
sphere & 3.480159e-03 &  & 1.208178e+04 & exact & 1.20628e+04 \\
ellipsoid & 2.551262e-03 & 4.66806e+04 & 2.198986e+04 & exact & 2.21507e+04 \\
 & 2.551262e-03 & 4.66806e+04 & 2.198986e+04 & rough shell & 2.21714e+04 \\
dumbell & 2.41276e-03 & 1.42974e+04 &  & bead model & 7.12435e+04 \\
 & 2.41276e-03 & 1.42974e+04 &  & rough shell & 7.04765e+04 \\
banana & 2.9846e-03 & 6.38323e+04 &  & rough shell & 7.0945e+04 \\
lipid & 3.488661e-03 & 7.79595e+04 &  & rough shell & 7.78886e+04 \\
\end{tabular}
\label{tab:rotation}
\end{center}
\end{minipage}
\end{table*}

Langevin dynamics simulations are applied to study the formation of
the ripple phase of lipid membranes. The initial configuration is
taken from our molecular dynamics studies on lipid bilayers with
lennard-Jones sphere solvents. The solvent molecules are excluded from
the system, the experimental value of water viscosity is applied to
mimic the heat bath. Fig. XXX is the snapshot of the stable
configuration of the system, the ripple structure stayed stable after
100 ns run. The efficiency of the simulation is increased by one order
of magnitude.

\subsection{Langevin Dynamics of Banana Shaped Molecules}

In order to verify that Langevin dynamics can mimic the dynamics of
the systems absent of explicit solvents, we carried out two sets of
simulations and compare their dynamic properties.
Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation
made of 256 pentane molecules and two banana shaped molecules at
273~K. It has an equivalent implicit solvent system containing only
two banana shaped molecules with viscosity of 0.289 center poise. To
calculate the hydrodynamic properties of the banana shaped molecule,
we created a rough shell model (see Fig.~\ref{langevin:roughShell}),
in which the banana shaped molecule is represented as a ``shell''
made of 2266 small identical beads with size of 0.3 \AA on the
surface. Applying the procedure described in
Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
identified the center of resistance at (0 $\rm{\AA}$, 0.7482 $\rm{\AA}$,
-0.1988 $\rm{\AA}$), as well as the resistance tensor,
\[
\left( {\begin{array}{*{20}c}
0.9261 & 0 & 0&0&0.08585&0.2057\\
0& 0.9270&-0.007063& 0.08585&0&0\\
0&-0.007063&0.7494&0.2057&0&0\\
0&0.0858&0.2057& 58.64& 0&0\\
0.08585&0&0&0&48.30&3.219&\\
0.2057&0&0&0&3.219&10.7373\\
\end{array}} \right).
\]
where the units for translational, translation-rotation coupling and rotational tensors are $\frac{kcal \cdot fs}{mol \cdot \rm{\AA}^2}$, $\frac{kcal \cdot fs}{mol \cdot \rm{\AA} \cdot rad}$ and $\frac{kcal \cdot fs}{mol \cdot rad^2}$ respectively.
Curves of the velocity auto-correlation functions in
Fig.~\ref{langevin:vacf} were shown to match each other very well.
However, because of the stochastic nature, simulation using Langevin
dynamics was shown to decay slightly faster than MD. In order to
study the rotational motion of the molecules, we also calculated the
auto-correlation function of the principle axis of the second GB
particle, $u$. The discrepancy shown in Fig.~\ref{langevin:uacf} was
probably due to the reason that we used the experimental viscosity directly instead of calculating bulk viscosity from simulation.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{roughShell.pdf}
\caption[Rough shell model for banana shaped molecule]{Rough shell
model for banana shaped molecule.} \label{langevin:roughShell}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{twoBanana.pdf}
\caption[Snapshot from Simulation of Two Banana Shaped Molecules and
256 Pentane Molecules]{Snapshot from simulation of two Banana shaped
molecules and 256 pentane molecules.} \label{langevin:twoBanana}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{vacf.pdf}
\caption[Plots of Velocity Auto-correlation Functions]{Velocity
auto-correlation functions of NVE (explicit solvent) in blue and
Langevin dynamics (implicit solvent) in red.} \label{langevin:vacf}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{uacf.pdf}
\caption[Auto-correlation functions of the principle axis of the
middle GB particle]{Auto-correlation functions of the principle axis
of the middle GB particle of NVE (blue) and Langevin dynamics
(red).} \label{langevin:uacf}
\end{figure}

\section{Conclusions}

We have presented a new Langevin algorithm by incorporating the
hydrodynamics properties of arbitrary shaped molecules into an
advanced symplectic integration scheme. The temperature control
ability of this algorithm was demonstrated by a set of simulations
with different viscosities. It was also shown to have significant
advantage of producing rapid thermal equilibration over
Nos\'{e}-Hoover method. Further studies in systems involving banana
shaped molecules illustrated that the dynamic properties could be
preserved by using this new algorithm as an implicit solvent model.


\section{Acknowledgments}
Support for this project was provided by the National Science
Foundation under grant CHE-0134881. T.L. also acknowledges the
financial support from center of applied mathematics at University
of Notre Dame.
\newpage

\bibliographystyle{jcp2}
\bibliography{langevin}

\end{document}
Revision:	3298
Committed:	Wed Jan 2 21:06:31 2008 UTC (16 years, 7 months ago) by xsun
Content type:	application/x-tex
File size:	40320 byte(s)
Log Message:	Add the text about langevin simulation on different shaped molecules.