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 referenc note

\begin{document}

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

\author{Xiuquan Sun, Teng Lin 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
Langevin dynamics, which mimics a simple heat bath with stochastic and
dissipative forces, has been applied in a variety of situations as an
alternative to molecular dynamics with explicit solvent molecules.
The stochastic treatment of the solvent allows the use of simulations
with substantially longer time and length scales.  In general, the
dynamic and structural properties obtained from Langevin simulations
agree quite well with similar properties obtained from explicit
solvent simulations.

Recent examples of the usefulness of Langevin simulations include a
study of met-enkephalin in which Langevin simulations predicted
dynamical properties that were largely in agreement with explicit
solvent simulations.\cite{Shen2002} By applying Langevin dynamics with
the UNRES model, Liow and his coworkers suggest that protein folding
pathways can be explored within a reasonable amount of
time.\cite{Liwo2005}

The stochastic nature of Langevin dynamics also enhances the sampling
of the system and increases the probability of crossing energy
barriers.\cite{Cui2003,Banerjee2004} 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 the implicit solvent model, Kaya
incorporated a desolvation free energy barrier into protein
folding/unfolding studies and discovered a higher free energy barrier
between the native and denatured states.\cite{HuseyinKaya07012005}

Because of its stability against noise, Langevin dynamics has also
proven useful for studying remagnetization processes in various
systems.\cite{Palacios1998,Berkov2002,Denisov2003} [Check: For
instance, the oscillation power spectrum of nanoparticles from
Langevin dynamics has the same peak frequencies for different wave
vectors, which recovers the property of magnetic excitations in small
finite structures.\cite{Berkov2005a}]

In typical LD simulations, the friction and random forces on
individual atoms are taken from the Stokes-Einstein hydrodynamic
approximation,
\begin{eqnarray}
m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + R(t) \\
\langle R(t) \rangle & = & 0 \\
\langle R(t) R(t') \rangle & = & 2 k_B T \xi m \delta(t - t')
\end{eqnarray}
where $\xi \approx 6 \pi \eta a$.  Here $\eta$ is the viscosity of the
implicit solvent, and $a$ is the hydrodynamic radius of the atom.

The use of rigid substructures,\cite{???} 
coarse-graining,\cite{Ayton,Sun,Zannoni} and ellipsoidal
representations of protein side chains~\cite{Schulten} has made the
use of the Stokes-Einstein approximation problematic.  A rigid
substructure moves as a single unit with orientational as well as
translational degrees of freedom.  This requires a more general
treatment of the hydrodynamics than the spherical approximation
provides.  The atoms involved in a rigid or coarse-grained structure
should properly have solvent-mediated interactions with each
other. The theory of interactions {\it between} bodies moving through
a fluid has been developed over the past century and has been applied
to simulations of Brownian
motion.\cite{MarshallNewton,GarciaDeLaTorre} There a need to have a
more thorough treatment of hydrodynamics included in the library of
methods available for performing Langevin simulations.

\subsection{Rigid Body Dynamics}
Rigid bodies are frequently involved in the modeling of large
collections of particles that move as a single unit.  In molecular
simulations, rigid bodies have been used to simplify protein-protein
docking,\cite{Gray2003} and lipid bilayer simulations.\cite{Sun2008}
Many of the water models in common use are also rigid-body
models,\cite{TIPs,SPC/E} although they are typically evolved using
constraints rather than rigid body equations of motion.

Euler angles are a natural choice to describe the rotational
degrees of freedom.  However, due to $1 \over \sin \theta$
singularities, the numerical integration of corresponding equations of
these motion can become inaccurate (and inefficient).  Although an
alternative integrator using multiple sets of Euler angles can
overcome this problem,\cite{Barojas1973} the computational penalty and
the loss of angular momentum conservation remain. A singularity-free
representation utilizing quaternions was developed by Evans in
1977.\cite{Evans1977} Unfortunately, this approach uses a nonseparable
Hamiltonian resulting from the quaternion representation, which
prevented symplectic algorithms from being utilized until very
recently.\cite{Miller2002} Another approach is the application of
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 breakthrough 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 suffer from some related inefficiencies. A symplectic
Lie-Poisson integrator for rigid bodies developed by Dullweber {\it et
al.}\cite{Dullweber1997} gets around most of the limitations mentioned
above and has become the basis for our Langevin integrator.


\subsection{The Hydrodynamic tensor and Brownian dynamics}
Combining Brownian dynamics with rigid substructures, one can study
slow processes in biomolecular systems.  Modeling DNA as a chain of
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 Noguchi and others.\cite{Noguchi2001,Noguchi2002,Shillcock2005}

Due to the difficulty of numerically integrating anisotropic
rotational motion, most of the coarse-grained rigid body models are
treated as spheres, cylinders, ellipsoids or other regular shapes in
stochastic simulations.  In an effort to account for the diffusion
anisotropy of arbitrarily-shaped particles, Fernandes and Garc\'{i}a
de la Torre improved the original Brownian dynamics simulation
algorithm~\cite{Ermak1978,Allison1991} by incorporating a generalized
$6\times6$ diffusion tensor and introducing a rotational evolution
scheme consisting of three consecutive rotations.\cite{Fernandes2002}
Unfortunately, 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 proteins. 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 macromolecular assemblies.

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
a symplectic rigid body dynamics propagator.  In the sections below,
we review some of the theory of hydrodynamic tensors developed for
Brownian simulations of rigid multi-particle systems, we then present
our integration method for a set of generalized Langevin equations of
motion, and we compare the behavior of the new Langevin integrator to
dynamical quantities obtained via explicit solvent molecular dynamics.

\subsection{\label{introSection:frictionTensor}The Friction Tensor}
Theoretically, a complete friction kernel can be determined using the
velocity autocorrelation function. However, this approach becomes
impractical when the solute becomes complex.  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
\begin{equation}
\Xi  = \left( {\begin{array}{*{20}c}
   {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
   {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
\end{array}} \right).
\end{equation}
Here, $\Xi^{tt}$ and $\Xi^{rr}$ are $3 \times 3$ translational and
rotational resistance (friction) tensors 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 ($\mathbf{F}_f$) and torque
($\mathbf{\tau}_f$) in opposition to the directions of the velocity
($\mathbf{v}$) and body-fixed angular velocity ($\mathbf{\omega}$),
\begin{equation}
\left( \begin{array}{l}
 \mathbf{F}_f  \\
 \mathbf{\tau}_f  \\
 \end{array} \right) =  - \left( {\begin{array}{*{20}c}
   \Xi ^{tt} & \Xi ^{rt}  \\
   \Xi ^{tr} & \Xi ^{rr}  \\
\end{array}} \right)\left( \begin{array}{l}
 \mathbf{v} \\
 \mathbf{\omega} \\
 \end{array} \right).
\end{equation}

\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}}
For a spherical particle under ``stick'' boundary conditions, the
translational and rotational friction tensors can be calculated from
Stoke's law,
\begin{equation}
\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)
\end{equation}
and
\begin{equation}
\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)
\end{equation}
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 theories,
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}
\begin{equation}
\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1
\end{equation}
where the semi-axes are of lengths $a$, $b$, and $c$. Due to the
complexity of the elliptic integral, only uniaxial ellipsoids,
{\it 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 :
\begin{equation}
S = \frac{2}{\sqrt{a^2  - b^2}} \ln \frac{a + \sqrt{a^2  - b^2}}{b},
\end{equation}
and oblate ellipsoids:
\begin{equation}
S = \frac{2}{\sqrt {b^2  - a^2 }} \arctan \frac{\sqrt {b^2  - a^2}}{a},
\end{equation}
one can write down the translational and rotational resistance
tensors for oblate,
\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 prolate,
\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*}
ellipsoids. For both spherical and ellipsoidal particles, the
translation-rotation and rotation-translation coupling tensors are
zero.

\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, the mapping from all
possible ellipsoidal spaces, $r$-space, to all possible combination of
rotational diffusion coefficients, $D$-space, is not
unique.\cite{Wegener1979} Additionally, because there is 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  + 6 \eta V {\bf I}. \notag
\label{introEquation:ResistanceTensorArbitraryOrigin}
\end{eqnarray}
The final term in the expression for $\Xi^{rr}$ is correction that
accounts for errors in the rotational motion of certain kinds of bead
models. The additive correction uses the solvent viscosity ($\eta$)
as well as the total volume of the beads that contribute to the
hydrodynamic model,
\begin{equation}
V = \frac{4 \pi}{3} \sum_{i=1}^{N} \sigma_i^3,
\end{equation}
where $\sigma_i$ is the radius of bead $i$.  This correction term was
rigorously tested and compared with the analytical results for
two-sphere and ellipsoidal systems by Garcia de la Torre and
Rodes.\cite{Torre:1983lr}


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


\section{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{Validating the Method\label{sec:validating}}
In order to validate our Langevin integrator for arbitrarily-shaped
rigid bodies, we implemented the algorithm in {\sc
oopse}\cite{Meineke2005} and  compared the results of this algorithm
with the known 
hydrodynamic limiting behavior for a few model systems, and to
microcanonical molecular dynamics simulations for some more
complicated bodies. The model systems and their analytical behavior
(if known) are summarized below. Parameters for the primary particles
comprising our model systems are given in table \ref{tab:parameters},
and a sketch of the arrangement of these primary particles into the
model rigid bodies is shown in figure \ref{fig:models}. In table
\ref{tab:parameters}, $d$ and $l$ are the physical dimensions of
ellipsoidal (Gay-Berne) particles.  For spherical particles, the value
of the Lennard-Jones $\sigma$ parameter is the particle diameter
($d$).  Gay-Berne ellipsoids have an energy scaling parameter,
$\epsilon^s$, which describes the well depth for two identical
ellipsoids in a {\it side-by-side} configuration.  Additionally, a
well depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$,
describes the ratio between the well depths in the {\it end-to-end}
and side-by-side configurations.  For spheres, $\epsilon^r \equiv 1$.
Moments of inertia are also required to describe the motion of primary
particles with orientational degrees of freedom.

\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{Parameters for the primary particles in use by the rigid body
models in figure \ref{fig:models}.}
\begin{tabular}{lrcccccccc}
\hline 
 & & & & & & & \multicolumn{3}c{$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$)} \\
 & & $d$ (\AA) & $l$ (\AA) & $\epsilon^s$ (kcal/mol) & $\epsilon^r$ &
