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 2876 by tim, Wed Jun 21 19:46:01 2006 UTC

# Line 1 | Line 1
1  
2 < \chapter{\label{chapt:methodology}Langevin Dynamics for Rigid Bodies of Arbitrary Shape}
2 > \chapter{\label{chapt:methodology}LANGEVIN DYNAMICS FOR RIGID BODIES OF ARBITRARY SHAPE}
3  
4   \section{Introduction}
5  
# Line 47 | Line 47 | efficiently\cite{Sandu1999}.
47   has been shown to be able to damp out the resonance artifact more
48   efficiently\cite{Sandu1999}.
49  
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 –
50   %review langevin/browninan dynamics for arbitrarily shaped rigid body
51   Combining Langevin or Brownian dynamics with rigid body dynamics,
52   one can study the slow processes in biomolecular systems. Modeling
# Line 151 | Line 98 | sophisticated rigid body dynamics.
98   estimation of friction tensor from hydrodynamics theory into the
99   sophisticated rigid body dynamics.
100  
101 < \section{Method{\label{methodSec}}}
101 > \section{Computational Methods{\label{methodSec}}}
102  
103 < \subsection{\label{introSection:frictionTensor} Friction Tensor}
104 < Theoretically, the friction kernel can be determined using velocity
105 < autocorrelation function. However, this approach become impractical
106 < when the system become more and more complicate. Instead, various
107 < approaches based on hydrodynamics have been developed to calculate
108 < the friction coefficients. The friction effect is isotropic in
109 < Equation, $\zeta$ can be taken as a scalar. In general, friction
163 < tensor $\Xi$ is a $6\times 6$ matrix given by
103 > \subsection{\label{introSection:frictionTensor}Friction Tensor}
104 > Theoretically, the friction kernel can be determined using the
105 > velocity autocorrelation function. However, this approach become
106 > impractical when the system become more and more complicate.
107 > Instead, various approaches based on hydrodynamics have been
108 > developed to calculate the friction coefficients. In general,
109 > friction tensor $\Xi$ is a $6\times 6$ matrix given by
110   \[
111   \Xi  = \left( {\begin{array}{*{20}c}
112     {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
# Line 191 | Line 137 | For a spherical particle, the translational and rotati
137  
138   \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}}
139  
140 < For a spherical particle, the translational and rotational friction
141 < constant can be calculated from Stoke's law,
140 > For a spherical particle with slip boundary conditions, the
141 > translational and rotational friction constant can be calculated
142 > from Stoke's law,
143   \[
144   \Xi ^{tt}  = \left( {\begin{array}{*{20}c}
145     {6\pi \eta R} & 0 & 0  \\
# Line 214 | Line 161 | exactly. In 1936, Perrin extended Stokes's law to gene
161   Other non-spherical shape, such as cylinder and ellipsoid
162   \textit{etc}, are widely used as reference for developing new
163   hydrodynamics theory, because their properties can be calculated
164 < exactly. In 1936, Perrin extended Stokes's law to general ellipsoid,
165 < also called a triaxial ellipsoid, which is given in Cartesian
166 < coordinates by\cite{Perrin1934, Perrin1936}
164 > exactly. In 1936, Perrin extended Stokes's law to general
165 > ellipsoids, also called a triaxial ellipsoid, which is given in
166 > Cartesian coordinates by\cite{Perrin1934, Perrin1936}
167   \[
168   \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
169   }} = 1
# Line 251 | Line 198 | and
198   \end{array}.
199   \]
200  
201 < \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shape}}
201 > \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}}
202  
203 < Unlike spherical and other regular shaped molecules, there is not
204 < analytical solution for friction tensor of any arbitrary shaped
203 > Unlike spherical and other simply shaped molecules, there is no
204 > analytical solution for the friction tensor for arbitrarily shaped
205   rigid molecules. The ellipsoid of revolution model and general
206   triaxial ellipsoid model have been used to approximate the
207   hydrodynamic properties of rigid bodies. However, since the mapping
208 < from all possible ellipsoidal space, $r$-space, to all possible
209 < combination of rotational diffusion coefficients, $D$-space is not
208 > from all possible ellipsoidal spaces, $r$-space, to all possible
209 > combination of rotational diffusion coefficients, $D$-space, is not
210   unique\cite{Wegener1979} as well as the intrinsic coupling between
211 < translational and rotational motion of rigid body, general ellipsoid
212 < is not always suitable for modeling arbitrarily shaped rigid
213 < molecule. A number of studies have been devoted to determine the
214 < friction tensor for irregularly shaped rigid bodies using more
211 > translational and rotational motion of rigid body, general
212 > ellipsoids are not always suitable for modeling arbitrarily shaped
213 > rigid molecule. A number of studies have been devoted to determine
214 > the friction tensor for irregularly shaped rigid bodies using more
215   advanced method where the molecule of interest was modeled by
216   combinations of spheres(beads)\cite{Carrasco1999} and the
217   hydrodynamics properties of the molecule can be calculated using the
# Line 326 | Line 273 | where $\delta _{ij}$ is Kronecker delta function. Inve
273   B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
274   )T_{ij}
275   \]
276 < where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
277 < $B$, we obtain
276 > where $\delta _{ij}$ is Kronecker delta function. Inverting the $B$
277 > matrix, we obtain
278  
279   \[
280   C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
# Line 359 | Line 306 | resistance (reaction), at which the trace of rotationa
306  
307   The resistance tensor depends on the origin to which they refer. The
308   proper location for applying friction force is the center of
309 < resistance (reaction), at which the trace of rotational resistance
310 < tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
311 < resistance is defined as an unique point of the rigid body at which
312 < the translation-rotation coupling tensor are symmetric,
309 > resistance (or center of reaction), at which the trace of rotational
310 > resistance tensor, $ \Xi ^{rr}$ reaches a minimum value.
311 > Mathematically, the center of resistance is defined as an unique
312 > point of the rigid body at which the translation-rotation coupling
313 > tensor are symmetric,
314   \begin{equation}
315   \Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
316   \label{introEquation:definitionCR}
317   \end{equation}
318 < Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
318 > From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
319   we can easily find out that the translational resistance tensor is
320   origin independent, while the rotational resistance tensor and
321   translation-rotation coupling resistance tensor depend on the
322 < origin. Given resistance tensor at an arbitrary origin $O$, and a
323 < vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
322 > origin. Given the resistance tensor at an arbitrary origin $O$, and
323 > a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
324   obtain the resistance tensor at $P$ by
325   \begin{equation}
326   \begin{array}{l}
# Line 413 | Line 361 | joining center of resistance $R$ and origin $O$.
361   where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
362   joining center of resistance $R$ and origin $O$.
363  
364 < \subsection{Langevin dynamics for rigid particles of arbitrary shape}
364 > \subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}
365  
366   Consider a Langevin equation of motions in generalized coordinates
367   \begin{equation}
# Line 422 | Line 370 | $V_i = V_i(v_i,\omega _i)$. The right side of Eq.
370   \end{equation}
371   where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
372   and moment of inertial) matrix and $V_i$ is a generalized velocity,
373 < $V_i = V_i(v_i,\omega _i)$. The right side of Eq.
373 > $V_i = V_i(v_i,\omega _i)$. The right side of Eq.~
374   (\ref{LDGeneralizedForm}) consists of three generalized forces in
375   lab-fixed frame, systematic force $F_{s,i}$, dissipative force
376   $F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
# Line 456 | Line 404 | with zero mean and variance
404   \begin{equation}
405   \left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
406   \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
407 < 2k_B T\Xi _R \delta (t - t').
407 > 2k_B T\Xi _R \delta (t - t'). \label{randomForce}
408   \end{equation}
409  
410   The equation of motion for $v_i$ can be written as
# Line 568 | Line 516 | be advanced to the same time value.
516      + \frac{h}{2} {\bf \tau}^b(t + h) .
517   \end{align*}
518  
519 < \section{Results and discussion}
519 > \section{Results and Discussion}
520  
521 + The Langevin algorithm described in previous section has been
522 + implemented in {\sc oopse}\cite{Meineke2005} and applied to the
523 + studies of kinetics and thermodynamic properties in several systems.
524 +
525 + \subsection{Temperature Control}
526 +
527 + As shown in Eq.~\ref{randomForce}, random collisions associated with
528 + the solvent's thermal motions is controlled by the external
529 + temperature. The capability to maintain the temperature of the whole
530 + system was usually used to measure the stability and efficiency of
531 + the algorithm. In order to verify the stability of this new
532 + algorithm, a series of simulations are performed on system
533 + consisiting of 256 SSD water molecules with different viscosities.
534 + Initial configuration for the simulations is taken from a 1ns NVT
535 + simulation with a cubic box of 19.7166~\AA. All simulation are
536 + carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns
537 + with reference temperature at 300~K. Average temperature as a
538 + function of $\eta$ is shown in Table \ref{langevin:viscosity} where
539 + the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 -
540 + 1$ poise. The better temperature control at higher viscosity can be
541 + explained by the finite size effect and relative slow relaxation
542 + rate at lower viscosity regime.
543 + \begin{table}
544 + \caption{Average temperatures from Langevin dynamics simulations of
545 + SSD water molecules with reference temperature at 300~K.}
546 + \label{langevin:viscosity}
547 + \begin{center}
548 + \begin{tabular}{|l|l|l|}
549 +  \hline
550 +  $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\
551 +  1    & 300.47 & 10.99 \\
552 +  0.1  & 301.19 & 11.136 \\
553 +  0.01 & 303.04 & 11.796 \\
554 +  \hline
555 + \end{tabular}
556 + \end{center}
557 + \end{table}
558 +
559 + Another set of calculation were performed to study the efficiency of
560 + temperature control using different temperature coupling schemes.
561 + The starting configuration is cooled to 173~K and evolved using NVE,
562 + NVT, and Langevin dynamic with time step of 2 fs.
563 + Fig.~\ref{langevin:temperature} shows the heating curve obtained as
564 + the systems reach equilibrium. The orange curve in
565 + Fig.~\ref{langevin:temperature} represents the simulation using
566 + Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps
567 + which gives reasonable tight coupling, while the blue one from
568 + Langevin dynamics with viscosity of 0.1 poise demonstrates a faster
569 + scaling to the desire temperature. In extremely lower friction
570 + regime (when $ \eta  \approx 0$), Langevin dynamics becomes normal
571 + NVE (see green curve in Fig.~\ref{langevin:temperature}) which loses
572 + the temperature control ability.
573 +
574 + \begin{figure}
575 + \centering
576 + \includegraphics[width=\linewidth]{temperature.eps}
577 + \caption[Plot of Temperature Fluctuation Versus Time]{Plot of
578 + temperature fluctuation versus time.} \label{langevin:temperature}
579 + \end{figure}
580 +
581 + \subsection{Langevin Dynamics of Banana Shaped Molecule}
582 +
583 + In order to verify that Langevin dynamics can mimic the dynamics of
584 + the systems absent of explicit solvents, we carried out two sets of
585 + simulations and compare their dynamic properties.
586 +
587 + \subsubsection{Simulations Contain One Banana Shaped Molecule}
588 +
589 + Fig.~\ref{langevin:oneBanana} shows a snapshot of the simulation
590 + made of 256 pentane molecules and one banana shaped molecule at
591 + 273~K. It has an equivalent implicit solvent model containing only
592 + one banana shaped molecule with viscosity of 0.289 center poise. To
593 + calculate the hydrodynamic properties of the banana shaped molecule,
594 + we create a rough shell model (see Fig.~\ref{langevin:roughShell}),
595 + in which the banana shaped molecule is represented as a ``shell''
596 + made of 2266 small identical beads with size of 0.3 $\AA$ on the
597 + surface. Applying the procedure described in
598 + Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
599 + identified the center of resistance at $(0, 0.7482, -0.1988)$, as
600 + well as the resistance tensor,
601 + \[
602 + \left( {\begin{array}{*{20}c}
603 + 0.9261 & 1.310e-14 & -7.292e-15&5.067e-14&0.08585&0.2057\\
604 + 3.968e-14& 0.9270&-0.007063& 0.08585&6.764e-14&4.846e-14\\
605 + -6.561e-16&-0.007063&0.7494&0.2057&4.846e-14&1.5036e-14\\
606 + 5.067e-14&0.0858&0.2057& 58.64& 8.563e-13&-8.5736\\
607 + 0.08585&6.764e-14&4.846e-14&1.555e-12&48.30&3.219&\\
608 + 0.2057&4.846e-14&1.5036e-14&-3.904e-13&3.219&10.7373\\
609 + \end{array}} \right).
610 + \]
611 +
612 + \begin{figure}
613 + \centering
614 + \includegraphics[width=\linewidth]{roughShell.eps}
615 + \caption[Rough shell model for banana shaped molecule]{Rough shell
616 + model for banana shaped molecule.} \label{langevin:roughShell}
617 + \end{figure}
618 +
619 + \begin{figure}
620 + \centering
621 + \includegraphics[width=\linewidth]{oneBanana.eps}
622 + \caption[Snapshot from Simulation of One Banana Shaped Molecules and
623 + 256 Pentane Molecules]{Snapshot from simulation of one Banana shaped
624 + molecules and 256 pentane molecules.} \label{langevin:oneBanana}
625 + \end{figure}
626 +
627 + \subsubsection{Simulations Contain Two Banana Shaped Molecules}
628 +
629 + \begin{figure}
630 + \centering
631 + \includegraphics[width=\linewidth]{twoBanana.eps}
632 + \caption[Snapshot from Simulation of Two Banana Shaped Molecules and
633 + 256 Pentane Molecules]{Snapshot from simulation of two Banana shaped
634 + molecules and 256 pentane molecules.} \label{langevin:twoBanana}
635 + \end{figure}
636 +
637   \section{Conclusions}

Diff Legend

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