ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/langevin/langevin.tex
(Generate patch)

Comparing trunk/langevin/langevin.tex (file contents):
Revision 3310 by gezelter, Fri Jan 11 22:08:57 2008 UTC vs.
Revision 3337 by gezelter, Thu Jan 31 14:03:06 2008 UTC

# Line 2 | Line 2
2   %\documentclass[aps,prb,preprint]{revtex4}
3   \documentclass[11pt]{article}
4   \usepackage{endfloat}
5 < \usepackage{amsmath,bm}
5 > \usepackage{amsmath}
6   \usepackage{amssymb}
7   \usepackage{times}
8   \usepackage{mathptmx}
# Line 48 | Line 48 | As alternative to Newtonian dynamics, Langevin dynamic
48   \section{Introduction}
49  
50   %applications of langevin dynamics
51 < As alternative to Newtonian dynamics, Langevin dynamics, which
52 < mimics a simple heat bath with stochastic and dissipative forces,
53 < has been applied in a variety of studies. The stochastic treatment
54 < of the solvent enables us to carry out substantially longer time
55 < simulations. Implicit solvent Langevin dynamics simulations of
56 < met-enkephalin not only outperform explicit solvent simulations for
57 < computational efficiency, but also agrees very well with explicit
58 < solvent simulations for dynamical properties.\cite{Shen2002}
59 < Recently, applying Langevin dynamics with the UNRES model, Liow and
60 < his coworkers suggest that protein folding pathways can be possibly
61 < explored within a reasonable amount of time.\cite{Liwo2005} The
62 < stochastic nature of the Langevin dynamics also enhances the
63 < sampling of the system and increases the probability of crossing
64 < energy barriers.\cite{Banerjee2004, Cui2003} Combining Langevin
65 < dynamics with Kramers's theory, Klimov and Thirumalai identified
66 < free-energy barriers by studying the viscosity dependence of the
67 < protein folding rates.\cite{Klimov1997} In order to account for
68 < solvent induced interactions missing from implicit solvent model,
69 < Kaya incorporated desolvation free energy barrier into implicit
70 < coarse-grained solvent model in protein folding/unfolding studies
71 < and discovered a higher free energy barrier between the native and
72 < denatured states. Because of its stability against noise, Langevin
73 < dynamics is very suitable for studying remagnetization processes in
74 < various systems.\cite{Palacios1998,Berkov2002,Denisov2003} For
51 > Langevin dynamics, which mimics a simple heat bath with stochastic and
52 > dissipative forces, has been applied in a variety of situations as an
53 > alternative to molecular dynamics with explicit solvent molecules.
54 > The stochastic treatment of the solvent allows the use of simulations
55 > with substantially longer time and length scales.  In general, the
56 > dynamic and structural properties obtained from Langevin simulations
57 > agree quite well with similar properties obtained from explicit
58 > solvent simulations.
59 >
60 > Recent examples of the usefulness of Langevin simulations include a
61 > study of met-enkephalin in which Langevin simulations predicted
62 > dynamical properties that were largely in agreement with explicit
63 > solvent simulations.\cite{Shen2002} By applying Langevin dynamics with
64 > the UNRES model, Liow and his coworkers suggest that protein folding
65 > pathways can be explored within a reasonable amount of
66 > time.\cite{Liwo2005}
67 >
68 > The stochastic nature of Langevin dynamics also enhances the sampling
69 > of the system and increases the probability of crossing energy
70 > barriers.\cite{Cui2003,Banerjee2004} Combining Langevin dynamics with
71 > Kramers' theory, Klimov and Thirumalai identified free-energy
72 > barriers by studying the viscosity dependence of the protein folding
73 > rates.\cite{Klimov1997} In order to account for solvent induced
74 > interactions missing from the implicit solvent model, Kaya
75 > incorporated a desolvation free energy barrier into protein
76 > folding/unfolding studies and discovered a higher free energy barrier
77 > between the native and denatured states.\cite{HuseyinKaya07012005}
78 >
79 > Because of its stability against noise, Langevin dynamics has also
80 > proven useful for studying remagnetization processes in various
81 > systems.\cite{Palacios1998,Berkov2002,Denisov2003} [Check: For
82   instance, the oscillation power spectrum of nanoparticles from
83 < Langevin dynamics simulation has the same peak frequencies for
84 < different wave vectors, which recovers the property of magnetic
85 < excitations in small finite structures.\cite{Berkov2005a}
83 > Langevin dynamics has the same peak frequencies for different wave
84 > vectors, which recovers the property of magnetic excitations in small
85 > finite structures.\cite{Berkov2005a}]
86  
87 < %review rigid body dynamics
88 < Rigid bodies are frequently involved in the modeling of different
89 < areas, from engineering, physics, to chemistry. For example,
90 < missiles and vehicle are usually modeled by rigid bodies.  The
91 < movement of the objects in 3D gaming engine or other physics
92 < simulator is governed by the rigid body dynamics. In molecular
93 < simulation, rigid body is used to simplify the model in
94 < protein-protein docking study{\cite{Gray2003}}.
87 > In typical LD simulations, the friction and random forces on
88 > individual atoms are taken from Stokes' law,
89 > \begin{eqnarray}
90 > m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + R(t) \\
91 > \langle R(t) \rangle & = & 0 \\
92 > \langle R(t) R(t') \rangle & = & 2 k_B T \xi m \delta(t - t')
93 > \end{eqnarray}
94 > where $\xi \approx 6 \pi \eta a$.  Here $\eta$ is the viscosity of the
95 > implicit solvent, and $a$ is the hydrodynamic radius of the atom.
96  
97 < It is very important to develop stable and efficient methods to
98 < integrate the equations of motion for orientational degrees of
99 < freedom. Euler angles are the natural choice to describe the
100 < rotational degrees of freedom. However, due to $\frac {1}{sin
101 < \theta}$ singularities, the numerical integration of corresponding
102 < equations of these motion is very inefficient and inaccurate.
103 < Although an alternative integrator using multiple sets of Euler
104 < angles can overcome this difficulty\cite{Barojas1973}, the
97 > The use of rigid substructures,\cite{Chun:2000fj}
98 > coarse-graining,\cite{Ayton01,Golubkov06,Orlandi:2006fk,SunGezelter08}
99 > and ellipsoidal representations of protein side chains~\cite{Fogolari:1996lr}
100 > has made the use of the Stokes-Einstein approximation problematic.  A
101 > rigid substructure moves as a single unit with orientational as well
102 > as translational degrees of freedom.  This requires a more general
103 > treatment of the hydrodynamics than the spherical approximation
104 > provides.  The atoms involved in a rigid or coarse-grained structure
105 > should properly have solvent-mediated interactions with each
106 > other. The theory of interactions {\it between} bodies moving through
107 > a fluid has been developed over the past century and has been applied
108 > to simulations of Brownian
109 > motion.\cite{FIXMAN:1986lr,Ramachandran1996}
110 >
111 > In order to account for the diffusion anisotropy of arbitrarily-shaped
112 > particles, Fernandes and Garc\'{i}a de la Torre improved the original
113 > Brownian dynamics simulation algorithm~\cite{Ermak1978,Allison1991} by
114 > incorporating a generalized $6\times6$ diffusion tensor and
115 > introducing a rotational evolution scheme consisting of three
116 > consecutive rotations.\cite{Fernandes2002} Unfortunately, biases are
117 > introduced into the system due to the arbitrary order of applying the
118 > noncommuting rotation operators.\cite{Beard2003} Based on the
119 > observation the momentum relaxation time is much less than the time
120 > step, one may ignore the inertia in Brownian dynamics.  However, the
121 > assumption of zero average acceleration is not always true for
122 > cooperative motion which is common in proteins. An inertial Brownian
123 > dynamics (IBD) was proposed to address this issue by adding an
124 > inertial correction term.\cite{Beard2000} As a complement to IBD which
125 > has a lower bound in time step because of the inertial relaxation
126 > time, long-time-step inertial dynamics (LTID) can be used to
127 > investigate the inertial behavior of linked polymer segments in a low
128 > friction regime.\cite{Beard2000} LTID can also deal with the
129 > rotational dynamics for nonskew bodies without translation-rotation
130 > coupling by separating the translation and rotation motion and taking
131 > advantage of the analytical solution of hydrodynamics
132 > properties. However, typical nonskew bodies like cylinders and
133 > ellipsoids are inadequate to represent most complex macromolecular
134 > assemblies.  There is therefore a need for incorporating the
135 > hydrodynamics of complex (and potentially skew) rigid bodies in the
136 > library of methods available for performing Langevin simulations.
137 >
138 > \subsection{Rigid Body Dynamics}
139 > Rigid bodies are frequently involved in the modeling of large
140 > collections of particles that move as a single unit.  In molecular
141 > simulations, rigid bodies have been used to simplify protein-protein
142 > docking,\cite{Gray2003} and lipid bilayer
143 > simulations.\cite{SunGezelter08} Many of the water models in common
144 > use are also rigid-body
145 > models,\cite{Jorgensen83,Berendsen81,Berendsen87} although they are
146 > typically evolved using constraints rather than rigid body equations
147 > of motion.
148 >
149 > Euler angles are a natural choice to describe the rotational degrees
150 > of freedom.  However, due to $\frac{1}{\sin \theta}$ singularities, the
151 > numerical integration of corresponding equations of these motion can
152 > become inaccurate (and inefficient).  Although the use of multiple
153 > sets of Euler angles can overcome this problem,\cite{Barojas1973} the
154   computational penalty and the loss of angular momentum conservation
155 < still remain. A singularity-free representation utilizing
156 < quaternions was developed by Evans in 1977.\cite{Evans1977}
157 < Unfortunately, this approach used a nonseparable Hamiltonian
158 < resulting from the quaternion representation, which prevented the
159 < symplectic algorithm from being utilized. Another different approach
103 < is to apply holonomic constraints to the atoms belonging to the
104 < rigid body. Each atom moves independently under the normal forces
105 < deriving from potential energy and constraint forces which are used
106 < to guarantee the rigidness. However, due to their iterative nature,
107 < the SHAKE and Rattle algorithms also converge very slowly when the
108 < number of constraints increases.\cite{Ryckaert1977, Andersen1983}
155 > remain.  A singularity-free representation utilizing quaternions was
156 > developed by Evans in 1977.\cite{Evans1977} The Evans quaternion
157 > approach uses a nonseparable Hamiltonian, and this has prevented
158 > symplectic algorithms from being utilized until very
159 > recently.\cite{Miller2002}
160  
161 < A break-through in geometric literature suggests that, in order to
162 < develop a long-term integration scheme, one should preserve the
163 < symplectic structure of the propagator. By introducing a conjugate
164 < momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's
165 < equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was
166 < proposed to evolve the Hamiltonian system in a constraint manifold
167 < by iteratively satisfying the orthogonality constraint $Q^T Q = 1$.
117 < An alternative method using the quaternion representation was
118 < developed by Omelyan.\cite{Omelyan1998} However, both of these
119 < methods are iterative and inefficient. In this section, we descibe a
120 < symplectic Lie-Poisson integrator for rigid bodies developed by
121 < Dullweber and his coworkers\cite{Dullweber1997} in depth.
161 > Another approach is the application of holonomic constraints to the
162 > atoms belonging to the rigid body.  Each atom moves independently
163 > under the normal forces deriving from potential energy and constraints
164 > are used to guarantee rigidity. However, due to their iterative
165 > nature, the SHAKE and RATTLE algorithms converge very slowly when the
166 > number of constraints (and the number of particles that belong to the
167 > rigid body) increases.\cite{Ryckaert1977,Andersen1983}
168  
169 < %review langevin/browninan dynamics for arbitrarily shaped rigid body
170 < Combining Langevin or Brownian dynamics with rigid body dynamics,
171 < one can study slow processes in biomolecular systems. Modeling DNA
172 < as a chain of rigid beads, which are subject to harmonic potentials
173 < as well as excluded volume potentials, Mielke and his coworkers
174 < discovered rapid superhelical stress generations from the stochastic
175 < simulation of twin supercoiling DNA with response to induced
176 < torques.\cite{Mielke2004} Membrane fusion is another key biological
177 < process which controls a variety of physiological functions, such as
178 < release of neurotransmitters \textit{etc}. A typical fusion event
179 < happens on the time scale of a millisecond, which is impractical to
180 < study using atomistic models with newtonian mechanics. With the help
181 < of coarse-grained rigid body model and stochastic dynamics, the
136 < fusion pathways were explored by many
137 < researchers.\cite{Noguchi2001,Noguchi2002,Shillcock2005} Due to the
138 < difficulty of numerical integration of anisotropic rotation, most of
139 < the rigid body models are simply modeled using spheres, cylinders,
140 < ellipsoids or other regular shapes in stochastic simulations. In an
141 < effort to account for the diffusion anisotropy of arbitrary
142 < particles, Fernandes and de la Torre improved the original Brownian
143 < dynamics simulation algorithm\cite{Ermak1978,Allison1991} by
144 < incorporating a generalized $6\times6$ diffusion tensor and
145 < introducing a simple rotation evolution scheme consisting of three
146 < consecutive rotations.\cite{Fernandes2002} Unfortunately, unexpected
147 < errors and biases are introduced into the system due to the
148 < arbitrary order of applying the noncommuting rotation
149 < operators.\cite{Beard2003} Based on the observation the momentum
150 < relaxation time is much less than the time step, one may ignore the
151 < inertia in Brownian dynamics. However, the assumption of zero
152 < average acceleration is not always true for cooperative motion which
153 < is common in protein motion. An inertial Brownian dynamics (IBD) was
154 < proposed to address this issue by adding an inertial correction
155 < term.\cite{Beard2000} As a complement to IBD which has a lower bound
156 < in time step because of the inertial relaxation time, long-time-step
157 < inertial dynamics (LTID) can be used to investigate the inertial
158 < behavior of the polymer segments in low friction
159 < regime.\cite{Beard2000} LTID can also deal with the rotational
160 < dynamics for nonskew bodies without translation-rotation coupling by
161 < separating the translation and rotation motion and taking advantage
162 < of the analytical solution of hydrodynamics properties. However,
163 < typical nonskew bodies like cylinders and ellipsoids are inadequate
164 < to represent most complex macromolecule assemblies. These intricate
165 < molecules have been represented by a set of beads and their
166 < hydrodynamic properties can be calculated using variants on the
167 < standard hydrodynamic interaction tensors.
169 > In order to develop a stable and efficient integration scheme that
170 > preserves most constants of the motion, symplectic propagators are
171 > necessary.  By introducing a conjugate momentum to the rotation matrix
172 > $Q$ and re-formulating Hamilton's equations, a symplectic
173 > orientational integrator, RSHAKE,\cite{Kol1997} was proposed to evolve
174 > rigid bodies on a constraint manifold by iteratively satisfying the
175 > orthogonality constraint $Q^T Q = 1$.  An alternative method using the
176 > quaternion representation was developed by Omelyan.\cite{Omelyan1998}
177 > However, both of these methods are iterative and suffer from some
178 > related inefficiencies. A symplectic Lie-Poisson integrator for rigid
179 > bodies developed by Dullweber {\it et al.}\cite{Dullweber1997} removes
180 > most of the limitations mentioned above and is therefore the basis for
181 > our Langevin integrator.
182  
183   The goal of the present work is to develop a Langevin dynamics
184   algorithm for arbitrary-shaped rigid particles by integrating the
185 < accurate estimation of friction tensor from hydrodynamics theory
186 < into the sophisticated rigid body dynamics algorithms.
185 > accurate estimation of friction tensor from hydrodynamics theory into
186 > a symplectic rigid body dynamics propagator.  In the sections below,
187 > we review some of the theory of hydrodynamic tensors developed
188 > primarily for Brownian simulations of multi-particle systems, we then
189 > present our integration method for a set of generalized Langevin
190 > equations of motion, and we compare the behavior of the new Langevin
191 > integrator to dynamical quantities obtained via explicit solvent
192 > molecular dynamics.
193  
194 < \subsection{\label{introSection:frictionTensor}Friction Tensor}
195 < Theoretically, the friction kernel can be determined using the
194 > \subsection{\label{introSection:frictionTensor}The Friction Tensor}
195 > Theoretically, a complete friction kernel can be determined using the
196   velocity autocorrelation function. However, this approach becomes
197 < impractical when the system becomes more and more complicated.
198 < Instead, various approaches based on hydrodynamics have been
199 < developed to calculate the friction coefficients. In general, the
200 < friction tensor $\Xi$ is a $6\times 6$ matrix given by
201 < \[
202 < \Xi  = \left( {\begin{array}{*{20}c}
203 <   {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
204 <   {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
205 < \end{array}} \right).
206 < \]
207 < Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are $3 \times 3$
208 < translational friction tensor and rotational resistance (friction)
209 < tensor respectively, while ${\Xi^{tr} }$ is translation-rotation
210 < coupling tensor and $ {\Xi^{rt} }$ is rotation-translation coupling
211 < tensor. When a particle moves in a fluid, it may experience friction
212 < force or torque along the opposite direction of the velocity or
213 < angular velocity,
214 < \[
197 > impractical when the solute becomes complex. Instead, various
198 > approaches based on hydrodynamics have been developed to calculate the
199 > friction coefficients. In general, the friction tensor $\Xi$ is a
200 > $6\times 6$ matrix given by
201 > \begin{equation}
202 > \Xi  = \left( \begin{array}{*{20}c}
203 >   \Xi^{tt} & \Xi^{rt}  \\
204 >   \Xi^{tr} & \Xi^{rr}  \\
205 > \end{array} \right).
206 > \end{equation}
207 > Here, $\Xi^{tt}$ and $\Xi^{rr}$ are $3 \times 3$ translational and
208 > rotational resistance (friction) tensors respectively, while
209 > $\Xi^{tr}$ is translation-rotation coupling tensor and $\Xi^{rt}$ is
210 > rotation-translation coupling tensor. When a particle moves in a
211 > fluid, it may experience friction force ($\mathbf{F}_f$) and torque
212 > ($\mathbf{\tau}_f$) in opposition to the directions of the velocity
213 > ($\mathbf{v}$) and body-fixed angular velocity ($\mathbf{\omega}$),
214 > \begin{equation}
215   \left( \begin{array}{l}
216 < F_R  \\
217 < \tau _R  \\
218 < \end{array} \right) =  - \left( {\begin{array}{*{20}c}
219 <   {\Xi ^{tt} } & {\Xi ^{rt} }  \\
220 <   {\Xi ^{tr} } & {\Xi ^{rr} }  \\
221 < \end{array}} \right)\left( \begin{array}{l}
222 < v \\
223 < w \\
224 < \end{array} \right)
225 < \]
206 < where $F_r$ is the friction force and $\tau _R$ is the friction
207 < torque.
216 > \mathbf{F}_f  \\
217 > \mathbf{\tau}_f  \\
218 > \end{array} \right) =  - \left( \begin{array}{*{20}c}
219 >   \Xi^{tt} & \Xi^{rt}  \\
220 >   \Xi^{tr} & \Xi^{rr}  \\
221 > \end{array} \right)\left( \begin{array}{l}
222 > \mathbf{v} \\
223 > \mathbf{\omega} \\
224 > \end{array} \right).
225 > \end{equation}
226  
227   \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}}
228 <
229 < For a spherical particle with slip boundary conditions, the
230 < translational and rotational friction constant can be calculated
231 < from Stoke's law,
232 < \[
215 < \Xi ^{tt}  = \left( {\begin{array}{*{20}c}
228 > For a spherical particle under ``stick'' boundary conditions, the
229 > translational and rotational friction tensors can be calculated from
230 > Stokes' law,
231 > \begin{equation}
232 > \Xi^{tt}  = \left( \begin{array}{*{20}c}
233     {6\pi \eta R} & 0 & 0  \\
234     0 & {6\pi \eta R} & 0  \\
235     0 & 0 & {6\pi \eta R}  \\
236 < \end{array}} \right)
237 < \]
236 > \end{array} \right)
237 > \end{equation}
238   and
239 < \[
240 < \Xi ^{rr}  = \left( {\begin{array}{*{20}c}
239 > \begin{equation}
240 > \Xi^{rr}  = \left( \begin{array}{*{20}c}
241     {8\pi \eta R^3 } & 0 & 0  \\
242     0 & {8\pi \eta R^3 } & 0  \\
243     0 & 0 & {8\pi \eta R^3 }  \\
244 < \end{array}} \right)
245 < \]
244 > \end{array} \right)
245 > \end{equation}
246   where $\eta$ is the viscosity of the solvent and $R$ is the
247   hydrodynamic radius.
248  
249   Other non-spherical shapes, such as cylinders and ellipsoids, are
250 < widely used as references for developing new hydrodynamics theory,
250 > widely used as references for developing new hydrodynamics theories,
251   because their properties can be calculated exactly. In 1936, Perrin
252 < extended Stokes's law to general ellipsoids, also called a triaxial
253 < ellipsoid, which is given in Cartesian coordinates
254 < by\cite{Perrin1934, Perrin1936}
255 < \[
256 < \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
257 < }} = 1
258 < \]
259 < where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
260 < due to the complexity of the elliptic integral, only the ellipsoid
261 < with the restriction of two axes being equal, \textit{i.e.}
262 < prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
263 < exactly. Introducing an elliptic integral parameter $S$ for prolate
264 < ellipsoids :
265 < \[
266 < S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
267 < } }}{b},
268 < \]
269 < and oblate ellipsoids:
253 < \[
254 < S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
255 < }}{a},
256 < \]
257 < one can write down the translational and rotational resistance
258 < tensors
252 > extended Stokes' law to general ellipsoids which are given in
253 > Cartesian coordinates by~\cite{Perrin1934,Perrin1936}
254 > \begin{equation}
255 > \frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1.
256 > \end{equation}
257 > Here, the semi-axes are of lengths $a$, $b$, and $c$. Due to the
258 > complexity of the elliptic integral, only uniaxial ellipsoids, either
259 > prolate ($a \ge b = c$) or oblate ($a < b = c$), can be solved
260 > exactly. Introducing an elliptic integral parameter $S$ for prolate,
261 > \begin{equation}
262 > S = \frac{2}{\sqrt{a^2  - b^2}} \ln \frac{a + \sqrt{a^2  - b^2}}{b},
263 > \end{equation}
264 > and oblate,
265 > \begin{equation}
266 > S = \frac{2}{\sqrt {b^2  - a^2 }} \arctan \frac{\sqrt {b^2  - a^2}}{a},
267 > \end{equation}
268 > ellipsoids, one can write down the translational and rotational
269 > resistance tensors:
270   \begin{eqnarray*}
271 < \Xi _a^{tt}  & = & 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}}. \\
272 < \Xi _b^{tt}  & = & \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S +
262 < 2a}},
271 > \Xi_a^{tt}  & = & 16\pi \eta \frac{a^2  - b^2}{(2a^2  - b^2 )S - 2a}. \\
272 > \Xi_b^{tt} =  \Xi_c^{tt} & = & 32\pi \eta \frac{a^2  - b^2 }{(2a^2 - 3b^2 )S + 2a},
273   \end{eqnarray*}
274 < and
274 > for oblate, and
275   \begin{eqnarray*}
276 < \Xi _a^{rr} & = & \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}}, \\
277 < \Xi _b^{rr} & = & \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}}.
276 > \Xi_a^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^2  - b^2 )b^2}{2a - b^2 S}, \\
277 > \Xi_b^{rr} = \Xi_c^{rr} & = & \frac{32\pi}{3} \eta \frac{(a^4  - b^4)}{(2a^2  - b^2 )S - 2a}
278   \end{eqnarray*}
279 + for prolate ellipsoids. For both spherical and ellipsoidal particles,
280 + the translation-rotation and rotation-translation coupling tensors are
281 + zero.
282  
283   \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}}
271
284   Unlike spherical and other simply shaped molecules, there is no
285   analytical solution for the friction tensor for arbitrarily shaped
286   rigid molecules. The ellipsoid of revolution model and general
287   triaxial ellipsoid model have been used to approximate the
288 < hydrodynamic properties of rigid bodies. However, since the mapping
289 < from all possible ellipsoidal spaces, $r$-space, to all possible
290 < combination of rotational diffusion coefficients, $D$-space, is not
291 < unique\cite{Wegener1979} as well as the intrinsic coupling between
292 < translational and rotational motion of rigid bodies, general
293 < ellipsoids are not always suitable for modeling arbitrarily shaped
294 < rigid molecules. A number of studies have been devoted to
288 > hydrodynamic properties of rigid bodies. However, the mapping from all
289 > possible ellipsoidal spaces, $r$-space, to all possible combination of
290 > rotational diffusion coefficients, $D$-space, is not
291 > unique.\cite{Wegener1979} Additionally, because there is intrinsic
292 > coupling between translational and rotational motion of rigid bodies,
293 > general ellipsoids are not always suitable for modeling arbitrarily
294 > shaped rigid molecules.  A number of studies have been devoted to
295   determining the friction tensor for irregularly shaped rigid bodies
296 < using more advanced methods where the molecule of interest was
297 < modeled by a combinations of spheres\cite{Carrasco1999} and the
298 < hydrodynamics properties of the molecule can be calculated using the
299 < hydrodynamic interaction tensor. Let us consider a rigid assembly of
300 < $N$ beads immersed in a continuous medium. Due to hydrodynamic
301 < interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different
302 < than its unperturbed velocity $v_i$,
303 < \[
296 > using more advanced methods where the molecule of interest was modeled
297 > by a combinations of spheres\cite{Carrasco1999} and the hydrodynamics
298 > properties of the molecule can be calculated using the hydrodynamic
299 > interaction tensor.
300 >
301 > Consider a rigid assembly of $N$ beads immersed in a continuous
302 > medium. Due to hydrodynamic interaction, the ``net'' velocity of $i$th
303 > bead, $v'_i$ is different than its unperturbed velocity $v_i$,
304 > \begin{equation}
305   v'_i  = v_i  - \sum\limits_{j \ne i} {T_{ij} F_j }
306 < \]
307 < where $F_i$ is the frictional force, and $T_{ij}$ is the
308 < hydrodynamic interaction tensor. The friction force of $i$th bead is
309 < proportional to its ``net'' velocity
306 > \end{equation}
307 > where $F_i$ is the frictional force, and $T_{ij}$ is the hydrodynamic
308 > interaction tensor. The frictional force on the $i^\mathrm{th}$ bead
309 > is proportional to its ``net'' velocity
310   \begin{equation}
311   F_i  = \zeta _i v_i  - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
312   \label{introEquation:tensorExpression}
# Line 330 | Line 343 | B = \left( {\begin{array}{*{20}c}
343   construct a $3N \times 3N$ matrix consisting of $N \times N$
344   $B_{ij}$ blocks
345   \begin{equation}
346 < B = \left( {\begin{array}{*{20}c}
347 <   {B_{11} } &  \ldots  & {B_{1N} }  \\
346 > B = \left( \begin{array}{*{20}c}
347 >   B_{11} &  \ldots  & B_{1N}   \\
348      \vdots  &  \ddots  &  \vdots   \\
349 <   {B_{N1} } &  \cdots  & {B_{NN} }  \\
350 < \end{array}} \right),
349 >   B_{N1} &  \cdots  & B_{NN} \\
350 > \end{array} \right),
351   \end{equation}
352   where $B_{ij}$ is given by
353 < \[
353 > \begin{equation}
354   B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
355   )T_{ij}
356 < \]
356 > \end{equation}
357   where $\delta _{ij}$ is the Kronecker delta function. Inverting the
358   $B$ matrix, we obtain
359   \[
# Line 425 | Line 438 | locate the position of center of resistance,
438   x_{OR}  \\
439   y_{OR}  \\
440   z_{OR}  \\
441 < \end{array} \right) & = &\left( {\begin{array}{*{20}c}
441 > \end{array} \right) & = &\left( \begin{array}{*{20}c}
442     {(\Xi _O^{rr} )_{yy}  + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} }  \\
443     { - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz}  + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} }  \\
444     { - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx}  + (\Xi _O^{rr} )_{yy} }  \\
445 < \end{array}} \right)^{ - 1}  \\
445 > \end{array} \right)^{ - 1}  \\
446    & & \left( \begin{array}{l}
447   (\Xi _O^{tr} )_{yz}  - (\Xi _O^{tr} )_{zy}  \\
448   (\Xi _O^{tr} )_{zx}  - (\Xi _O^{tr} )_{xz}  \\
# Line 441 | Line 454 | Consider the Langevin equations of motion in generaliz
454  
455  
456   \section{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}
457 +
458   Consider the Langevin equations of motion in generalized coordinates
459   \begin{equation}
460 < M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)}  + F_{r,i} (t)
460 > \mathbf{M} \dot{\mathbf{V}}(t) = \mathbf{F}_{s}(t) +
461 > \mathbf{F}_{f}(t)  + \mathbf{R}(t)
462   \label{LDGeneralizedForm}
463   \end{equation}
464 < where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
465 < and moment of inertial) matrix and $V_i$ is a generalized velocity,
466 < $V_i = V_i(v_i,\omega _i)$. The right side of
467 < Eq.~\ref{LDGeneralizedForm} consists of three generalized forces in
468 < lab-fixed frame, systematic force $F_{s,i}$, dissipative force
469 < $F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
470 < system in Newtownian mechanics typically refers to lab-fixed frame,
471 < it is also convenient to handle the rotation of rigid body in
472 < body-fixed frame. Thus the friction and random forces are calculated
473 < in body-fixed frame and converted back to lab-fixed frame by:
474 < \[
464 > where $\mathbf{M}$ is a $6 \times 6$ diagonal mass matrix (which
465 > includes the mass of the rigid body as well as the moments of inertia
466 > in the body-fixed frame) and $\mathbf{V}$ is a generalized velocity,
467 > $\mathbf{V} =
468 > \left\{\mathbf{v},\mathbf{\omega}\right\}$. The right side of
469 > Eq.~\ref{LDGeneralizedForm} consists of three generalized forces: a
470 > system force $\mathbf{F}_{s}$, a frictional or dissipative force
471 > $\mathbf{F}_{f}$ and stochastic force $\mathbf{R}$. While the
472 > evolution of the system in Newtownian mechanics is typically done in the
473 > lab-fixed frame, it is convenient to handle the rotation of rigid
474 > bodies in the body-fixed frame. Thus the friction and random forces are
475 > calculated in body-fixed frame and converted back to lab-fixed frame
476 > using the rigid body's rotation matrix ($Q$):
477 > \begin{equation}
478   \begin{array}{l}
479 < F_{f,i}^l (t) = Q^T F_{f,i}^b (t), \\
480 < F_{r,i}^l (t) = Q^T F_{r,i}^b (t). \\
479 > \mathbf{F}_{f}(t) = Q^{T} \mathbf{F}_{f}^b (t), \\
480 > \mathbf{R}(t) = Q^{T} \mathbf{R}^b (t). \\
481   \end{array}
482 < \]
483 < Here, the body-fixed friction force $F_{r,i}^b$ is proportional to
484 < the body-fixed velocity at center of resistance $v_{R,i}^b$ and
485 < angular velocity $\omega _i$
482 > \end{equation}
483 > Here, the body-fixed friction force $\mathbf{F}_{f,i}^b$ is proportional to
484 > the body-fixed velocity at the center of resistance $\mathbf{v}_{R,i}^b$ and
485 > angular velocity $\mathbf{\omega}_i$
486   \begin{equation}
487 < F_{r,i}^b (t) = \left( \begin{array}{l}
488 < f_{r,i}^b (t) \\
489 < \tau _{r,i}^b (t) \\
490 < \end{array} \right) =  - \left( {\begin{array}{*{20}c}
491 <   {\Xi _{R,t} } & {\Xi _{R,c}^T }  \\
492 <   {\Xi _{R,c} } & {\Xi _{R,r} }  \\
493 < \end{array}} \right)\left( \begin{array}{l}
494 < v_{R,i}^b (t) \\
495 < \omega _i (t) \\
487 > \mathbf{F}_{f}^b (t) = \left( \begin{array}{l}
488 > \mathbf{f}_{f}^b (t) \\
489 > \mathbf{\tau}_{f}^b (t) \\
490 > \end{array} \right) =  - \left( \begin{array}{*{20}c}
491 >   \Xi_{R,t} & \Xi_{R,c}^T  \\
492 >   \Xi_{R,c} & \Xi_{R,r}    \\
493 > \end{array} \right)\left( \begin{array}{l}
494 > \mathbf{v}_{R}^b (t) \\
495 > \mathbf{\omega} (t) \\
496   \end{array} \right),
497   \end{equation}
498 < while the random force $F_{r,i}^l$ is a Gaussian stochastic variable
498 > while the random force $\mathbf{R}^l$ is a Gaussian stochastic variable
499   with zero mean and variance
500   \begin{equation}
501 < \left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
502 < \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
503 < 2k_B T\Xi _R \delta (t - t'). \label{randomForce}
501 > \left\langle {\mathbf{R}^l (t) (\mathbf{R}^l (t'))^T } \right\rangle  =
502 > \left\langle {\mathbf{R}^b (t) (\mathbf{R}^b (t'))^T } \right\rangle  =
503 > 2 k_B T \Xi_R \delta(t - t'). \label{randomForce}
504   \end{equation}
505 < The equation of motion for $v_i$ can be written as
505 > Once the $6\times6$ resistance tensor at the center of resistance
506 > ($\Xi_R$) is known, obtaining a stochastic vector that has the
507 > properties in Eq. (\ref{eq:randomForce}) can be done efficiently by
508 > carrying out a one-time Cholesky decomposition to obtain the square
509 > root matrix of $\Xi_R$.\cite{SchlickBook} Each time a random force
510 > vector is needed, a gaussian random vector is generated and then the
511 > square root matrix is multiplied onto this vector.
512 >
513 > The equation of motion for $\mathbf{v}$ can be written as
514   \begin{equation}
515 < m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
516 < f_{r,i}^l (t)
515 > m \dot{\mathbf{v}} (t) =  \mathbf{f}_{s} (t) + \mathbf{f}_{f}^l (t) +
516 > \mathbf{R}^l (t)
517   \end{equation}
518   Since the frictional force is applied at the center of resistance
519   which generally does not coincide with the center of mass, an extra
520   torque is exerted at the center of mass. Thus, the net body-fixed
521 < frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is
521 > frictional torque at the center of mass, $\tau_{f}^b (t)$, is
522   given by
523   \begin{equation}
524 < \tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b
524 > \tau_{f}^b \leftarrow \tau_{f}^b + \mathbf{r}_{MR} \times \mathbf{f}_{r}^b
525   \end{equation}
526   where $r_{MR}$ is the vector from the center of mass to the center
527   of the resistance. Instead of integrating the angular velocity in
528   lab-fixed frame, we consider the equation of angular momentum in
529   body-fixed frame
530   \begin{equation}
531 < \dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t)
506 < + \tau _{r,i}^b(t)
531 > \dot j(t) = \tau_{s} (t) + \tau_{f}^b (t) + \mathbf{R}^b(t)
532   \end{equation}
533 < Embedding the friction terms into force and torque, one can
534 < integrate the langevin equations of motion for rigid body of
535 < arbitrary shape in a velocity-Verlet style 2-part algorithm, where
511 < $h= \delta t$:
533 > Embedding the friction terms into force and torque, one can integrate
534 > the Langevin equations of motion for rigid body of arbitrary shape in
535 > a velocity-Verlet style 2-part algorithm, where $h= \delta t$:
536  
537   {\tt moveA:}
538   \begin{align*}
# Line 1061 | Line 1085 | used recently as models for lipid molecules.\cite{SunG
1085  
1086   \subsection{Composite sphero-ellipsoids}
1087   Spherical heads perched on the ends of Gay-Berne ellipsoids have been
1088 < used recently as models for lipid molecules.\cite{SunGezelter08,Ayton01}
1089 <
1088 > used recently as models for lipid
1089 > molecules.\cite{SunGezelter08,Ayton01}
1090   MORE DETAILS
1091  
1092 + A reference system composed of a single lipid rigid body embedded in a
1093 + sea of 1929 solvent particles was created and run under standard
1094 + (microcanonical) molecular dynamics.  The resulting viscosity of this
1095 + mixture was 0.349 centipoise (as estimated using
1096 + Eq. (\ref{eq:shear})).  To calculate the hydrodynamic properties of
1097 + the lipid rigid body model, we created a rough shell (see
1098 + Fig.~\ref{fig:roughShell}), in which the lipid is represented as a
1099 + ``shell'' made of 3550 identical beads (0.25 \AA\ in diameter)
1100 + distributed on the surface.  Applying the procedure described in
1101 + Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
1102 + identified the center of resistance, ${\bf r} = $(0 \AA, 0 \AA, 1.46
1103 + \AA).
1104  
1105 +
1106   \subsection{Summary}
1107   According to our simulations, the langevin dynamics is a reliable
1108   theory to apply to replace the explicit solvents, especially for the
1109   translation properties. For large molecules, the rotation properties
1110   are also mimiced reasonablly well.
1111  
1112 + \begin{figure}
1113 + \centering
1114 + \includegraphics[width=\linewidth]{graph}
1115 + \caption[Mean squared displacements and orientational
1116 + correlation functions for each of the model rigid bodies.]{The
1117 + mean-squared displacements ($\langle r^2(t) \rangle$) and
1118 + orientational correlation functions ($C_2(t)$) for each of the model
1119 + rigid bodies studied.  The circles are the results for microcanonical
1120 + simulations with explicit solvent molecules, while the other data sets
1121 + are results for Langevin dynamics using the different hydrodynamic
1122 + tensor approximations.  The Perrin model for the ellipsoids is
1123 + considered the ``exact'' hydrodynamic behavior (this can also be said
1124 + for the translational motion of the dumbbell operating under the bead
1125 + model). In most cases, the various hydrodynamics models reproduce
1126 + each other quantitatively.}
1127 + \label{fig:results}
1128 + \end{figure}
1129 +
1130   \begin{table*}
1131   \begin{minipage}{\linewidth}
1132   \begin{center}
# Line 1090 | Line 1145 | sphere    & 0.261  & ?    & & 2.59 & exact       & 2.5
1145   \cline{2-3} \cline{5-7}
1146   model & $\eta$ (centipoise)  & D & & Analytical & method & Hydrodynamics & simulation \\
1147   \hline
1148 < sphere    & 0.261  & ?    & & 2.59 & exact       & 2.59 & 2.56 \\
1148 > sphere    & 0.279  & 3.06 & & 2.42 & exact       & 2.42 & 2.33 \\
1149   ellipsoid & 0.255  & 2.44 & & 2.34 & exact       & 2.34 & 2.37 \\
1150            & 0.255  & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\
1151 < dumbbell  & 0.322  & ?    & & 1.57 & bead model  & 1.57 & 1.57 \\
1152 <          & 0.322  & ?    & & 1.57 & rough shell & ?    & ?    \\
1151 > dumbbell  & 0.308  & 2.06 & & 1.64 & bead model  & 1.65 & 1.62 \\
1152 >          & 0.308  & 2.06 & & 1.64 & rough shell & 1.59 & 1.62 \\
1153   banana    & 0.298  & 1.53 & &      & rough shell & 1.56 & 1.55 \\
1154   lipid     & 0.349  & 0.96 & &      & rough shell & 1.33 & 1.32 \\
1155   \end{tabular}
# Line 1119 | Line 1174 | sphere    & 0.261  &      & & 9.06 & exact       & 9.0
1174   \cline{2-3} \cline{5-7}
1175   model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic  & simulation \\
1176   \hline
1177 < sphere    & 0.261  &      & & 9.06 & exact       & 9.06 & 9.11 \\
1177 > sphere    & 0.279  &      & & 9.69 & exact       & 9.69 & 9.64 \\
1178   ellipsoid & 0.255  & 46.7 & & 22.0 & exact       & 22.0 & 22.2 \\
1179            & 0.255  & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\
1180 < dumbbell  & 0.322  & 14.0 & &      & bead model  & 52.3 & 52.8 \\
1181 <          & 0.322  & 14.0 & &      & rough shell & ?    & ?    \\
1180 > dumbbell  & 0.308  & 14.1 & &      & bead model  & 50.0 & 50.1 \\
1181 >          & 0.308  & 14.1 & &      & rough shell & 41.5 & 41.3 \\
1182   banana    & 0.298  & 63.8 & &      & rough shell & 70.9 & 70.9 \\
1183   lipid     & 0.349  & 78.0 & &      & rough shell & 76.9 & 77.9 \\
1184   \hline

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines