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 2876 by tim, Wed Jun 21 19:46:01 2006 UTC vs.
Revision 2905 by tim, Wed Jun 28 19:09:25 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 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}
184 < \],
183 > }}{a}.
184 > \]
185   one can write down the translational and rotational resistance
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   \]
# Line 230 | Line 231 | This equation is the basis for deriving the hydrodynam
231   \label{introEquation:tensorExpression}
232   \end{equation}
233   This equation is the basis for deriving the hydrodynamic tensor. In
234 < 1930, Oseen and Burgers gave a simple solution to Equation
235 < \ref{introEquation:tensorExpression}
234 > 1930, Oseen and Burgers gave a simple solution to
235 > Eq.~\ref{introEquation:tensorExpression}
236   \begin{equation}
237   T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
238   R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor}
# Line 247 | Line 248 | Both of the Equation \ref{introEquation:oseenTensor} a
248   \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
249   \label{introEquation:RPTensorNonOverlapped}
250   \end{equation}
251 < Both of the Equation \ref{introEquation:oseenTensor} and Equation
252 < \ref{introEquation:RPTensorNonOverlapped} have an assumption $R_{ij}
253 < \ge \sigma _i  + \sigma _j$. An alternative expression for
251 > Both of the Eq.~\ref{introEquation:oseenTensor} and
252 > Eq.~\ref{introEquation:RPTensorNonOverlapped} have an assumption
253 > $R_{ij} \ge \sigma _i  + \sigma _j$. An alternative expression for
254   overlapping beads with the same radius, $\sigma$, is given by
255   \begin{equation}
256   T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
# Line 257 | Line 258 | T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left(
258   \frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
259   \label{introEquation:RPTensorOverlapped}
260   \end{equation}
260 –
261   To calculate the resistance tensor at an arbitrary origin $O$, we
262   construct a $3N \times 3N$ matrix consisting of $N \times N$
263   $B_{ij}$ blocks
# Line 275 | Line 275 | matrix, we obtain
275   \]
276   where $\delta _{ij}$ is Kronecker delta function. Inverting the $B$
277   matrix, we obtain
278 –
278   \[
279   C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
280     {C_{11} } &  \ldots  & {C_{1N} }  \\
281      \vdots  &  \ddots  &  \vdots   \\
282     {C_{N1} } &  \cdots  & {C_{NN} }  \\
283 < \end{array}} \right)
283 > \end{array}} \right),
284   \]
285 < , which can be partitioned into $N \times N$ $3 \times 3$ block
285 > which can be partitioned into $N \times N$ $3 \times 3$ block
286   $C_{ij}$. With the help of $C_{ij}$ and skew matrix $U_i$
287   \[
288   U_i  = \left( {\begin{array}{*{20}c}
# Line 295 | Line 294 | arbitrary origin $O$ can be written as
294   where $x_i$, $y_i$, $z_i$ are the components of the vector joining
295   bead $i$ and origin $O$. Hence, the elements of resistance tensor at
296   arbitrary origin $O$ can be written as
297 < \begin{equation}
298 < \begin{array}{l}
299 < \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\
300 < \Xi _{}^{tr}  = \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
302 < \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j  \\
303 < \end{array}
297 > \begin{eqnarray}
298 > \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\
299 > \Xi _{}^{tr}  = \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
300 > \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j. \notag \\
301   \label{introEquation:ResistanceTensorArbitraryOrigin}
302 < \end{equation}
302 > \end{eqnarray}
303  
304   The resistance tensor depends on the origin to which they refer. The
305   proper location for applying friction force is the center of
# Line 338 | Line 335 | Using Equations \ref{introEquation:definitionCR} and
335     { - y_{OP} } & {x_{OP} } & 0  \\
336   \end{array}} \right)
337   \]
338 < Using Equations \ref{introEquation:definitionCR} and
339 < \ref{introEquation:resistanceTensorTransformation}, one can locate
340 < the position of center of resistance,
338 > Using Eq.~\ref{introEquation:definitionCR} and
339 > Eq.~\ref{introEquation:resistanceTensorTransformation}, one can
340 > locate the position of center of resistance,
341   \begin{eqnarray*}
342   \left( \begin{array}{l}
343   x_{OR}  \\
# Line 357 | Line 354 | the position of center of resistance,
354   (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
355   \end{array} \right) \\
356   \end{eqnarray*}
360 –
357   where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
358   joining center of resistance $R$ and origin $O$.
359  
# Line 370 | Line 366 | $V_i = V_i(v_i,\omega _i)$. The right side of Eq.~
366   \end{equation}
367   where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
368   and moment of inertial) matrix and $V_i$ is a generalized velocity,
369 < $V_i = V_i(v_i,\omega _i)$. The right side of Eq.~
370 < (\ref{LDGeneralizedForm}) consists of three generalized forces in
369 > $V_i = V_i(v_i,\omega _i)$. The right side of
370 > Eq.~\ref{LDGeneralizedForm} consists of three generalized forces in
371   lab-fixed frame, systematic force $F_{s,i}$, dissipative force
372   $F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
373   system in Newtownian mechanics typically refers to lab-fixed frame,
# Line 381 | Line 377 | in body-fixed frame and converted back to lab-fixed fr
377   \[
378   \begin{array}{l}
379   F_{f,i}^l (t) = A^T F_{f,i}^b (t), \\
380 < F_{r,i}^l (t) = A^T F_{r,i}^b (t) \\
381 < \end{array}.
380 > F_{r,i}^l (t) = A^T F_{r,i}^b (t). \\
381 > \end{array}
382   \]
383   Here, the body-fixed friction force $F_{r,i}^b$ is proportional to
384   the body-fixed velocity at center of resistance $v_{R,i}^b$ and
385 < angular velocity $\omega _i$,
385 > angular velocity $\omega _i$
386   \begin{equation}
387   F_{r,i}^b (t) = \left( \begin{array}{l}
388   f_{r,i}^b (t) \\
# Line 406 | Line 402 | with zero mean and variance
402   \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
403   2k_B T\Xi _R \delta (t - t'). \label{randomForce}
404   \end{equation}
409 –
405   The equation of motion for $v_i$ can be written as
406   \begin{equation}
407   m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
# Line 428 | Line 423 | + \tau _{r,i}^b(t)
423   \dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t)
424   + \tau _{r,i}^b(t)
425   \end{equation}
431 –
426   Embedding the friction terms into force and torque, one can
427   integrate the langevin equations of motion for rigid body of
428   arbitrary shape in a velocity-Verlet style 2-part algorithm, where
# Line 448 | Line 442 | $h= \delta t$:
442   \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
443      (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
444   \end{align*}
451 –
445   In this context, the $\mathrm{rotate}$ function is the reversible
446   product of the three body-fixed rotations,
447   \begin{equation}
# Line 483 | Line 476 | All other rotations follow in a straightforward manner
476   \end{array}
477   \right).
478   \end{equation}
479 < All other rotations follow in a straightforward manner.
479 > All other rotations follow in a straightforward manner. After the
480 > first part of the propagation, the forces and body-fixed torques are
481 > calculated at the new positions and orientations
482  
488 – After the first part of the propagation, the forces and body-fixed
489 – torques are calculated at the new positions and orientations
490 –
483   {\tt doForces:}
484   \begin{align*}
485   {\bf f}(t + h) &\leftarrow
# Line 499 | Line 491 | torques are calculated at the new positions and orient
491   {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
492      \cdot {\bf \tau}^s(t + h).
493   \end{align*}
494 + Once the forces and torques have been obtained at the new time step,
495 + the velocities can be advanced to the same time value.
496  
503 – {\sc oopse} automatically updates ${\bf u}$ when the rotation matrix
504 – $\mathsf{A}$ is calculated in {\tt moveA}.  Once the forces and
505 – torques have been obtained at the new time step, the velocities can
506 – be advanced to the same time value.
507 –
497   {\tt moveB:}
498   \begin{align*}
499   {\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2
# Line 541 | Line 530 | rate at lower viscosity regime.
530   explained by the finite size effect and relative slow relaxation
531   rate at lower viscosity regime.
532   \begin{table}
533 < \caption{Average temperatures from Langevin dynamics simulations of
534 < SSD water molecules with reference temperature at 300~K.}
533 > \caption{AVERAGE TEMPERATURES FROM LANGEVIN DYNAMICS SIMULATIONS OF
534 > SSD WATER MOLECULES WITH REFERENCE TEMPERATURE AT 300~K.}
535   \label{langevin:viscosity}
536   \begin{center}
537 < \begin{tabular}{|l|l|l|}
537 > \begin{tabular}{lll}
538    \hline
539    $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\
540 +  \hline
541    1    & 300.47 & 10.99 \\
542    0.1  & 301.19 & 11.136 \\
543    0.01 & 303.04 & 11.796 \\
# Line 578 | Line 568 | temperature fluctuation versus time.} \label{langevin:
568   temperature fluctuation versus time.} \label{langevin:temperature}
569   \end{figure}
570  
571 < \subsection{Langevin Dynamics of Banana Shaped Molecule}
571 > \subsection{Langevin Dynamics of Banana Shaped Molecules}
572  
573   In order to verify that Langevin dynamics can mimic the dynamics of
574   the systems absent of explicit solvents, we carried out two sets of
575   simulations and compare their dynamic properties.
576 <
577 < \subsubsection{Simulations Contain One Banana Shaped Molecule}
578 <
579 < 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
576 > Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation
577 > made of 256 pentane molecules and two banana shaped molecules at
578 > 273~K. It has an equivalent implicit solvent system containing only
579 > two banana shaped molecules with viscosity of 0.289 center poise. To
580   calculate the hydrodynamic properties of the banana shaped molecule,
581   we create a rough shell model (see Fig.~\ref{langevin:roughShell}),
582   in which the banana shaped molecule is represented as a ``shell''
583 < made of 2266 small identical beads with size of 0.3 $\AA$ on the
583 > made of 2266 small identical beads with size of 0.3 \AA on the
584   surface. Applying the procedure described in
585   Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
586   identified the center of resistance at $(0, 0.7482, -0.1988)$, as
587   well as the resistance tensor,
588   \[
589   \left( {\begin{array}{*{20}c}
590 < 0.9261 & 1.310e-14 & -7.292e-15&5.067e-14&0.08585&0.2057\\
591 < 3.968e-14& 0.9270&-0.007063& 0.08585&6.764e-14&4.846e-14\\
592 < -6.561e-16&-0.007063&0.7494&0.2057&4.846e-14&1.5036e-14\\
593 < 5.067e-14&0.0858&0.2057& 58.64& 8.563e-13&-8.5736\\
594 < 0.08585&6.764e-14&4.846e-14&1.555e-12&48.30&3.219&\\
595 < 0.2057&4.846e-14&1.5036e-14&-3.904e-13&3.219&10.7373\\
590 > 0.9261 & 0 & 0&0&0.08585&0.2057\\
591 > 0& 0.9270&-0.007063& 0.08585&0&0\\
592 > 0&0.007063&0.7494&0.2057&0&0\\
593 > 0&0.0858&0.2057& 58.64& 0&-8.5736\\
594 > 0.08585&0&0&0&48.30&3.219&\\
595 > 0.2057&0&0&0&3.219&10.7373\\
596   \end{array}} \right).
597   \]
598 + %\[
599 + %\left( {\begin{array}{*{20}c}
600 + %0.9261 & 1.310e-14 & -7.292e-15&5.067e-14&0.08585&0.2057\\
601 + %3.968e-14& 0.9270&-0.007063& 0.08585&6.764e-14&4.846e-14\\
602 + %-6.561e-16&-0.007063&0.7494&0.2057&4.846e-14&1.5036e-14\\
603 + %5.067e-14&0.0858&0.2057& 58.64& 8.563e-13&-8.5736\\
604 + %0.08585&6.764e-14&4.846e-14&1.555e-12&48.30&3.219&\\
605 + %0.2057&4.846e-14&1.5036e-14&-3.904e-13&3.219&10.7373\\
606 + %\end{array}} \right).
607 + %\]
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]{roughShell.eps}
# Line 618 | Line 623 | model for banana shaped molecule.} \label{langevin:rou
623  
624   \begin{figure}
625   \centering
626 < \includegraphics[width=\linewidth]{oneBanana.eps}
627 < \caption[Snapshot from Simulation of One Banana Shaped Molecules and
628 < 256 Pentane Molecules]{Snapshot from simulation of one Banana shaped
629 < molecules and 256 pentane molecules.} \label{langevin:oneBanana}
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 < \subsubsection{Simulations Contain Two Banana Shaped Molecules}
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]{twoBanana.eps}
643 < \caption[Snapshot from Simulation of Two Banana Shaped Molecules and
644 < 256 Pentane Molecules]{Snapshot from simulation of two Banana shaped
645 < molecules and 256 pentane molecules.} \label{langevin:twoBanana}
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