--- trunk/langevin/langevin.tex	2008/01/02 21:06:31	3298
+++ trunk/langevin/langevin.tex	2008/02/29 22:02:46	3352
@@ -1,503 +1,612 @@
 %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
-%\documentclass[aps,prb,preprint]{revtex4}
+%\documentclass[aps,prb,preprint,amsmath,amssymb,endfloats]{revtex4}
 \documentclass[11pt]{article}
-\usepackage{endfloat}
-\usepackage{amsmath,bm}
+\usepackage{amsmath}
 \usepackage{amssymb}
 \usepackage{times}
 \usepackage{mathptmx}
 \usepackage{setspace}
-\usepackage{tabularx}
+\usepackage{endfloat}
+%\usepackage{tabularx}
 \usepackage{graphicx}
-\usepackage{booktabs}
-\usepackage{bibentry}
-\usepackage{mathrsfs}
+%\usepackage{booktabs}
+%\usepackage{bibentry}
+%\usepackage{mathrsfs}
 \usepackage[ref]{overcite}
 \pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
 \evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
 9.0in \textwidth 6.5in \brokenpenalty=10000
-\renewcommand{\baselinestretch}{1.2}
-\renewcommand\citemid{\ } % no comma in optional reference note
 
+\renewcommand\citemid{\ } % no comma in optional referenc note
+
 \begin{document}
 
-\title{An algorithm for performing Langevin dynamics on rigid bodies of arbitrary shape }
+\title{Langevin dynamics for rigid bodies of arbitrary shape}
 
-\author{Teng Lin, Xiuquan Sun and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail:
-gezelter@nd.edu} \\
-Department of Chemistry and Biochemistry\\
+\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}
+\maketitle 
 
+
+
+\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 symplectic 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
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
+\begin{doublespace}
+
 \section{Introduction}
 
 %applications of langevin dynamics
-As alternative to Newtonian dynamics, Langevin dynamics, which
-mimics a simple heat bath with stochastic and dissipative forces,
-has been applied in a variety of studies. The stochastic treatment
-of the solvent enables us to carry out substantially longer time
-simulations. Implicit solvent Langevin dynamics simulations of
-met-enkephalin not only outperform explicit solvent simulations for
-computational efficiency, but also agrees very well with explicit
-solvent simulations for dynamical properties.\cite{Shen2002}
-Recently, applying Langevin dynamics with the UNRES model, Liow and
-his coworkers suggest that protein folding pathways can be possibly
-explored within a reasonable amount of time.\cite{Liwo2005} The
-stochastic nature of the Langevin dynamics also enhances the
-sampling of the system and increases the probability of crossing
-energy barriers.\cite{Banerjee2004, Cui2003} Combining Langevin
-dynamics with Kramers's theory, Klimov and Thirumalai identified
-free-energy barriers by studying the viscosity dependence of the
-protein folding rates.\cite{Klimov1997} In order to account for
-solvent induced interactions missing from implicit solvent model,
-Kaya incorporated desolvation free energy barrier into implicit
-coarse-grained solvent model in protein folding/unfolding studies
-and discovered a higher free energy barrier between the native and
-denatured states. Because of its stability against noise, Langevin
-dynamics is very suitable for studying remagnetization processes in
-various systems.\cite{Palacios1998,Berkov2002,Denisov2003} For
-instance, the oscillation power spectrum of nanoparticles from
-Langevin dynamics simulation has the same peak frequencies for
-different wave vectors, which recovers the property of magnetic
-excitations in small finite structures.\cite{Berkov2005a}
+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.
 
-%review rigid body dynamics
-Rigid bodies are frequently involved in the modeling of different
-areas, from engineering, physics, to chemistry. For example,
-missiles and vehicle are usually modeled by rigid bodies.  The
-movement of the objects in 3D gaming engine or other physics
-simulator is governed by the rigid body dynamics. In molecular
-simulation, rigid body is used to simplify the model in
-protein-protein docking study{\cite{Gray2003}}.
+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}
 
