--- 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}