--- trunk/tengDissertation/Langevin.tex 2006/06/25 19:43:39 2887 +++ trunk/tengDissertation/Langevin.tex 2006/07/18 14:19:49 2945 @@ -4,107 +4,94 @@ As an excellent alternative to newtonian dynamics, Lan \section{Introduction} %applications of langevin dynamics -As an excellent 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 simulation. Implicit solvent Langevin -dynamics simulation of met-enkephalin not only outperforms explicit -solvent simulation on computation efficiency, but also agrees very -well with explicit solvent simulation on dynamics -properties\cite{Shen2002}. Recently, applying Langevin dynamics with -UNRES model, Liow and his coworkers suggest that protein folding -pathways can be possibly exploited 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 barrier\cite{Banerjee2004, Cui2003}. -Combining Langevin dynamics with Kramers's theory, Klimov and -Thirumalai identified the free-energy barrier 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 study 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}. In an attempt to reduce the -computational cost of simulation, multiple time stepping (MTS) -methods have been introduced and have been of great interest to -macromolecule and protein community\cite{Tuckerman1992}. Relying on -the observation that forces between distant atoms generally -demonstrate slower fluctuations than forces between close atoms, MTS -method are generally implemented by evaluating the slowly -fluctuating forces less frequently than the fast ones. -Unfortunately, nonlinear instability resulting from increasing -timestep in MTS simulation have became a critical obstruction -preventing the long time simulation. Due to the coupling to the heat -bath, Langevin dynamics has been shown to be able to damp out the -resonance artifact more efficiently\cite{Sandu1999}. +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} %review langevin/browninan dynamics for arbitrarily shaped rigid body Combining Langevin or Brownian dynamics with rigid body dynamics, -one can study the slow processes in biomolecular systems. Modeling -the DNA as a chain of rigid spheres beads, which subject to harmonic -potentials as well as excluded volume potentials, Mielke and his -coworkers discover 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 millisecond, which -is impracticable to study using all atomistic model with newtonian -mechanics. With the help of coarse-grained rigid body model and -stochastic dynamics, the fusion pathways were exploited by many -researchers\cite{Noguchi2001,Noguchi2002,Shillcock2005}. Due to the -difficulty of numerical integration of anisotropy rotation, most of -the rigid body models are simply modeled by sphere, cylinder, -ellipsoid or other regular shapes in stochastic simulations. In an -effort to account for the diffusion anisotropy of the arbitrary +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 incorporating a generalized $6\times6$ diffusion tensor and introducing a simple rotation evolution scheme consisting of three -consecutive rotations\cite{Fernandes2002}. Unfortunately, unexpected -error and bias are introduced into the system due to the arbitrary -order of applying the noncommuting rotation -operators\cite{Beard2003}. Based on the observation the momentum +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, assumption of the zero +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 +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 +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 cylinder and ellipsoid are inadequate to -represent most complex macromolecule assemblies. These intricate +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 -hydrodynamics properties can be calculated using variant -hydrodynamic interaction tensors. +hydrodynamic properties can be calculated using variants on the +standard hydrodynamic interaction tensors. The goal of the present work is to develop a Langevin dynamics -algorithm for arbitrary rigid particles by integrating the accurate -estimation of friction tensor from hydrodynamics theory into the -sophisticated rigid body 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. \section{Computational Methods{\label{methodSec}}} \subsection{\label{introSection:frictionTensor}Friction Tensor} Theoretically, the friction kernel can be determined using the -velocity autocorrelation function. However, this approach become -impractical when the system become more and more complicate. +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, +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} @@ -112,13 +99,13 @@ Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are translatio {\Xi _{}^{tr} } & {\Xi _{}^{rr} } \\ \end{array}} \right). \] -Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are 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, +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, \[ \left( \begin{array}{l} F_R \\ @@ -132,9 +119,9 @@ toque. \end{array} \right) \] where $F_r$ is the friction force and $\tau _R$ is the friction -toque. +torque. -\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}} +\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 @@ -155,49 +142,45 @@ hydrodynamics radius. \end{array}} \right) \] where $\eta$ is the viscosity of the solvent and $R$ is the -hydrodynamics radius. +hydrodynamic radius. -Other non-spherical shape, such as cylinder and ellipsoid -\textit{etc}, are widely used as reference for developing new -hydrodynamics theory, 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} +Other non-spherical shapes, such as cylinders and ellipsoids, are +widely used as references for developing new hydrodynamics theory, +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 having to be equal, \textit{i.e.} +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 : +and oblate ellipsoids: \[ S = \frac{2}{{\sqrt {b^2 - a^2 } }}arctg\frac{{\sqrt {b^2 - a^2 } -}}{a} -\], +}}{a}, +\] one can write down the translational and rotational resistance tensors -\[ -\begin{array}{l} - \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 + +\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{array} -\] +\end{eqnarray*} and -\[ -\begin{array}{l} - \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{array}. -\] +\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*} \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}} @@ -209,15 +192,15 @@ translational and rotational motion of rigid body, gen 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 body, general +translational and rotational motion of rigid bodies, general ellipsoids are not always suitable for modeling arbitrarily shaped -rigid molecule. A number of studies have been devoted to determine -the friction tensor for irregularly shaped rigid bodies using more -advanced method where the molecule of interest was modeled by -combinations of spheres(beads)\cite{Carrasco1999} and the +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 hydrodynamics +$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$, \[ @@ -231,8 +214,8 @@ This equation is the basis for deriving the hydrodynam \label{introEquation:tensorExpression} \end{equation} This equation is the basis for deriving the hydrodynamic tensor. In -1930, Oseen and Burgers gave a simple solution to Equation -\ref{introEquation:tensorExpression} +1930, Oseen and Burgers gave a simple solution to +Eq.~\ref{introEquation:tensorExpression} \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} @@ -240,7 +223,7 @@ Garc\'{i}a de la Torre and Bloomfield\cite{Torre1977}, 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 -Garc\'{i}a de la Torre and Bloomfield\cite{Torre1977}, +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 @@ -248,9 +231,9 @@ Both of the Equation \ref{introEquation:oseenTensor} a \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right]. \label{introEquation:RPTensorNonOverlapped} \end{equation} -Both of the Equation \ref{introEquation:oseenTensor} and Equation -\ref{introEquation:RPTensorNonOverlapped} have an assumption $R_{ij} -\ge \sigma _i + \sigma _j$. An alternative expression for +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 \begin{equation} T_{ij} = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 - @@ -258,7 +241,6 @@ T_{ij} = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( \frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right] \label{introEquation:RPTensorOverlapped} \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 @@ -274,18 +256,17 @@ where $\delta _{ij}$ is Kronecker delta function. Inve B_{ij} = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij} )T_{ij} \] -where $\delta _{ij}$ is Kronecker delta function. Inverting the $B$ -matrix, we obtain - +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) +\end{array}} \right), \] -, which can be partitioned into $N \times N$ $3 \times 3$ block -$C_{ij}$. With the help of $C_{ij}$ and skew matrix $U_i$ +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 } \\ @@ -294,30 +275,27 @@ bead $i$ and origin $O$. Hence, the elements of resist \end{array}} \right) \] where $x_i$, $y_i$, $z_i$ are the components of the vector joining -bead $i$ and origin $O$. Hence, the elements of resistance tensor at +bead $i$ and origin $O$, the elements of resistance tensor at arbitrary origin $O$ can be written as -\begin{equation} -\begin{array}{l} - \Xi _{}^{tt} = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\ - \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 \\ - \end{array} +\begin{eqnarray} + \Xi _{}^{tt} & = & \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\ + \Xi _{}^{tr} & = & \Xi _{}^{rt} = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\ + \Xi _{}^{rr} & = & - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j. \notag \\ \label{introEquation:ResistanceTensorArbitraryOrigin} -\end{equation} - +\end{eqnarray} The resistance tensor depends on the origin to which they refer. The -proper location for applying friction force is the center of +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 -tensor are symmetric, +tensors are symmetric, \begin{equation} \Xi^{tr} = \left( {\Xi^{tr} } \right)^T \label{introEquation:definitionCR} \end{equation} From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin}, -we can easily find out that the translational resistance tensor is +we can easily derive that the translational resistance tensor is origin independent, while the rotational resistance tensor and translation-rotation coupling resistance tensor depend on the origin. Given the resistance tensor at an arbitrary origin $O$, and @@ -339,9 +317,9 @@ Using Equations \ref{introEquation:definitionCR} and { - y_{OP} } & {x_{OP} } & 0 \\ \end{array}} \right) \] -Using Equations \ref{introEquation:definitionCR} and -\ref{introEquation:resistanceTensorTransformation}, one can locate -the position of center of resistance, +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} \\ @@ -358,21 +336,20 @@ the position of center of resistance, (\Xi _O^{tr} )_{xy} - (\Xi _O^{tr} )_{yx} \\ \end{array} \right) \\ \end{eqnarray*} - where $x_OR$, $y_OR$, $z_OR$ are the components of the vector joining center of resistance $R$ and origin $O$. \subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}} -Consider a Langevin equation of motions in generalized coordinates +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) \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 +$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, @@ -381,13 +358,13 @@ in body-fixed frame and converted back to lab-fixed fr in body-fixed frame and converted back to lab-fixed frame by: \[ \begin{array}{l} - F_{f,i}^l (t) = A^T F_{f,i}^b (t), \\ - F_{r,i}^l (t) = A^T F_{r,i}^b (t) \\ - \end{array}. + 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$, +angular velocity $\omega _i$ \begin{equation} F_{r,i}^b (t) = \left( \begin{array}{l} f_{r,i}^b (t) \\ @@ -407,7 +384,6 @@ with zero mean and variance \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 \begin{equation} m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) + @@ -422,14 +398,13 @@ of the resistance. Instead of integrating angular velo \tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times 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 angular velocity in -lab-fixed frame, we consider the equation of motion of angular -momentum in body-fixed frame +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) \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 @@ -446,10 +421,9 @@ $h= \delta t$: {\bf j}\left(t + h / 2 \right) &\leftarrow {\bf j}(t) + \frac{h}{2} {\bf \tau}^b(t), \\ % -\mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j} +\mathsf{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, \begin{equation} @@ -458,13 +432,13 @@ rotates both the rotation matrix ($\mathsf{A}$) and th / 2) \cdot \mathsf{G}_x(a_x /2), \end{equation} where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, -rotates both the rotation matrix ($\mathsf{A}$) and the body-fixed +rotates both the rotation matrix ($\mathsf{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{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\ +\mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\ {\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0). \end{array} @@ -484,11 +458,10 @@ All other rotations follow in a straightforward manner \end{array} \right). \end{equation} -All other rotations follow in a straightforward manner. +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 -After the first part of the propagation, the forces and body-fixed -torques are calculated at the new positions and orientations - {\tt doForces:} \begin{align*} {\bf f}(t + h) &\leftarrow @@ -497,14 +470,11 @@ torques are calculated at the new positions and orient {\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h) \times \frac{\partial V}{\partial {\bf u}}, \\ % -{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h) +{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h) \cdot {\bf \tau}^s(t + h). \end{align*} - -{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix -$\mathsf{A}$ is calculated in {\tt moveA}. Once the forces and -torques have been obtained at the new time step, the velocities can -be advanced to the same time value. +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:} \begin{align*} @@ -520,8 +490,8 @@ implemented in {\sc oopse}\cite{Meineke2005} and appli \section{Results and Discussion} The Langevin algorithm described in previous section has been -implemented in {\sc oopse}\cite{Meineke2005} and applied to the -studies of kinetics and thermodynamic properties in several systems. +implemented in {\sc oopse}\cite{Meineke2005} and applied to studies +of the static and dynamic properties in several systems. \subsection{Temperature Control} @@ -532,23 +502,24 @@ Initial configuration for the simulations is taken fro the algorithm. In order to verify the stability of this new algorithm, a series of simulations are performed on system consisiting of 256 SSD water molecules with different viscosities. -Initial configuration for the simulations is taken from a 1ns NVT -simulation with a cubic box of 19.7166~\AA. All simulation are +The initial configuration for the simulations is taken from a 1ns +NVT simulation with a cubic box of 19.7166~\AA. All simulation are carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns -with reference temperature at 300~K. Average temperature as a +with reference temperature at 300~K. The average temperature as a function of $\eta$ is shown in Table \ref{langevin:viscosity} where the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 - 1$ poise. The better temperature control at higher viscosity can be explained by the finite size effect and relative slow relaxation rate at lower viscosity regime. \begin{table} -\caption{Average temperatures from Langevin dynamics simulations of -SSD water molecules with reference temperature at 300~K.} +\caption{AVERAGE TEMPERATURES FROM LANGEVIN DYNAMICS SIMULATIONS OF +SSD WATER MOLECULES WITH REFERENCE TEMPERATURE AT 300~K.} \label{langevin:viscosity} \begin{center} -\begin{tabular}{|l|l|l|} +\begin{tabular}{lll} \hline $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\ + \hline 1 & 300.47 & 10.99 \\ 0.1 & 301.19 & 11.136 \\ 0.01 & 303.04 & 11.796 \\ @@ -557,7 +528,7 @@ Another set of calculation were performed to study the \end{center} \end{table} -Another set of calculation were performed to study the efficiency of +Another set of calculations were performed to study the efficiency of temperature control using different temperature coupling schemes. The starting configuration is cooled to 173~K and evolved using NVE, NVT, and Langevin dynamic with time step of 2 fs. @@ -567,10 +538,9 @@ scaling to the desire temperature. In extremely lower Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps which gives reasonable tight coupling, while the blue one from Langevin dynamics with viscosity of 0.1 poise demonstrates a faster -scaling to the desire temperature. In extremely lower friction -regime (when $ \eta \approx 0$), Langevin dynamics becomes normal -NVE (see green curve in Fig.~\ref{langevin:temperature}) which loses -the temperature control ability. +scaling to the desire temperature. When $ \eta = 0$, Langevin dynamics becomes normal +NVE (see orange curve in Fig.~\ref{langevin:temperature}) which +loses the temperature control ability. \begin{figure} \centering @@ -579,37 +549,47 @@ temperature fluctuation versus time.} \label{langevin: temperature fluctuation versus time.} \label{langevin:temperature} \end{figure} -\subsection{Langevin Dynamics of Banana Shaped Molecule} +\subsection{Langevin Dynamics of Banana Shaped Molecules} In order to verify that Langevin dynamics can mimic the dynamics of the systems absent of explicit solvents, we carried out two sets of simulations and compare their dynamic properties. - Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation made of 256 pentane molecules and two banana shaped molecules at 273~K. It has an equivalent implicit solvent system containing only two banana shaped molecules with viscosity of 0.289 center poise. To calculate the hydrodynamic properties of the banana shaped molecule, -we create a rough shell model (see Fig.~\ref{langevin:roughShell}), +we created a rough shell model (see Fig.~\ref{langevin:roughShell}), in which the banana shaped molecule is represented as a ``shell'' -made of 2266 small identical beads with size of 0.3 $\AA$ on the +made of 2266 small identical beads with size of 0.3 \AA on the surface. Applying the procedure described in Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we -identified the center of resistance at $(0, 0.7482, -0.1988)$, as -well as the resistance tensor, +identified the center of resistance at (0 $\rm{\AA}$, 0.7482 $\rm{\AA}$, +-0.1988 $\rm{\AA}$), as well as the resistance tensor, \[ \left( {\begin{array}{*{20}c} -0.9261 & 1.310e-14 & -7.292e-15&5.067e-14&0.08585&0.2057\\ -3.968e-14& 0.9270&-0.007063& 0.08585&6.764e-14&4.846e-14\\ --6.561e-16&-0.007063&0.7494&0.2057&4.846e-14&1.5036e-14\\ -5.067e-14&0.0858&0.2057& 58.64& 8.563e-13&-8.5736\\ -0.08585&6.764e-14&4.846e-14&1.555e-12&48.30&3.219&\\ -0.2057&4.846e-14&1.5036e-14&-3.904e-13&3.219&10.7373\\ +0.9261 & 0 & 0&0&0.08585&0.2057\\ +0& 0.9270&-0.007063& 0.08585&0&0\\ +0&-0.007063&0.7494&0.2057&0&0\\ +0&0.0858&0.2057& 58.64& 0&0\\ +0.08585&0&0&0&48.30&3.219&\\ +0.2057&0&0&0&3.219&10.7373\\ \end{array}} \right). \] +where the units for translational, translation-rotation coupling and +rotational tensors are $\frac{kcal \cdot fs}{mol \cdot \rm{\AA}^2}$, +$\frac{kcal \cdot fs}{mol \cdot \rm{\AA} \cdot rad}$ and $\frac{kcal +\cdot fs}{mol \cdot rad^2}$ respectively. Curves of the velocity +auto-correlation functions in Fig.~\ref{langevin:vacf} were shown to +match each other very well. However, because of the stochastic +nature, simulation using Langevin dynamics was shown to decay +slightly faster than MD. In order to study the rotational motion of +the molecules, we also calculated the auto-correlation function of +the principle axis of the second GB particle, $u$. The discrepancy +shown in Fig.~\ref{langevin:uacf} was probably due to the reason +that we used the experimental viscosity directly instead of +calculating bulk viscosity from simulation. - - \begin{figure} \centering \includegraphics[width=\linewidth]{roughShell.eps} @@ -628,17 +608,28 @@ molecules and 256 pentane molecules.} \label{langevin: \begin{figure} \centering \includegraphics[width=\linewidth]{vacf.eps} -\caption[Plots of Velocity Auto-correlation functions]{Velocity -Auto-correlation function of NVE (blue) and Langevin dynamics -(red).} \label{langevin:twoBanana} +\caption[Plots of Velocity Auto-correlation Functions]{Velocity +auto-correlation functions of NVE (explicit solvent) in blue and +Langevin dynamics (implicit solvent) in red.} \label{langevin:vacf} \end{figure} \begin{figure} \centering \includegraphics[width=\linewidth]{uacf.eps} -\caption[Snapshot from Simulation of Two Banana Shaped Molecules and -256 Pentane Molecules]{Snapshot from simulation of two Banana shaped -molecules and 256 pentane molecules.} \label{langevin:twoBanana} +\caption[Auto-correlation functions of the principle axis of the +middle GB particle]{Auto-correlation functions of the principle axis +of the middle GB particle of NVE (blue) and Langevin dynamics +(red).} \label{langevin:uacf} \end{figure} \section{Conclusions} + +We have presented a new Langevin algorithm by incorporating the +hydrodynamics properties of arbitrary shaped molecules into an +advanced symplectic integration scheme. The temperature control +ability of this algorithm was demonstrated by a set of simulations +with different viscosities. It was also shown to have significant +advantage of producing rapid thermal equilibration over +Nos\'{e}-Hoover method. Further studies in systems involving banana +shaped molecules illustrated that the dynamic properties could be +preserved by using this new algorithm as an implicit solvent model.