-It is very important to develop stable and efficient methods to
-integrate the equations of motion for orientational degrees of
-freedom. Euler angles are the natural choice to describe the
-rotational degrees of freedom. However, due to $\frac {1}{sin
-\theta}$ singularities, the numerical integration of corresponding
-equations of these motion is very inefficient and inaccurate.
-Although an alternative integrator using multiple sets of Euler
-angles can overcome this difficulty\cite{Barojas1973}, the
-computational penalty and the loss of angular momentum conservation
-still remain. A singularity-free representation utilizing
-quaternions was developed by Evans in 1977.\cite{Evans1977}
-Unfortunately, this approach used a nonseparable Hamiltonian
-resulting from the quaternion representation, which prevented the
-symplectic algorithm from being utilized. Another different approach
-is to apply holonomic constraints to the atoms belonging to the
-rigid body. Each atom moves independently under the normal forces
-deriving from potential energy and constraint forces which are used
-to guarantee the rigidness. However, due to their iterative nature,
-the SHAKE and Rattle algorithms also converge very slowly when the
-number of constraints increases.\cite{Ryckaert1977, Andersen1983}
+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}
 
-A break-through in geometric literature suggests that, in order to
-develop a long-term integration scheme, one should preserve the
-symplectic structure of the propagator. By introducing a conjugate
-momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's
-equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was
-proposed to evolve the Hamiltonian system in a constraint manifold
-by iteratively satisfying the orthogonality constraint $Q^T Q = 1$.
-An alternative method using the quaternion representation was
-developed by Omelyan.\cite{Omelyan1998} However, both of these
-methods are iterative and inefficient. In this section, we descibe a
-symplectic Lie-Poisson integrator for rigid bodies developed by
-Dullweber and his coworkers\cite{Dullweber1997} in depth.
+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.
 
-%review langevin/browninan dynamics for arbitrarily shaped rigid body
-Combining Langevin or Brownian dynamics with rigid body dynamics,
-one can study slow processes in biomolecular systems. Modeling DNA
-as a chain of rigid beads, which are subject to harmonic potentials
-as well as excluded volume potentials, Mielke and his coworkers
-discovered rapid superhelical stress generations from the stochastic
-simulation of twin supercoiling DNA with response to induced
-torques.\cite{Mielke2004} Membrane fusion is another key biological
-process which controls a variety of physiological functions, such as
-release of neurotransmitters \textit{etc}. A typical fusion event
-happens on the time scale of a millisecond, which is impractical to
-study using atomistic models with newtonian mechanics. With the help
-of coarse-grained rigid body model and stochastic dynamics, the
-fusion pathways were explored by many
-researchers.\cite{Noguchi2001,Noguchi2002,Shillcock2005} Due to the
-difficulty of numerical integration of anisotropic rotation, most of
-the rigid body models are simply modeled using spheres, cylinders,
-ellipsoids or other regular shapes in stochastic simulations. In an
-effort to account for the diffusion anisotropy of arbitrary
-particles, Fernandes and de la Torre improved the original Brownian
-dynamics simulation algorithm\cite{Ermak1978,Allison1991} by
+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 simple rotation evolution scheme consisting of three
-consecutive rotations.\cite{Fernandes2002} Unfortunately, unexpected
-errors and biases are introduced into the system due to the
-arbitrary order of applying the noncommuting rotation
-operators.\cite{Beard2003} Based on the observation the momentum
-relaxation time is much less than the time step, one may ignore the
-inertia in Brownian dynamics. However, the assumption of zero
-average acceleration is not always true for cooperative motion which
-is common in protein motion. An inertial Brownian dynamics (IBD) was
-proposed to address this issue by adding an inertial correction
-term.\cite{Beard2000} As a complement to IBD which has a lower bound
-in time step because of the inertial relaxation time, long-time-step
-inertial dynamics (LTID) can be used to investigate the inertial
-behavior of the polymer segments in low friction
-regime.\cite{Beard2000} LTID can also deal with the rotational
-dynamics for nonskew bodies without translation-rotation coupling by
-separating the translation and rotation motion and taking advantage
-of the analytical solution of hydrodynamics properties. However,
-typical nonskew bodies like cylinders and ellipsoids are inadequate
-to represent most complex macromolecule assemblies. These intricate
-molecules have been represented by a set of beads and their
-hydrodynamic properties can be calculated using variants on the
-standard hydrodynamic interaction tensors.
+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 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.
 
-The goal of the present work is to develop a Langevin dynamics
-algorithm for arbitrary-shaped rigid particles by integrating the
-accurate estimation of friction tensor from hydrodynamics theory
-into the sophisticated rigid body dynamics algorithms.
+\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.
 
-\section{Computational Methods{\label{methodSec}}}
+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} 
 
-\subsection{\label{introSection:frictionTensor}Friction Tensor}
-Theoretically, the friction kernel can be determined using the
-velocity autocorrelation function. However, this approach becomes
-impractical when the system becomes more and more complicated.
-Instead, various approaches based on hydrodynamics have been
-developed to calculate the friction coefficients. In general, the
-friction tensor $\Xi$ is a $6\times 6$ matrix given by
-\[
-\Xi  = \left( {\begin{array}{*{20}c}
-   {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
-   {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
-\end{array}} \right).
-\]
-Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are $3 \times 3$
-translational friction tensor and rotational resistance (friction)
-tensor respectively, while ${\Xi^{tr} }$ is translation-rotation
-coupling tensor and $ {\Xi^{rt} }$ is rotation-translation coupling
-tensor. When a particle moves in a fluid, it may experience friction
-force or torque along the opposite direction of the velocity or
-angular velocity,
-\[
+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, 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 an
+accurate estimate of the 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
+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}
- F_R  \\
- \tau _R  \\
- \end{array} \right) =  - \left( {\begin{array}{*{20}c}
-   {\Xi ^{tt} } & {\Xi ^{rt} }  \\
-   {\Xi ^{tr} } & {\Xi ^{rr} }  \\
-\end{array}} \right)\left( \begin{array}{l}
- v \\
- w \\
- \end{array} \right)
-\]
-where $F_r$ is the friction force and $\tau _R$ is the friction
-torque.
+ \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 with slip boundary conditions, the
-translational and rotational friction constant can be calculated
-from Stoke's law,
-\[
-\Xi ^{tt}  = \left( {\begin{array}{*{20}c}
-   {6\pi \eta R} & 0 & 0  \\
-   0 & {6\pi \eta R} & 0  \\
-   0 & 0 & {6\pi \eta R}  \\
-\end{array}} \right)
-\]
+For a  spherical body under ``stick'' boundary conditions,
+the translational and rotational friction tensors can be calculated
+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
-\[
-\Xi ^{rr}  = \left( {\begin{array}{*{20}c}
-   {8\pi \eta R^3 } & 0 & 0  \\
-   0 & {8\pi \eta R^3 } & 0  \\
-   0 & 0 & {8\pi \eta R^3 }  \\
-\end{array}} \right)
-\]
-where $\eta$ is the viscosity of the solvent and $R$ is the
-hydrodynamic radius.
+\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 hydrodynamics theory,
+widely used as references for developing new hydrodynamic 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}
-\[
-\frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
-}} = 1
-\]
-where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
-due to the complexity of the elliptic integral, only the ellipsoid
-with the restriction of two axes being equal, \textit{i.e.}
-prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
-exactly. Introducing an elliptic integral parameter $S$ for prolate
-ellipsoids :
-\[
-S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
-} }}{b},
-\]
-and oblate ellipsoids:
-\[
-S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
-}}{a},
-\]
-one can write down the translational and rotational resistance
-tensors
-\begin{eqnarray*}
- \Xi _a^{tt}  & = & 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}}. \\
- \Xi _b^{tt}  & = & \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S +
- 2a}},
-\end{eqnarray*}
-and
-\begin{eqnarray*}
- \Xi _a^{rr} & = & \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}}, \\
- \Xi _b^{rr} & = & \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}}.
-\end{eqnarray*}
+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}}
+There is no analytical solution 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}.
 
