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 2858 by tim, Tue Jun 13 02:08:23 2006 UTC vs.
Revision 2890 by tim, Mon Jun 26 13:43:22 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  
# Line 28 | Line 28 | systems\cite{Garcia-Palacios1998,Berkov2002,Denisov200
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
32 < instance, the oscillation power spectrum of nanoparticles from
33 < Langevin dynamics simulation has the same peak frequencies for
34 < different wave vectors,which recovers the property of magnetic
35 < 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}.
31 > systems\cite{Palacios1998,Berkov2002,Denisov2003}. For instance, the
32 > oscillation power spectrum of nanoparticles from Langevin dynamics
33 > simulation has the same peak frequencies for different wave
34 > vectors,which recovers the property of magnetic excitations in small
35 > finite structures\cite{Berkov2005a}. In an attempt to reduce the
36 > computational cost of simulation, multiple time stepping (MTS)
37 > methods have been introduced and have been of great interest to
38 > macromolecule and protein community\cite{Tuckerman1992}. Relying on
39 > the observation that forces between distant atoms generally
40 > demonstrate slower fluctuations than forces between close atoms, MTS
41 > method are generally implemented by evaluating the slowly
42 > fluctuating forces less frequently than the fast ones.
43 > Unfortunately, nonlinear instability resulting from increasing
44 > timestep in MTS simulation have became a critical obstruction
45 > preventing the long time simulation. Due to the coupling to the heat
46 > bath, Langevin dynamics has been shown to be able to damp out the
47 > resonance artifact more efficiently\cite{Sandu1999}.
48  
49   %review langevin/browninan dynamics for arbitrarily shaped rigid body
50   Combining Langevin or Brownian dynamics with rigid body dynamics,
# Line 79 | Line 78 | term\cite{Beard2001}. As a complement to IBD which has
78   average acceleration is not always true for cooperative motion which
79   is common in protein motion. An inertial Brownian dynamics (IBD) was
80   proposed to address this issue by adding an inertial correction
81 < term\cite{Beard2001}. As a complement to IBD which has a lower bound
81 > term\cite{Beard2000}. As a complement to IBD which has a lower bound
82   in time step because of the inertial relaxation time, long-time-step
83   inertial dynamics (LTID) can be used to investigate the inertial
84   behavior of the polymer segments in low friction
85 < regime\cite{Beard2001}. LTID can also deal with the rotational
85 > regime\cite{Beard2000}. LTID can also deal with the rotational
86   dynamics for nonskew bodies without translation-rotation coupling by
87   separating the translation and rotation motion and taking advantage
88   of the analytical solution of hydrodynamics properties. However,
# Line 98 | Line 97 | sophisticated rigid body dynamics.
97   estimation of friction tensor from hydrodynamics theory into the
98   sophisticated rigid body dynamics.
99  
100 < \section{Computational methods{\label{methodSec}}}
100 > \section{Computational Methods{\label{methodSec}}}
101  
102   \subsection{\label{introSection:frictionTensor}Friction Tensor}
103 < Theoretically, the friction kernel can be determined using velocity
104 < autocorrelation function. However, this approach become impractical
105 < when the system become more and more complicate. Instead, various
106 < approaches based on hydrodynamics have been developed to calculate
107 < the friction coefficients. The friction effect is isotropic in
108 < Equation, $\zeta$ can be taken as a scalar. In general, friction
110 < tensor $\Xi$ is a $6\times 6$ matrix given by
103 > Theoretically, the friction kernel can be determined using the
104 > velocity autocorrelation function. However, this approach become
105 > impractical when the system become more and more complicate.
106 > Instead, various approaches based on hydrodynamics have been
107 > developed to calculate the friction coefficients. In general,
108 > friction tensor $\Xi$ is a $6\times 6$ matrix given by
109   \[
110   \Xi  = \left( {\begin{array}{*{20}c}
111     {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
# Line 138 | Line 136 | For a spherical particle, the translational and rotati
136  
137   \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}}
138  
139 < For a spherical particle, the translational and rotational friction
140 < constant can be calculated from Stoke's law,
139 > For a spherical particle with slip boundary conditions, the
140 > translational and rotational friction constant can be calculated
141 > from Stoke's law,
142   \[
143   \Xi ^{tt}  = \left( {\begin{array}{*{20}c}
144     {6\pi \eta R} & 0 & 0  \\
# Line 161 | Line 160 | exactly. In 1936, Perrin extended Stokes's law to gene
160   Other non-spherical shape, such as cylinder and ellipsoid
161   \textit{etc}, are widely used as reference for developing new
162   hydrodynamics theory, because their properties can be calculated
163 < exactly. In 1936, Perrin extended Stokes's law to general ellipsoid,
164 < also called a triaxial ellipsoid, which is given in Cartesian
165 < coordinates by\cite{Perrin1934, Perrin1936}
163 > exactly. In 1936, Perrin extended Stokes's law to general
164 > ellipsoids, also called a triaxial ellipsoid, which is given in
165 > Cartesian coordinates by\cite{Perrin1934, Perrin1936}
166   \[
167   \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
168   }} = 1
# Line 172 | Line 171 | exactly. Introducing an elliptic integral parameter $S
171   due to the complexity of the elliptic integral, only the ellipsoid
172   with the restriction of two axes having to be equal, \textit{i.e.}
173   prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
174 < exactly. Introducing an elliptic integral parameter $S$ for prolate,
174 > exactly. Introducing an elliptic integral parameter $S$ for prolate
175 > :
176   \[
177   S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
178   } }}{b},
179   \]
180 < and oblate,
180 > and oblate :
181   \[
182   S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
183   }}{a}
# Line 186 | Line 186 | tensors
186   tensors
187   \[
188   \begin{array}{l}
189 < \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}} \\
190 < \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S + 2a}} \\
191 < \end{array},
189 > \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}}. \\
190 > \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S +
191 > 2a}},
192 > \end{array}
193   \]
194   and
195   \[
196   \begin{array}{l}
197 < \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}} \\
197 > \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}}, \\
198   \Xi _b^{rr}  = \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}} \\
199   \end{array}.
200   \]
201  
202 < \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shape}}
202 > \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}}
203  
204 < Unlike spherical and other regular shaped molecules, there is not
205 < analytical solution for friction tensor of any arbitrary shaped
204 > Unlike spherical and other simply shaped molecules, there is no
205 > analytical solution for the friction tensor for arbitrarily shaped
206   rigid molecules. The ellipsoid of revolution model and general
207   triaxial ellipsoid model have been used to approximate the
208   hydrodynamic properties of rigid bodies. However, since the mapping
209 < from all possible ellipsoidal space, $r$-space, to all possible
210 < combination of rotational diffusion coefficients, $D$-space is not
209 > from all possible ellipsoidal spaces, $r$-space, to all possible
210 > combination of rotational diffusion coefficients, $D$-space, is not
211   unique\cite{Wegener1979} as well as the intrinsic coupling between
212 < translational and rotational motion of rigid body, general ellipsoid
213 < is not always suitable for modeling arbitrarily shaped rigid
214 < molecule. A number of studies have been devoted to determine the
215 < friction tensor for irregularly shaped rigid bodies using more
212 > translational and rotational motion of rigid body, general
213 > ellipsoids are not always suitable for modeling arbitrarily shaped
214 > rigid molecule. A number of studies have been devoted to determine
215 > the friction tensor for irregularly shaped rigid bodies using more
216   advanced method where the molecule of interest was modeled by
217   combinations of spheres(beads)\cite{Carrasco1999} and the
218   hydrodynamics properties of the molecule can be calculated using the
# Line 273 | Line 274 | where $\delta _{ij}$ is Kronecker delta function. Inve
274   B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
275   )T_{ij}
276   \]
277 < where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
278 < $B$, we obtain
277 > where $\delta _{ij}$ is Kronecker delta function. Inverting the $B$
278 > matrix, we obtain
279  
280   \[
281   C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
# Line 306 | Line 307 | resistance (reaction), at which the trace of rotationa
307  
308   The resistance tensor depends on the origin to which they refer. The
309   proper location for applying friction force is the center of
310 < resistance (reaction), at which the trace of rotational resistance
311 < tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
312 < resistance is defined as an unique point of the rigid body at which
313 < the translation-rotation coupling tensor are symmetric,
310 > resistance (or center of reaction), at which the trace of rotational
311 > resistance tensor, $ \Xi ^{rr}$ reaches a minimum value.
312 > Mathematically, the center of resistance is defined as an unique
313 > point of the rigid body at which the translation-rotation coupling
314 > tensor are symmetric,
315   \begin{equation}
316   \Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
317   \label{introEquation:definitionCR}
318   \end{equation}
319 < Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
319 > From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
320   we can easily find out that the translational resistance tensor is
321   origin independent, while the rotational resistance tensor and
322   translation-rotation coupling resistance tensor depend on the
323 < origin. Given resistance tensor at an arbitrary origin $O$, and a
324 < vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
323 > origin. Given the resistance tensor at an arbitrary origin $O$, and
324 > a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
325   obtain the resistance tensor at $P$ by
326   \begin{equation}
327   \begin{array}{l}
# Line 360 | Line 362 | joining center of resistance $R$ and origin $O$.
362   where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
363   joining center of resistance $R$ and origin $O$.
364  
365 < \subsection{Langevin dynamics for rigid particles of arbitrary shape\label{LDRB}}
365 > \subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}
366  
367   Consider a Langevin equation of motions in generalized coordinates
368   \begin{equation}
# Line 403 | Line 405 | with zero mean and variance
405   \begin{equation}
406   \left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
407   \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
408 < 2k_B T\Xi _R \delta (t - t').
408 > 2k_B T\Xi _R \delta (t - t'). \label{randomForce}
409   \end{equation}
410  
411   The equation of motion for $v_i$ can be written as
# Line 515 | Line 517 | be advanced to the same time value.
517      + \frac{h}{2} {\bf \tau}^b(t + h) .
518   \end{align*}
519  
520 < \section{Results and discussion}
520 > \section{Results and Discussion}
521  
522 < The Langevin algorithm described in Sec.~\ref{LDRB} has been
523 < implemented in {\sc oopse}\cite{Meineke2005} and applied to several
524 < test systems.
523 <
524 < \subsection{Langevin dynamics of}
522 > The Langevin algorithm described in previous section has been
523 > implemented in {\sc oopse}\cite{Meineke2005} and applied to the
524 > studies of kinetics and thermodynamic properties in several systems.
525  
526 + \subsection{Temperature Control}
527 +
528 + As shown in Eq.~\ref{randomForce}, random collisions associated with
529 + the solvent's thermal motions is controlled by the external
530 + temperature. The capability to maintain the temperature of the whole
531 + system was usually used to measure the stability and efficiency of
532 + the algorithm. In order to verify the stability of this new
533 + algorithm, a series of simulations are performed on system
534 + consisiting of 256 SSD water molecules with different viscosities.
535 + Initial configuration for the simulations is taken from a 1ns NVT
536 + simulation with a cubic box of 19.7166~\AA. All simulation are
537 + carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns
538 + with reference temperature at 300~K. Average temperature as a
539 + function of $\eta$ is shown in Table \ref{langevin:viscosity} where
540 + the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 -
541 + 1$ poise. The better temperature control at higher viscosity can be
542 + explained by the finite size effect and relative slow relaxation
543 + rate at lower viscosity regime.
544 + \begin{table}
545 + \caption{AVERAGE TEMPERATURES FROM LANGEVIN DYNAMICS SIMULATIONS OF
546 + SSD WATER MOLECULES WITH REFERENCE TEMPERATURE AT 300~K.}
547 + \label{langevin:viscosity}
548 + \begin{center}
549 + \begin{tabular}{|l|l|l|}
550 +  \hline
551 +  $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\
552 +  1    & 300.47 & 10.99 \\
553 +  0.1  & 301.19 & 11.136 \\
554 +  0.01 & 303.04 & 11.796 \\
555 +  \hline
556 + \end{tabular}
557 + \end{center}
558 + \end{table}
559 +
560 + Another set of calculation were performed to study the efficiency of
561 + temperature control using different temperature coupling schemes.
562 + The starting configuration is cooled to 173~K and evolved using NVE,
563 + NVT, and Langevin dynamic with time step of 2 fs.
564 + Fig.~\ref{langevin:temperature} shows the heating curve obtained as
565 + the systems reach equilibrium. The orange curve in
566 + Fig.~\ref{langevin:temperature} represents the simulation using
567 + Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps
568 + which gives reasonable tight coupling, while the blue one from
569 + Langevin dynamics with viscosity of 0.1 poise demonstrates a faster
570 + scaling to the desire temperature. In extremely lower friction
571 + regime (when $ \eta  \approx 0$), Langevin dynamics becomes normal
572 + NVE (see green curve in Fig.~\ref{langevin:temperature}) which loses
573 + the temperature control ability.
574 +
575   \begin{figure}
576   \centering
577   \includegraphics[width=\linewidth]{temperature.eps}
578 < \caption[]{.} \label{langevin:temperature}
578 > \caption[Plot of Temperature Fluctuation Versus Time]{Plot of
579 > temperature fluctuation versus time.} \label{langevin:temperature}
580   \end{figure}
581  
582 < \subsection{LD of banana-shaped molecule}
582 > \subsection{Langevin Dynamics of Banana Shaped Molecules}
583  
584 + In order to verify that Langevin dynamics can mimic the dynamics of
585 + the systems absent of explicit solvents, we carried out two sets of
586 + simulations and compare their dynamic properties.
587 + Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation
588 + made of 256 pentane molecules and two banana shaped molecules at
589 + 273~K. It has an equivalent implicit solvent system containing only
590 + two banana shaped molecules with viscosity of 0.289 center poise. To
591 + calculate the hydrodynamic properties of the banana shaped molecule,
592 + we create a rough shell model (see Fig.~\ref{langevin:roughShell}),
593 + in which the banana shaped molecule is represented as a ``shell''
594 + made of 2266 small identical beads with size of 0.3 \AA on the
595 + surface. Applying the procedure described in
596 + Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
597 + identified the center of resistance at $(0, 0.7482, -0.1988)$, as
598 + well as the resistance tensor,
599 + \[
600 + \left( {\begin{array}{*{20}c}
601 + 0.9261 & 1.310e-14 & -7.292e-15&5.067e-14&0.08585&0.2057\\
602 + 3.968e-14& 0.9270&-0.007063& 0.08585&6.764e-14&4.846e-14\\
603 + -6.561e-16&-0.007063&0.7494&0.2057&4.846e-14&1.5036e-14\\
604 + 5.067e-14&0.0858&0.2057& 58.64& 8.563e-13&-8.5736\\
605 + 0.08585&6.764e-14&4.846e-14&1.555e-12&48.30&3.219&\\
606 + 0.2057&4.846e-14&1.5036e-14&-3.904e-13&3.219&10.7373\\
607 + \end{array}} \right).
608 + \]
609 + Curves of velocity auto-correlation functions in
610 + Fig.~\ref{langevin:vacf} were shown to match each other very well.
611 + However, because of the stochastic nature, simulation using Langevin
612 + dynamics was shown to decay slightly fast. In order to study the
613 + rotational motion of the molecules, we also calculated the auto-
614 + correlation function of the principle axis of the second GB
615 + particle, $u$.
616 +
617   \begin{figure}
618   \centering
619 < \includegraphics[width=\linewidth]{one_banana.eps}
620 < \caption[]{.} \label{langevin:banana}
619 > \includegraphics[width=\linewidth]{roughShell.eps}
620 > \caption[Rough shell model for banana shaped molecule]{Rough shell
621 > model for banana shaped molecule.} \label{langevin:roughShell}
622   \end{figure}
623  
624 + \begin{figure}
625 + \centering
626 + \includegraphics[width=\linewidth]{twoBanana.eps}
627 + \caption[Snapshot from Simulation of Two Banana Shaped Molecules and
628 + 256 Pentane Molecules]{Snapshot from simulation of two Banana shaped
629 + molecules and 256 pentane molecules.} \label{langevin:twoBanana}
630 + \end{figure}
631 +
632 + \begin{figure}
633 + \centering
634 + \includegraphics[width=\linewidth]{vacf.eps}
635 + \caption[Plots of Velocity Auto-correlation Functions]{Velocity
636 + auto-correlation functions in NVE (blue) and Langevin dynamics
637 + (red).} \label{langevin:vacf}
638 + \end{figure}
639 +
640 + \begin{figure}
641 + \centering
642 + \includegraphics[width=\linewidth]{uacf.eps}
643 + \caption[Auto-correlation functions of the principle axis of the
644 + middle GB particle]{Auto-correlation functions of the principle axis
645 + of the middle GB particle in NVE (blue) and Langevin dynamics
646 + (red).} \label{langevin:twoBanana}
647 + \end{figure}
648 +
649   \section{Conclusions}
650 +
651 + We have presented a new Langevin algorithm by incorporating the
652 + hydrodynamics properties of arbitrary shaped molecules into an
653 + advanced symplectic integration scheme. The temperature control
654 + ability of this algorithm was demonstrated by a set of simulations
655 + with different viscosities. It was also shown to have significant
656 + advantage of producing rapid thermal equilibration over
657 + Nos\'{e}-Hoover method. Further studies in systems involving banana
658 + shaped molecules illustrated that the dynamic properties could be
659 + preserved by using this new algorithm as an implicit solvent model.

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines