trunk/langevin/langevin.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
%\documentclass[aps,prb,preprint,amsmath,amssymb,endfloats]{revtex4}
\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{times}
\usepackage{mathptmx}
\usepackage{setspace}
\usepackage{endfloat}
\usepackage{caption}
%\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

% double space list of tables and figures
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{1.66}}
\setlength{\abovecaptionskip}{20 pt}
\setlength{\belowcaptionskip}{30 pt}

\renewcommand\citemid{\ } % no comma in optional referenc note

\begin{document}

\title{Langevin dynamics for 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 


\begin{doublespace}

\begin{abstract}
We present an algorithm for carrying out Langevin dynamics simulations
on complex rigid bodies by incorporating the hydrodynamic resistance
tensors for arbitrary shapes into an advanced rotational integration
scheme.  The integrator gives quantitative agreement with both
analytic and approximate hydrodynamic theories for a number of model
rigid bodies, and works well at reproducing the solute dynamical
properties (diffusion constants, and orientational relaxation times)
obtained from explicitly-solvated simulations.
\end{abstract}

\newpage

%\narrowtext

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


\section{Introduction}

%applications of langevin dynamics
Langevin dynamics, which mimics a heat bath using both 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, Liwo 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' 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}

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

The use of rigid substructures,\cite{Chun:2000fj}
coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunX._jp0762020}
and ellipsoidal representations of protein side
chains~\cite{Fogolari:1996lr} 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.  Also, the atoms involved in a rigid
or coarse-grained structure have solvent-mediated interactions with
each other, and these interactions are ignored if all atoms are
treated as separate spherical particles.  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{FIXMAN:1986lr,Ramachandran1996} 

In order to account for the diffusion anisotropy of complex shapes,
Fernandes and Garc\'{i}a de la Torre improved an earlier 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} The effects of these
biases can be reduced by decreasing the size of the time step. 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 linked polymer segments in a
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 hydrodynamic
properties. However, typical nonskew bodies like cylinders and
ellipsoids are inadequate to represent most complex macromolecular
assemblies. Therefore, the goal of this work is to adapt some of the
hydrodynamic methodologies developed to treat Brownian motion of
complex assemblies into a Langevin integrator for rigid bodies with
arbitrary shapes.

\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{SunX._jp0762020} Many of the water models in common
use are also rigid-body
models,\cite{Jorgensen83,Berendsen81,Berendsen87} although they are
typically evolved in molecular dynamics simulations 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 $\frac{1}{\sin \theta}$ singularities, the
numerical integration of corresponding equations of these motion can
become inaccurate (and inefficient).  Although the use of 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} The Evans quaternion
approach uses a nonseparable Hamiltonian, and this has 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 constraints
are used to guarantee rigidity. However, due to their iterative
nature, the SHAKE and RATTLE algorithms converge very slowly when the
number of constraints (and the number of particles that belong to the
rigid body) increases.\cite{Ryckaert1977,Andersen1983}

In order to develop a stable and efficient integration scheme that
preserves most constants of the motion in microcanonical simulations,
symplectic propagators are necessary.  By introducing a conjugate
momentum to the rotation matrix ${\bf Q}$ and re-formulating
Hamilton's equations, a symplectic orientational integrator,
RSHAKE,\cite{Kol1997} was proposed to evolve rigid bodies on a
constraint manifold by iteratively satisfying the orthogonality
constraint ${\bf Q}^T {\bf 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} removes
most of the limitations mentioned above and is therefore the basis for
our Langevin integrator.

The goal of the present work is to develop a Langevin dynamics
algorithm for arbitrary-shaped rigid particles by integrating accurate
estimates of the friction tensor that can be obtaind from classical
hydrodynamics theory (or other recently-developed hydrodynamic
procedures) into a stable and efficient rigid body dynamics
propagator.  In the sections below, we review some of the theory of
hydrodynamic tensors developed primarily for Brownian simulations of
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 for a solute particle can be
determined using the velocity autocorrelation function from a
simulation with explicit solvent molecules. However, this approach
becomes impractical when the solute becomes complex.  Instead, various
approaches based on hydrodynamics have been developed to calculate
static 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 a friction force ($\mathbf{f}_f$) and torque
($\mathbf{\tau}_f$) in opposition to 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}
For an arbitrary body moving in a fluid, Peters has derived a set of
fluctuation-dissipation relations for the friction
tensors,\cite{Peters:1999qy,Peters:1999uq,Peters:2000fk}
\begin{eqnarray}
\Xi^{tt} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
F}(0) {\bf F}(-s) \rangle_{eq} - \langle {\bf F} \rangle_{eq}^2
\right] ds \\ 
\notag \\ 
\Xi^{tr} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
F}(0) {\bf \tau}(-s) \rangle_{eq} - \langle {\bf F} \rangle_{eq}
\langle {\bf \tau} \rangle_{eq} \right] ds \\ 
\notag \\
\Xi^{rt} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
\tau}(0) {\bf F}(-s) \rangle_{eq} - \langle {\bf \tau} \rangle_{eq}
\langle {\bf F} \rangle_{eq} \right] ds \\ 
\notag \\
\Xi^{rr} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
\tau}(0) {\bf \tau}(-s) \rangle_{eq} - \langle {\bf \tau} \rangle_{eq}^2
\right] ds 
\end{eqnarray}
In these expressions, the forces (${\bf F}$) and torques (${\bf
\tau}$) are those that arise solely from the interactions of the body with
the surrounding fluid. For a single solute body in an isotropic fluid,
the average forces and torques in these expressions ($\langle {\bf F}
\rangle_{eq}$ and $\langle {\bf \tau} \rangle_{eq}$) 
vanish, and one obtains the simpler force-torque correlation formulae
of Nienhuis.\cite{Nienhuis:1970lr} Molecular dynamics simulations with
explicit solvent molecules can be used to obtain estimates of the
friction tensors with these formulae. In practice, however, one needs
relatively long simulations with frequently-stored force and torque
information to compute friction tensors, and this becomes
prohibitively expensive when there are large numbers of large solute
particles.  For bodies with simple shapes, there are a number of
approximate expressions that allow computation of these tensors
without the need for expensive simulations that utilize explicit
solvent particles.

\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}}
For a spherical body under ``stick'' boundary conditions, the
translational and rotational friction tensors can be estimated from
Stokes' law,
\begin{equation}
\label{eq:StokesTranslation}
\Xi^{tt}  = \left( \begin{array}{*{20}c}
   {6\pi \eta \rho} & 0 & 0  \\
   0 & {6\pi \eta \rho} & 0  \\
   0 & 0 & {6\pi \eta \rho}  \\
\end{array} \right)
\end{equation}
and
\begin{equation}
\label{eq:StokesRotation}
\Xi^{rr}  = \left( \begin{array}{*{20}c}
   {8\pi \eta \rho^3 } & 0 & 0  \\
   0 & {8\pi \eta \rho^3 } & 0  \\
   0 & 0 & {8\pi \eta \rho^3 }  \\
\end{array} \right)
\end{equation}
where $\eta$ is the viscosity of the solvent and $\rho$ is the
hydrodynamic radius.  The presence of the rotational resistance tensor
implies that the spherical body has internal structure and
orientational degrees of freedom that must be propagated in time.  For
non-structured spherical bodies (i.e. the atoms in a traditional
molecular dynamics simulation) these degrees of freedom do not exist.

Other non-spherical shapes, such as cylinders and ellipsoids, are
widely used as references for developing new hydrodynamic theories,
because their properties can be calculated exactly. In 1936, Perrin
extended Stokes' law to general
ellipsoids,\cite{Perrin1934,Perrin1936} described in Cartesian
coordinates as
\begin{equation}
\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1.
\end{equation}
Here, the semi-axes are of lengths $a$, $b$, and $c$. Due to the
complexity of the elliptic integral, only uniaxial ellipsoids, either
prolate ($a \ge b = c$) or oblate ($a < b = c$), were solved
exactly. Introducing an elliptic integral parameter $S$ for prolate,
\begin{equation}
S = \frac{2}{\sqrt{a^2  - b^2}} \ln \frac{a + \sqrt{a^2  - b^2}}{b},
\end{equation}
and oblate,
\begin{equation}
S = \frac{2}{\sqrt {b^2  - a^2 }} \arctan \frac{\sqrt {b^2  - a^2}}{a},
\end{equation}
ellipsoids, it is possible to write down exact solutions for the
resistance tensors.  As is the case for spherical bodies, the translational,
\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 rotational,
\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}
resistance tensors are diagonal $3 \times 3$ matrices. 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}}
Other than the fluctuation dissipation formulae given by
Peters,\cite{Peters:1999qy,Peters:1999uq,Peters:2000fk} there are no
analytic solutions for the friction tensor for rigid molecules of
arbitrary shape. The ellipsoid of revolution and general triaxial
ellipsoid models have been widely used to approximate the hydrodynamic
properties of rigid bodies. However, the mapping from all possible
ellipsoidal spaces ($r$-space) to all possible combinations of
rotational diffusion coefficients ($D$-space) is not
unique.\cite{Wegener1979} Additionally, because there is intrinsic
coupling between translational and rotational motion of {\it skew}
rigid bodies, general ellipsoids are not always suitable for modeling
rigid molecules.  A number of studies have been devoted to determining
the friction tensor for irregular shapes using methods in which the
molecule of interest is modeled with a combination of
spheres\cite{Carrasco1999} and the hydrodynamic properties of the
molecule are then calculated using a set of two-point interaction
tensors.  We have found the {\it bead} and {\it rough shell} models of
Carrasco and Garc\'{i}a de la Torre to be the most useful of these
methods,\cite{Carrasco1999} and we review the basic outline of the
rough shell approach here.  A more thorough explanation can be found
in Ref. \citen{Carrasco1999}.

Consider a rigid assembly of $N$ small beads moving through a
continuous medium.  Due to hydrodynamic interactions between the
beads, the net velocity of the $i^\mathrm{th}$ bead relative to the
medium, ${\bf v}'_i$, is different than its unperturbed velocity ${\bf
v}_i$,
\begin{equation}
{\bf v}'_i  = {\bf v}_i  - \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }
\end{equation}
where ${\bf F}_j$ is the frictional force on the medium due to bead $j$, and
${\bf T}_{ij}$ is the hydrodynamic interaction tensor between the two beads.
The frictional force felt by the $i^\mathrm{th}$ bead is proportional to
its net velocity
\begin{equation}
{\bf F}_i  = \xi_i {\bf v}_i  - \xi_i \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }.
\label{introEquation:tensorExpression}
\end{equation}
Eq. (\ref{introEquation:tensorExpression}) defines the two-point
hydrodynamic tensor, ${\bf T}_{ij}$.  There have been many proposed
solutions to this equation, including the simple solution given by
Oseen and Burgers in 1930 for two beads of identical radius.  A second
order expression for beads of different hydrodynamic radii was
introduced by Rotne and Prager,\cite{Rotne1969} and improved by
Garc\'{i}a de la Torre and Bloomfield,\cite{Torre1977}
\begin{equation}
{\bf T}_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {{\bf I} +
\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right) + \frac{{\rho
_i^2  + \rho_j^2 }}{{R_{ij}^2 }}\left( {\frac{{\bf I}}{3} -
\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
\label{introEquation:RPTensorNonOverlapped}
\end{equation}
Here ${\bf R}_{ij}$ is the distance vector between beads $i$ and $j$.  Both
the Oseen-Burgers tensor and
Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption that
the beads do not overlap ($R_{ij} \ge \rho_i + \rho_j$). 

To calculate the resistance tensor for a body represented as the union
of many non-overlapping beads, we first pick an arbitrary origin $O$
and then construct a $3N \times 3N$ supermatrix consisting of $N
\times N$ ${\bf B}_{ij}$ blocks
\begin{equation}
{\bf B} = \left( \begin{array}{*{20}c}
 {\bf B}_{11} &  \ldots  & {\bf B}_{1N}   \\
\vdots  &  \ddots  &  \vdots   \\
 {\bf B}_{N1} &  \cdots  & {\bf B}_{NN}
\end{array} \right)
\end{equation}
${\bf B}_{ij}$ is a version of the hydrodynamic tensor which includes the
self-contributions for spheres,
\begin{equation}
{\bf B}_{ij}  = \delta _{ij} \frac{{\bf I}}{{6\pi \eta R_{ij}}} + (1 - \delta_{ij}
){\bf T}_{ij}
\end{equation}
where $\delta_{ij}$ is the Kronecker delta function. Inverting the
${\bf B}$ matrix, we obtain
\begin{equation}
{\bf C} = {\bf B}^{ - 1}  = \left(\begin{array}{*{20}c}
  {\bf C}_{11} &  \ldots  & {\bf C}_{1N} \\
 \vdots  &  \ddots  &  \vdots   \\
  {\bf C}_{N1} &  \cdots  & {\bf C}_{NN} 
\end{array} \right),
\end{equation}
which can be partitioned into $N \times N$ blocks labeled ${\bf C}_{ij}$.
(Each of the ${\bf C}_{ij}$ blocks is a $3 \times 3$ matrix.)  Using the
skew matrix,
\begin{equation}
{\bf 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)
\label{eq:skewMatrix}
\end{equation}
where $x_i$, $y_i$, $z_i$ are the components of the vector joining
bead $i$ and origin $O$, the elements of the resistance tensor (at the
arbitrary origin $O$) can be written as
\begin{eqnarray}
\label{introEquation:ResistanceTensorArbitraryOrigin}
 \Xi^{tt}  & = & \sum\limits_i {\sum\limits_j {{\bf C}_{ij} } } \notag , \\
 \Xi^{tr}  = \Xi _{}^{rt}  & = & \sum\limits_i {\sum\limits_j {{\bf U}_i {\bf C}_{ij} } } , \\
 \Xi^{rr}  & = &  -\sum\limits_i \sum\limits_j {\bf U}_i {\bf C}_{ij} {\bf U}_j + 6 \eta V {\bf I}. \notag
\end{eqnarray}
The final term in the expression for $\Xi^{rr}$ is a correction that
accounts for errors in the rotational motion of the 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} \rho_i^3,
\end{equation}
where $\rho_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 Garc\'{i}a de la Torre and
Rodes.\cite{Torre:1983lr}

In general, resistance tensors depend on the origin at which they were
computed.  However, the proper location for applying the friction
force is the center of resistance, the special point at which the
trace of rotational resistance tensor, $\Xi^{rr}$ reaches a minimum
value.  Mathematically, the center of resistance can also be defined
as the unique point for a 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 Eq. \ref{introEquation:ResistanceTensorArbitraryOrigin}, we can
easily derive that the {\it 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 ,${\bf r}_{OP} = (x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we
can obtain the resistance tensor at $P$ by
\begin{eqnarray}
\label{introEquation:resistanceTensorTransformation}
 \Xi_P^{tt}  & = & \Xi_O^{tt}  \notag \\
 \Xi_P^{tr}  = \Xi_P^{rt}  & = & \Xi_O^{tr}  - {\bf U}_{OP} \Xi _O^{tt}  \\
 \Xi_P^{rr}  & = &\Xi_O^{rr}  - {\bf U}_{OP} \Xi_O^{tt} {\bf U}_{OP}
+ \Xi_O^{tr} {\bf U}_{OP}  - {\bf U}_{OP} \left( \Xi_O^{tr}
\right)^{^T} \notag
\end{eqnarray}
where ${\bf U}_{OP}$ is the skew matrix (Eq. (\ref{eq:skewMatrix}))
for the vector between the origin $O$ and the point $P$. Using
Eqs.~\ref{introEquation:definitionCR}~and~\ref{introEquation:resistanceTensorTransformation},
one can locate the position of center of resistance,
\begin{equation*}
\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{equation*}
where $x_{OR}$, $y_{OR}$, $z_{OR}$ are the components of the vector
joining center of resistance $R$ and origin $O$.

For a general rigid molecular substructure, finding the $6 \times 6$
resistance tensor can be a computationally demanding task.  First, a
lattice of small beads that extends well beyond the boundaries of the
rigid substructure is created.  The lattice is typically composed of
0.25 \AA\ beads on a dense FCC lattice.  The lattice constant is taken
to be the bead diameter, so that adjacent beads are touching, but do
not overlap. To make a shape corresponding to the rigid structure,
beads that sit on lattice sites that are outside the van der Waals
radii of all of the atoms comprising the rigid body are excluded from
the calculation.

For large structures, most of the beads will be deep within the rigid
body and will not contribute to the hydrodynamic tensor.  In the {\it
rough shell} approach, beads which have all of their lattice neighbors
inside the structure are considered interior beads, and are removed
from the calculation.  After following this procedure, only those
beads in direct contact with the van der Waals surface of the rigid
body are retained.  For reasonably large molecular structures, this
truncation can still produce bead assemblies with thousands of
members.

If all of the {\it atoms} comprising the rigid substructure are
spherical and non-overlapping, the tensor in
Eq.~(\ref{introEquation:RPTensorNonOverlapped}) may be used directly
using the atoms themselves as the hydrodynamic beads.  This is a
variant of the {\it bead model} approach of Carrasco and Garc\'{i}a de
la Torre.\cite{Carrasco1999} In this case, the size of the ${\bf B}$
matrix can be quite small, and the calculation of the hydrodynamic
tensor is straightforward.

In general, the inversion of the ${\bf B}$ matrix is the most
computationally demanding task.  This inversion is done only once for
each type of rigid structure.  We have used straightforward
LU-decomposition to solve the linear system and to obtain the elements
of ${\bf C}$. Once ${\bf C}$ has been obtained, the location of the
center of resistance ($R$) is found and the resistance tensor at this
point is calculated.  The $3 \times 1$ vector giving the location of
the rigid body's center of resistance and the $6 \times 6$ resistance
tensor are then stored for use in the Langevin dynamics calculation.
These quantities depend on solvent viscosity and temperature and must
be recomputed if different simulation conditions are required.

\section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}

Consider the Langevin equations of motion in generalized coordinates
\begin{equation}
{\bf M} \dot{{\bf V}}(t) = {\bf F}_{s}(t) +
{\bf F}_{f}(t)  + {\bf F}_{r}(t) 
\label{LDGeneralizedForm}
\end{equation}
where ${\bf M}$ is a $6 \times 6$ diagonal mass matrix (which
includes the mass of the rigid body as well as the moments of inertia
in the body-fixed frame) and ${\bf V}$ is a generalized velocity,
${\bf V} =
\left\{{\bf v},{\bf \omega}\right\}$. The right side of
Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a
system force (${\bf F}_{s}$), a frictional or dissipative force (${\bf
F}_{f}$) and a stochastic force (${\bf F}_{r}$). While the evolution
of the system in Newtonian mechanics is typically done in the lab
frame, it is convenient to handle the dynamics of rigid bodies in
body-fixed frames. Thus the friction and random forces on each
substructure are calculated in a body-fixed frame and may converted
back to the lab frame using that substructure's rotation matrix (${\bf
Q}$):
\begin{equation}
{\bf F}_{f,r} = 
\left( \begin{array}{c}
{\bf f}_{f,r} \\
{\bf \tau}_{f,r}
\end{array} \right)
=
\left( \begin{array}{c}
{\bf Q}^{T} {\bf f}^{~b}_{f,r} \\
{\bf Q}^{T} {\bf \tau}^{~b}_{f,r}
\end{array} \right)
\end{equation}
The body-fixed friction force, ${\bf F}_{f}^{~b}$, is proportional to
the (body-fixed) velocity at the center of resistance
${\bf v}_{R}^{~b}$ and the angular velocity ${\bf \omega}$
\begin{equation}
{\bf F}_{f}^{~b}(t) = \left( \begin{array}{l}
 {\bf f}_{f}^{~b}(t) \\
 {\bf \tau}_{f}^{~b}(t) \\
 \end{array} \right) =  - \left( \begin{array}{*{20}c}
   \Xi_{R}^{tt} & \Xi_{R}^{rt} \\
   \Xi_{R}^{tr} & \Xi_{R}^{rr}    \\
\end{array} \right)\left( \begin{array}{l}
 {\bf v}_{R}^{~b}(t) \\
 {\bf \omega}(t) \\
 \end{array} \right),