$m$ (amu) & $I_{xx}$ & $I_{yy}$ & $I_{zz}$ \\ \hline
Sphere   & & 6.5 & $= d$ & 0.8 & 1 & 190 & 802.75 & 802.75  & 802.75  \\
Ellipsoid & & 4.6 & 13.8  & 0.8 & 0.2 & 200 & 2105 & 2105 & 421 \\
Dumbbell &(2 identical spheres) & 6.5 & $= d$ & 0.8 & 1   & 190 & 802.75 & 802.75 & 802.75 \\
Banana  &(3 identical ellipsoids)& 4.2 & 11.2  & 0.8 & 0.2 & 240 & 10000 & 10000 & 0 \\
Lipid: & Spherical Head & 6.5 & $= d$ & 0.185 & 1 & 196 & & & \\
       & Ellipsoidal Tail & 4.6 & 13.8  & 0.8   & 0.2 & 760 & 45000 & 45000 & 9000 \\
Solvent &  & 4.7 & $= d$ & 0.8 & 1   & 72.06 & & & \\
\hline
\end{tabular}
\label{tab:parameters}
\end{center}
\end{minipage}
\end{table*}

\begin{figure}
\centering
\includegraphics[width=3in]{sketch}
\caption[Sketch of the model systems]{A sketch of the model systems
used in evaluating the behavior of the rigid body Langevin
integrator.} \label{fig:models}
\end{figure}

\subsection{Simulation Methodology}
We performed reference microcanonical simulations with explicit
solvents for each of the different model system.  In each case there
was one solute model and 1929 solvent molecules present in the
simulation box.  All simulations were equilibrated using a
constant-pressure and temperature integrator with target values of 300
K for the temperature and 1 atm for pressure.  Following this stage,
further equilibration and sampling was done in a microcanonical
ensemble.  Since the model bodies are typically quite massive, we were
able to use a time step of 25 fs. 

The model systems studied used both Lennard-Jones spheres as well as
uniaxial Gay-Berne ellipoids. In its original form, the Gay-Berne
potential was a single site model for the interactions of rigid
ellipsoidal molecules.\cite{Gay81} It can be thought of as a
modification of the Gaussian overlap model originally described by
Berne and Pechukas.\cite{Berne72} The potential is constructed in the
familiar form of the Lennard-Jones function using
orientation-dependent $\sigma$ and $\epsilon$ parameters,
\begin{equation*}
V_{ij}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat
r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
{\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u
}_i},
{\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12}
-\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j},
{\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right]
\label{eq:gb}
\end{equation*}

The range $(\sigma({\bf \hat{u}}_{i},{\bf \hat{u}}_{j},{\bf
\hat{r}}_{ij}))$, and strength $(\epsilon({\bf \hat{u}}_{i},{\bf
\hat{u}}_{j},{\bf \hat{r}}_{ij}))$ parameters 
are dependent on the relative orientations of the two ellipsoids (${\bf
\hat{u}}_{i},{\bf \hat{u}}_{j}$) as well as the direction of the
inter-ellipsoid separation (${\bf \hat{r}}_{ij}$).  The shape and
attractiveness of each ellipsoid is governed by a relatively small set
of parameters: $l$ and $d$ describe the length and width of each
uniaxial ellipsoid, while $\epsilon^s$, which describes the well depth
for two identical ellipsoids in a {\it side-by-side} configuration.
Additionally, a well depth aspect ratio, $\epsilon^r = \epsilon^e /
\epsilon^s$, describes the ratio between the well depths in the {\it
end-to-end} and side-by-side configurations.  Details of the potential
are given elsewhere,\cite{Luckhurst90,Golubkov06,SunGezelter08} and an
excellent overview of the computational methods that can be used to
efficiently compute forces and torques for this potential can be found
in Ref. \citen{Golubkov06}

