--- trunk/langevin/langevin.tex 2008/01/11 22:08:57 3310 +++ trunk/langevin/langevin.tex 2008/02/29 22:02:46 3352 @@ -1,516 +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 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{Xiuquan Sun, Teng Lin 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. +\subsection{Rigid Body Dynamics} +Rigid bodies are frequently involved in the modeling of large +collections of particles that move as a single unit. In molecular +simulations, rigid bodies have been used to simplify protein-protein +docking,\cite{Gray2003} and lipid bilayer +simulations.\cite{SunX._jp0762020} Many of the water models in common +use are also rigid-body +models,\cite{Jorgensen83,Berendsen81,Berendsen87} although they are +typically evolved in molecular dynamics simulations using constraints +rather than rigid body equations of motion. + +Euler angles are a natural choice to describe the rotational degrees +of freedom. However, due to $\frac{1}{\sin \theta}$ singularities, the +numerical integration of corresponding equations of these motion can +become inaccurate (and inefficient). Although the use of multiple +sets of Euler angles can overcome this problem,\cite{Barojas1973} the +computational penalty and the loss of angular momentum conservation +remain. A singularity-free representation utilizing quaternions was +developed by Evans in 1977.\cite{Evans1977} The Evans quaternion +approach uses a nonseparable Hamiltonian, and this has prevented +symplectic algorithms from being utilized until very +recently.\cite{Miller2002} + +Another approach is the application of holonomic constraints to the +atoms belonging to the rigid body. Each atom moves independently +under the normal forces deriving from potential energy and constraints +are used to guarantee rigidity. However, due to their iterative +nature, the SHAKE and RATTLE algorithms converge very slowly when the +number of constraints (and the number of particles that belong to the +rigid body) increases.\cite{Ryckaert1977,Andersen1983} + +In order to develop a stable and efficient integration scheme that +preserves most constants of the motion, 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 the -accurate estimation of friction tensor from hydrodynamics theory -into the sophisticated rigid body dynamics algorithms. +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}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, -\[ +\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 + 6 \eta V {\bf I}. \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 final term in the expression for $\Xi^{rr}$ is correction that -accounts for errors in the rotational motion of certain kinds of bead -models. The additive correction uses the solvent viscosity ($\eta$) -as well as the total volume of the beads that contribute to the -hydrodynamic model, +The final term in the expression for $\Xi^{rr}$ is a correction that +accounts for errors in the rotational motion of the bead models. The +additive correction uses the solvent viscosity ($\eta$) as well as the +total volume of the beads that contribute to the hydrodynamic model, \begin{equation} -V = \frac{4 \pi}{3} \sum_{i=1}^{N} \sigma_i^3, +V = \frac{4 \pi}{3} \sum_{i=1}^{N} \rho_i^3, \end{equation} -where $\sigma_i$ is the radius of bead $i$. This correction term was +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 Garcia de la Torre and +two-sphere and ellipsoidal systems by Garc\'{i}a de la Torre and Rodes.\cite{Torre:1983lr} - -The resistance tensor depends on the origin to which they refer. The -proper location for applying the friction force is the center of -resistance (or center of reaction), at which the trace of rotational -resistance tensor, $ \Xi ^{rr}$ reaches a minimum value. -Mathematically, the center of resistance is defined as an unique -point of the rigid body at which the translation-rotation coupling -tensors are symmetric, +In general, resistance tensors depend on the origin at which they were +computed. However, the proper location for applying the friction +force is the center of resistance, the special point at which the +trace of rotational resistance tensor, $\Xi^{rr}$ reaches a minimum +value. Mathematically, the center of resistance can also be defined +as the unique point for a rigid body at which the translation-rotation +coupling tensors are symmetric, \begin{equation} -\Xi^{tr} = \left( {\Xi^{tr} } \right)^T +\Xi^{tr} = \left(\Xi^{tr}\right)^T \label{introEquation:definitionCR} \end{equation} -From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin}, -we can easily derive that the translational resistance tensor is -origin independent, while the rotational resistance tensor and +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 ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can -obtain the resistance tensor at $P$ by -\begin{equation} -\begin{array}{l} - \Xi _P^{tt} = \Xi _O^{tt} \\ - \Xi _P^{tr} = \Xi _P^{rt} = \Xi _O^{tr} - U_{OP} \Xi _O^{tt} \\ - \Xi _P^{rr} = \Xi _O^{rr} - U_{OP} \Xi _O^{tt} U_{OP} + \Xi _O^{tr} U_{OP} - U_{OP} \Xi _O^{{tr} ^{^T }} \\ - \end{array} - \label{introEquation:resistanceTensorTransformation} -\end{equation} -where -\[ -U_{OP} = \left( {\begin{array}{*{20}c} - 0 & { - z_{OP} } & {y_{OP} } \\ - {z_i } & 0 & { - x_{OP} } \\ - { - y_{OP} } & {x_{OP} } & 0 \\ -\end{array}} \right) -\] -Using Eq.~\ref{introEquation:definitionCR} and -Eq.~\ref{introEquation:resistanceTensorTransformation}, one can -locate the position of center of resistance, -\begin{eqnarray*} - \left( \begin{array}{l} - x_{OR} \\ - y_{OR} \\ - z_{OR} \\ - \end{array} \right) & = &\left( {\begin{array}{*{20}c} - {(\Xi _O^{rr} )_{yy} + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} } \\ - { - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz} + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} } \\ - { - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx} + (\Xi _O^{rr} )_{yy} } \\ -\end{array}} \right)^{ - 1} \\ - & & \left( \begin{array}{l} - (\Xi _O^{tr} )_{yz} - (\Xi _O^{tr} )_{zy} \\ - (\Xi _O^{tr} )_{zx} - (\Xi _O^{tr} )_{xz} \\ - (\Xi _O^{tr} )_{xy} - (\Xi _O^{tr} )_{yx} \\ - \end{array} \right) \\ -\end{eqnarray*} -where $x_OR$, $y_OR$, $z_OR$ are the components of the vector +origin. Given the resistance tensor at an arbitrary origin $O$, and a +vector ,${\bf r}_{OP} = (x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we +can obtain the resistance tensor at $P$ by +\begin{eqnarray} +\label{introEquation:resistanceTensorTransformation} + \Xi_P^{tt} & = & \Xi_O^{tt} \notag \\ + \Xi_P^{tr} = \Xi_P^{rt} & = & \Xi_O^{tr} - {\bf U}_{OP} \Xi _O^{tt} \\ + \Xi_P^{rr} & = &\Xi_O^{rr} - {\bf U}_{OP} \Xi_O^{tt} {\bf U}_{OP} ++ \Xi_O^{tr} {\bf U}_{OP} - {\bf U}_{OP} \left( \Xi_O^{tr} +\right)^{^T} \notag +\end{eqnarray} +where ${\bf U}_{OP}$ is the skew matrix (Eq. (\ref{eq:skewMatrix})) +for the vector between the origin $O$ and the point $P$. Using +Eqs.~\ref{introEquation:definitionCR}~and~\ref{introEquation:resistanceTensorTransformation}, +one can locate the position of center of resistance, +\begin{equation*} +\left(\begin{array}{l} +x_{OR} \\ +y_{OR} \\ +z_{OR} +\end{array}\right) = +\left(\begin{array}{*{20}c} +(\Xi_O^{rr})_{yy} + (\Xi_O^{rr})_{zz} & -(\Xi_O^{rr})_{xy} & -(\Xi_O^{rr})_{xz} \\ +-(\Xi_O^{rr})_{xy} & (\Xi_O^{rr})_{zz} + (\Xi_O^{rr})_{xx} & -(\Xi_O^{rr})_{yz} \\ +-(\Xi_O^{rr})_{xz} & -(\Xi_O^{rr})_{yz} & (\Xi_O^{rr})_{xx} + (\Xi_O^{rr})_{yy} \\ +\end{array}\right)^{-1} +\left(\begin{array}{l} +(\Xi_O^{tr})_{yz} - (\Xi_O^{tr})_{zy} \\ +(\Xi_O^{tr})_{zx} - (\Xi_O^{tr})_{xz} \\ +(\Xi_O^{tr})_{xy} - (\Xi_O^{tr})_{yx} +\end{array}\right) +\end{equation*} +where $x_{OR}$, $y_{OR}$, $z_{OR}$ are the components of the vector joining center of resistance $R$ and origin $O$. +For a general rigid molecular substructure, finding the $6 \times 6$ +resistance tensor can be a computationally demanding task. First, a +lattice of small beads that extends well beyond the boundaries of the +rigid substructure is created. The lattice is typically composed of +0.25 \AA\ beads on a dense FCC lattice. The lattice constant is taken +to be the bead diameter, so that adjacent beads are touching, but do +not overlap. To make a shape corresponding to the rigid structure, +beads that sit on lattice sites that are outside the van der Waals +radii of all of the atoms comprising the rigid body are excluded from +the calculation. +For large structures, most of the beads will be deep within the rigid +body and will not contribute to the hydrodynamic tensor. In the {\it +rough shell} approach, beads which have all of their lattice neighbors +inside the structure are considered interior beads, and are removed +from the calculation. After following this procedure, only those +beads in direct contact with the van der Waals surface of the rigid +body are retained. For reasonably large molecular structures, this +truncation can still produce bead assemblies with thousands of +members. + +If all of the {\it atoms} comprising the rigid substructure are +spherical and non-overlapping, the tensor in +Eq.~(\ref{introEquation:RPTensorNonOverlapped}) may be used directly +using the atoms themselves as the hydrodynamic beads. This is a +variant of the {\it bead model} approach of Carrasco and Garc\'{i}a de +la Torre.\cite{Carrasco1999} In this case, the size of the ${\bf B}$ +matrix can be quite small, and the calculation of the hydrodynamic +tensor is straightforward. + +In general, the inversion of the ${\bf B}$ matrix is the most +computationally demanding task. This inversion is done only once for +each type of rigid structure. We have used straightforward +LU-decomposition to solve the linear system and to obtain the elements +of ${\bf C}$. Once ${\bf C}$ has been obtained, the location of the +center of resistance ($R$) is found and the resistance tensor at this +point is calculated. The $3 \times 1$ vector giving the location of +the rigid body's center of resistance and the $6 \times 6$ resistance +tensor are then stored for use in the Langevin dynamics calculation. +These quantities depend on solvent viscosity and temperature and must +be recomputed if different simulation conditions are required. + \section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} + Consider the Langevin equations of motion in generalized coordinates \begin{equation} -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) \\ - \end{array} \right), +{\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} -while the random force $F_{r,i}^l$ is a Gaussian stochastic variable -with zero mean and variance +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} -\left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle = -\left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle = -2k_B T\Xi _R \delta (t - t'). \label{randomForce} -\end{equation} -The equation of motion for $v_i$ can be written as +{\bf F}_{f}^{~b}(t) = \left( \begin{array}{l} + {\bf f}_{f}^{~b}(t) \\ + {\bf \tau}_{f}^{~b}(t) \\ + \end{array} \right) = - \left( \begin{array}{*{20}c} + \Xi_{R}^{tt} & \Xi_{R}^{rt} \\ + \Xi_{R}^{tr} & \Xi_{R}^{rr} \\ +\end{array} \right)\left( \begin{array}{l} + {\bf v}_{R}^{~b}(t) \\ + {\bf \omega}(t) \\ + \end{array} \right), +\end{equation} +while the random force, ${\bf F}_{r}$, is a Gaussian stochastic +variable with zero mean and variance, \begin{equation} -m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) + -f_{r,i}^l (t) +\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} -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 +$\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} -\tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b +\langle {\bf Z}_i \rangle = 0, \hspace{1in} \langle {\bf Z}_i \cdot +{\bf Z}_j \rangle = \frac{2 k_B T}{\delta t} \delta_{ij}, \end{equation} -where $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 +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} -\dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t) -+ \tau _{r,i}^b(t) +m \dot{{\bf v}} (t) = {\bf f}_{s}(t) + {\bf f}_{f}(t) + +{\bf f}_{r}(t) \end{equation} -Embedding the friction terms into force and torque, one can -integrate the langevin equations of motion for rigid body of -arbitrary shape in a velocity-Verlet style 2-part algorithm, where -$h= \delta t$: +Since the frictional and random forces are applied at the center of +resistance, which generally does not coincide with the center of mass, +extra torques are exerted at the center of mass. Thus, the net +body-fixed torque at the center of mass, $\tau^{~b}(t)$, +is given by +\begin{equation} +\tau^{~b} \leftarrow \tau_{s}^{~b} + \tau_{f}^{~b} + \tau_{r}^{~b} + {\bf r}_{MR} \times \left( {\bf f}_{f}^{~b} + {\bf f}_{r}^{~b} \right) +\end{equation} +where ${\bf r}_{MR}$ is the vector from the center of mass to the center of +resistance. Instead of integrating the angular velocity in lab-fixed +frame, we consider the equation of motion for the angular momentum +(${\bf j}$) in the body-fixed frame +\begin{equation} +\frac{\partial}{\partial t}{\bf j}(t) = \tau^{~b}(t) +\end{equation} +Embedding the friction and random forces into the the total force and +torque, one can integrate the Langevin equations of motion for a rigid +body of arbitrary shape in a velocity-Verlet style 2-part algorithm, +where $h = \delta t$: -{\tt 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), \\ @@ -519,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} @@ -560,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) @@ -584,7 +706,7 @@ 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{Validating the Method\label{sec:validating}} @@ -647,29 +769,29 @@ simulation box. All simulations were equilibrated usi We performed reference microcanonical simulations with explicit solvents for each of the different model system. In each case there was one solute model and 1929 solvent molecules present in the -simulation box. All simulations were equilibrated using a +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 and sampling was done in a microcanonical -ensemble. Since the model bodies are typically quite massive, we were -able to use a time step of 25 fs. +further equilibration (5 ns) and sampling (10 ns) was done in a +microcanonical ensemble. Since the model bodies are typically quite +massive, we were able to use a time step of 25 fs. The model systems studied used both Lennard-Jones spheres as well as uniaxial Gay-Berne ellipoids. In its original form, the Gay-Berne potential was a single site model for the interactions of rigid -ellipsoidal molecules.\cite{Gay81} It can be thought of as a +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}({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, {\mathbf{\hat -r}_{ij}}) = 4\epsilon ({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, -{\mathbf{\hat r}_{ij}})\left[\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u +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}, -{\mathbf{\hat u}_j}, {\mathbf{\hat r}_{ij}})+\sigma_0}\right)^{12} --\left(\frac{\sigma_0}{r_{ij}-\sigma({\mathbf{\hat u}_i}, {\mathbf{\hat u}_j}, -{\mathbf{\hat r}_{ij}})+\sigma_0}\right)^6\right] +{{\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*} @@ -686,7 +808,7 @@ are given elsewhere,\cite{Luckhurst90,Golubkov06,SunGe Additionally, a well depth aspect ratio, $\epsilon^r = \epsilon^e / \epsilon^s$, describes the ratio between the well depths in the {\it end-to-end} and side-by-side configurations. Details of the potential -are given elsewhere,\cite{Luckhurst90,Golubkov06,SunGezelter08} and an +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} @@ -721,21 +843,8 @@ A similar form exists for the bulk viscosity \int_{t_0}^{t_0 + t} P_{xz}(t') dt' \right)^2 \right\rangle_{t_0}. \label{eq:shear} \end{equation} -A similar form exists for the bulk viscosity -\begin{equation} -\kappa = \lim_{t->\infty} \frac{V}{2 k_B T} \frac{d}{dt} \left\langle \left( -\int_{t_0}^{t_0 + t} -\left(P\left(t'\right)-\left\langle P \right\rangle \right)dt' -\right)^2 \right\rangle_{t_0}. -\end{equation} -Alternatively, the shear viscosity can also be calculated using a -Green-Kubo formula with the off-diagonal pressure tensor correlation function, -\begin{equation} -\eta = \frac{V}{k_B T} \int_0^{\infty} \left\langle P_{xz}(t_0) P_{xz}(t_0 -+ t) \right\rangle_{t_0} dt, -\end{equation} -although this method converges extremely slowly and is not practical -for obtaining viscosities from molecular dynamics simulations. +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 @@ -761,10 +870,10 @@ have used in this study, there are differences between compute the diffusive behavior for motion parallel to each body-fixed axis by projecting the displacement of the particle onto the body-fixed reference frame at $t=0$. With an isotropic solvent, as we -have used in this study, there are differences between the three -diffusion constants, but these must converge to the same value at -longer times. Translational diffusion constants for the different -shaped models are shown in table \ref{tab:translation}. +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 @@ -844,14 +953,14 @@ hydrodynamic flows is well known, giving translational an arbitrary value of 0.8 kcal/mol. The Stokes-Einstein behavior of large spherical particles in -hydrodynamic flows is well known, giving translational friction -coefficients of $6 \pi \eta R$ (stick boundary conditions) and -rotational friction coefficients of $8 \pi \eta R^3$. Recently, -Schmidt and Skinner have computed the behavior of spherical tag -particles in molecular dynamics simulations, and have shown that {\it -slip} boundary conditions ($\Xi_{tt} = 4 \pi \eta R$) may be more -appropriate for molecule-sized spheres embedded in a sea of spherical -solvent particles.\cite{Schmidt:2004fj,Schmidt:2003kx} +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 @@ -876,20 +985,19 @@ D = \frac{k_B T}{6 \pi \eta a} G(\rho), can be combined to give a single translational diffusion constant,\cite{Berne90} \begin{equation} -D = \frac{k_B T}{6 \pi \eta a} G(\rho), +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 - \rho^2) -G(\rho) - 1}{1 - \rho^4} \right\}. +\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(\rho)$ is a function of the axial ratio -($\rho = b / a$), which for prolate ellipsoids, is +In these expressions, $G(s)$ is a function of the axial ratio +($s = b / a$), which for prolate ellipsoids, is \begin{equation} -G(\rho) = (1- \rho^2)^{-1/2} \ln \left\{ \frac{1 + (1 - -\rho^2)^{1/2}}{\rho} \right\} +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 @@ -914,17 +1022,18 @@ The translational diffusion constants from the microca exact treatment of the diffusion tensor as well as the rough-shell model for the ellipsoid. -The translational diffusion constants from the microcanonical simulations -agree well with the predictions of the Perrin model, although the rotational -correlation times are a factor of 2 shorter than expected from hydrodynamic -theory. One explanation for the slower rotation -of explicitly-solvated ellipsoids is the possibility that solute-solvent -collisions happen at both ends of the solute whenever the principal -axis of the ellipsoid is turning. In the upper portion of figure -\ref{fig:explanation} we sketch a physical picture of this explanation. -Since our Langevin integrator is providing nearly quantitative agreement with -the Perrin model, it also predicts orientational diffusion for ellipsoids that -exceed explicitly solvated correlation times by a factor of two. +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 @@ -932,10 +1041,9 @@ model. Equation (\ref{introEquation:oseenTensor}) abov two-sphere dumbbell model. This model consists of two non-overlapping spheres held by a rigid bond connecting their centers. There are competing expressions for the 6x6 resistance tensor for this -model. Equation (\ref{introEquation:oseenTensor}) above gives the -original Oseen tensor, while the second order expression introduced by -Rotne and Prager,\cite{Rotne1969} and improved by Garc\'{i}a de la -Torre and Bloomfield,\cite{Torre1977} is given above as +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 @@ -948,7 +1056,7 @@ those derived from the 6 x 6 tensors mentioned above). motion in a flow {\it perpendicular} to the inter-sphere axis.\cite{Davis:1969uq} We know of no analytic solutions for the {\it orientational} correlation times for this model system (other than -those derived from the 6 x 6 tensors mentioned above). +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 @@ -1016,8 +1124,6 @@ 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 @@ -1029,60 +1135,91 @@ A reference system composed of a single banana rigid b behavior of this model, we have left out the dipolar interactions of the original Orlandi model. -A reference system composed of a single banana rigid body embedded in a -sea of 1929 solvent particles was created and run under standard -(microcanonical) molecular dynamics. The resulting viscosity of this -mixture was 0.298 centipoise (as estimated using Eq. (\ref{eq:shear})). -To calculate the hydrodynamic properties of the banana rigid body model, -we created a rough shell (see Fig.~\ref{fig:roughShell}), in which -the banana is represented as a ``shell'' made of 3321 identical beads -(0.25 \AA\ in diameter) distributed on the surface. Applying the -procedure described in Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we -identified the center of resistance, ${\bf r} = $(0 \AA, 0.81 \AA, 0 \AA), as -well as the resistance tensor, -\begin{equation*} -\Xi = -\left( {\begin{array}{*{20}c} -0.9261 & 0 & 0&0&0.08585&0.2057\\ -0& 0.9270&-0.007063& 0.08585&0&0\\ -0&-0.007063&0.7494&0.2057&0&0\\ -0&0.0858&0.2057& 58.64& 0&0\\0.08585&0&0&0&48.30&3.219&\\0.2057&0&0&0&3.219&10.7373\\\end{array}} \right), -\end{equation*} -where the units for translational, translation-rotation coupling and -rotational tensors are (kcal fs / mol \AA$^2$), (kcal fs / mol \AA\ rad), -and (kcal fs / mol rad$^2$), respectively. +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. +The Langevin rigid-body integrator (and the hydrodynamic diffusion +tensor) are essentially quantitative for translational diffusion of +this model. Orientational correlation times under the Langevin +rigid-body integrator are within 11\% of the values obtained from +explicit solvent, but these models also exhibit some solvent +inaccessible surface area in the explicitly-solvated case. \subsection{Composite sphero-ellipsoids} -Spherical heads perched on the ends of Gay-Berne ellipsoids have been -used recently as models for lipid molecules.\cite{SunGezelter08,Ayton01} -MORE DETAILS +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} -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. +\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{Translational diffusion constants (D) for the model systems calculated using microcanonical simulations (with explicit solvent), theoretical predictions, and Langevin simulations (with implicit solvent). -Analytical solutions for the exactly-solved hydrodynamics models are -from Refs. \citen{Einstein05} (sphere), \citen{Perrin1934} and \citen{Perrin1936} +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 (= +All diffusion constants are reported in units of $10^{-3}$ cm$^2$ / ps (= $10^{-4}$ \AA$^2$ / fs). } \begin{tabular}{lccccccc} \hline @@ -1090,13 +1227,13 @@ sphere & 0.261 & ? & & 2.59 & exact & 2.5 \cline{2-3} \cline{5-7} model & $\eta$ (centipoise) & D & & Analytical & method & Hydrodynamics & simulation \\ \hline -sphere & 0.261 & ? & & 2.59 & exact & 2.59 & 2.56 \\ +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.322 & ? & & 1.57 & bead model & 1.57 & 1.57 \\ - & 0.322 & ? & & 1.57 & rough shell & ? & ? \\ +dumbbell & 0.308 & 2.06 & & 1.64 & bead model & 1.65 & 1.62 \\ + & 0.308 & 2.06 & & 1.64 & rough shell & 1.59 & 1.62 \\ banana & 0.298 & 1.53 & & & rough shell & 1.56 & 1.55 \\ -lipid & 0.349 & 0.96 & & & rough shell & 1.33 & 1.32 \\ +lipid & 0.349 & 1.41 & & & rough shell & 1.33 & 1.32 \\ \end{tabular} \label{tab:translation} \end{center} @@ -1119,11 +1256,11 @@ sphere & 0.261 & & & 9.06 & exact & 9.0 \cline{2-3} \cline{5-7} model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic & simulation \\ \hline -sphere & 0.261 & & & 9.06 & exact & 9.06 & 9.11 \\ +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.322 & 14.0 & & & bead model & 52.3 & 52.8 \\ - & 0.322 & 14.0 & & & rough shell & ? & ? \\ +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 @@ -1135,49 +1272,86 @@ The Langevin dynamics integrator was applied to study \section{Application: A rigid-body lipid bilayer} -The Langevin dynamics integrator was applied to study the formation of -corrugated structures emerging from simulations of the coarse grained -lipid molecular models presented above. The initial configuration is -taken from our molecular dynamics studies on lipid bilayers with -lennard-Jones sphere solvents. The solvent molecules were excluded -from the system, and the experimental value for the viscosity of water -at 20C ($\eta = 1.00$ cp) was used to mimic the hydrodynamic effects -of the solvent. The absence of explicit solvent molecules and the -stability of the integrator allowed us to take timesteps of 50 fs. A -total simulation run time of 100 ns was sampled. -Fig. \ref{fig:bilayer} shows the configuration of the system after 100 -ns, and the ripple structure remains stable during the entire -trajectory. Compared with using explicit bead-model solvent -molecules, the efficiency of the simulation has increased by an order -of magnitude. +To test the accuracy and efficiency of the new integrator, we applied +it to study the formation of corrugated structures emerging from +simulations of the coarse grained lipid molecular models presented +above. The initial configuration is taken from earlier molecular +dynamics studies on lipid bilayers which had used spherical +(Lennard-Jones) solvent particles and moderate (480 solvated lipid +molecules) system sizes.\cite{SunX._jp0762020} the solvent molecules +were excluded from the system and the box was replicated three times +in the x- and y- axes (giving a single simulation cell containing 4320 +lipids). The experimental value for the viscosity of water at 20C +($\eta = 1.00$ cp) was used with the Langevin integrator to mimic the +hydrodynamic effects of the solvent. The absence of explicit solvent +molecules and the stability of the integrator allowed us to take +timesteps of 50 fs. A simulation run time of 30 ns was sampled to +calculate structural properties. Fig. \ref{fig:bilayer} shows the +configuration of the system after 30 ns. Structural properties of the +bilayer (e.g. the head and body $P_2$ order parameters) are nearly +identical to those obtained via solvated molecular dynamics. The +ripple structure remained stable during the entire trajectory. +Compared with using explicit bead-model solvent molecules, the 30 ns +trajectory for 4320 lipids with the Langevin integrator is now {\it +faster} on the same hardware than the same length trajectory was for +the 480-lipid system previously studied. \begin{figure} \centering \includegraphics[width=\linewidth]{bilayer} \caption[Snapshot of a bilayer of rigid-body models for lipids]{A -snapshot of a bilayer composed of rigid-body models for lipid +snapshot of a bilayer composed of 4320 rigid-body models for lipid molecules evolving using the Langevin integrator described in this work.} \label{fig:bilayer} \end{figure} \section{Conclusions} -We have presented a new Langevin algorithm by incorporating the -hydrodynamics properties of arbitrary shaped molecules into an -advanced symplectic integration scheme. Further studies in systems -involving banana shaped molecules illustrated that the dynamic -properties could be preserved by using this new algorithm as an -implicit solvent model. +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{jcp} +\bibliographystyle{jcp2} \bibliography{langevin} - \end{document}