\end{equation}
while the random force, ${\bf F}_{r}$, is a Gaussian stochastic
variable with zero mean and variance,
\begin{equation}
\left\langle {{\bf F}_{r}(t) ({\bf F}_{r}(t'))^T } \right\rangle  =
\left\langle {{\bf F}_{r}^{~b} (t) ({\bf F}_{r}^{~b} (t'))^T } \right\rangle  =
2 k_B T \Xi_R \delta(t - t'). \label{eq:randomForce}
\end{equation}
$\Xi_R$ is the $6\times6$ resistance tensor at the center of
resistance.  Once this tensor is known for a given rigid body (as
described in the previous section) obtaining a stochastic vector that
has the properties in Eq. (\ref{eq:randomForce}) can be done
efficiently by carrying out a one-time Cholesky decomposition to
obtain the square root matrix of the resistance tensor,
\begin{equation} 
\Xi_R = {\bf S} {\bf S}^{T},
\label{eq:Cholesky}
\end{equation} 
where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A
vector with the statistics required for the random force can then be
obtained by multiplying ${\bf S}$ onto a random 6-vector ${\bf Z}$ which
has elements chosen from a Gaussian distribution, such that:
\begin{equation}
\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot
{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij},
\end{equation}
where $\delta t$ is the timestep in use during the simulation. The
random force, ${\bf F}_{r}^{~b} = {\bf S} {\bf Z}$, can be shown to have the
correct properties required by Eq. (\ref{eq:randomForce}).

The equation of motion for the translational velocity of the center of
mass (${\bf v}$) can be written as
\begin{equation}
m \dot{{\bf v}} (t) =  {\bf f}_{s}(t) + {\bf f}_{f}(t) +
{\bf f}_{r}(t)
\end{equation}
Since the frictional and random forces are applied at the center of
resistance, which generally does not coincide with the center of mass,
extra torques are exerted at the center of mass. Thus, the net
body-fixed torque at the center of mass, $\tau^{~b}(t)$,
is given by
\begin{equation}
\tau^{~b} \leftarrow \tau_{s}^{~b} + \tau_{f}^{~b} + \tau_{r}^{~b} + {\bf r}_{MR} \times \left( {\bf f}_{f}^{~b} + {\bf f}_{r}^{~b} \right)
\end{equation}
where ${\bf r}_{MR}$ is the vector from the center of mass to the center of
resistance. Instead of integrating the angular velocity in lab-fixed
frame, we consider the equation of motion for the angular momentum
(${\bf j}$) in the body-fixed frame
\begin{equation}
\frac{\partial}{\partial t}{\bf j}(t) = \tau^{~b}(t)
\end{equation}
Embedding the friction and random forces into the the total force and
torque, one can integrate the Langevin equations of motion for a rigid
body of arbitrary shape in a velocity-Verlet style 2-part algorithm,
where $h = \delta t$:

{\tt move A:}
\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), \\
%
{\bf Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
    (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
\end{align*}
In this context, $\overleftrightarrow{\mathsf{I}}$ is the diagonal
moment of inertia tensor, and 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 ($\mathbf{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}
\mathbf{Q}(t) & \leftarrow & \mathbf{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.  The system forces
and torques are derivatives of the total potential energy function
($U$) with respect to the rigid body positions (${\bf r}$) and the
columns of the transposed rotation matrix ${\bf Q}^T = \left({\bf
u}_x, {\bf u}_y, {\bf u}_z \right)$:

{\tt Forces:}
\begin{align*}
{\bf f}_{s}(t + h) & \leftarrow
    - \left(\frac{\partial U}{\partial {\bf r}}\right)_{{\bf r}(t + h)} \\
%
{\bf \tau}_{s}(t + h) &\leftarrow {\bf u}(t + h)
    \times \frac{\partial U}{\partial {\bf u}} \\
%
{\bf v}^{b}_{R}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \left({\bf v}(t+h) + {\bf \omega}(t+h) \times {\bf r}_{MR} \right) \\
%
{\bf f}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tt} \cdot
{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rt} \cdot {\bf \omega}(t+h) \\
%
{\bf \tau}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tr} \cdot
{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rr} \cdot {\bf \omega}(t+h) \\
%
Z & \leftarrow  {\tt GaussianNormal}(2 k_B T / h, 6) \\
{\bf F}_{R,r}^{b}(t+h) & \leftarrow {\bf S} \cdot Z \\
%
{\bf f}(t+h) & \leftarrow {\bf f}_{s}(t+h) + \mathbf{Q}^{T}(t+h)
\cdot \left({\bf f}_{R,f}^{~b} + {\bf f}_{R,r}^{~b} \right) \\
%
\tau(t+h) & \leftarrow \tau_{s}(t+h) + \mathbf{Q}^{T}(t+h) \cdot \left(\tau_{R,f}^{~b} + \tau_{R,r}^{~b} \right) + {\bf r}_{MR} \times \left({\bf f}_{f}(t+h) + {\bf f}_{r}(t+h) \right) \\
\tau^{~b}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \tau(t+h) \\
\end{align*}
Frictional (and random) forces and torques must be computed at the
center of resistance, so there are additional steps required to find
the body-fixed velocity (${\bf v}_{R}^{~b}$) at this location.  Mapping
the frictional and random forces at the center of resistance back to
the center of mass also introduces an additional term in the torque
one obtains at the center of mass.

Once the forces and torques have been obtained at the new time step,
the velocities can be advanced to the same time value.

{\tt move B:}
\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 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. Code
to estimate the hydrodynamic tensors using the bead and rough shell
models was also implemented in OOPSE, although a number of other
programs exist that can perform these
tasks.\cite{JoseGarciadelaTorre02012000,Garcia-de-la-Torre:2001wd,GarciadelaTorreJ2002}

 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[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 for 5 ns 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 (5 ns) and sampling (10 ns) 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{Gay1981} 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}({{\bf \hat u}_i}, {{\bf \hat u}_j}, {{\bf \hat
r}_{ij}}) = 4\epsilon ({{\bf \hat u}_i}, {{\bf \hat u}_j},
{{\bf \hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u
}_i},
{{\bf \hat u}_j}, {{\bf \hat r}_{ij}})+\sigma_0}\right)^{12}
-\left(\frac{\sigma_0}{r_{ij}-\sigma({{\bf \hat u}_i}, {{\bf \hat u}_j},
{{\bf \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,SunX._jp0762020} 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}
which converges much more rapidly in molecular dynamics simulations
than the traditional Green-Kubo formula.

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 may be differences between the three
diffusion constants at short times, 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 with ``stick'' boundary conditions is well known,
and is given in Eqs. (\ref{eq:StokesTranslation}) and
(\ref{eq:StokesRotation}).  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 \rho$) 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(s),
\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 - s^2)
G(s) - 1}{1 - s^4} \right\}.
\label{ThetaPerrin}
\end{equation}
In these expressions, $G(s)$ is a function of the axial ratio
($s = b / a$), which for prolate ellipsoids, is
\begin{equation}
G(s) = (1- s^2)^{-1/2} \ln \left\{ \frac{1 + (1 - s^2)^{1/2}}{s} \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 {\it 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. 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 tensor 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]{RoughShellReduced}
\caption[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 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).  

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{SunX._jp0762020,Ayton01} A reference system composed
of a single lipid rigid body embedded in a sea of 1929 solvent
particles was created and run under a microcanonical ensemble.  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 by
Eq. (\ref{introEquation:ResistanceTensorArbitraryOrigin}), we
identified the center of resistance, ${\bf r} = $(0 \AA, 0 \AA, 1.46
\AA).

The translational diffusion constants and rotational correlation times
obtained using the Langevin rigid-body integrator (and the
hydrodynamic tensor) are essentially quantitative when compared with
the explicit solvent simulations for this model system.  

\subsection{Summary of comparisons with explicit solvent simulations}
The Langevin rigid-body integrator we have developed is a reliable way
to replace explicit solvent simulations in cases where the detailed
solute-solvent interactions do not greatly impact the behavior of the
solute.  As such, it has the potential to greatly increase the length
and time scales of coarse grained simulations of large solvated
molecules.  In cases where the dielectric screening of the solvent, or
specific solute-solvent interactions become important for structural
or dynamic features of the solute molecule, this integrator may be
less useful.  However, for the kinds of coarse-grained modeling that
have become popular in recent years (ellipsoidal side chains, rigid
bodies, and molecular-scale models), this integrator may prove itself
to be quite valuable.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{graph}
\caption[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 obtained
from: Stokes' law (sphere), and Refs. \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  & 1.41 & &      & 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}

To test the accuracy and efficiency of the new integrator, we applied
it 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 earlier molecular
dynamics studies on lipid bilayers which had used spherical
(Lennard-Jones) solvent particles and moderate (480 solvated lipid
molecules) system sizes.\cite{SunX._jp0762020} the solvent molecules
were excluded from the system and the box was replicated three times
in the x- and y- axes (giving a single simulation cell containing 4320
lipids).  The experimental value for the viscosity of water at 20C
($\eta = 1.00$ cp) was used with the Langevin integrator 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 simulation run time of 30 ns was sampled to
calculate structural properties.  Fig. \ref{fig:bilayer} shows the
configuration of the system after 30 ns.  Structural properties of the
bilayer (e.g. the head and body $P_2$ order parameters) are nearly
identical to those obtained via solvated molecular dynamics. The
ripple structure remained stable during the entire trajectory.
Compared with using explicit bead-model solvent molecules, the 30 ns
trajectory for 4320 lipids with the Langevin integrator is now {\it
faster} on the same hardware than the same length trajectory was for
the 480-lipid system previously studied.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{bilayer}
\caption[A snapshot of a bilayer composed of 4320 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 algorithm for carrying out Langevin dynamics
simulations on complex rigid bodies by incorporating the hydrodynamic
resistance tensors for arbitrary shapes into a stable and efficient
integration scheme.  The integrator gives quantitative agreement with
both analytic and approximate hydrodynamic theories, and works
reasonably well at reproducing the solute dynamical properties
(diffusion constants, and orientational relaxation times) from
explicitly-solvated simulations.  For the cases where there are
discrepancies between our Langevin integrator and the explicit solvent
simulations, two features of molecular simulations help explain the
differences.

First, the use of ``stick'' boundary conditions for molecular-sized
solutes in a sea of similarly-sized solvent particles may be
problematic.  We are certainly not the first group to notice this
difference between hydrodynamic theories and explicitly-solvated
molecular
simulations.\cite{Schmidt:2004fj,Schmidt:2003kx,Ravichandran:1999fk,TANG:1993lr}
The problem becomes particularly noticable in both the translational
diffusion of the spherical particles and the rotational diffusion of
the ellipsoids.  In both of these cases it is clear that the
approximations that go into hydrodynamics are the source of the error,
and not the integrator itself.

Second, in the case of structures which have substantial surface area
that is inaccessible to solvent particles, the hydrodynamic theories
(and the Langevin integrator) may overestimate the effects of solvent
friction because they overestimate the exposed surface area of the
rigid body.  This is particularly noticable in the rotational
diffusion of the dumbbell model.  We believe that given a solvent of
known radius, it may be possible to modify the rough shell approach to
place beads on solvent-accessible surface, instead of on the geometric
surface defined by the van der Waals radii of the components of the
rigid body.  Further work to confirm the behavior of this new
approximation is ongoing.

\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{doublespace}
\end{document}
Revision:	3391
Committed:	Wed Apr 30 16:16:25 2008 UTC (17 years, 6 months ago) by gezelter
Content type:	application/x-tex
File size:	67647 byte(s)
Log Message:	Nearing the end