[ViewVC] Diff of: group/trunk/xDissertation/ld.tex

Comparing trunk/xDissertation/ld.tex (file contents):
Revision 3336 by xsun, Wed Jan 30 16:01:02 2008 UTC vs.
Revision 3383 by xsun, Wed Apr 16 21:56:34 2008 UTC

+\chapter{\label{chap:ld}LANGEVIN DYNAMICS}
-+
-+
+Recent examples of the usefulness of Langevin simulations include a
-+
+study of met-enkephalin in which Langevin simulations predicted
-+
+dynamical properties that were large\-ly in agreement with explicit
-+
+solvent simulations.\cite{Shen2002} By applying Langevin dynamics with
-+
+the UNRES model, Liwo and his coworkers suggest that protein folding
-+
+pathways can be explored within a reasonable amount of
-+
+time.\cite{Liwo2005}
-+
-+
+The stochastic nature of Langevin dynamics also enhances the sampling
-+
+of the system and increases the probability of crossing energy
-+
+barriers.\cite{Cui2003,Banerjee2004} Combining Langevin dynamics with
-+
+Kramers' theory, Klimov and Thirumalai identified free-energy
-+
+barriers by studying the viscosity dependence of the protein folding
-+
+rates.\cite{Klimov1997} In order to account for solvent induced
-+
+interactions missing from the implicit solvent model, Kaya
-+
+incorporated a desolvation free energy barrier into protein
-+
+folding/unfolding studies and discovered a higher free energy barrier
-+
+between the native and denatured states.\cite{HuseyinKaya07012005}
-+
-+
+In typical LD simulations, the friction and random ($f_r$) forces on
-+
+individual atoms are taken from Stokes' law,
-+
+\begin{eqnarray}
-+
+m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + f_r(t) \notag \\
-+
+\langle f_r(t) \rangle & = & 0 \\
-+
+\langle f_r(t) f_r(t') \rangle & = & 2 k_B T \xi m \delta(t - t') \notag
-+
+\end{eqnarray}
-+
+where $\xi \approx 6 \pi \eta \rho$.  Here $\eta$ is the viscosity of the
-+
+implicit solvent, and $\rho$ is the hydrodynamic radius of the atom.
-+
-+
+The use of rigid substructures,\cite{Chun:2000fj}
-+
+coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunX._jp0762020}
-+
+and ellipsoidal representations of protein side
-+
+chains~\cite{Fogolari:1996lr} has made the use of the Stokes-Einstein
-+
+approximation problematic.  A rigid substructure moves as a single
-+
+unit with orientational as well as translational degrees of freedom.
-+
+This requires a more general treatment of the hydrodynamics than the
-+
+spherical approximation provides.  Also, the atoms involved in a rigid
-+
+or coarse-grained structure have solvent-mediated interactions with
-+
+each other, and these interactions are ignored if all atoms are
-+
+treated as separate spherical particles.  The theory of interactions
-+
+{\it between} bodies moving through a fluid has been developed over
-+
+the past century and has been applied to simulations of Brownian
-+
+motion.\cite{FIXMAN:1986lr,Ramachandran1996}
-+
-+
+In order to account for the diffusion anisotropy of complex shapes,
-+
+Fernandes and Garc\'{i}a de la Torre improved an earlier Brownian
-+
+dynamics simulation algorithm~\cite{Ermak1978,Allison1991} by
-+
+incorporating a generalized $6\times6$ diffusion tensor and
-+
+introducing a rotational evolution scheme consisting of three
-+
+consecutive rotations.\cite{Fernandes2002} Unfortunately, biases are
-+
+introduced into the system due to the arbitrary order of applying the
-+
+noncommuting rotation operators.\cite{Beard2003} Based on the
-+
+observation the momentum relaxation time is much less than the time
-+
+step, one may ignore the inertia in Brownian dynamics.  However, the
-+
+assumption of zero average acceleration is not always true for
-+
+cooperative motion which is common in proteins. An inertial Brownian
-+
+dynamics (IBD) was proposed to address this issue by adding an
-+
+inertial correction term.\cite{Beard2000} As a complement to IBD,
-+
+which has a lower bound in time step because of the inertial
-+
+relaxation time, long-time-step inertial dynamics (LTID) can be used
-+
+to investigate the inertial behavior of linked polymer segments in a
-+
+low friction regime.\cite{Beard2000} LTID can also deal with the
-+
+rotational dynamics for nonskew bodies without translation-rotation
-+
+coupling by separating the translation and rotation motion and taking
-+
+advantage of the analytical solution of hydrodynamic
-+
+properties. However, typical nonskew bodies like cylinders and
-+
+ellipsoids are inadequate to represent most complex macromolecular
-+
+assemblies. Therefore, the goal of this work is to adapt some of the
-+
+hydrodynamic methodologies developed to treat Brownian motion of
-+
+complex assemblies into a Langevin integrator for rigid bodies with
-+
+arbitrary shapes.
-+
-+
+\subsection{Rigid Body Dynamics}
-+
+Rigid bodies are frequently involved in the modeling of large
-+
+collections of particles that move as a single unit.  In molecular
-+
+simulations, rigid bodies have been used to simplify protein-protein
-+
+docking,\cite{Gray2003} and lipid bilayer
-+
+simulations.\cite{SunX._jp0762020} Many of the water models in common
-+
+use are also rigid-body
-+
+models,\cite{Jorgensen83,Berendsen81,Berendsen87} although they are
-+
+typically evolved in molecular dynamics simulations using constraints
-+
+rather than rigid body equations of motion.
-+
-+
+Euler angles are a natural choice to describe the rotational degrees
-+
+of freedom.  However, due to $\frac{1}{\sin \theta}$ singularities, the
-+
+numerical integration of corresponding equations of these motion can
-+
+become inaccurate (and inefficient).  Although the use of multiple
-+
+sets of Euler angles can overcome this problem,\cite{Barojas1973} the
-+
+computational penalty and the loss of angular momentum conservation
-+
+remain.  A singularity-free representation utilizing quaternions was
-+
+developed by Evans in 1977.\cite{Evans1977} The Evans quaternion
-+
+approach uses a nonseparable Hamiltonian, and this has prevented
-+
+symplectic algorithms from being utilized until very
-+
+recently.\cite{Miller2002}
-+
-+
+Another approach is the application of holonomic constraints to the
-+
+atoms belonging to the rigid body.  Each atom moves independently
-+
+under the normal forces deriving from potential energy and constraints
-+
+are used to guarantee rigidity. However, due to their iterative
-+
+nature, the SHAKE and RATTLE algorithms converge very slowly when the
-+
+number of constraints (and the number of particles that belong to the
-+
+rigid body) increases.\cite{Ryckaert1977,Andersen1983}
-+
-+
+In order to develop a stable and efficient integration scheme that
-+
+preserves most constants of the motion in microcanonical simulations,
-+
+symplectic propagators are necessary.  By introducing a conjugate
-+
+momentum to the rotation matrix ${\bf Q}$ and re-formulating
-+
+Hamilton's equations, a symplectic orientational integrator,
-+
+RSHAKE,\cite{Kol1997} was proposed to evolve rigid bodies on a
-+
+constraint manifold by iteratively satisfying the orthogonality
-+
+constraint ${\bf Q}^T {\bf Q} = 1$.  An alternative method using the
-+
+quaternion representation was developed by Omelyan.\cite{Omelyan1998}
-+
+However, both of these methods are iterative and suffer from some
-+
+related inefficiencies. A symplectic Lie-Poisson integrator for rigid
-+
+bodies developed by Dullweber {\it et al.}\cite{Dullweber1997} removes
-+
+most of the limitations mentioned above and is therefore the basis for
-+
+our Langevin integrator.
-+
-+
+The goal of the present work is to develop a Langevin dynamics
-+
+algorithm for ar\-bi\-trary-shaped rigid particles by integrating an
-+
+accurate estimate of the friction tensor from hydrodynamics theory
-+
+into a stable and efficient rigid body dynamics propagator.  In the
-+
+sections below, we review some of the theory of hydrodynamic tensors
-+
+developed primarily for Brownian simulations of multi-particle
-+
+systems, we then present our integration method for a set of
-+
+generalized Langevin equations of motion, and we compare the behavior
-+
+of the new Langevin integrator to dynamical quantities obtained via
-+
+explicit solvent molecular dynamics.
-+
-+
+\subsection{\label{ldintroSection:frictionTensor}The Friction Tensor}
-+
+Theoretically, a complete friction kernel for a solute particle can be
-+
+determined using the velocity autocorrelation function from a
-+
+simulation with explicit solvent molecules. However, this approach
-+
+becomes impractical when the solute becomes complex.  Instead, various
-+
+approaches based on hydrodynamics have been developed to calculate
-+
+static friction coefficients. In general, the friction tensor $\Xi$ is
-+
+a $6\times 6$ matrix given by
-+
+\begin{equation}
-+
+\Xi  = \left( \begin{array}{*{20}c}
-+
+   \Xi^{tt} & \Xi^{rt}  \\
-+
+   \Xi^{tr} & \Xi^{rr}  \\
-+
+\end{array} \right).
-+
+\end{equation}
-+
+Here, $\Xi^{tt}$ and $\Xi^{rr}$ are $3 \times 3$ translational and
-+
+rotational resistance (friction) tensors respectively, while
-+
+$\Xi^{tr}$ is translation-rotation coupling tensor and $\Xi^{rt}$ is
-+
+rotation-translation coupling tensor. When a particle moves in a
-+
+fluid, it may experience a friction force ($\mathbf{f}_f$) and torque
-+
+($\mathbf{\tau}_f$) in opposition to the velocity ($\mathbf{v}$) and
-+
+body-fixed angular velocity ($\mathbf{\omega}$),
-+
+\begin{equation}
-+
+\left( \begin{array}{l}
-+
+ \mathbf{f}_f  \\
-+
+ \mathbf{\tau}_f  \\
-+
+ \end{array} \right) =  - \left( \begin{array}{*{20}c}
-+
+   \Xi^{tt} & \Xi^{rt}  \\
-+
+   \Xi^{tr} & \Xi^{rr}  \\
-+
+\end{array} \right)\left( \begin{array}{l}
-+
+ \mathbf{v} \\
-+
+ \mathbf{\omega} \\
-+
+ \end{array} \right).
-+
+\end{equation}
-+
+For an arbitrary body moving in a fluid, Peters has derived a set of
-+
+fluctuation-dissipation relations for the friction
-+
+tensors,\cite{Peters:1999qy,Peters:1999uq,Peters:2000fk}
-+
+\begin{eqnarray}
-+
+\Xi^{tt} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
-+
+F}(0) {\bf F}(-s) \rangle_{eq} - \langle {\bf F} \rangle_{eq}^2
-+
+\right] ds \\
-+
+\notag \\
-+
+\Xi^{tr} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
-+
+F}(0) {\bf \tau}(-s) \rangle_{eq} - \langle {\bf F} \rangle_{eq}
-+
+\langle {\bf \tau} \rangle_{eq} \right] ds \\
-+
+\notag \\
-+
+\Xi^{rt} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
-+
+\tau}(0) {\bf F}(-s) \rangle_{eq} - \langle {\bf \tau} \rangle_{eq}
-+
+\langle {\bf F} \rangle_{eq} \right] ds \\
-+
+\notag \\
-+
+\Xi^{rr} & = & \frac{1}{k_B T} \int_0^\infty \left[ \langle {\bf
-+
+\tau}(0) {\bf \tau}(-s) \rangle_{eq} - \langle {\bf \tau} \rangle_{eq}^2
-+
+\right] ds
-+
+\end{eqnarray}
-+
+In these expressions, the forces (${\bf F}$) and torques (${\bf
-+
+\tau}$) are those that arise solely from the interactions of the body with
-+
+the surrounding fluid. For a single solute body in an isotropic fluid,
-+
+the average forces and torques in these expressions ($\langle {\bf F}
-+
+\rangle_{eq}$ and $\langle {\bf \tau} \rangle_{eq}$)
-+
+vanish, and one obtains the simpler force-torque correlation formulae
-+
+of Nienhuis.\cite{Nienhuis:1970lr} Molecular dynamics simulations with
-+
+explicit solvent molecules can be used to obtain estimates of the
-+
+friction tensors with these formulae. In practice, however, one needs
-+
+relatively long simulations with frequently-stored force and torque
-+
+information to compute friction tensors, and this becomes
-+
+prohibitively expensive when there are large numbers of large solute
-+
+particles.  For bodies with simple shapes, there are a number of
-+
+approximate expressions that allow computation of these tensors
-+
+without the need for expensive simulations that utilize explicit
-+
+solvent particles.
-+
-+
+\subsubsection{\label{ldintroSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}}
-+
+For a spherical body under ``stick'' boundary conditions, the
-+
+translational and rotational friction tensors can be estimated from
-+
+Stokes' law,
-+
+\begin{equation}
-+
+\label{ldeq:StokesTranslation}
-+
+\Xi^{tt}  = \left( \begin{array}{*{20}c}
-+
+   {6\pi \eta \rho} & 0 & 0  \\
-+
+& {6\pi \eta \rho} & 0  \\
-+
+& 0 & {6\pi \eta \rho}  \\
-+
+\end{array} \right)
-+
+\end{equation}
-+
+and
-+
+\begin{equation}
-+
+\label{ldeq:StokesRotation}
-+
+\Xi^{rr}  = \left( \begin{array}{*{20}c}
-+
+   {8\pi \eta \rho^3 } & 0 & 0  \\
-+
+& {8\pi \eta \rho^3 } & 0  \\
-+
+& 0 & {8\pi \eta \rho^3 }  \\
-+
+\end{array} \right)
-+
+\end{equation}
-+
+where $\eta$ is the viscosity of the solvent and $\rho$ is the
-+
+hydrodynamic radius.  The presence of the rotational resistance tensor
-+
+implies that the spherical body has internal structure and
-+
+orientational degrees of freedom that must be propagated in time.  For
-+
+non-structured spherical bodies (i.e. the atoms in a traditional
-+
+molecular dynamics simulation) these degrees of freedom do not exist.
-+
-+
+Other non-spherical shapes, such as cylinders and ellipsoids, are
-+
+widely used as references for developing new hydrodynamic theories,
-+
+because their properties can be calculated exactly. In 1936, Perrin
-+
+extended Stokes' law to general
-+
+ellipsoids,\cite{Perrin1934,Perrin1936} described in Cartesian
-+
+coordinates as
-+
+\begin{equation}
-+
+\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1.
-+
+\end{equation}
-+
+Here, the semi-axes are of lengths $a$, $b$, and $c$. Due to the
-+
+complexity of the elliptic integral, only uniaxial ellipsoids, either
-+
+prolate ($a \ge b = c$) or oblate ($a < b = c$), were solved
-+
+exactly. Introducing an elliptic integral parameter $S$ for prolate,
-+
+\begin{equation}
-+
+S = \frac{2}{\sqrt{a^2  - b^2}} \ln \frac{a + \sqrt{a^2  - b^2}}{b},
-+
+\end{equation}
-+
+and oblate,
-+
+\begin{equation}
-+
+S = \frac{2}{\sqrt {b^2  - a^2 }} \arctan \frac{\sqrt {b^2  - a^2}}{a},
-+
+\end{equation}
-+
+ellipsoids, it is possible to write down exact solutions for the
-+
+resistance tensors.  As is the case for spherical bodies, the translational,
-+
+\begin{eqnarray}
-+
+ \Xi_a^{tt}  & = & 16\pi \eta \frac{a^2  - b^2}{(2a^2  - b^2 )S - 2a}. \\
-+
+ \Xi_b^{tt} =  \Xi_c^{tt} & = & 32\pi \eta \frac{a^2  - b^2 }{(2a^2 - 3b^2 )S + 2a},
-+
+\end{eqnarray}
-+
+and rotational,
-+
+\begin{eqnarray}
-+
+ \Xi_a^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^2  - b^2 )b^2}{2a - b^2 S}, \\
-+
+ \Xi_b^{rr} = \Xi_c^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^4  - b^4)}{(2a^2  - b^2 )S - 2a}
-+
+\end{eqnarray}
-+
+resistance tensors are diagonal $3 \times 3$ matrices. For both
-+
+spherical and ellipsoidal particles, the translation-rotation and
-+
+rotation-translation coupling tensors are zero.
-+
-+
+\subsubsection{\label{ldintroSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}}
-+
+Other than the fluctuation dissipation formulae given by
-+
+Peters,\cite{Peters:1999qy,Peters:1999uq,Peters:2000fk} there are no
-+
+analytic solutions for the friction tensor for rigid molecules of
-+
+arbitrary shape. The ellipsoid of revolution and general triaxial
-+
+ellipsoid models have been widely used to approximate the hydrodynamic
-+
+properties of rigid bodies. However, the mapping from all possible
-+
+ellipsoidal spaces ($r$-space) to all possible combinations of
-+
+rotational diffusion coefficients ($D$-space) is not
-+
+unique.\cite{Wegener1979} Additionally, because there is intrinsic
-+
+coupling between translational and rotational motion of {\it skew}
-+
+rigid bodies, general ellipsoids are not always suitable for modeling
-+
+rigid molecules.  A number of studies have been devoted to determining
-+
+the friction tensor for irregular shapes using methods in which the
-+
+molecule of interest is modeled with a combination of
-+
+spheres\cite{Carrasco1999} and the hydrodynamic properties of the
-+
+molecule are then calculated using a set of two-point interaction
-+
+tensors.  We have found the {\it bead} and {\it rough shell} models of
-+
+Carrasco and Garc\'{i}a de la Torre to be the most useful of these
-+
+methods,\cite{Carrasco1999} and we review the basic outline of the
-+
+rough shell approach here.  A more thorough explanation can be found
-+
+in Ref. \citen{Carrasco1999}.
-+
-+
+Consider a rigid assembly of $N$ small beads moving through a
-+
+continuous medium.  Due to hydrodynamic interactions between the
-+
+beads, the net velocity of the $i^\mathrm{th}$ bead relative to the
-+
+medium, ${\bf v}'_i$, is different than its unperturbed velocity ${\bf
-+
+v}_i$,
-+
+\begin{equation}
-+
+{\bf v}'_i  = {\bf v}_i  - \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }
-+
+\end{equation}
-+
+where ${\bf F}_j$ is the frictional force on the medium due to bead $j$, and
-+
+${\bf T}_{ij}$ is the hydrodynamic interaction tensor between the two beads.
-+
+The frictional force felt by the $i^\mathrm{th}$ bead is proportional to
-+
+its net velocity
-+
+\begin{equation}
-+
+{\bf F}_i  = \xi_i {\bf v}_i  - \xi_i \sum\limits_{j \ne i} {{\bf T}_{ij} {\bf F}_j }.
-+
+\label{ldintroEquation:tensorExpression}
-+
+\end{equation}
-+
+Eq. (\ref{ldintroEquation: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,
-+
+\begin{equation}
-+
+{\bf T}_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left( {{\bf I} +
-+
+\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right).
-+
+\label{ldintroEquation:oseenTensor}
-+
+\end{equation}
-+
+Here ${\bf R}_{ij}$ is the distance vector between beads $i$ and
-+
+$j$. A second order expression for beads of different hydrodynamic
-+
+radii was introduced by Rotne and Prager,\cite{Rotne1969} and improved
-+
+by Garc\'{i}a de la Torre and Bloomfield,\cite{Torre1977}
-+
+\begin{equation}
-+
+{\bf T}_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {{\bf I} +
-+
+\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right) + \frac{{\rho
-+
+_i^2  + \rho_j^2 }}{{R_{ij}^2 }}\left( {\frac{{\bf I}}{3} -
-+
+\frac{{{\bf R}_{ij} {\bf R}_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
-+
+\label{ldintroEquation:RPTensorNonOverlapped}
-+
+\end{equation}
-+
+ Both the Oseen-Burgers tensor and
-+
+ Eq.~\ref{ldintroEquation:RPTensorNonOverlapped} have an assumption
-+
+ that the beads do not overlap ($R_{ij} \ge \rho_i + \rho_j$).
-+
-+
+To calculate the resistance tensor for a body represented as the union
-+
+of many non-overlapping beads, we first pick an arbitrary origin $O$
-+
+and then construct a $3N \times 3N$ supermatrix consisting of $N
-+
+\times N$ ${\bf B}_{ij}$ blocks
-+
+\begin{equation}
-+
+{\bf B} = \left( \begin{array}{*{20}c}
-+
+ {\bf B}_{11} &  \ldots  & {\bf B}_{1N}   \\
-+
+\vdots  &  \ddots  &  \vdots   \\
-+
+ {\bf B}_{N1} &  \cdots  & {\bf B}_{NN}
-+
+\end{array} \right)
-+
+\end{equation}
-+
+${\bf B}_{ij}$ is a version of the hydrodynamic tensor which includes the
-+
+self-contributions for spheres,
-+
+\begin{equation}
-+
+{\bf B}_{ij}  = \delta _{ij} \frac{{\bf I}}{{6\pi \eta R_{ij}}} + (1 - \delta_{ij}
-+
+){\bf T}_{ij}
-+
+\end{equation}
-+
+where $\delta_{ij}$ is the Kronecker delta function. Inverting the
-+
+${\bf B}$ matrix, we obtain
-+
+\begin{equation}
-+
+{\bf C} = {\bf B}^{ - 1}  = \left(\begin{array}{*{20}c}
-+
+  {\bf C}_{11} &  \ldots  & {\bf C}_{1N} \\
-+
+ \vdots  &  \ddots  &  \vdots   \\
-+
+  {\bf C}_{N1} &  \cdots  & {\bf C}_{NN}
-+
+\end{array} \right),
-+
+\end{equation}
-+
+which can be partitioned into $N \times N$ blocks labeled ${\bf C}_{ij}$.
-+
+(Each of the ${\bf C}_{ij}$ blocks is a $3 \times 3$ matrix.)  Using the
-+
+skew matrix,
-+
+\begin{equation}
-+
+{\bf U}_i  = \left(\begin{array}{*{20}c}
-+
+& -z_i & y_i \\
-+
+z_i &  0   & - x_i \\
-+
+-y_i & x_i & 0
-+
+\end{array}\right)
-+
+\label{ldeq:skewMatrix}
-+
+\end{equation}
-+
+where $x_i$, $y_i$, $z_i$ are the components of the vector joining
-+
+bead $i$ and origin $O$, the elements of the resistance tensor (at the
-+
+arbitrary origin $O$) can be written as
-+
+\begin{eqnarray}
-+
+\label{ldintroEquation:ResistanceTensorArbitraryOrigin}
-+
+ \Xi^{tt}  & = & \sum\limits_i {\sum\limits_j {{\bf C}_{ij} } } \notag , \\
-+
+ \Xi^{tr}  = \Xi _{}^{rt}  & = & \sum\limits_i {\sum\limits_j {{\bf U}_i {\bf C}_{ij} } } , \\
-+
+ \Xi^{rr}  & = &  -\sum\limits_i \sum\limits_j {\bf U}_i {\bf C}_{ij} {\bf U}_j + 6 \eta V {\bf I}. \notag
-+
+\end{eqnarray}
-+
+The final term in the expression for $\Xi^{rr}$ is a correction that
-+
+accounts for errors in the rotational motion of the bead models. The
-+
+additive correction uses the solvent viscosity ($\eta$) as well as the
-+
+total volume of the beads that contribute to the hydrodynamic model,
-+
+\begin{equation}
-+
+V = \frac{4 \pi}{3} \sum_{i=1}^{N} \rho_i^3,
-+
+\end{equation}
-+
+where $\rho_i$ is the radius of bead $i$.  This correction term was
-+
+rigorously tested and compared with the analytical results for
-+
+two-sphere and ellipsoidal systems by Garc\'{i}a de la Torre and
-+
+Rodes.\cite{Torre:1983lr}
-+
-+
+In general, resistance tensors depend on the origin at which they were
-+
+computed.  However, the proper location for applying the friction
-+
+force is the center of resistance, the special point at which the
-+
+trace of rotational resistance tensor, $\Xi^{rr}$ reaches a minimum
-+
+value.  Mathematically, the center of resistance can also be defined
-+
+as the unique point for a rigid body at which the translation-rotation
-+
+coupling tensors are symmetric,
-+
+\begin{equation}
-+
+\Xi^{tr}  = \left(\Xi^{tr}\right)^T
-+
+\label{ldintroEquation:definitionCR}
-+
+\end{equation}
-+
+From Eq. \ref{ldintroEquation: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{ldintroEquation: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{ldeq:skewMatrix}))
-+
+for the vector between the origin $O$ and the point $P$. Using
-+
+Eqs.~\ref{ldintroEquation:definitionCR}~and~\ref{ldintroEquation:resistanceTensorTransformation},
-+
+one can locate the position of center of resistance,
-+
+\begin{eqnarray*}
-+
+\left(\begin{array}{l}
-+
+x_{OR} \\
-+
+y_{OR} \\
-+
+z_{OR} \\
-+
+\end{array}\right) & = &
-+
+\left(\begin{array}{*{20}c}
-+
+(\Xi_O^{rr})_{yy} + (\Xi_O^{rr})_{zz} & -(\Xi_O^{rr})_{xy} & -(\Xi_O^{rr})_{xz} \\
-+
+-(\Xi_O^{rr})_{xy} & (\Xi_O^{rr})_{zz} + (\Xi_O^{rr})_{xx} & -(\Xi_O^{rr})_{yz} \\
-+
+-(\Xi_O^{rr})_{xz} & -(\Xi_O^{rr})_{yz} & (\Xi_O^{rr})_{xx} + (\Xi_O^{rr})_{yy} \\
-+
+\end{array}\right)^{-1} \\
-+
+ & & \left(\begin{array}{l}
-+
+(\Xi_O^{tr})_{yz} - (\Xi_O^{tr})_{zy} \\
-+
+(\Xi_O^{tr})_{zx} - (\Xi_O^{tr})_{xz} \\
-+
+(\Xi_O^{tr})_{xy} - (\Xi_O^{tr})_{yx} \\
-+
+\end{array}\right) \\
-+
+\end{eqnarray*}
-+
+where $x_{OR}$, $y_{OR}$, $z_{OR}$ are the components of the vector
-+
+joining center of resistance $R$ and origin $O$.
-+
-+
+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
-+
+.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{ldintroEquation: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{ldLDRB}}
-+
-+
+Consider the Langevin equations of motion in generalized coordinates
-+
+\begin{equation}
-+
+{\bf M} \dot{{\bf V}}(t) = {\bf F}_{s}(t) +
-+
+{\bf F}_{f}(t)  + {\bf F}_{r}(t)
-+
+\label{ldLDGeneralizedForm}
-+
+\end{equation}
-+
+where ${\bf M}$ is a $6 \times 6$ diagonal mass matrix (which
-+
+includes the mass of the rigid body as well as the moments of inertia
-+
+in the body-fixed frame) and ${\bf V}$ is a generalized velocity,
-+
+${\bf V} =
-+
+\left\{{\bf v},{\bf \omega}\right\}$. The right side of
-+
+Eq.~\ref{ldLDGeneralizedForm} consists of three generalized forces: a
-+
+system force (${\bf F}_{s}$), a frictional or dissipative force (${\bf
-+
+F}_{f}$) and a stochastic force (${\bf F}_{r}$). While the evolution
-+
+of the system in Newtonian mechanics is typically done in the lab
-+
+frame, it is convenient to handle the dynamics of rigid bodies in
-+
+body-fixed frames. Thus the friction and random forces on each
-+
+substructure are calculated in a body-fixed frame and may converted
-+
+back to the lab frame using that substructure's rotation matrix (${\bf
-+
+Q}$):
-+
+\begin{equation}
-+
+{\bf F}_{f,r} =
-+
+\left( \begin{array}{c}
-+
+{\bf f}_{f,r} \\
-+
+{\bf \tau}_{f,r}
-+
+\end{array} \right)
-+
+=
-+
+\left( \begin{array}{c}
-+
+{\bf Q}^{T} {\bf f}^{~b}_{f,r} \\
-+
+{\bf Q}^{T} {\bf \tau}^{~b}_{f,r}
-+
+\end{array} \right)
-+
+\end{equation}
-+
+The body-fixed friction force, ${\bf F}_{f}^{~b}$, is proportional to
-+
+the (body-fixed) velocity at the center of resistance
-+
+${\bf v}_{R}^{~b}$ and the angular velocity ${\bf \omega}$
-+
+\begin{equation}
-+
+{\bf F}_{f}^{~b}(t) = \left( \begin{array}{l}
-+
+ {\bf f}_{f}^{~b}(t) \\
-+
+ {\bf \tau}_{f}^{~b}(t) \\
-+
+ \end{array} \right) =  - \left( \begin{array}{*{20}c}
-+
+   \Xi_{R}^{tt} & \Xi_{R}^{rt} \\
-+
+   \Xi_{R}^{tr} & \Xi_{R}^{rr}    \\
-+
+\end{array} \right)\left( \begin{array}{l}
-+
+ {\bf v}_{R}^{~b}(t) \\
-+
+ {\bf \omega}(t) \\
-+
+ \end{array} \right),
-+
+\end{equation}
-+
+while the random force, ${\bf F}_{r}$, is a Gaussian stochastic
-+
+variable with zero mean and variance,
-+
+\begin{equation}
-+
+\left\langle {{\bf F}_{r}(t) ({\bf F}_{r}(t'))^T } \right\rangle  =
-+
+\left\langle {{\bf F}_{r}^{~b} (t) ({\bf F}_{r}^{~b} (t'))^T } \right\rangle  =
-+
+k_B T \Xi_R \delta(t - t'). \label{ldeq:randomForce}
-+
+\end{equation}
-+
+$\Xi_R$ is the $6\times6$ resistance tensor at the center of
-+
+resistance.  Once this tensor is known for a given rigid body (as
-+
+described in the previous section) obtaining a stochastic vector that
-+
+has the properties in Eq. (\ref{ldeq: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{ldeq:Cholesky}
-+
+\end{equation}
-+
+where ${\bf S}$ is a lower triangular matrix.\cite{Schlick2002} A
-+
+vector with the statistics required for the random force can then be
-+
+obtained by multiplying ${\bf S}$ onto a random 6-vector ${\bf Z}$ which
-+
+has elements chosen from a Gaussian distribution, such that:
-+
+\begin{equation}
-+
+\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot
-+
+{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij},
-+
+\end{equation}
-+
+where $\delta t$ is the timestep in use during the simulation. The
-+
+random force, ${\bf F}_{r}^{~b} = {\bf S} {\bf Z}$, can be shown to have the
-+
+correct properties required by Eq. (\ref{ldeq:randomForce}).
-+
-+
+The equation of motion for the translational velocity of the center of
-+
+mass (${\bf v}$) can be written as
-+
+\begin{equation}
-+
+m \dot{{\bf v}} (t) =  {\bf f}_{s}(t) + {\bf f}_{f}(t) +
-+
+{\bf f}_{r}(t)
-+
+\end{equation}
-+
+Since the frictional and random forces are applied at the center of
-+
+resistance, which generally does not coincide with the center of mass,
-+
+extra torques are exerted at the center of mass. Thus, the net
-+
+body-fixed torque at the center of mass, $\tau^{~b}(t)$,
-+
+is given by
-+
+\begin{equation}
-+
+\tau^{~b} \leftarrow \tau_{s}^{~b} + \tau_{f}^{~b} + \tau_{r}^{~b} + {\bf r}_{MR} \times \left( {\bf f}_{f}^{~b} + {\bf f}_{r}^{~b} \right)
-+
+\end{equation}
-+
+where ${\bf r}_{MR}$ is the vector from the center of mass to the center of
-+
+resistance. Instead of integrating the angular velocity in lab-fixed
-+
+frame, we consider the equation of motion for the angular momentum
-+
+(${\bf j}$) in the body-fixed frame
-+
+\begin{equation}
-+
+\frac{\partial}{\partial t}{\bf j}(t) = \tau^{~b}(t)
-+
+\end{equation}
-+
+Embedding the friction and random forces into the the total force and
-+
+torque, one can integrate the Langevin equations of motion for a rigid
-+
+body of arbitrary shape in a velocity-Verlet style 2-part algorithm,
-+
+where $h = \delta t$:
-+
-+
+{\tt move A:}
-+
+\begin{align*}
-+
+{\bf v}\left(t + h / 2\right)  &\leftarrow  {\bf v}(t)
-+
+    + \frac{h}{2} \left( {\bf f}(t) / m \right), \\
-+
+%
-+
+{\bf r}(t + h) &\leftarrow {\bf r}(t)
-+
+    + h  {\bf v}\left(t + h / 2 \right), \\
-+
+%
-+
+{\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
-+
+    + \frac{h}{2} {\bf \tau}^{~b}(t), \\
-+
+%
-+
+{\bf Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
-+
+    (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
-+
+\end{align*}
-+
+In this context, $\overleftrightarrow{\mathsf{I}}$ is the diagonal
-+
+moment of inertia tensor, and the $\mathrm{rotate}$ function is the
-+
+reversible product of the three body-fixed rotations,
-+
+\begin{equation}
-+
+\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
-+
+\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y
-+
+/ 2) \cdot \mathsf{G}_x(a_x /2),
-+
+\end{equation}
-+
+where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
-+
+rotates both the rotation matrix ($\mathbf{Q}$) and the body-fixed
-+
+angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed
-+
+axis $\alpha$,
-+
+\begin{equation}
-+
+\mathsf{G}_\alpha( \theta ) = \left\{
-+
+\begin{array}{lcl}
-+
+\mathbf{Q}(t) & \leftarrow & \mathbf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
-+
+{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
-+
+j}(0).
-+
+\end{array}
-+
+\right.
-+
+\end{equation}
-+
+$\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis
-+
+rotation matrix.  For example, in the small-angle limit, the
-+
+rotation matrix around the body-fixed x-axis can be approximated as
-+
+\begin{equation}
-+
+\mathsf{R}_x(\theta) \approx \left(
-+
+\begin{array}{ccc}
-+
+& 0 & 0 \\
-+
+& \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
-+
+\theta^2 / 4} \\
-+
+& \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 +
-+
+\theta^2 / 4}
-+
+\end{array}
-+
+\right).
-+
+\end{equation}
-+
+All other rotations follow in a straightforward manner. After the
-+
+first part of the propagation, the forces and body-fixed torques are
-+
+calculated at the new positions and orientations.  The system forces
-+
+and torques are derivatives of the total potential energy function
-+
+($U$) with respect to the rigid body positions (${\bf r}$) and the
-+
+columns of the transposed rotation matrix ${\bf Q}^T = \left({\bf
-+
+u}_x, {\bf u}_y, {\bf u}_z \right)$:
-+
-+
+{\tt Forces:}
-+
+\begin{align*}
-+
+{\bf f}_{s}(t + h) & \leftarrow
-+
+    - \left(\frac{\partial U}{\partial {\bf r}}\right)_{{\bf r}(t + h)} \\
-+
+%
-+
+{\bf \tau}_{s}(t + h) &\leftarrow {\bf u}(t + h)
-+
+    \times \frac{\partial U}{\partial {\bf u}} \\
-+
+%
-+
+{\bf v}^{b}_{R}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \left({\bf v}(t+h) + {\bf \omega}(t+h) \times {\bf r}_{MR} \right) \\
-+
+%
-+
+{\bf f}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tt} \cdot
-+
+{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rt} \cdot {\bf \omega}(t+h) \\
-+
+%
-+
+{\bf \tau}_{R,f}^{b}(t+h) & \leftarrow - \Xi_{R}^{tr} \cdot
-+
+{\bf v}^{b}_{R}(t+h) - \Xi_{R}^{rr} \cdot {\bf \omega}(t+h) \\
-+
+%
-+
+Z & \leftarrow  {\tt GaussianNormal}(2 k_B T / h, 6) \\
-+
+{\bf F}_{R,r}^{b}(t+h) & \leftarrow {\bf S} \cdot Z \\
-+
+%
-+
+{\bf f}(t+h) & \leftarrow {\bf f}_{s}(t+h) + \mathbf{Q}^{T}(t+h)
-+
+\cdot \left({\bf f}_{R,f}^{~b} + {\bf f}_{R,r}^{~b} \right) \\
-+
+%
-+
+\tau(t+h) & \leftarrow \tau_{s}(t+h) + \mathbf{Q}^{T}(t+h) \cdot \left(\tau_{R,f}^{~b} + \tau_{R,r}^{~b} \right) + {\bf r}_{MR} \times \left({\bf f}_{f}(t+h) + {\bf f}_{r}(t+h) \right) \\
-+
+\tau^{~b}(t+h) & \leftarrow \mathbf{Q}(t+h) \cdot \tau(t+h) \\
-+
+\end{align*}
-+
+Frictional (and random) forces and torques must be computed at the
-+
+center of resistance, so there are additional steps required to find
-+
+the body-fixed velocity (${\bf v}_{R}^{~b}$) at this location.  Mapping
-+
+the frictional and random forces at the center of resistance back to
-+
+the center of mass also introduces an additional term in the torque
-+
+one obtains at the center of mass.
-+
-+
+Once the forces and torques have been obtained at the new time step,
-+
+the velocities can be advanced to the same time value.
-+
-+
+{\tt move B:}
-+
+\begin{align*}
-+
+{\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2
-+
+\right)
-+
+    + \frac{h}{2} \left( {\bf f}(t + h) / m \right), \\
-+
+%
-+
+{\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t + h / 2
-+
+\right)
-+
+    + \frac{h}{2} {\bf \tau}^{~b}(t + h) .
-+
+\end{align*}
-+
-+
+\section{Validating the Method\label{ldsec: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{ldtab:parameters},
-+
+and a sketch of the arrangement of these primary particles into the
-+
+model rigid bodies is shown in figure \ref{ldfig:models}. In table
-+
+\ref{ldtab:parameters}, $d$ and $l$ are the physical dimensions of
-+
+ellipsoidal (Gay-Berne) particles.  For spherical particles, the value
-+
+of the Lennard-Jones $\sigma$ parameter is the particle diameter
-+
+($d$).  Gay-Berne ellipsoids have an energy scaling parameter,
-+
+$\epsilon^s$, which describes the well depth for two identical
-+
+ellipsoids in a {\it side-by-side} configuration.  Additionally, a
-+
+well depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$,
-+
+describes the ratio between the well depths in the {\it end-to-end}
-+
+and side-by-side configurations.  For spheres, $\epsilon^r \equiv 1$.
-+
+Moments of inertia are also required to describe the motion of primary
-+
+particles with orientational degrees of freedom.
-+
-+
+\begin{table*}
-+
+\begin{minipage}{\linewidth}
-+
+\begin{center}
-+
+\caption{PARAMETERS FOR THE PRIMARY PARTICLES IN USE BY THE RIGID BODY
-+
+MODELS IN FIGURE \ref{ldfig:models}}
-+
+\begin{tabular}{lrcccccccc}
-+
+\hline
-+
+ & & & & & & \multicolumn{3}c{$\overleftrightarrow{\mathsf I}$ (amu \AA$^2$)} \\
-+
+ & $d$ (\AA) & $l$ (\AA) & $\epsilon^s$ ($\frac{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: & & & & & & & & \\
-+
+\quad {\it 2 spheres} & 6.5 & $= d$ & 0.8 & 1   & 190 & 802.75 & 802.75 & 802.75 \\
-+
+Banana: & & & & & & & & \\
-+
+\quad {\it 3 ellipsoids} & 4.2 & 11.2  & 0.8 & 0.2 & 240 & 10000 & 10000 & 0 \\
-+
+Lipid: & & & & & & & & \\
-+
+\quad {\it head} & 6.5 & $= d$ & 0.185 & 1 & 196 & & & \\
-+
+\quad {\it 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{ldtab:parameters}
-+
+\end{center}
-+
+\end{minipage}
-+
+\end{table*}
-+
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=3in]{./figures/ldSketch}
-+
+\caption[Sketch of the model systems]{A sketch of the model systems
-+
+used in evaluating the behavior of the rigid body Langevin
-+
+integrator.} \label{ldfig:models}
-+
+\end{figure}
-+
-+
+\subsection{Simulation Methodology}
-+
+We performed reference microcanonical simulations with explicit
-+
+solvents for each of the different model system.  In each case there
-+
+was one solute model and 1929 solvent molecules present in the
-+
+simulation box.  All simulations were equilibrated for 5 ns using a
-+
+constant-pressure and temperature integrator with target values of 300
-+
+K for the temperature and 1 atm for pressure.  Following this stage,
-+
+further equilibration (5 ns) and sampling (10 ns) was done in a
-+
+microcanonical ensemble.  Since the model bodies are typically quite
-+
+massive, we were able to use a time step of 25 fs.
-+
-+
+The model systems studied used both Lennard-Jones spheres as well as
-+
+uniaxial Gay-Berne ellipoids. The Gay-Berne potential is given by
-+
+equation~\ref{mdeq:gb}.  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$).
-+
+A switching function (Eq.~\ref{mdeq:dipoleSwitching}) was applied to
-+
+all potentials to smoothly turn off the interactions between a range
-+
+of $22$ and $25$ \AA.
-+
-+
+To measure shear viscosities from our microcanonical simulations, we
-+
+used the Einstein form of the pressure correlation function,\cite{hess:209}
-+
+\begin{equation}
-+
+\eta = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left(
-+
+\int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \right\rangle_{t_0}.
-+
+\label{ldeq:shear}
-+
+\end{equation}
-+
+which converges much more rapidly in molecular dynamics simulations
-+
+than the traditional Green-Kubo formula.
-+
-+
+The Langevin dynamics for the different model systems were performed
-+
+at the same temperature as the average temperature of the
-+
+microcanonical simulations and with a solvent viscosity taken from
-+
+Eq. (\ref{ldeq: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 (Eq.~\ref{mdeq:msdisplacement}) 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{ldtab: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}
-+
+/ \tau_1 & = & 6 D_r + 2 \Delta \\
-+
+/ \tau_2 & = & 6 D_r - 2 \Delta \\
-+
+/ \tau_3 & = & 3 (D_r + D_1)  \\
-+
+/ \tau_4 & = & 3 (D_r + D_2) \\
-+
+/ \tau_5 & = & 3 (D_r + D_3)
-+
+\end{eqnarray}
-+
+where
-+
+\begin{equation}
-+
+D_r = \frac{1}{3} \left(D_1 + D_2 + D_3 \right)
-+
+\end{equation}
-+
+and
-+
+\begin{equation}
-+
+\Delta = \left( (D_1 - D_2)^2 + (D_3 - D_1 )(D_3 - D_2)\right)^{1/2}
-+
+\end{equation}
-+
+Each of these characteristic times can be used to predict the decay of
-+
+part of the rotational correlation function when $\ell = 2$,
-+
+\begin{equation}
-+
+C_2(t)  =  \frac{a^2}{N^2} e^{-t/\tau_1} + \frac{b^2}{N^2} e^{-t/\tau_2}.
-+
+\end{equation}
-+
+This is the same as the $F^2_{0,0}(t)$ correlation function that
-+
+appears in Ref. \citen{Berne90}.  The amplitudes of the two decay
-+
+terms are expressed in terms of three dimensionless functions of the
-+
+eigenvalues: $a = \sqrt{3} (D_1 - D_2)$, $b = (2D_3 - D_1 - D_2 +
-+
+\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}
-+
+/ \tau_0 = \frac{1}{5} \sum_{i=1}^5 \tau_{i}^{-1},
-+
+\end{equation}
-+
+although for the cases reported here, this averaging is not necessary
-+
+and only one of the five relaxation times is relevant.
-+
-+
+To test the Langevin integrator's behavior for rotational relaxation,
-+
+we have compared the analytical orientational relaxation times (if
-+
+they are known) with the general result from the diffusion tensor and
-+
+with the results from both the explicitly solvated molecular dynamics
-+
+and Langevin simulations.  Relaxation times from simulations (both
-+
+microcanonical and Langevin), were computed using Legendre polynomial
-+
+correlation functions for a unit vector (${\bf u}$) fixed along one or
-+
+more of the body-fixed axes of the model.
-+
+\begin{equation}
-+
+C_{\ell}(t)  =  \left\langle P_{\ell}\left({\bf u}_{i}(t) \cdot {\bf
-+
+u}_{i}(0) \right) \right\rangle
-+
+\end{equation}
-+
+For simulations in the high-friction limit, orientational correlation
-+
+times can then be obtained from exponential fits of this function, or by
-+
+integrating,
-+
+\begin{equation}
-+
+\tau = \ell (\ell + 1) \int_0^{\infty} C_{\ell}(t) dt.
-+
+\end{equation}
-+
+In lower-friction solvents, the Legendre correlation functions often
-+
+exhibit non-ex\-po\-nen\-tial decay, and may not be characterized by a
-+
+single decay constant.
-+
-+
+In table \ref{ldtab: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$ =
-+
+.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{ldeq:StokesTranslation}) and
-+
+(\ref{ldeq: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{ldDperrin}
-+
+\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{ldThetaPerrin}
-+
+\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{ldGPerrin}
-+
+\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{ldfig: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{ldintroEquation: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]{./figures/ldRoughShell}
-+
+\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{ldfig: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{ldfig:explanation}.
-+
-+
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=4in]{./figures/ldExplanation}
-+
+\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{ldfig: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{ldeq:shear})).  To calculate the hydrodynamic properties of
-+
+the banana rigid body model, we created a rough shell (see
-+
+Fig.~\ref{ldfig: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{ldintroEquation: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{ldeq:shear})).  To calculate the hydrodynamic properties
-+
+of the lipid rigid body model, we created a rough shell (see
-+
+Fig.~\ref{ldfig:roughShell}), in which the lipid is represented as a
-+
+``shell'' made of 3550 identical beads (0.25 \AA\ in diameter)
-+
+distributed on the surface.  Applying the procedure described by
-+
+Eq. (\ref{ldintroEquation: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]{./figures/ldGraph}
-+
+\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{ldfig:results}
-+
+\end{figure}
-+
-+
+\begin{table*}
-+
+\begin{center}
-+
+\caption{TRANSLATIONAL DIFFUSION CONSTANTS (D) FOR THE MODEL SYSTEMS
-+
+CALCULATED USING MICROCANONICAL SIM\-U\-LA\-TIONS (WITH EXPLICIT
-+
+SOLVENT), THEORETICAL PREDICTIONS, AND LANGEVIN SIMULATIONS (WITH
-+
+IMPLICIT SOLVENT)}
-+
+\begin{tabular}{lccccccc}
-+
+\hline
-+
+ & \multicolumn{2}c{microcanonical} & & \multicolumn{3}c{Theoretical} & Langevin \\
-+
+\cline{2-3} \cline{5-7}
-+
+model & $\eta$ (cP)  & D & & Analytical & method & Hydro &  \\
-+
+\hline
-+
+sphere    & 0.279  & 3.06 & & 2.42 & exact       & 2.42 & 2.33 \\
-+
+ellipsoid & 0.255  & 2.44 & & 2.34 & exact       & 2.34 & 2.37 \\
-+
+          & 0.255  & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\
-+
+dumbbell  & 0.308  & 2.06 & & 1.64 & bead model  & 1.65 & 1.62 \\
-+
+          & 0.308  & 2.06 & & 1.64 & rough shell & 1.59 & 1.62 \\
-+
+banana    & 0.298  & 1.53 & &      & rough shell & 1.56 & 1.55 \\
-+
+lipid     & 0.349  & 1.41 & &      & rough shell & 1.33 & 1.32 \\
-+
+\hline
-+
+\end{tabular}
-+
+\begin{minipage}{\linewidth}
-+
+%\centering
-+
+\vspace{2mm}
-+
+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).
-+
+\label{ldtab:translation}
-+
+\end{minipage}
-+
+\end{center}
-+
+\end{table*}
-+
-+
+\begin{table*}
-+
+\begin{center}
-+
+\caption{ORIENTATIONAL RELAXATION TIMES ($\tau$) FOR THE MODEL SYSTEMS USING
-+
+MICROCANONICAL SIMULATION (WITH EXPLICIT SOLVENT), THEORETICAL
-+
+PREDICTIONS, AND LANGEVIN SIMULATIONS (WITH IMPLICIT SOLVENT)}
-+
+\begin{tabular}{lccccccc}
-+
+\hline
-+
+ & \multicolumn{2}c{microcanonical} & & \multicolumn{3}c{Theoretical} & Langevin \\
-+
+\cline{2-3} \cline{5-7}
-+
+model & $\eta$ (cP) & $\tau$ & & Analytical & method & Hydro &  \\
-+
+\hline
-+
+sphere    & 0.279  &      & & 9.69 & exact       & 9.69 & 9.64 \\
-+
+ellipsoid & 0.255  & 46.7 & & 22.0 & exact       & 22.0 & 22.2 \\
-+
+          & 0.255  & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\
-+
+dumbbell  & 0.308  & 14.1 & &      & bead model  & 50.0 & 50.1 \\
-+
+          & 0.308  & 14.1 & &      & rough shell & 41.5 & 41.3 \\
-+
+banana    & 0.298  & 63.8 & &      & rough shell & 70.9 & 70.9 \\
-+
+lipid     & 0.349  & 78.0 & &      & rough shell & 76.9 & 77.9 \\
-+
+\hline
-+
+\end{tabular}
-+
+\begin{minipage}{\linewidth}
-+
+%\centering
-+
+\vspace{2mm}
-+
+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.
-+
+\label{ldtab:rotation}
-+
+\end{minipage}
-+
+\end{center}
-+
+\end{table*}
-+
-+
+\section{Application: A rigid-body lipid bilayer}
-+
-+
+To test the accuracy and efficiency of the new integrator, we applied
-+
+it to study the formation of corrugated structures emerging from
-+
+simulations of the coarse grained lipid molecular models presented
-+
+above.  The initial configuration is taken from earlier molecular
-+
+dynamics studies on lipid bilayers which had used spherical
-+
+(Lennard-Jones) solvent particles and moderate (480 solvated lipid
-+
+molecules) system sizes.\cite{SunX._jp0762020} the solvent molecules
-+
+were excluded from the system and the box was replicated three times
-+
+in the x- and y- axes (giving a single simulation cell containing 4320
-+
+lipids).  The experimental value for the viscosity of water at 20C
-+
+($\eta = 1.00$ cp) was used with the Langevin integrator to mimic the
-+
+hydrodynamic effects of the solvent.  The absence of explicit solvent
-+
+molecules and the stability of the integrator allowed us to take
-+
+timesteps of 50 fs. A simulation run time of 30 ns was sampled to
-+
+calculate structural properties.  Fig. \ref{ldfig:bilayer} shows the
-+
+configuration of the system after 30 ns.  Structural properties of the
-+
+bilayer (e.g. the head and body $P_2$ order parameters) are nearly
-+
+identical to those obtained via solvated molecular dynamics. The
-+
+ripple structure remained stable during the entire trajectory.
-+
+Compared with using explicit bead-model solvent molecules, the 30 ns
-+
+trajectory for 4320 lipids with the Langevin integrator is now {\it
-+
+faster} on the same hardware than the same length trajectory was for
-+
+the 480-lipid system previously studied.
-+
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=\linewidth]{./figures/ldBilayer}
-+
+\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{ldfig:bilayer}
-+
+\end{figure}
-+
-+
+\section{Conclusions}
-+
-+
+We have presented a new algorithm for carrying out Langevin dynamics
-+
+simulations on complex rigid bodies by incorporating the hydrodynamic
-+
+resistance tensors for arbitrary shapes into a stable and efficient
-+
+integration scheme.  The integrator gives quantitative agreement with
-+
+both analytic and approximate hydrodynamic theories, and works
-+
+reasonably well at reproducing the solute dynamical properties
-+
+(diffusion constants, and orientational relaxation times) from
-+
+explicitly-solvated simulations.  For the cases where there are
-+
+discrepancies between our Langevin integrator and the explicit solvent
-+
+simulations, two features of molecular simulations help explain the
-+
+differences.
-+
-+
+First, the use of ``stick'' boundary conditions for molecular-sized
-+
+solutes in a sea of similarly-sized solvent particles may be
-+
+problematic.  We are certainly not the first group to notice this
-+
+difference between hydrodynamic theories and explicitly-solvated
-+
+molecular
-+
+simulations.\cite{Schmidt:2004fj,Schmidt:2003kx,Ravichandran:1999fk,TANG:1993lr}
-+
+The problem becomes particularly noticable in both the translational
-+
+diffusion of the spherical particles and the rotational diffusion of
-+
+the ellipsoids.  In both of these cases it is clear that the
-+
+approximations that go into hydrodynamics are the source of the error,
-+
+and not the integrator itself.
-+
-+
+Second, in the case of structures which have substantial surface area
-+
+that is inaccessible to solvent particles, the hydrodynamic theories
-+
+(and the Langevin integrator) may overestimate the effects of solvent
-+
+friction because they overestimate the exposed surface area of the
-+
+rigid body.  This is particularly noticable in the rotational
-+
+diffusion of the dumbbell model.  We believe that given a solvent of
-+
+known radius, it may be possible to modify the rough shell approach to
-+
+place beads on solvent-accessible surface, instead of on the geometric
-+
+surface defined by the van der Waals radii of the components of the
-+
+rigid body.  Further work to confirm the behavior of this new
-+
+approximation is ongoing.

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing trunk/xDissertation/ld.tex (file contents): Revision 3336 by xsun, Wed Jan 30 16:01:02 2008 UTC vs. Revision 3383 by xsun, Wed Apr 16 21:56:34 2008 UTC

Diff Legend

Comparing trunk/xDissertation/ld.tex (file contents):
Revision 3336 by xsun, Wed Jan 30 16:01:02 2008 UTC vs.
Revision 3383 by xsun, Wed Apr 16 21:56:34 2008 UTC