--- trunk/langevin/langevin.tex 2008/01/23 22:36:36 3332 +++ trunk/langevin/langevin.tex 2008/01/24 14:16:07 3333 @@ -68,7 +68,7 @@ Kramers's theory, Klimov and Thirumalai identified fre 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's theory, Klimov and Thirumalai identified free-energy +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 @@ -85,8 +85,7 @@ individual atoms are taken from the Stokes-Einstein hy finite structures.\cite{Berkov2005a}] In typical LD simulations, the friction and random forces on -individual atoms are taken from the Stokes-Einstein hydrodynamic -approximation, +individual atoms are taken from Stokes' law, \begin{eqnarray} m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + R(t) \\ \langle R(t) \rangle & = & 0 \\ @@ -95,130 +94,115 @@ The use of rigid substructures,\cite{???} where $\xi \approx 6 \pi \eta a$. Here $\eta$ is the viscosity of the implicit solvent, and $a$ is the hydrodynamic radius of the atom. -The use of rigid substructures,\cite{???} -coarse-graining,\cite{Ayton,Sun,Zannoni} and ellipsoidal -representations of protein side chains~\cite{Schulten} 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 +The use of rigid substructures,\cite{Chun:2000fj} +coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunGezelter08} +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. The atoms involved in a rigid or coarse-grained structure should properly have solvent-mediated interactions with each other. 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{MarshallNewton,GarciaDeLaTorre} There a need to have a -more thorough treatment of hydrodynamics included in the library of -methods available for performing Langevin simulations. +motion.\cite{FIXMAN:1986lr,Ramachandran1996} +In order to account for the diffusion anisotropy of arbitrarily-shaped +particles, Fernandes and Garc\'{i}a de la Torre improved the original +Brownian dynamics simulation algorithm~\cite{Ermak1978,Allison1991} by +incorporating a generalized $6\times6$ diffusion tensor and +introducing a rotational evolution scheme consisting of three +consecutive rotations.\cite{Fernandes2002} Unfortunately, biases are +introduced into the system due to the arbitrary order of applying the +noncommuting rotation operators.\cite{Beard2003} Based on the +observation the momentum relaxation time is much less than the time +step, one may ignore the inertia in Brownian dynamics. However, the +assumption of zero average acceleration is not always true for +cooperative motion which is common in proteins. An inertial Brownian +dynamics (IBD) was proposed to address this issue by adding an +inertial correction term.\cite{Beard2000} As a complement to IBD which +has a lower bound in time step because of the inertial relaxation +time, long-time-step inertial dynamics (LTID) can be used to +investigate the inertial behavior of linked polymer segments in a low +friction regime.\cite{Beard2000} LTID can also deal with the +rotational dynamics for nonskew bodies without translation-rotation +coupling by separating the translation and rotation motion and taking +advantage of the analytical solution of hydrodynamics +properties. However, typical nonskew bodies like cylinders and +ellipsoids are inadequate to represent most complex macromolecular +assemblies. There is therefore a need for incorporating the +hydrodynamics of complex (and potentially skew) rigid bodies in the +library of methods available for performing Langevin simulations. + \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{Sun2008} -Many of the water models in common use are also rigid-body -models,\cite{TIPs,SPC/E} although they are typically evolved using -constraints rather than rigid body equations of motion. +docking,\cite{Gray2003} and lipid bilayer +simulations.\cite{SunGezelter08} Many of the water models in common +use are also rigid-body +models,\cite{Jorgensen83,Berendsen81,Berendsen87} although they are +typically evolved 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 $1 \over \sin \theta$ -singularities, the numerical integration of corresponding equations of -these motion can become inaccurate (and inefficient). Although an -alternative integrator using 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} Unfortunately, this approach uses a nonseparable -Hamiltonian resulting from the quaternion representation, which -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 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} +Euler angles are a natural choice to describe the rotational degrees +of freedom. However, due to $1 \over \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} -A breakthrough 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 suffer from some related inefficiencies. A symplectic -Lie-Poisson integrator for rigid bodies developed by Dullweber {\it et -al.}\cite{Dullweber1997} gets around most of the limitations mentioned -above and has become the basis for our Langevin integrator. +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 +$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 $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 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. -\subsection{The Hydrodynamic tensor and Brownian dynamics} -Combining Brownian dynamics with rigid substructures, one can study -slow processes in biomolecular systems. Modeling DNA as a chain of -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 Noguchi and others.\cite{Noguchi2001,Noguchi2002,Shillcock2005} - -Due to the difficulty of numerically integrating anisotropic -rotational motion, most of the coarse-grained rigid body models are -treated as spheres, cylinders, ellipsoids or other regular shapes in -stochastic simulations. In an effort to account for the diffusion -anisotropy of arbitrarily-shaped particles, Fernandes and Garc\'{i}a -de la Torre improved the original Brownian dynamics simulation -algorithm~\cite{Ermak1978,Allison1991} by incorporating a generalized -$6\times6$ diffusion tensor and introducing a rotational evolution -scheme consisting of three consecutive rotations.\cite{Fernandes2002} -Unfortunately, biases are introduced into the system due to the -arbitrary order of applying the noncommuting rotation -operators.\cite{Beard2003} Based on the observation the momentum -relaxation time is much less than the time step, one may ignore the -inertia in Brownian dynamics. However, the assumption of zero average -acceleration is not always true for cooperative motion which is common -in proteins. An inertial Brownian dynamics (IBD) was proposed to -address this issue by adding an inertial correction -term.\cite{Beard2000} As a complement to IBD which has a lower bound -in time step because of the inertial relaxation time, long-time-step -inertial dynamics (LTID) can be used to investigate the inertial -behavior of 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 macromolecular assemblies. - 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 a symplectic rigid body dynamics propagator. In the sections below, -we review some of the theory of hydrodynamic tensors developed for -Brownian simulations of rigid 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. +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 can be determined using the velocity autocorrelation function. However, this approach becomes -impractical when the solute becomes complex. Instead, various +impractical when the solute becomes complex. 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 \begin{equation} -\Xi = \left( {\begin{array}{*{20}c} - {\Xi _{}^{tt} } & {\Xi _{}^{rt} } \\ - {\Xi _{}^{tr} } & {\Xi _{}^{rr} } \\ -\end{array}} \right). +\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 @@ -231,10 +215,10 @@ fluid, it may experience friction force ($\mathbf{F}_f \left( \begin{array}{l} \mathbf{F}_f \\ \mathbf{\tau}_f \\ - \end{array} \right) = - \left( {\begin{array}{*{20}c} - \Xi ^{tt} & \Xi ^{rt} \\ - \Xi ^{tr} & \Xi ^{rr} \\ -\end{array}} \right)\left( \begin{array}{l} + \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). @@ -243,21 +227,21 @@ Stoke's law, \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}} For a spherical particle under ``stick'' boundary conditions, the translational and rotational friction tensors can be calculated from -Stoke's law, +Stokes' law, \begin{equation} -\Xi^{tt} = \left( {\begin{array}{*{20}c} +\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) +\end{array} \right) \end{equation} and \begin{equation} -\Xi ^{rr} = \left( {\begin{array}{*{20}c} +\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) +\end{array} \right) \end{equation} where $\eta$ is the viscosity of the solvent and $R$ is the hydrodynamic radius. @@ -265,41 +249,38 @@ extended Stokes's law to general ellipsoids, also call Other non-spherical shapes, such as cylinders and ellipsoids, are widely used as references for developing new hydrodynamics 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} +extended Stokes' law to general ellipsoids which are given in +Cartesian coordinates by~\cite{Perrin1934,Perrin1936} \begin{equation} -\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1 +\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1. \end{equation} -where the semi-axes are of lengths $a$, $b$, and $c$. Due to the -complexity of the elliptic integral, only uniaxial ellipsoids, -{\it 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 : +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$), can be 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 ellipsoids: +and oblate, \begin{equation} S = \frac{2}{\sqrt {b^2 - a^2 }} \arctan \frac{\sqrt {b^2 - a^2}}{a}, \end{equation} -one can write down the translational and rotational resistance -tensors for oblate, +ellipsoids, 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 prolate, +for oblate, 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*} -ellipsoids. For both spherical and ellipsoidal particles, the -translation-rotation and rotation-translation coupling tensors are +for prolate ellipsoids. 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}} - 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 @@ -315,16 +296,17 @@ interaction tensor. Let us consider a rigid assembly o 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$, -\[ +interaction tensor. + +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$, +\begin{equation} 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 +\end{equation} +where $F_i$ is the frictional force, and $T_{ij}$ is the hydrodynamic +interaction tensor. The frictional force on the $i^\mathrm{th}$ bead +is proportional to its ``net'' velocity \begin{equation} F_i = \zeta _i v_i - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }. \label{introEquation:tensorExpression} @@ -361,17 +343,17 @@ B = \left( {\begin{array}{*{20}c} construct a $3N \times 3N$ matrix consisting of $N \times N$ $B_{ij}$ blocks \begin{equation} -B = \left( {\begin{array}{*{20}c} - {B_{11} } & \ldots & {B_{1N} } \\ +B = \left( \begin{array}{*{20}c} + B_{11} & \ldots & B_{1N} \\ \vdots & \ddots & \vdots \\ - {B_{N1} } & \cdots & {B_{NN} } \\ -\end{array}} \right), + B_{N1} & \cdots & B_{NN} \\ +\end{array} \right), \end{equation} where $B_{ij}$ is given by -\[ +\begin{equation} B_{ij} = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij} )T_{ij} -\] +\end{equation} where $\delta _{ij}$ is the Kronecker delta function. Inverting the $B$ matrix, we obtain \[ @@ -456,11 +438,11 @@ locate the position of center of resistance, x_{OR} \\ y_{OR} \\ z_{OR} \\ - \end{array} \right) & = &\left( {\begin{array}{*{20}c} + \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} \\ +\end{array} \right)^{ - 1} \\ & & \left( \begin{array}{l} (\Xi _O^{tr} )_{yz} - (\Xi _O^{tr} )_{zy} \\ (\Xi _O^{tr} )_{zx} - (\Xi _O^{tr} )_{xz} \\ @@ -474,72 +456,80 @@ M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)} + F_{r,i} \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) +\mathbf{M}_i \dot \mathbf{V}_i(t) = \mathbf{F}_{s,i}(t) + \mathbf{F}_{f,i}(t) + \mathbf{R}_{i}(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: -\[ +where $\mathbf{M}_i$ is a $6\times6$ diagonal mass matrix (which +includes the rigid body mass and moments of inertia) and $\mathbf{V}_i$ is a +generalized velocity, $\mathbf{V}_i = +\left\{\mathbf{v}_i,\mathbf{\omega}_i \right\}$. The right side of +Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a +system force $\mathbf{F}_{s,i}$, a frictional or dissipative force +$\mathbf{F}_{f,i}$ and stochastic force $\mathbf{R}_{i}$. While the +evolution of the system in Newtownian mechanics is typically done in the +lab-fixed frame, it is convenient to handle the rotation of rigid +bodies in the body-fixed frame. Thus the friction and random forces are +calculated in body-fixed frame and converted back to lab-fixed frame +using the rigid body's rotation matrix ($Q_i$): +\begin{equation} \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). \\ + \mathbf{F}_{f,i}(t) = Q_{i}^{T} \mathbf{F}_{f,i}^b (t), \\ + \mathbf{R}_{i}(t) = Q_{i}^{T} \mathbf{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$ +\end{equation} +Here, the body-fixed friction force $\mathbf{F}_{f,i}^b$ is proportional to +the body-fixed velocity at the center of resistance $\mathbf{v}_{R,i}^b$ and +angular velocity $\mathbf{\omega}_i$ \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) \\ +\mathbf{F}_{f,i}^b (t) = \left( \begin{array}{l} + \mathbf{f}_{f,i}^b (t) \\ + \mathbf{\tau}_{f,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} + \mathbf{v}_{R,i}^b (t) \\ + \mathbf{\omega}_i (t) \\ \end{array} \right), \end{equation} -while the random force $F_{r,i}^l$ is a Gaussian stochastic variable +while the random force $\mathbf{R}_{i}^l$ 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 {\mathbf{R}_{i}^l (t) (\mathbf{R}_{i}^l (t'))^T } \right\rangle = +\left\langle {\mathbf{R}_{i}^b (t) (\mathbf{R}_{i}^b (t'))^T } \right\rangle = +2 k_B T \Xi_R \delta(t - t'). \label{randomForce} \end{equation} -The equation of motion for $v_i$ can be written as +Once the $6\times6$ resistance tensor at the center of resistance +($\Xi_R$) is known, 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 $\Xi_R$.\cite{SchlickBook} Each time a random force +vector is needed, a gaussian random vector is generated and then the +square root matrix is multiplied onto this vector. + +The equation of motion for $\mathbf{v}_i$ can be written as \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) +m\dot \mathbf{v}_i (t) = \mathbf{f}_{s,i} (t) + \mathbf{f}_{f,i}^l (t) + +\mathbf{R}_{i}^l (t) \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 +frictional torque at the center of mass, $\tau_{f,i}^b (t)$, is given by \begin{equation} -\tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b +\tau_{f,i}^b \leftarrow \tau_{f,i}^b + \mathbf{r}_{MR} \times \mathbf{f}_{r,i}^b \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 \begin{equation} -\dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t) -+ \tau _{r,i}^b(t) +\dot j_i (t) = \tau_{s,i} (t) + \tau_{f,i}^b (t) + \mathbf{R}_{i}^b(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$: +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$: {\tt moveA:} \begin{align*}