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

Comparing trunk/tengDissertation/Langevin.tex (file contents):
Revision 2851 by tim, Sun Jun 11 02:06:01 2006 UTC vs.
Revision 2909 by tim, Thu Jun 29 23:00:35 2006 UTC

# Line 1 | Line 1
1  
2 < \chapter{\label{chapt:methodology}Langevin Dynamics for Rigid Bodies of Arbitrary Shape}
2 > \chapter{\label{chapt:langevin}LANGEVIN DYNAMICS FOR RIGID BODIES OF ARBITRARY SHAPE}
3  
4   \section{Introduction}
5  
6   %applications of langevin dynamics
7 < As an excellent alternative to newtonian dynamics, Langevin
8 < dynamics, which mimics a simple heat bath with stochastic and
9 < dissipative forces, has been applied in a variety of studies. The
10 < stochastic treatment of the solvent enables us to carry out
11 < substantially longer time simulation. Implicit solvent Langevin
12 < dynamics simulation of met-enkephalin not only outperforms explicit
13 < solvent simulation on computation efficiency, but also agrees very
14 < well with explicit solvent simulation on dynamics
15 < properties\cite{Shen2002}. Recently, applying Langevin dynamics with
16 < UNRES model, Liow and his coworkers suggest that protein folding
17 < pathways can be possibly exploited within a reasonable amount of
18 < time\cite{Liwo2005}. The stochastic nature of the Langevin dynamics
19 < also enhances the sampling of the system and increases the
20 < probability of crossing energy barrier\cite{Banerjee2004, Cui2003}.
21 < Combining Langevin dynamics with Kramers's theory, Klimov and
22 < Thirumalai identified the free-energy barrier by studying the
23 < viscosity dependence of the protein folding rates\cite{Klimov1997}.
24 < In order to account for solvent induced interactions missing from
25 < implicit solvent model, Kaya incorporated desolvation free energy
26 < barrier into implicit coarse-grained solvent model in protein
27 < folding/unfolding study and discovered a higher free energy barrier
28 < between the native and denatured states. Because of its stability
29 < against noise, Langevin dynamics is very suitable for studying
30 < remagnetization processes in various
31 < systems\cite{Garcia-Palacios1998,Berkov2002,Denisov2003}. For
7 > As alternative to Newtonian dynamics, Langevin dynamics, which
8 > mimics a simple heat bath with stochastic and dissipative forces,
9 > has been applied in a variety of studies. The stochastic treatment
10 > of the solvent enables us to carry out substantially longer time
11 > simulations. Implicit solvent Langevin dynamics simulations of
12 > met-enkephalin not only outperform explicit solvent simulations for
13 > computational efficiency, but also agrees very well with explicit
14 > solvent simulations for dynamical properties\cite{Shen2002}.
15 > Recently, applying Langevin dynamics with the UNRES model, Liow and
16 > his coworkers suggest that protein folding pathways can be possibly
17 > explored within a reasonable amount of time\cite{Liwo2005}. The
18 > stochastic nature of the Langevin dynamics also enhances the
19 > sampling of the system and increases the probability of crossing
20 > energy barriers\cite{Banerjee2004, Cui2003}. Combining Langevin
21 > dynamics with Kramers's theory, Klimov and Thirumalai identified
22 > free-energy barriers by studying the viscosity dependence of the
23 > protein folding rates\cite{Klimov1997}. In order to account for
24 > solvent induced interactions missing from implicit solvent model,
25 > Kaya incorporated desolvation free energy barrier into implicit
26 > coarse-grained solvent model in protein folding/unfolding studies
27 > and discovered a higher free energy barrier between the native and
28 > denatured states. Because of its stability against noise, Langevin
29 > dynamics is very suitable for studying remagnetization processes in
30 > various systems\cite{Palacios1998,Berkov2002,Denisov2003}. For
31   instance, the oscillation power spectrum of nanoparticles from
32   Langevin dynamics simulation has the same peak frequencies for
33 < different wave vectors,which recovers the property of magnetic
34 < excitations in small finite structures\cite{Berkov2005a}. In an
36 < attempt to reduce the computational cost of simulation, multiple
37 < time stepping (MTS) methods have been introduced and have been of
38 < great interest to macromolecule and protein
39 < community\cite{Tuckerman1992}. Relying on the observation that
40 < forces between distant atoms generally demonstrate slower
41 < fluctuations than forces between close atoms, MTS method are
42 < generally implemented by evaluating the slowly fluctuating forces
43 < less frequently than the fast ones. Unfortunately, nonlinear
44 < instability resulting from increasing timestep in MTS simulation
45 < have became a critical obstruction preventing the long time
46 < simulation. Due to the coupling to the heat bath, Langevin dynamics
47 < has been shown to be able to damp out the resonance artifact more
48 < efficiently\cite{Sandu1999}.
33 > different wave vectors, which recovers the property of magnetic
34 > excitations in small finite structures\cite{Berkov2005a}.
35  
50 – %review rigid body dynamics
51 – Rigid bodies are frequently involved in the modeling of different
52 – areas, from engineering, physics, to chemistry. For example,
53 – missiles and vehicle are usually modeled by rigid bodies.  The
54 – movement of the objects in 3D gaming engine or other physics
55 – simulator is governed by the rigid body dynamics. In molecular
56 – simulation, rigid body is used to simplify the model in
57 – protein-protein docking study{\cite{Gray2003}}.
58 –
59 – It is very important to develop stable and efficient methods to
60 – integrate the equations of motion of orientational degrees of
61 – freedom. Euler angles are the nature choice to describe the
62 – rotational degrees of freedom. However, due to its singularity, the
63 – numerical integration of corresponding equations of motion is very
64 – inefficient and inaccurate. Although an alternative integrator using
65 – different sets of Euler angles can overcome this difficulty\cite{},
66 – the computational penalty and the lost of angular momentum
67 – conservation still remain. In 1977, a singularity free
68 – representation utilizing quaternions was developed by
69 – Evans\cite{Evans1977}. Unfortunately, this approach suffer from the
70 – nonseparable Hamiltonian resulted from quaternion representation,
71 – which prevents the symplectic algorithm to be utilized. Another
72 – different approach is to apply holonomic constraints to the atoms
73 – belonging to the rigid body\cite{}. Each atom moves independently
74 – under the normal forces deriving from potential energy and
75 – constraint forces which are used to guarantee the rigidness.
76 – However, due to their iterative nature, SHAKE and Rattle algorithm
77 – converge very slowly when the number of constraint increases.
78 –
79 – The break through in geometric literature suggests that, in order to
80 – develop a long-term integration scheme, one should preserve the
81 – geometric structure of the flow. Matubayasi and Nakahara developed a
82 – time-reversible integrator for rigid bodies in quaternion
83 – representation. Although it is not symplectic, this integrator still
84 – demonstrates a better long-time energy conservation than traditional
85 – methods because of the time-reversible nature. Extending
86 – Trotter-Suzuki to general system with a flat phase space, Miller and
87 – his colleagues devised an novel symplectic, time-reversible and
88 – volume-preserving integrator in quaternion representation. However,
89 – all of the integrators in quaternion representation suffer from the
90 – computational penalty of constructing a rotation matrix from
91 – quaternions to evolve coordinates and velocities at every time step.
92 – An alternative integration scheme utilizing rotation matrix directly
93 – is RSHAKE , in which a conjugate momentum to rotation matrix is
94 – introduced to re-formulate the Hamiltonian's equation and the
95 – Hamiltonian is evolved in a constraint manifold by iteratively
96 – satisfying the orthogonality constraint. However, RSHAKE is
97 – inefficient because of the iterative procedure. An extremely
98 – efficient integration scheme in rotation matrix representation,
99 – which also preserves the same structural properties of the
100 – Hamiltonian flow as Miller's integrator, is proposed by Dullweber,
101 – Leimkuhler and McLachlan (DLM).
102 –
36   %review langevin/browninan dynamics for arbitrarily shaped rigid body
37   Combining Langevin or Brownian dynamics with rigid body dynamics,
38 < one can study the slow processes in biomolecular systems. Modeling
39 < the DNA as a chain of rigid spheres beads, which subject to harmonic
40 < potentials as well as excluded volume potentials, Mielke and his
41 < coworkers discover rapid superhelical stress generations from the
42 < stochastic simulation of twin supercoiling DNA with response to
43 < induced torques\cite{Mielke2004}. Membrane fusion is another key
44 < biological process which controls a variety of physiological
45 < functions, such as release of neurotransmitters \textit{etc}. A
46 < typical fusion event happens on the time scale of millisecond, which
47 < is impracticable to study using all atomistic model with newtonian
48 < mechanics. With the help of coarse-grained rigid body model and
49 < stochastic dynamics, the fusion pathways were exploited by many
38 > one can study slow processes in biomolecular systems. Modeling DNA
39 > as a chain of rigid beads, which are subject to harmonic potentials
40 > as well as excluded volume potentials, Mielke and his coworkers
41 > discovered rapid superhelical stress generations from the stochastic
42 > simulation of twin supercoiling DNA with response to induced
43 > torques\cite{Mielke2004}. Membrane fusion is another key biological
44 > process which controls a variety of physiological functions, such as
45 > release of neurotransmitters \textit{etc}. A typical fusion event
46 > happens on the time scale of millisecond, which is impractical to
47 > study using atomistic models with newtonian mechanics. With the help
48 > of coarse-grained rigid body model and stochastic dynamics, the
49 > fusion pathways were explored by many
50   researchers\cite{Noguchi2001,Noguchi2002,Shillcock2005}. Due to the
51 < difficulty of numerical integration of anisotropy rotation, most of
52 < the rigid body models are simply modeled by sphere, cylinder,
53 < ellipsoid or other regular shapes in stochastic simulations. In an
54 < effort to account for the diffusion anisotropy of the arbitrary
51 > difficulty of numerical integration of anisotropic rotation, most of
52 > the rigid body models are simply modeled using spheres, cylinders,
53 > ellipsoids or other regular shapes in stochastic simulations. In an
54 > effort to account for the diffusion anisotropy of arbitrary
55   particles, Fernandes and de la Torre improved the original Brownian
56   dynamics simulation algorithm\cite{Ermak1978,Allison1991} by
57   incorporating a generalized $6\times6$ diffusion tensor and
58   introducing a simple rotation evolution scheme consisting of three
59   consecutive rotations\cite{Fernandes2002}. Unfortunately, unexpected
60 < error and bias are introduced into the system due to the arbitrary
61 < order of applying the noncommuting rotation
60 > errors and biases are introduced into the system due to the
61 > arbitrary order of applying the noncommuting rotation
62   operators\cite{Beard2003}. Based on the observation the momentum
63   relaxation time is much less than the time step, one may ignore the
64 < inertia in Brownian dynamics. However, assumption of the zero
64 > inertia in Brownian dynamics. However, the assumption of zero
65   average acceleration is not always true for cooperative motion which
66   is common in protein motion. An inertial Brownian dynamics (IBD) was
67   proposed to address this issue by adding an inertial correction
68 < term\cite{Beard2001}. As a complement to IBD which has a lower bound
68 > term\cite{Beard2000}. As a complement to IBD which has a lower bound
69   in time step because of the inertial relaxation time, long-time-step
70   inertial dynamics (LTID) can be used to investigate the inertial
71   behavior of the polymer segments in low friction
72 < regime\cite{Beard2001}. LTID can also deal with the rotational
72 > regime\cite{Beard2000}. LTID can also deal with the rotational
73   dynamics for nonskew bodies without translation-rotation coupling by
74   separating the translation and rotation motion and taking advantage
75   of the analytical solution of hydrodynamics properties. However,
76 < typical nonskew bodies like cylinder and ellipsoid are inadequate to
77 < represent most complex macromolecule assemblies. These intricate
76 > typical nonskew bodies like cylinders and ellipsoids are inadequate
77 > to represent most complex macromolecule assemblies. These intricate
78   molecules have been represented by a set of beads and their
79 < hydrodynamics properties can be calculated using variant
80 < hydrodynamic interaction tensors.
79 > hydrodynamic properties can be calculated using variants on the
80 > standard hydrodynamic interaction tensors.
81  
82   The goal of the present work is to develop a Langevin dynamics
83 < algorithm for arbitrary rigid particles by integrating the accurate
84 < estimation of friction tensor from hydrodynamics theory into the
85 < sophisticated rigid body dynamics.
83 > algorithm for arbitrary-shaped rigid particles by integrating the
84 > accurate estimation of friction tensor from hydrodynamics theory
85 > into the sophisticated rigid body dynamics algorithms.
86  
87 < \section{Method{\label{methodSec}}}
87 > \section{Computational Methods{\label{methodSec}}}
88  
89 < \subsection{\label{introSection:frictionTensor} Friction Tensor}
90 < Theoretically, the friction kernel can be determined using velocity
91 < autocorrelation function. However, this approach become impractical
92 < when the system become more and more complicate. Instead, various
93 < approaches based on hydrodynamics have been developed to calculate
94 < the friction coefficients. The friction effect is isotropic in
95 < Equation, $\zeta$ can be taken as a scalar. In general, friction
163 < tensor $\Xi$ is a $6\times 6$ matrix given by
89 > \subsection{\label{introSection:frictionTensor}Friction Tensor}
90 > Theoretically, the friction kernel can be determined using the
91 > velocity autocorrelation function. However, this approach becomes
92 > impractical when the system becomes more and more complicated.
93 > Instead, various approaches based on hydrodynamics have been
94 > developed to calculate the friction coefficients. In general, the
95 > friction tensor $\Xi$ is a $6\times 6$ matrix given by
96   \[
97   \Xi  = \left( {\begin{array}{*{20}c}
98     {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
99     {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
100   \end{array}} \right).
101   \]
102 < Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are translational friction
103 < tensor and rotational resistance (friction) tensor respectively,
104 < while ${\Xi^{tr} }$ is translation-rotation coupling tensor and $
105 < {\Xi^{rt} }$ is rotation-translation coupling tensor. When a
106 < particle moves in a fluid, it may experience friction force or
107 < torque along the opposite direction of the velocity or angular
108 < velocity,
102 > Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are $3 \times 3$
103 > translational friction tensor and rotational resistance (friction)
104 > tensor respectively, while ${\Xi^{tr} }$ is translation-rotation
105 > coupling tensor and $ {\Xi^{rt} }$ is rotation-translation coupling
106 > tensor. When a particle moves in a fluid, it may experience friction
107 > force or torque along the opposite direction of the velocity or
108 > angular velocity,
109   \[
110   \left( \begin{array}{l}
111   F_R  \\
# Line 187 | Line 119 | toque.
119   \end{array} \right)
120   \]
121   where $F_r$ is the friction force and $\tau _R$ is the friction
122 < toque.
122 > torque.
123  
124 < \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}}
124 > \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}}
125  
126 < For a spherical particle, the translational and rotational friction
127 < constant can be calculated from Stoke's law,
126 > For a spherical particle with slip boundary conditions, the
127 > translational and rotational friction constant can be calculated
128 > from Stoke's law,
129   \[
130   \Xi ^{tt}  = \left( {\begin{array}{*{20}c}
131     {6\pi \eta R} & 0 & 0  \\
# Line 209 | Line 142 | hydrodynamics radius.
142   \end{array}} \right)
143   \]
144   where $\eta$ is the viscosity of the solvent and $R$ is the
145 < hydrodynamics radius.
145 > hydrodynamic radius.
146  
147 < Other non-spherical shape, such as cylinder and ellipsoid
148 < \textit{etc}, are widely used as reference for developing new
149 < hydrodynamics theory, because their properties can be calculated
150 < exactly. In 1936, Perrin extended Stokes's law to general ellipsoid,
151 < also called a triaxial ellipsoid, which is given in Cartesian
152 < coordinates by\cite{Perrin1934, Perrin1936}
147 > Other non-spherical shapes, such as cylinders and ellipsoids, are
148 > widely used as references for developing new hydrodynamics theory,
149 > because their properties can be calculated exactly. In 1936, Perrin
150 > extended Stokes's law to general ellipsoids, also called a triaxial
151 > ellipsoid, which is given in Cartesian coordinates
152 > by\cite{Perrin1934, Perrin1936}
153   \[
154   \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
155   }} = 1
156   \]
157   where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
158   due to the complexity of the elliptic integral, only the ellipsoid
159 < with the restriction of two axes having to be equal, \textit{i.e.}
159 > with the restriction of two axes being equal, \textit{i.e.}
160   prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
161 < exactly. Introducing an elliptic integral parameter $S$ for prolate,
161 > exactly. Introducing an elliptic integral parameter $S$ for prolate
162 > ellipsoids :
163   \[
164   S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
165   } }}{b},
166   \]
167 < and oblate,
167 > and oblate ellipsoids:
168   \[
169   S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
170 < }}{a}
171 < \],
170 > }}{a},
171 > \]
172   one can write down the translational and rotational resistance
173   tensors
174 < \[
175 < \begin{array}{l}
176 < \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}} \\
177 < \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S + 2a}} \\
178 < \end{array},
245 < \]
174 > \begin{eqnarray*}
175 > \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}}. \\
176 > \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S +
177 > 2a}},
178 > \end{eqnarray*}
179   and
180 < \[
181 < \begin{array}{l}
182 < \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}} \\
183 < \Xi _b^{rr}  = \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}} \\
251 < \end{array}.
252 < \]
180 > \begin{eqnarray*}
181 > \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}}, \\
182 > \Xi _b^{rr}  = \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}}.
183 > \end{eqnarray}
184  
185 < \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shape}}
185 > \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}}
186  
187 < Unlike spherical and other regular shaped molecules, there is not
188 < analytical solution for friction tensor of any arbitrary shaped
187 > Unlike spherical and other simply shaped molecules, there is no
188 > analytical solution for the friction tensor for arbitrarily shaped
189   rigid molecules. The ellipsoid of revolution model and general
190   triaxial ellipsoid model have been used to approximate the
191   hydrodynamic properties of rigid bodies. However, since the mapping
192 < from all possible ellipsoidal space, $r$-space, to all possible
193 < combination of rotational diffusion coefficients, $D$-space is not
192 > from all possible ellipsoidal spaces, $r$-space, to all possible
193 > combination of rotational diffusion coefficients, $D$-space, is not
194   unique\cite{Wegener1979} as well as the intrinsic coupling between
195 < translational and rotational motion of rigid body, general ellipsoid
196 < is not always suitable for modeling arbitrarily shaped rigid
197 < molecule. A number of studies have been devoted to determine the
198 < friction tensor for irregularly shaped rigid bodies using more
199 < advanced method where the molecule of interest was modeled by
200 < combinations of spheres(beads)\cite{Carrasco1999} and the
195 > translational and rotational motion of rigid bodies, general
196 > ellipsoids are not always suitable for modeling arbitrarily shaped
197 > rigid molecules. A number of studies have been devoted to
198 > determining the friction tensor for irregularly shaped rigid bodies
199 > using more advanced methods where the molecule of interest was
200 > modeled by a combinations of spheres\cite{Carrasco1999} and the
201   hydrodynamics properties of the molecule can be calculated using the
202   hydrodynamic interaction tensor. Let us consider a rigid assembly of
203 < $N$ beads immersed in a continuous medium. Due to hydrodynamics
203 > $N$ beads immersed in a continuous medium. Due to hydrodynamic
204   interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different
205   than its unperturbed velocity $v_i$,
206   \[
# Line 283 | Line 214 | This equation is the basis for deriving the hydrodynam
214   \label{introEquation:tensorExpression}
215   \end{equation}
216   This equation is the basis for deriving the hydrodynamic tensor. In
217 < 1930, Oseen and Burgers gave a simple solution to Equation
218 < \ref{introEquation:tensorExpression}
217 > 1930, Oseen and Burgers gave a simple solution to
218 > Eq.~\ref{introEquation:tensorExpression}
219   \begin{equation}
220   T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
221   R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor}
# Line 300 | Line 231 | Both of the Equation \ref{introEquation:oseenTensor} a
231   \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
232   \label{introEquation:RPTensorNonOverlapped}
233   \end{equation}
234 < Both of the Equation \ref{introEquation:oseenTensor} and Equation
235 < \ref{introEquation:RPTensorNonOverlapped} have an assumption $R_{ij}
236 < \ge \sigma _i  + \sigma _j$. An alternative expression for
234 > Both of the Eq.~\ref{introEquation:oseenTensor} and
235 > Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption
236 > $R_{ij} \ge \sigma _i  + \sigma _j$. An alternative expression for
237   overlapping beads with the same radius, $\sigma$, is given by
238   \begin{equation}
239   T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
# Line 310 | Line 241 | T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left(
241   \frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
242   \label{introEquation:RPTensorOverlapped}
243   \end{equation}
313 –
244   To calculate the resistance tensor at an arbitrary origin $O$, we
245   construct a $3N \times 3N$ matrix consisting of $N \times N$
246   $B_{ij}$ blocks
# Line 326 | Line 256 | where $\delta _{ij}$ is Kronecker delta function. Inve
256   B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
257   )T_{ij}
258   \]
259 < where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
260 < $B$, we obtain
331 <
259 > where $\delta _{ij}$ is the Kronecker delta function. Inverting the
260 > $B$ matrix, we obtain
261   \[
262   C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
263     {C_{11} } &  \ldots  & {C_{1N} }  \\
264      \vdots  &  \ddots  &  \vdots   \\
265     {C_{N1} } &  \cdots  & {C_{NN} }  \\
266 < \end{array}} \right)
266 > \end{array}} \right),
267   \]
268 < , which can be partitioned into $N \times N$ $3 \times 3$ block
269 < $C_{ij}$. With the help of $C_{ij}$ and skew matrix $U_i$
268 > which can be partitioned into $N \times N$ $3 \times 3$ block
269 > $C_{ij}$. With the help of $C_{ij}$ and the skew matrix $U_i$
270   \[
271   U_i  = \left( {\begin{array}{*{20}c}
272     0 & { - z_i } & {y_i }  \\
# Line 346 | Line 275 | bead $i$ and origin $O$. Hence, the elements of resist
275   \end{array}} \right)
276   \]
277   where $x_i$, $y_i$, $z_i$ are the components of the vector joining
278 < bead $i$ and origin $O$. Hence, the elements of resistance tensor at
278 > bead $i$ and origin $O$, the elements of resistance tensor at
279   arbitrary origin $O$ can be written as
280 < \begin{equation}
281 < \begin{array}{l}
353 < \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\
280 > \begin{eqnarray}
281 > \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\
282   \Xi _{}^{tr}  = \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
283 < \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j  \\
356 < \end{array}
283 > \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j. \notag \\
284   \label{introEquation:ResistanceTensorArbitraryOrigin}
285 < \end{equation}
359 <
285 > \end{eqnarray}
286   The resistance tensor depends on the origin to which they refer. The
287 < proper location for applying friction force is the center of
288 < resistance (reaction), at which the trace of rotational resistance
289 < tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
290 < resistance is defined as an unique point of the rigid body at which
291 < the translation-rotation coupling tensor are symmetric,
287 > proper location for applying the friction force is the center of
288 > resistance (or center of reaction), at which the trace of rotational
289 > resistance tensor, $ \Xi ^{rr}$ reaches a minimum value.
290 > Mathematically, the center of resistance is defined as an unique
291 > point of the rigid body at which the translation-rotation coupling
292 > tensors are symmetric,
293   \begin{equation}
294   \Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
295   \label{introEquation:definitionCR}
296   \end{equation}
297 < Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
298 < we can easily find out that the translational resistance tensor is
297 > From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
298 > we can easily derive that the translational resistance tensor is
299   origin independent, while the rotational resistance tensor and
300   translation-rotation coupling resistance tensor depend on the
301 < origin. Given resistance tensor at an arbitrary origin $O$, and a
302 < vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
301 > origin. Given the resistance tensor at an arbitrary origin $O$, and
302 > a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
303   obtain the resistance tensor at $P$ by
304   \begin{equation}
305   \begin{array}{l}
# Line 390 | Line 317 | Using Equations \ref{introEquation:definitionCR} and
317     { - y_{OP} } & {x_{OP} } & 0  \\
318   \end{array}} \right)
319   \]
320 < Using Equations \ref{introEquation:definitionCR} and
321 < \ref{introEquation:resistanceTensorTransformation}, one can locate
322 < the position of center of resistance,
320 > Using Eq.~\ref{introEquation:definitionCR} and
321 > Eq.~\ref{introEquation:resistanceTensorTransformation}, one can
322 > locate the position of center of resistance,
323   \begin{eqnarray*}
324   \left( \begin{array}{l}
325   x_{OR}  \\
# Line 409 | Line 336 | the position of center of resistance,
336   (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
337   \end{array} \right) \\
338   \end{eqnarray*}
412 –
339   where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
340   joining center of resistance $R$ and origin $O$.
341  
342 < \subsection{Langevin dynamics for rigid particles of arbitrary shape}
342 > \subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}
343  
344 < Consider a Langevin equation of motions in generalized coordinates
344 > Consider the Langevin equations of motion in generalized coordinates
345   \begin{equation}
346   M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)}  + F_{r,i} (t)
347   \label{LDGeneralizedForm}
348   \end{equation}
349   where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
350   and moment of inertial) matrix and $V_i$ is a generalized velocity,
351 < $V_i = V_i(v_i,\omega _i)$. The right side of Eq.
352 < (\ref{LDGeneralizedForm}) consists of three generalized forces in
351 > $V_i = V_i(v_i,\omega _i)$. The right side of
352 > Eq.~\ref{LDGeneralizedForm} consists of three generalized forces in
353   lab-fixed frame, systematic force $F_{s,i}$, dissipative force
354   $F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
355   system in Newtownian mechanics typically refers to lab-fixed frame,
# Line 432 | Line 358 | in body-fixed frame and converted back to lab-fixed fr
358   in body-fixed frame and converted back to lab-fixed frame by:
359   \[
360   \begin{array}{l}
361 < F_{f,i}^l (t) = A^T F_{f,i}^b (t), \\
362 < F_{r,i}^l (t) = A^T F_{r,i}^b (t) \\
363 < \end{array}.
361 > F_{f,i}^l (t) = Q^T F_{f,i}^b (t), \\
362 > F_{r,i}^l (t) = Q^T F_{r,i}^b (t). \\
363 > \end{array}
364   \]
365   Here, the body-fixed friction force $F_{r,i}^b$ is proportional to
366   the body-fixed velocity at center of resistance $v_{R,i}^b$ and
367 < angular velocity $\omega _i$,
367 > angular velocity $\omega _i$
368   \begin{equation}
369   F_{r,i}^b (t) = \left( \begin{array}{l}
370   f_{r,i}^b (t) \\
# Line 456 | Line 382 | with zero mean and variance
382   \begin{equation}
383   \left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
384   \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
385 < 2k_B T\Xi _R \delta (t - t').
385 > 2k_B T\Xi _R \delta (t - t'). \label{randomForce}
386   \end{equation}
461 –
387   The equation of motion for $v_i$ can be written as
388   \begin{equation}
389   m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
# Line 473 | Line 398 | of the resistance. Instead of integrating angular velo
398   \tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b
399   \end{equation}
400   where $r_{MR}$ is the vector from the center of mass to the center
401 < of the resistance. Instead of integrating angular velocity in
402 < lab-fixed frame, we consider the equation of motion of angular
403 < momentum in body-fixed frame
401 > of the resistance. Instead of integrating the angular velocity in
402 > lab-fixed frame, we consider the equation of angular momentum in
403 > body-fixed frame
404   \begin{equation}
405   \dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t)
406   + \tau _{r,i}^b(t)
407   \end{equation}
483 –
408   Embedding the friction terms into force and torque, one can
409   integrate the langevin equations of motion for rigid body of
410   arbitrary shape in a velocity-Verlet style 2-part algorithm, where
# Line 497 | Line 421 | $h= \delta t$:
421   {\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
422      + \frac{h}{2} {\bf \tau}^b(t), \\
423   %
424 < \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
424 > \mathsf{Q}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
425      (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
426   \end{align*}
503 –
427   In this context, the $\mathrm{rotate}$ function is the reversible
428   product of the three body-fixed rotations,
429   \begin{equation}
# Line 509 | Line 432 | rotates both the rotation matrix ($\mathsf{A}$) and th
432   / 2) \cdot \mathsf{G}_x(a_x /2),
433   \end{equation}
434   where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
435 < rotates both the rotation matrix ($\mathsf{A}$) and the body-fixed
435 > rotates both the rotation matrix ($\mathsf{Q}$) and the body-fixed
436   angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed
437   axis $\alpha$,
438   \begin{equation}
439   \mathsf{G}_\alpha( \theta ) = \left\{
440   \begin{array}{lcl}
441 < \mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
441 > \mathsf{Q}(t) & \leftarrow & \mathsf{Q}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
442   {\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
443   j}(0).
444   \end{array}
# Line 535 | Line 458 | All other rotations follow in a straightforward manner
458   \end{array}
459   \right).
460   \end{equation}
461 < All other rotations follow in a straightforward manner.
461 > All other rotations follow in a straightforward manner. After the
462 > first part of the propagation, the forces and body-fixed torques are
463 > calculated at the new positions and orientations
464  
540 – After the first part of the propagation, the forces and body-fixed
541 – torques are calculated at the new positions and orientations
542 –
465   {\tt doForces:}
466   \begin{align*}
467   {\bf f}(t + h) &\leftarrow
# Line 548 | Line 470 | torques are calculated at the new positions and orient
470   {\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h)
471      \times \frac{\partial V}{\partial {\bf u}}, \\
472   %
473 < {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
473 > {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{Q}(t + h)
474      \cdot {\bf \tau}^s(t + h).
475   \end{align*}
476 + Once the forces and torques have been obtained at the new time step,
477 + the velocities can be advanced to the same time value.
478  
555 – {\sc oopse} automatically updates ${\bf u}$ when the rotation matrix
556 – $\mathsf{A}$ is calculated in {\tt moveA}.  Once the forces and
557 – torques have been obtained at the new time step, the velocities can
558 – be advanced to the same time value.
559 –
479   {\tt moveB:}
480   \begin{align*}
481   {\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2
# Line 568 | Line 487 | be advanced to the same time value.
487      + \frac{h}{2} {\bf \tau}^b(t + h) .
488   \end{align*}
489  
490 < \section{Results and discussion}
490 > \section{Results and Discussion}
491 >
492 > The Langevin algorithm described in previous section has been
493 > implemented in {\sc oopse}\cite{Meineke2005} and applied to studies
494 > of the static and dynamic properties in several systems.
495 >
496 > \subsection{Temperature Control}
497 >
498 > As shown in Eq.~\ref{randomForce}, random collisions associated with
499 > the solvent's thermal motions is controlled by the external
500 > temperature. The capability to maintain the temperature of the whole
501 > system was usually used to measure the stability and efficiency of
502 > the algorithm. In order to verify the stability of this new
503 > algorithm, a series of simulations are performed on system
504 > consisiting of 256 SSD water molecules with different viscosities.
505 > The initial configuration for the simulations is taken from a 1ns
506 > NVT simulation with a cubic box of 19.7166~\AA. All simulation are
507 > carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns
508 > with reference temperature at 300~K. The average temperature as a
509 > function of $\eta$ is shown in Table \ref{langevin:viscosity} where
510 > the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 -
511 > 1$ poise. The better temperature control at higher viscosity can be
512 > explained by the finite size effect and relative slow relaxation
513 > rate at lower viscosity regime.
514 > \begin{table}
515 > \caption{AVERAGE TEMPERATURES FROM LANGEVIN DYNAMICS SIMULATIONS OF
516 > SSD WATER MOLECULES WITH REFERENCE TEMPERATURE AT 300~K.}
517 > \label{langevin:viscosity}
518 > \begin{center}
519 > \begin{tabular}{lll}
520 >  \hline
521 >  $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\
522 >  \hline
523 >  1    & 300.47 & 10.99 \\
524 >  0.1  & 301.19 & 11.136 \\
525 >  0.01 & 303.04 & 11.796 \\
526 >  \hline
527 > \end{tabular}
528 > \end{center}
529 > \end{table}
530  
531 + Another set of calculation were performed to study the efficiency of
532 + temperature control using different temperature coupling schemes.
533 + The starting configuration is cooled to 173~K and evolved using NVE,
534 + NVT, and Langevin dynamic with time step of 2 fs.
535 + Fig.~\ref{langevin:temperature} shows the heating curve obtained as
536 + the systems reach equilibrium. The orange curve in
537 + Fig.~\ref{langevin:temperature} represents the simulation using
538 + Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps
539 + which gives reasonable tight coupling, while the blue one from
540 + Langevin dynamics with viscosity of 0.1 poise demonstrates a faster
541 + scaling to the desire temperature. In extremely lower friction
542 + regime (when $ \eta  \approx 0$), Langevin dynamics becomes normal
543 + NVE (see orange curve in Fig.~\ref{langevin:temperature}) which
544 + loses the temperature control ability.
545 +
546 + \begin{figure}
547 + \centering
548 + \includegraphics[width=\linewidth]{temperature.eps}
549 + \caption[Plot of Temperature Fluctuation Versus Time]{Plot of
550 + temperature fluctuation versus time.} \label{langevin:temperature}
551 + \end{figure}
552 +
553 + \subsection{Langevin Dynamics of Banana Shaped Molecules}
554 +
555 + In order to verify that Langevin dynamics can mimic the dynamics of
556 + the systems absent of explicit solvents, we carried out two sets of
557 + simulations and compare their dynamic properties.
558 + Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation
559 + made of 256 pentane molecules and two banana shaped molecules at
560 + 273~K. It has an equivalent implicit solvent system containing only
561 + two banana shaped molecules with viscosity of 0.289 center poise. To
562 + calculate the hydrodynamic properties of the banana shaped molecule,
563 + we created a rough shell model (see Fig.~\ref{langevin:roughShell}),
564 + in which the banana shaped molecule is represented as a ``shell''
565 + made of 2266 small identical beads with size of 0.3 \AA on the
566 + surface. Applying the procedure described in
567 + Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
568 + identified the center of resistance at $(0\AA, 0.7482\AA,
569 + -0.1988\AA)$, as well as the resistance tensor,
570 + \[
571 + \left( {\begin{array}{*{20}c}
572 + 0.9261 & 0 & 0&0&0.08585&0.2057\\
573 + 0& 0.9270&-0.007063& 0.08585&0&0\\
574 + 0&0.007063&0.7494&0.2057&0&0\\
575 + 0&0.0858&0.2057& 58.64& 0&-8.5736\\
576 + 0.08585&0&0&0&48.30&3.219&\\
577 + 0.2057&0&0&0&3.219&10.7373\\
578 + \end{array}} \right).
579 + \]
580 + %\[
581 + %\left( {\begin{array}{*{20}c}
582 + %0.9261 & 1.310e-14 & -7.292e-15&5.067e-14&0.08585&0.2057\\
583 + %3.968e-14& 0.9270&-0.007063& 0.08585&6.764e-14&4.846e-14\\
584 + %-6.561e-16&-0.007063&0.7494&0.2057&4.846e-14&1.5036e-14\\
585 + %5.067e-14&0.0858&0.2057& 58.64& 8.563e-13&-8.5736\\
586 + %0.08585&6.764e-14&4.846e-14&1.555e-12&48.30&3.219&\\
587 + %0.2057&4.846e-14&1.5036e-14&-3.904e-13&3.219&10.7373\\
588 + %\end{array}} \right).
589 + %\]
590 +
591 + Curves of the velocity auto-correlation functions in
592 + Fig.~\ref{langevin:vacf} were shown to match each other very well.
593 + However, because of the stochastic nature, simulation using Langevin
594 + dynamics was shown to decay slightly faster than MD. In order to
595 + study the rotational motion of the molecules, we also calculated the
596 + auto- correlation function of the principle axis of the second GB
597 + particle, $u$. The discrepancy shown in Fig.~\ref{langevin:uacf} was
598 + probably due to the reason that the viscosity using in the
599 + simulations only partially preserved the dynamics of the system.
600 +
601 + \begin{figure}
602 + \centering
603 + \includegraphics[width=\linewidth]{roughShell.eps}
604 + \caption[Rough shell model for banana shaped molecule]{Rough shell
605 + model for banana shaped molecule.} \label{langevin:roughShell}
606 + \end{figure}
607 +
608 + \begin{figure}
609 + \centering
610 + \includegraphics[width=\linewidth]{twoBanana.eps}
611 + \caption[Snapshot from Simulation of Two Banana Shaped Molecules and
612 + 256 Pentane Molecules]{Snapshot from simulation of two Banana shaped
613 + molecules and 256 pentane molecules.} \label{langevin:twoBanana}
614 + \end{figure}
615 +
616 + \begin{figure}
617 + \centering
618 + \includegraphics[width=\linewidth]{vacf.eps}
619 + \caption[Plots of Velocity Auto-correlation Functions]{Velocity
620 + auto-correlation functions of NVE (explicit solvent) in blue) and
621 + Langevin dynamics (implicit solvent) in red.} \label{langevin:vacf}
622 + \end{figure}
623 +
624 + \begin{figure}
625 + \centering
626 + \includegraphics[width=\linewidth]{uacf.eps}
627 + \caption[Auto-correlation functions of the principle axis of the
628 + middle GB particle]{Auto-correlation functions of the principle axis
629 + of the middle GB particle of NVE (blue) and Langevin dynamics
630 + (red).} \label{langevin:uacf}
631 + \end{figure}
632 +
633   \section{Conclusions}
634 +
635 + We have presented a new Langevin algorithm by incorporating the
636 + hydrodynamics properties of arbitrary shaped molecules into an
637 + advanced symplectic integration scheme. The temperature control
638 + ability of this algorithm was demonstrated by a set of simulations
639 + with different viscosities. It was also shown to have significant
640 + advantage of producing rapid thermal equilibration over
641 + Nos\'{e}-Hoover method. Further studies in systems involving banana
642 + shaped molecules illustrated that the dynamic properties could be
643 + preserved by using this new algorithm as an implicit solvent model.

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines