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 2867 by tim, Sat Jun 17 18:43:58 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 101 | Line 100 | Theoretically, the friction kernel can be determined u
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}
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   \]
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 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 273 | 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
278 <
276 > where $\delta _{ij}$ is Kronecker delta function. Inverting the $B$
277 > matrix, we obtain
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}
300 < \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\
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  \\
303 < \end{array}
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
306 < resistance (reaction), at which the trace of rotational resistance
307 < tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
308 < resistance is defined as an unique point of the rigid body at which
309 < the translation-rotation coupling tensor are symmetric,
306 > resistance (or center of reaction), at which the trace of rotational
307 > resistance tensor, $ \Xi ^{rr}$ reaches a minimum value.
308 > Mathematically, the center of resistance is defined as an unique
309 > point of the rigid body at which the translation-rotation coupling
310 > tensor are symmetric,
311   \begin{equation}
312   \Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
313   \label{introEquation:definitionCR}
314   \end{equation}
315 < Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
315 > From Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
316   we can easily find out that the translational resistance tensor is
317   origin independent, while the rotational resistance tensor and
318   translation-rotation coupling resistance tensor depend on the
319 < origin. Given resistance tensor at an arbitrary origin $O$, and a
320 < vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
319 > origin. Given the resistance tensor at an arbitrary origin $O$, and
320 > a vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
321   obtain the resistance tensor at $P$ by
322   \begin{equation}
323   \begin{array}{l}
# Line 337 | 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 356 | Line 354 | the position of center of resistance,
354   (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
355   \end{array} \right) \\
356   \end{eqnarray*}
359 –
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 369 | 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 380 | 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 405 | 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}
408 –
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 427 | 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}
430 –
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 447 | 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*}
450 –
445   In this context, the $\mathrm{rotate}$ function is the reversible
446   product of the three body-fixed rotations,
447   \begin{equation}
# Line 482 | 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.
480 <
481 < After the first part of the propagation, the forces and body-fixed
488 < torques are calculated at the new positions and orientations
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  
483   {\tt doForces:}
484   \begin{align*}
# Line 498 | 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  
502 – {\sc oopse} automatically updates ${\bf u}$ when the rotation matrix
503 – $\mathsf{A}$ is calculated in {\tt moveA}.  Once the forces and
504 – torques have been obtained at the new time step, the velocities can
505 – be advanced to the same time value.
506 –
497   {\tt moveB:}
498   \begin{align*}
499   {\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2
# Line 540 | 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 577 | 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 <
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 > 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
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 & 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 < \begin{figure}
610 < \centering
611 < \includegraphics[width=\linewidth]{one_banana.eps}
612 < \caption[]{.} \label{langevin:banana}
613 < \end{figure}
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
# Line 602 | Line 629 | molecules and 256 pentane molecules.} \label{langevin:
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