For the interaction between nonequivalent uniaxial ellipsoids (or
between spheres and ellipsoids), the spheres are treated as ellipsoids
with an aspect ratio of 1 ($d = l$) and with an well depth ratio
($\epsilon^r$) of 1 ($\epsilon^e = \epsilon^s$).  The form of the
Gay-Berne potential we are using was generalized by Cleaver {\it et
al.} and is appropriate for dissimilar uniaxial
ellipsoids.\cite{Cleaver96}

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}

To measure shear viscosities from our microcanonical simulations, we
used the Einstein form of the pressure correlation function,\cite{hess:209}
\begin{equation}
\eta = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left(
\int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \right\rangle_{t_0}.
\label{eq:shear}
\end{equation}
A similar form exists for the bulk viscosity
\begin{equation}
\kappa = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left(
\int_{t_0}^{t_0 + t} 
\left(P\left(t'\right)-\left\langle P \right\rangle \right)dt'
\right)^2 \right\rangle_{t_0}.
\end{equation}
Alternatively, the shear viscosity can also be calculated using a
Green-Kubo formula with the off-diagonal pressure tensor correlation function,
\begin{equation}
\eta = \frac{V}{k_B T} \int_0^{\infty} \left\langle P_{xz}(t_0) P_{xz}(t_0
+ t) \right\rangle_{t_0} dt,
\end{equation}
although this method converges extremely slowly and is not practical
for obtaining viscosities from molecular dynamics simulations. 

The Langevin dynamics for the different model systems were performed
at the same temperature as the average temperature of the
microcanonical simulations and with a solvent viscosity taken from
Eq. (\ref{eq:shear}) applied to these simulations.  We used 1024
independent solute simulations to obtain statistics on our Langevin
integrator.

\subsection{Analysis}

The quantities of interest when comparing the Langevin integrator to
analytic hydrodynamic equations and to molecular dynamics simulations
are typically translational diffusion constants and orientational
relaxation times.  Translational diffusion constants for point
particles are computed easily from the long-time slope of the
mean-square displacement,
\begin{equation}
D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \left\langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \right\rangle,
\end{equation}
of the solute molecules.  For models in which the translational
diffusion tensor (${\bf D}_{tt}$) has non-degenerate eigenvalues
(i.e. any non-spherically-symmetric rigid body), it is possible to
compute the diffusive behavior for motion parallel to each body-fixed
axis by projecting the displacement of the particle onto the
body-fixed reference frame at $t=0$.  With an isotropic solvent, as we
have used in this study, there are differences between the three
diffusion constants, but these must converge to the same value at
longer times.  Translational diffusion constants for the different
shaped models are shown in table \ref{tab:translation}.

In general, the three eigenvalues ($D_1, D_2, D_3$) of the rotational
diffusion tensor (${\bf D}_{rr}$) measure the diffusion of an object
{\it around} a particular body-fixed axis and {\it not} the diffusion
of a vector pointing along the axis.  However, these eigenvalues can
be combined to find 5 characteristic rotational relaxation
times,\cite{PhysRev.119.53,Berne90}
\begin{eqnarray}
1 / \tau_1 & = & 6 D_r + 2 \Delta \\
1 / \tau_2 & = & 6 D_r - 2 \Delta \\
1 / \tau_3 & = & 3 (D_r + D_1)  \\
1 / \tau_4 & = & 3 (D_r + D_2) \\
1 / \tau_5 & = & 3 (D_r + D_3) 
\end{eqnarray}
where 
\begin{equation}
D_r = \frac{1}{3} \left(D_1 + D_2 + D_3 \right)
\end{equation}
and 
\begin{equation}
\Delta = \left( (D_1 - D_2)^2 + (D_3 - D_1 )(D_3 - D_2)\right)^{1/2}
\end{equation}
Each of these characteristic times can be used to predict the decay of
part of the rotational correlation function when $\ell = 2$,
\begin{equation}
C_2(t)  =  \frac{a^2}{N^2} e^{-t/\tau_1} + \frac{b^2}{N^2} e^{-t/\tau_2}.
\end{equation}
This is the same as the $F^2_{0,0}(t)$ correlation function that
appears in Ref. \citen{Berne90}.  The amplitudes of the two decay
terms are expressed in terms of three dimensionless functions of the
eigenvalues: $a = \sqrt{3} (D_1 - D_2)$, $b = (2D_3 - D_1 - D_2 +
2\Delta)$, and $N = 2 \sqrt{\Delta b}$.  Similar expressions can be
obtained for other angular momentum correlation
functions.\cite{PhysRev.119.53,Berne90} In all of the model systems we
studied, only one of the amplitudes of the two decay terms was
non-zero, so it was possible to derive a single relaxation time for
each of the hydrodynamic tensors. In many cases, these characteristic
times are averaged and reported in the literature as a single relaxation
time,\cite{Garcia-de-la-Torre:1997qy}
\begin{equation}
1 / \tau_0 = \frac{1}{5} \sum_{i=1}^5 \tau_{i}^{-1},
\end{equation}
although for the cases reported here, this averaging is not necessary
and only one of the five relaxation times is relevant.

To test the Langevin integrator's behavior for rotational relaxation,
we have compared the analytical orientational relaxation times (if
they are known) with the general result from the diffusion tensor and
with the results from both the explicitly solvated molecular dynamics
and Langevin simulations.  Relaxation times from simulations (both
microcanonical and Langevin), were computed using Legendre polynomial
correlation functions for a unit vector (${\bf u}$) fixed along one or
more of the body-fixed axes of the model.
\begin{equation}
C_{\ell}(t)  =  \left\langle P_{\ell}\left({\bf u}_{i}(t) \cdot {\bf
u}_{i}(0) \right) \right\rangle 
\end{equation}
For simulations in the high-friction limit, orientational correlation
times can then be obtained from exponential fits of this function, or by
integrating,
\begin{equation}
\tau = \ell (\ell + 1) \int_0^{\infty} C_{\ell}(t) dt.
\end{equation}
In lower-friction solvents, the Legendre correlation functions often
exhibit non-exponential decay, and may not be characterized by a
single decay constant.

In table \ref{tab:rotation} we show the characteristic rotational
relaxation times (based on the diffusion tensor) for each of the model
systems compared with the values obtained via microcanonical and Langevin
simulations.

\subsection{Spherical particles}
Our model system for spherical particles was a Lennard-Jones sphere of
diameter ($\sigma$) 6.5 \AA\ in a sea of smaller spheres ($\sigma$ =
4.7 \AA).  The well depth ($\epsilon$) for both particles was set to
an arbitrary value of 0.8 kcal/mol.  

The Stokes-Einstein behavior of large spherical particles in
hydrodynamic flows is well known, giving translational friction
coefficients of $6 \pi \eta R$ (stick boundary conditions) and
rotational friction coefficients of $8 \pi \eta R^3$.  Recently,
Schmidt and Skinner have computed the behavior of spherical tag
particles in molecular dynamics simulations, and have shown that {\it
slip} boundary conditions ($\Xi_{tt} = 4 \pi \eta R$) may be more
appropriate for molecule-sized spheres embedded in a sea of spherical
solvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx}

Our simulation results show similar behavior to the behavior observed
by Schmidt and Skinner.  The diffusion constant obtained from our
microcanonical molecular dynamics simulations lies between the slip
and stick boundary condition results obtained via Stokes-Einstein
behavior.  Since the Langevin integrator assumes Stokes-Einstein stick
boundary conditions in calculating the drag and random forces for
spherical particles, our Langevin routine obtains nearly quantitative
agreement with the hydrodynamic results for spherical particles.  One
avenue for improvement of the method would be to compute elements of
$\Xi_{tt}$ assuming behavior intermediate between the two boundary
conditions.

In the explicit solvent simulations, both our solute and solvent
particles were structureless, exerting no torques upon each other.
Therefore, there are not rotational correlation times available for
this model system.

\subsection{Ellipsoids} 
For uniaxial ellipsoids ($a > b = c$), Perrin's formulae for both
translational and rotational diffusion of each of the body-fixed axes
can be combined to give a single translational diffusion
constant,\cite{Berne90} 
\begin{equation}
D = \frac{k_B T}{6 \pi \eta a} G(\rho),
\label{Dperrin}
\end{equation}
as well as a single rotational diffusion coefficient,
\begin{equation}
\Theta = \frac{3 k_B T}{16 \pi \eta a^3} \left\{ \frac{(2 - \rho^2)
G(\rho) - 1}{1 - \rho^4} \right\}.
\label{ThetaPerrin}
\end{equation}
In these expressions, $G(\rho)$ is a function of the axial ratio
($\rho = b / a$), which for prolate ellipsoids, is
\begin{equation}
G(\rho) = (1- \rho^2)^{-1/2} \ln \left\{ \frac{1 + (1 -
\rho^2)^{1/2}}{\rho} \right\}
\label{GPerrin}
\end{equation}
Again, there is some uncertainty about the correct boundary conditions
to use for molecular-scale ellipsoids in a sea of similarly-sized
solvent particles.  Ravichandran and Bagchi found that {\it slip}
boundary conditions most closely resembled the simulation
results,\cite{Ravichandran:1999fk} in agreement with earlier work of
Tang and Evans.\cite{TANG:1993lr}

Even though there are analytic resistance tensors for ellipsoids, we
constructed a rough-shell model using 2135 beads (each with a diameter
of 0.25 \AA) to approximate the shape of the model ellipsoid.  We
compared the Langevin dynamics from both the simple ellipsoidal
resistance tensor and the rough shell approximation with
microcanonical simulations and the predictions of Perrin.  As in the
case of our spherical model system, the Langevin integrator reproduces
almost exactly the behavior of the Perrin formulae (which is
unsurprising given that the Perrin formulae were used to derive the
drag and random forces applied to the ellipsoid).  We obtain
translational diffusion constants and rotational correlation times
that are within a few percent of the analytic values for both the
exact treatment of the diffusion tensor as well as the rough-shell
model for the ellipsoid.

The translational diffusion constants from the microcanonical simulations 
agree well with the predictions of the Perrin model, although the rotational
correlation times are a factor of 2 shorter than expected from hydrodynamic 
theory.  One explanation for the slower rotation
of explicitly-solvated ellipsoids is the possibility that solute-solvent 
collisions happen at both ends of the solute whenever the principal 
axis of the ellipsoid is turning. In the upper portion of figure 
\ref{fig:explanation} we sketch a physical picture of this explanation.
Since our Langevin integrator is providing nearly quantitative agreement with 
the Perrin model, it also predicts orientational diffusion for ellipsoids that
exceed explicitly solvated correlation times by a factor of two.

\subsection{Rigid dumbbells}
Perhaps the only {\it composite} rigid body for which analytic
expressions for the hydrodynamic tensor are available is the
two-sphere dumbbell model.  This model consists of two non-overlapping
spheres held by a rigid bond connecting their centers. There are
competing expressions for the 6x6 resistance tensor for this
model. Equation (\ref{introEquation:oseenTensor}) above gives the
original Oseen tensor, while the second order expression introduced by
Rotne and Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la
Torre and Bloomfield,\cite{Torre1977} is given above as
Eq. (\ref{introEquation:RPTensorNonOverlapped}).  In our case, we use
a model dumbbell in which the two spheres are identical Lennard-Jones
particles ($\sigma$ = 6.5 \AA\ , $\epsilon$ = 0.8 kcal / mol) held at
a distance of 6.532 \AA.

The theoretical values for the translational diffusion constant of the
dumbbell are calculated from the work of Stimson and Jeffery, who
studied the motion of this system in a flow parallel to the
inter-sphere axis,\cite{Stimson:1926qy} and Davis, who studied the
motion in a flow {\it perpendicular} to the inter-sphere
axis.\cite{Davis:1969uq} We know of no analytic solutions for the {\it
orientational} correlation times for this model system (other than
those derived from the 6 x 6 tensors mentioned above).

The bead model for this model system comprises the two large spheres
by themselves, while the rough shell approximation used 3368 separate
beads (each with a diameter of 0.25 \AA) to approximate the shape of
the rigid body.  The hydrodynamics tensors computed from both the bead
and rough shell models are remarkably similar.  Computing the initial
hydrodynamic tensor for a rough shell model can be quite expensive (in
this case it requires inverting a 10104 x 10104 matrix), while the
bead model is typically easy to compute (in this case requiring
inversion of a 6 x 6 matrix).  

\begin{figure}
\centering
\includegraphics[width=2in]{RoughShell}
\caption[Model rigid bodies and their rough shell approximations]{The
model rigid bodies (left column) used to test this algorithm and their
rough-shell approximations (right-column) that were used to compute
the hydrodynamic tensors.  The top two models (ellipsoid and dumbbell)
have analytic solutions and were used to test the rough shell
approximation.  The lower two models (banana and lipid) were compared
with explicitly-solvated molecular dynamics simulations. } 
\label{fig:roughShell}
\end{figure}


Once the hydrodynamic tensor has been computed, there is no additional
penalty for carrying out a Langevin simulation with either of the two
different hydrodynamics models.  Our naive expectation is that since
the rigid body's surface is roughened under the various shell models,
the diffusion constants will be even farther from the ``slip''
boundary conditions than observed for the bead model (which uses a
Stokes-Einstein model to arrive at the hydrodynamic tensor).  For the
dumbbell, this prediction is correct although all of the Langevin
diffusion constants are within 6\% of the diffusion constant predicted
from the fully solvated system.

For rotational motion, Langevin integration (and the hydrodynamic tensor) 
yields rotational correlation times that are substantially shorter than those
obtained from explicitly-solvated simulations.  It is likely that this is due
to the large size of the explicit solvent spheres, a feature that prevents 
the solvent from coming in contact with a substantial fraction of the surface 
area of the dumbbell.  Therefore, the explicit solvent only provides drag 
over a substantially reduced surface area of this model, while the 
hydrodynamic theories utilize the entire surface area for estimating 
rotational diffusion.  A sketch of the free volume available in the explicit
solvent simulations is shown in figure \ref{fig:explanation}.


\begin{figure}
\centering
\includegraphics[width=6in]{explanation}
\caption[Explanations of the differences between orientational
correlation times for explicitly-solvated models and hydrodynamics
predictions]{Explanations of the differences between orientational
correlation times for explicitly-solvated models and hydrodynamic
predictions.   For the ellipsoids (upper figures), rotation of the
principal axis can involve correlated collisions at both sides of the
solute.  In the rigid dumbbell model (lower figures), the large size
of the explicit solvent spheres prevents them from coming in contact
with a substantial fraction of the surface area of the dumbbell.
Therefore, the explicit solvent only provides drag over a
substantially reduced surface area of this model, where the
hydrodynamic theories utilize the entire surface area for estimating
rotational diffusion. 
} \label{fig:explanation}
\end{figure}


\subsection{Composite banana-shaped molecules}
Banana-shaped rigid bodies composed of three Gay-Berne ellipsoids have
been used by Orlandi {\it et al.} to observe mesophases in
coarse-grained models for bent-core liquid crystalline
molecules.\cite{Orlandi:2006fk} We have used the same overlapping
ellipsoids as a way to test the behavior of our algorithm for a
structure of some interest to the materials science community,
although since we are interested in capturing only the hydrodynamic
behavior of this model, we have left out the dipolar interactions of
the original Orlandi model.
 
A reference system composed of a single banana rigid body embedded in a
sea of 1929 solvent particles was created and run under standard
(microcanonical) molecular dynamics.  The resulting viscosity of this 
mixture was 0.298 centipoise (as estimated using Eq. (\ref{eq:shear})).
To calculate the hydrodynamic properties of the banana rigid body model,
we created a rough shell (see Fig.~\ref{fig:roughShell}), in which 
the banana is represented as a ``shell'' made of 3321 identical beads 
(0.25 \AA\  in diameter) distributed on the surface.  Applying the 
procedure described in Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
identified the center of resistance, ${\bf r} = $(0 \AA, 0.81 \AA, 0 \AA), as
well as the resistance tensor, 
\begin{equation*}
\Xi = 
\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),
\end{equation*}
where the units for translational, translation-rotation coupling and
rotational tensors are (kcal fs / mol \AA$^2$), (kcal fs / mol \AA\ rad),
and (kcal fs / mol rad$^2$), respectively.

The Langevin rigid-body integrator (and the hydrodynamic diffusion tensor)
are essentially quantitative for translational diffusion of this model.  
Orientational correlation times under the Langevin rigid-body integrator 
are within 11\% of the values obtained from explicit solvent, but these 
models also exhibit some solvent inaccessible surface area in the 
explicitly-solvated case.  

\subsection{Composite sphero-ellipsoids}
Spherical heads perched on the ends of Gay-Berne ellipsoids have been
used recently as models for lipid
molecules.\cite{SunGezelter08,Ayton01}
MORE DETAILS

A reference system composed of a single lipid rigid body embedded in a
sea of 1929 solvent particles was created and run under standard
(microcanonical) molecular dynamics.  The resulting viscosity of this
mixture was 0.349 centipoise (as estimated using
Eq. (\ref{eq:shear})).  To calculate the hydrodynamic properties of
the lipid rigid body model, we created a rough shell (see
Fig.~\ref{fig:roughShell}), in which the lipid is represented as a
``shell'' made of 3550 identical beads (0.25 \AA\ in diameter)
distributed on the surface.  Applying the procedure described in
Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
identified the center of resistance, ${\bf r} = $(0 \AA, 0 \AA, 1.46
\AA).


\subsection{Summary}
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{figure}
\centering
\includegraphics[width=\linewidth]{graph}
\caption[Mean squared displacements and orientational
correlation functions for each of the model rigid bodies.]{The
mean-squared displacements ($\langle r^2(t) \rangle$) and
orientational correlation functions ($C_2(t)$) for each of the model
rigid bodies studied.  The circles are the results for microcanonical
simulations with explicit solvent molecules, while the other data sets
are results for Langevin dynamics using the different hydrodynamic
tensor approximations.  The Perrin model for the ellipsoids is
considered the ``exact'' hydrodynamic behavior (this can also be said
for the translational motion of the dumbbell operating under the bead
model). In most cases, the various hydrodynamics models reproduce
each other quantitatively.}
\label{fig:results} 
\end{figure}

\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{Translational diffusion constants (D) for the model systems 
calculated using microcanonical simulations (with explicit solvent), 
theoretical predictions, and Langevin simulations (with implicit solvent).
Analytical solutions for the exactly-solved hydrodynamics models are
from Refs. \citen{Einstein05} (sphere), \citen{Perrin1934} and \citen{Perrin1936}
(ellipsoid), \citen{Stimson:1926qy} and \citen{Davis:1969uq}
(dumbbell). The other model systems have no known analytic solution.
All  diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (=
$10^{-4}$ \AA$^2$  / fs). }
\begin{tabular}{lccccccc}
\hline
 & \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\
\cline{2-3} \cline{5-7}
model & $\eta$ (centipoise)  & D & & Analytical & method & Hydrodynamics & simulation \\ 
\hline
sphere    & 0.279  & 3.06 & & 2.42 & exact       & 2.42 & 2.33 \\
ellipsoid & 0.255  & 2.44 & & 2.34 & exact       & 2.34 & 2.37 \\
          & 0.255  & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\
dumbbell  & 0.308  & 2.06 & & 1.64 & bead model  & 1.65 & 1.62 \\
          & 0.308  & 2.06 & & 1.64 & rough shell & 1.59 & 1.62 \\
banana    & 0.298  & 1.53 & &      & rough shell & 1.56 & 1.55 \\
lipid     & 0.349  & 0.96 & &      & rough shell & 1.33 & 1.32 \\
\end{tabular}
\label{tab:translation}
\end{center}
\end{minipage}
\end{table*}

\begin{table*}
\begin{minipage}{\linewidth}
\begin{center}
\caption{Orientational relaxation times ($\tau$) for the model systems using
microcanonical simulation (with explicit solvent), theoretical
predictions, and Langevin simulations (with implicit solvent). All
relaxation times are for the rotational correlation function with
$\ell = 2$ and are reported in units of ps.  The ellipsoidal model has
an exact solution for the orientational correlation time due to
Perrin, but the other model systems have no known analytic solution.}
\begin{tabular}{lccccccc}
\hline
 & \multicolumn{2}c{microcanonical simulation} & & \multicolumn{3}c{Theoretical} & Langevin \\
\cline{2-3} \cline{5-7}
model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic  & simulation \\ 
\hline
sphere    & 0.279  &      & & 9.69 & exact       & 9.69 & 9.64 \\
ellipsoid & 0.255  & 46.7 & & 22.0 & exact       & 22.0 & 22.2 \\
          & 0.255  & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\
dumbbell  & 0.308  & 14.1 & &      & bead model  & 50.0 & 50.1 \\
          & 0.308  & 14.1 & &      & rough shell & 41.5 & 41.3 \\
banana    & 0.298  & 63.8 & &      & rough shell & 70.9 & 70.9 \\
lipid     & 0.349  & 78.0 & &      & rough shell & 76.9 & 77.9 \\
\hline
\end{tabular}
\label{tab:rotation}
\end{center}
\end{minipage}
\end{table*}

\section{Application: A rigid-body lipid bilayer}

The Langevin dynamics integrator was applied to study the formation of
corrugated structures emerging from simulations of the coarse grained
lipid molecular models presented above.  The initial configuration is
taken from our molecular dynamics studies on lipid bilayers with
lennard-Jones sphere solvents. The solvent molecules were excluded
from the system, and the experimental value for the viscosity of water
at 20C ($\eta = 1.00$ cp) was used to mimic the hydrodynamic effects
of the solvent.  The absence of explicit solvent molecules and the
stability of the integrator allowed us to take timesteps of 50 fs.  A
total simulation run time of 100 ns was sampled.
Fig. \ref{fig:bilayer} shows the configuration of the system after 100
ns, and the ripple structure remains stable during the entire
trajectory.  Compared with using explicit bead-model solvent
molecules, the efficiency of the simulation has increased by an order
of magnitude.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{bilayer}
\caption[Snapshot of a bilayer of rigid-body models for lipids]{A
snapshot of a bilayer composed of rigid-body models for lipid
molecules evolving using the Langevin integrator described in this
work.} \label{fig:bilayer}
\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. 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{jcp}
\bibliography{langevin}

\end{document}
Revision:	3317
Committed:	Fri Jan 18 22:07:35 2008 UTC (16 years, 7 months ago) by xsun
Content type:	application/x-tex
File size:	58027 byte(s)
Log Message:	add one ref