-Unlike spherical and other simply shaped molecules, there is no
-analytical solution for the friction tensor for arbitrarily shaped
-rigid molecules. The ellipsoid of revolution model and general
-triaxial ellipsoid model have been used to approximate the
-hydrodynamic properties of rigid bodies. However, since the mapping
-from all possible ellipsoidal spaces, $r$-space, to all possible
-combination of rotational diffusion coefficients, $D$-space, is not
-unique\cite{Wegener1979} as well as the intrinsic coupling between
-translational and rotational motion of rigid bodies, general
-ellipsoids are not always suitable for modeling arbitrarily shaped
-rigid molecules. A number of studies have been devoted to
-determining the friction tensor for irregularly shaped rigid bodies
-using more advanced methods where the molecule of interest was
-modeled by a combinations of spheres\cite{Carrasco1999} and the
-hydrodynamics properties of the molecule can be calculated using the
-hydrodynamic interaction tensor. Let us consider a rigid assembly of
-$N$ beads immersed in a continuous medium. Due to hydrodynamic
-interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different
-than its unperturbed velocity $v_i$,
-\[
-v'_i  = v_i  - \sum\limits_{j \ne i} {T_{ij} F_j }
-\]
-where $F_i$ is the frictional force, and $T_{ij}$ is the
-hydrodynamic interaction tensor. The friction force of $i$th bead is
-proportional to its ``net'' velocity
+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}
-F_i  = \zeta _i v_i  - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
-\label{introEquation:tensorExpression}
+{\bf v}'_i  = {\bf v}_i  - \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }
 \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}
+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}
-T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
-R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor}
+{\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}
-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
+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}
-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].
+{\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}
-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
+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}
-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}
+{\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}
-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
+${\bf B}_{ij}$ is a version of the hydrodynamic tensor which includes the
+self-contributions for spheres,
 \begin{equation}
-B = \left( {\begin{array}{*{20}c}
-   {B_{11} } &  \ldots  & {B_{1N} }  \\
-    \vdots  &  \ddots  &  \vdots   \\
-   {B_{N1} } &  \cdots  & {B_{NN} }  \\
-\end{array}} \right),
+{\bf B}_{ij}  = \delta _{ij} \frac{{\bf I}}{{6\pi \eta R_{ij}}} + (1 - \delta_{ij}
+){\bf T}_{ij}
 \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 $\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 resistance tensor at
-arbitrary origin $O$ can be written as
+bead $i$ and origin $O$, the elements of the resistance tensor (at the
+arbitrary origin $O$) can be written as
 \begin{eqnarray}
- \Xi _{}^{tt}  & = & \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\
- \Xi _{}^{tr}  & = & \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
- \Xi _{}^{rr}  & = &  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j. \notag \\
 \label{introEquation:ResistanceTensorArbitraryOrigin}
+ \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 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,
+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}
-\Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
-\label{introEquation:definitionCR}
+V = \frac{4 \pi}{3} \sum_{i=1}^{N} \rho_i^3,
 \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
+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}
-\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}
+\Xi^{tr}  = \left(\Xi^{tr}\right)^T
+\label{introEquation:definitionCR}
 \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
+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$.
 
-\subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}
+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}
-M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)}  + F_{r,i} (t)
+{\bf M} \dot{{\bf V}}(t) = {\bf F}_{s}(t) +
+{\bf F}_{f}(t)  + {\bf F}_{r}(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$
+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}
-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) \\
+{\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 $F_{r,i}^l$ is a Gaussian stochastic variable
-with zero mean and variance
+while the random force, ${\bf F}_{r}$, 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}
+\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}
-The equation of motion for $v_i$ can be written as
+$\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}
-m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
-f_{r,i}^l (t)
+\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}
-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
+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}
-\tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b
+m \dot{{\bf v}} (t) =  {\bf f}_{s}(t) + {\bf f}_{f}(t) +
+{\bf f}_{r}(t)
 \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
+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}
-\dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t)
-+ \tau _{r,i}^b(t)
+\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}
-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$:
+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 moveA:}
+{\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), \\
@@ -506,26 +615,27 @@ $h= \delta 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), \\
+    + \frac{h}{2} {\bf \tau}^{~b}(t), \\
 %
-\mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
+{\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, the $\mathrm{rotate}$ function is the reversible
-product of the three body-fixed rotations,
+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 ($\mathsf{Q}$) and the body-fixed
+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}
-\mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
+\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}
@@ -547,23 +657,48 @@ calculated at the new positions and orientations
 \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
+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 doForces:}
+{\tt Forces:}
 \begin{align*}
-{\bf f}(t + h) &\leftarrow
-    - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\
+{\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 V}{\partial {\bf u}}, \\
+{\bf \tau}_{s}(t + h) &\leftarrow {\bf u}(t + h)
+    \times \frac{\partial U}{\partial {\bf u}} \\
 %
-{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h)
-    \cdot {\bf \tau}^s(t + h).
+{\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 moveB:}
+{\tt move B:}
 \begin{align*}
 {\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2
 \right)
@@ -571,84 +706,124 @@ the velocities can be advanced to the same time value.
 %
 {\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t + h / 2
 \right)
-    + \frac{h}{2} {\bf \tau}^b(t + h) .
+    + \frac{h}{2} {\bf \tau}^{~b}(t + h) .
 \end{align*}
 
-\section{Results and Discussion}
+\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.
 
-The Langevin algorithm described in previous section has been
-implemented in {\sc oopse}\cite{Meineke2005} and applied to studies
-of the static and dynamic properties in several systems.
-
-\subsection{Temperature Control}
-
-As shown in Eq.~\ref{randomForce}, random collisions associated with
-the solvent's thermal motions is controlled by the external
-temperature. The capability to maintain the temperature of the whole
-system was usually used to measure the stability and efficiency of
-the algorithm. In order to verify the stability of this new
-algorithm, a series of simulations are performed on system
-consisiting of 256 SSD water molecules with different viscosities.
-The initial configuration for the simulations is taken from a 1ns
-NVT simulation with a cubic box of 19.7166~\AA. All simulation are
-carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns
-with reference temperature at 300~K. The average temperature as a
-function of $\eta$ is shown in Table \ref{langevin:viscosity} where
-the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 -
-1$ poise. The better temperature control at higher viscosity can be
-explained by the finite size effect and relative slow relaxation
-rate at lower viscosity regime.
-\begin{table}
-\caption{AVERAGE TEMPERATURES FROM LANGEVIN DYNAMICS SIMULATIONS OF
-SSD WATER MOLECULES WITH REFERENCE TEMPERATURE AT 300~K.}
-\label{langevin:viscosity}
+\begin{table*}
+\begin{minipage}{\linewidth}
 \begin{center}
-\begin{tabular}{lll}
-  \hline
-  $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\
-  \hline
-  1    & 300.47 & 10.99 \\
-  0.1  & 301.19 & 11.136 \\
-  0.01 & 303.04 & 11.796 \\
-  \hline
+\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{table}
+\end{minipage}
+\end{table*}
 
-Another set of calculations were performed to study the efficiency of
-temperature control using different temperature coupling schemes.
-The starting configuration is cooled to 173~K and evolved using NVE,
-NVT, and Langevin dynamic with time step of 2 fs.
-Fig.~\ref{langevin:temperature} shows the heating curve obtained as
-the systems reach equilibrium. The orange curve in
-Fig.~\ref{langevin:temperature} represents the simulation using
-Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps
-which gives reasonable tight coupling, while the blue one from
-Langevin dynamics with viscosity of 0.1 poise demonstrates a faster
-scaling to the desire temperature. When $ \eta = 0$, Langevin dynamics becomes normal
-NVE (see orange curve in Fig.~\ref{langevin:temperature}) which
-loses the temperature control ability.
-
 \begin{figure}
 \centering
-\includegraphics[width=\linewidth]{temperature.pdf}
-\caption[Plot of Temperature Fluctuation Versus Time]{Plot of
-temperature fluctuation versus time.} \label{langevin:temperature}
+\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{Langevin Dynamics simulation {\it vs} NVE simulations}
+\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.
 
-To validate our langevin dynamics simulation. We performed several NVE
-simulations with explicit solvents for different shaped
-molecules. There are one solute molecule and 1929 solvent molecules in
-NVE simulation. The parameters are shown in table
-\ref{tab:parameters}. The force field between spheres is standard
-Lennard-Jones, and ellipsoids interact with other ellipsoids and
-spheres with generalized Gay-Berne potential. All simulations are
-carried out at 300 K and 1 Atm. The time step is 25 ns, and a
-switching function was applied to all potentials to smoothly turn off
-the interactions between a range of $22$ and $25$ \AA.  The switching
-function was the standard (cubic) function,
+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}
@@ -660,127 +835,405 @@ We have computed translational diffusion constants for
         \end{cases}
 \label{eq:switchingFunc}
 \end{equation}
-We have computed translational diffusion constants for lipid molecules
-from the mean-square displacement,
+
+To measure shear viscosities from our microcanonical simulations, we
+used the Einstein form of the pressure correlation function,\cite{hess:209}
 \begin{equation}
-D = \lim_{t\rightarrow \infty} \frac{1}{6 t} \langle {|\left({\bf r}_{i}(t) - {\bf r}_{i}(0) \right)|}^2 \rangle,
+\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}
-of the solute molecules. Translational diffusion constants for the
-different shaped molecules are shown in table
-\ref{tab:translation}.  We have also computed orientational correlation
-times for different shaped molecules from fits of the second-order
-Legendre polynomial correlation function,
+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}
-C_{\ell}(t)  =  \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
-\mu}_{i}(0) \right) 
+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}
-the results are shown in table \ref{tab:rotation}. We used einstein
-format of the pressure correlation function,
+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}
-C_{\ell}(t)  =  \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
-\mu}_{i}(0) \right) 
+D_r = \frac{1}{3} \left(D_1 + D_2 + D_3 \right)
 \end{equation}
-to estimate the viscosity of the systems from NVE simulations. The
-viscosity can also be calculated by Green-Kubo pressure correlaton
-function,
+and 
 \begin{equation}
-C_{\ell}(t)  =  \langle P_{\ell}\left({\bf \mu}_{i}(t) \cdot {\bf
-\mu}_{i}(0) \right) 
+\Delta = \left( (D_1 - D_2)^2 + (D_3 - D_1 )(D_3 - D_2)\right)^{1/2}
 \end{equation}
-However, this method converges slowly, and the statistics are not good
-enough to give us a very accurate value. The langevin dynamics
-simulations for different shaped molecules are performed at the same
-conditions as the NVE simulations with viscosity estimated from NVE
-simulations. To get better statistics, 1024 non-interacting solute
-molecules are put into one simulation box for each langevin
-simulation, this is equal to 1024 simulations for single solute
-systems. The diffusion constants and rotation relaxation times for
-different shaped molecules are shown in table \ref{tab:translation}
-and \ref{tab:rotation} to compare to the results calculated from NVE
-simulations. The theoretical values for sphere is calculated from the
-Stokes-Einstein law, the theoretical values for ellipsoid is
-calculated from Perrin's fomula, the theoretical values for dumbbell
-is calculated from StinXX and Davis theory. The exact method is
-applied to the langevin dynamics simulations for sphere and ellipsoid,
-the bead model is applied to the simulation for dumbbell molecule, and
-the rough shell model is applied to ellipsoid, dumbbell, banana and
-lipid molecules. The results from all the langevin dynamics
-simulations, including exact, bead model and rough shell, match the
-theoretical values perfectly for all different shaped molecules. This
-indicates that our simulation package for langevin dynamics is working
-well. The approxiate methods ( bead model and rough shell model) are
-accurate enough for the current simulations. The goal of the langevin
-dynamics theory is to replace the explicit solvents by the friction
-forces. We compared the dynamic properties of different shaped
-molecules in langevin dynamics simulations with that in NVE
-simulations. The results are reasonable close. Overall, the
-translational diffusion constants calculated from langevin dynamics
-simulations are very close to the values from the NVE simulation. For
-sphere and lipid molecules, the diffusion constants are a little bit
-off from the NVE simulation results. One possible reason is that the
-calculation of the viscosity is very difficult to be accurate. Another
-possible reason is that although we save very frequently during the
-NVE simulations and run pretty long time simulations, there is only
-one solute molecule in the system which makes the calculation for the
-diffusion constant difficult. The sphere molecule behaves as a free
-rotor in the solvent, so there is no rotation relaxation time
-calculated from NVE simulations. The rotation relaxation time is not
-very close to the NVE simulations results. The banana and lipid
-molecules match the NVE simulations results pretty well. The mismatch
-between langevin dynamics and NVE simulation for ellipsoid is possibly
-caused by the slip boundary condition. For dumbbell, the mismatch is
-caused by the size of the solvent molecule is pretty large compared to
-dumbbell molecule in NVE simulations.
+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.
 
-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.
+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.
 
-\begin{table*}
-\begin{minipage}{\linewidth}
-\begin{center}
-\caption{}
-\begin{tabular}{llccccccc}
-\hline 
-  & & Sphere & Ellipsoid & Dumbbell(2 spheres) & Banana(3 ellpsoids) &
-Lipid(head) & lipid(tail) & Solvent \\
-\hline 
-$d$ (\AA) & & 6.5 & 4.6  & 6.5 &  4.2 & 6.5 & 4.6 & 4.7 \\
-$l$ (\AA) & & $= d$ & 13.8 & $=d$ & 11.2 & $=d$ & 13.8 & 4.7 \\
-$\epsilon^s$ (kcal/mol) & & 0.8 & 0.8 & 0.8 & 0.8 & 0.185 & 0.8 & 0.8 \\
-$\epsilon_r$ (well-depth aspect ratio)& & 1 & 0.2 & 1 & 0.2 & 1 & 0.2 & 1 \\
-$m$ (amu) & & 190 & 200 & 190 & 240 & 196 & 760 & 72.06 \\
-%$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$) & & & & \\
-%\multicolumn{2}c{$I_{xx}$} & 1125 & 45000 & N/A \\
-%\multicolumn{2}c{$I_{yy}$} & 1125 & 45000 & N/A \\
-%\multicolumn{2}c{$I_{zz}$} &  0 &    9000 & N/A \\
-%$\mu$ (Debye) & & varied & 0 & 0 \\
-\end{tabular}
-\label{tab:parameters}
-\end{center}
-\end{minipage}
-\end{table*}
+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]{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).  
+
+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 in
+Sec.~\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[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{}
-\begin{tabular}{lccccc}
+\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
- & & & & &Translation \\
-\hline 
- & NVE &  & Theoretical & Langevin & \\
-\hline 
- & $\eta$ & D & D & method & D \\
+ & \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 & 3.480159e-03 & 1.643135e-04 & 1.942779e-04 & exact & 1.982283e-04 \\
-ellipsoid & 2.551262e-03 & 2.437492e-04 & 2.335756e-04 & exact & 2.374905e-04 \\
- & 2.551262e-03  & 2.437492e-04 & 2.335756e-04 & rough shell & 2.284088e-04 \\
-dumbell & 2.41276e-03  & 2.129432e-04 & 2.090239e-04 & bead model & 2.148098e-04 \\
- & 2.41276e-03 & 2.129432e-04 & 2.090239e-04 & rough shell & 2.013219e-04 \\
-banana & 2.9846e-03 & 1.527819e-04 &  & rough shell & 1.54807e-04 \\
-lipid & 3.488661e-03 & 0.9562979e-04 &  & rough shell & 1.320987e-04 \\
+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}
@@ -790,128 +1243,115 @@ lipid & 3.488661e-03 & 0.9562979e-04 &  & rough shell 
 \begin{table*}
 \begin{minipage}{\linewidth}
 \begin{center}
-\caption{}
-\begin{tabular}{lccccc}
+\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
- & & & & &Rotation \\
-\hline 
- & NVE &  & Theoretical & Langevin & \\
-\hline 
- & $\eta$ & $\tau_0$ & $\tau_0$ & method & $\tau_0$ \\
+ & \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 & 3.480159e-03 &  & 1.208178e+04 & exact & 1.20628e+04 \\
-ellipsoid & 2.551262e-03 & 4.66806e+04 & 2.198986e+04 & exact & 2.21507e+04 \\
- & 2.551262e-03 & 4.66806e+04 & 2.198986e+04 & rough shell & 2.21714e+04 \\
-dumbell & 2.41276e-03 & 1.42974e+04 &  & bead model & 7.12435e+04 \\
- & 2.41276e-03 & 1.42974e+04 &  & rough shell & 7.04765e+04 \\
-banana & 2.9846e-03 & 6.38323e+04 &  & rough shell & 7.0945e+04 \\
-lipid & 3.488661e-03 & 7.79595e+04 &  & rough shell & 7.78886e+04 \\
+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*}
 
-Langevin dynamics simulations are applied to study the formation of
-the ripple phase of lipid membranes. The initial configuration is
-taken from our molecular dynamics studies on lipid bilayers with
-lennard-Jones sphere solvents. The solvent molecules are excluded from
-the system, the experimental value of water viscosity is applied to
-mimic the heat bath. Fig. XXX is the snapshot of the stable
-configuration of the system, the ripple structure stayed stable after
-100 ns run. The efficiency of the simulation is increased by one order
-of magnitude.
+\section{Application: A rigid-body lipid bilayer}
 
-\subsection{Langevin Dynamics of Banana Shaped Molecules}
+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.
 
-In order to verify that Langevin dynamics can mimic the dynamics of
-the systems absent of explicit solvents, we carried out two sets of
-simulations and compare their dynamic properties.
-Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation
-made of 256 pentane molecules and two banana shaped molecules at
-273~K. It has an equivalent implicit solvent system containing only
-two banana shaped molecules with viscosity of 0.289 center poise. To
-calculate the hydrodynamic properties of the banana shaped molecule,
-we created a rough shell model (see Fig.~\ref{langevin:roughShell}),
-in which the banana shaped molecule is represented as a ``shell''
-made of 2266 small identical beads with size of 0.3 \AA on the
-surface. Applying the procedure described in
-Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
-identified the center of resistance at (0 $\rm{\AA}$, 0.7482 $\rm{\AA}$,
--0.1988 $\rm{\AA}$), as well as the resistance tensor,
-\[
-\left( {\begin{array}{*{20}c}
-0.9261 & 0 & 0&0&0.08585&0.2057\\
-0& 0.9270&-0.007063& 0.08585&0&0\\
-0&-0.007063&0.7494&0.2057&0&0\\
-0&0.0858&0.2057& 58.64& 0&0\\
-0.08585&0&0&0&48.30&3.219&\\
-0.2057&0&0&0&3.219&10.7373\\
-\end{array}} \right).
-\]
-where the units for translational, translation-rotation coupling and rotational tensors are $\frac{kcal \cdot fs}{mol \cdot \rm{\AA}^2}$, $\frac{kcal \cdot fs}{mol \cdot \rm{\AA} \cdot rad}$ and $\frac{kcal \cdot fs}{mol \cdot rad^2}$ respectively.
-Curves of the velocity auto-correlation functions in
-Fig.~\ref{langevin:vacf} were shown to match each other very well.
-However, because of the stochastic nature, simulation using Langevin
-dynamics was shown to decay slightly faster than MD. In order to
-study the rotational motion of the molecules, we also calculated the
-auto-correlation function of the principle axis of the second GB
-particle, $u$. The discrepancy shown in Fig.~\ref{langevin:uacf} was
-probably due to the reason that we used the experimental viscosity directly instead of calculating bulk viscosity from simulation.
-
 \begin{figure}
 \centering
-\includegraphics[width=\linewidth]{roughShell.pdf}
-\caption[Rough shell model for banana shaped molecule]{Rough shell
-model for banana shaped molecule.} \label{langevin:roughShell}
+\includegraphics[width=\linewidth]{bilayer}
+\caption[Snapshot of a bilayer of rigid-body models for lipids]{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}
 
-\begin{figure}
-\centering
-\includegraphics[width=\linewidth]{twoBanana.pdf}
-\caption[Snapshot from Simulation of Two Banana Shaped Molecules and
-256 Pentane Molecules]{Snapshot from simulation of two Banana shaped
-molecules and 256 pentane molecules.} \label{langevin:twoBanana}
-\end{figure}
-
-\begin{figure}
-\centering
-\includegraphics[width=\linewidth]{vacf.pdf}
-\caption[Plots of Velocity Auto-correlation Functions]{Velocity
-auto-correlation functions of NVE (explicit solvent) in blue and
-Langevin dynamics (implicit solvent) in red.} \label{langevin:vacf}
-\end{figure}
-
-\begin{figure}
-\centering
-\includegraphics[width=\linewidth]{uacf.pdf}
-\caption[Auto-correlation functions of the principle axis of the
-middle GB particle]{Auto-correlation functions of the principle axis
-of the middle GB particle of NVE (blue) and Langevin dynamics
-(red).} \label{langevin:uacf}
-\end{figure}
-
 \section{Conclusions}
 
-We have presented a new Langevin algorithm by incorporating the
-hydrodynamics properties of arbitrary shaped molecules into an
-advanced symplectic integration scheme. The temperature control
-ability of this algorithm was demonstrated by a set of simulations
-with different viscosities. It was also shown to have significant
-advantage of producing rapid thermal equilibration over
-Nos\'{e}-Hoover method. Further studies in systems involving banana
-shaped molecules illustrated that the dynamic properties could be
-preserved by using this new algorithm as an implicit solvent model.
+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 an advanced symplectic
+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.
+financial support from Center of Applied Mathematics at University of
+Notre Dame.
+
+\end{doublespace}
 \newpage
 
 \bibliographystyle{jcp2}
 \bibliography{langevin}
-
 \end{document}
