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

Comparing trunk/tengDissertation/Introduction.tex (file contents):
Revision 2850 by tim, Sun Jun 11 01:55:48 2006 UTC vs.
Revision 2904 by tim, Wed Jun 28 17:36:32 2006 UTC

# Line 31 | Line 31 | Newton's third law states that
31   $F_{ji}$ be the force that particle $j$ exerts on particle $i$.
32   Newton's third law states that
33   \begin{equation}
34 < F_{ij} = -F_{ji}
34 > F_{ij} = -F_{ji}.
35   \label{introEquation:newtonThirdLaw}
36   \end{equation}
37 –
37   Conservation laws of Newtonian Mechanics play very important roles
38   in solving mechanics problems. The linear momentum of a particle is
39   conserved if it is free or it experiences no force. The second
# Line 63 | Line 62 | momentum of it is conserved. The last conservation the
62   \end{equation}
63   If there are no external torques acting on a body, the angular
64   momentum of it is conserved. The last conservation theorem state
65 < that if all forces are conservative, Energy
66 < \begin{equation}E = T + V \label{introEquation:energyConservation}
65 > that if all forces are conservative, energy is conserved,
66 > \begin{equation}E = T + V. \label{introEquation:energyConservation}
67   \end{equation}
68 < is conserved. All of these conserved quantities are
69 < important factors to determine the quality of numerical integration
70 < schemes for rigid bodies \cite{Dullweber1997}.
68 > All of these conserved quantities are important factors to determine
69 > the quality of numerical integration schemes for rigid bodies
70 > \cite{Dullweber1997}.
71  
72   \subsection{\label{introSection:lagrangian}Lagrangian Mechanics}
73  
74 < Newtonian Mechanics suffers from two important limitations: motions
75 < can only be described in cartesian coordinate systems. Moreover, It
76 < become impossible to predict analytically the properties of the
77 < system even if we know all of the details of the interaction. In
78 < order to overcome some of the practical difficulties which arise in
79 < attempts to apply Newton's equation to complex system, approximate
80 < numerical procedures may be developed.
74 > Newtonian Mechanics suffers from a important limitation: motions can
75 > only be described in cartesian coordinate systems which make it
76 > impossible to predict analytically the properties of the system even
77 > if we know all of the details of the interaction. In order to
78 > overcome some of the practical difficulties which arise in attempts
79 > to apply Newton's equation to complex system, approximate numerical
80 > procedures may be developed.
81  
82   \subsubsection{\label{introSection:halmiltonPrinciple}\textbf{Hamilton's
83   Principle}}
84  
85   Hamilton introduced the dynamical principle upon which it is
86   possible to base all of mechanics and most of classical physics.
87 < Hamilton's Principle may be stated as follows,
88 <
89 < The actual trajectory, along which a dynamical system may move from
90 < one point to another within a specified time, is derived by finding
91 < the path which minimizes the time integral of the difference between
93 < the kinetic, $K$, and potential energies, $U$.
87 > Hamilton's Principle may be stated as follows: the actual
88 > trajectory, along which a dynamical system may move from one point
89 > to another within a specified time, is derived by finding the path
90 > which minimizes the time integral of the difference between the
91 > kinetic, $K$, and potential energies, $U$,
92   \begin{equation}
93 < \delta \int_{t_1 }^{t_2 } {(K - U)dt = 0} ,
93 > \delta \int_{t_1 }^{t_2 } {(K - U)dt = 0}.
94   \label{introEquation:halmitonianPrinciple1}
95   \end{equation}
98 –
96   For simple mechanical systems, where the forces acting on the
97   different parts are derivable from a potential, the Lagrangian
98   function $L$ can be defined as the difference between the kinetic
99   energy of the system and its potential energy,
100   \begin{equation}
101 < L \equiv K - U = L(q_i ,\dot q_i ) ,
101 > L \equiv K - U = L(q_i ,\dot q_i ).
102   \label{introEquation:lagrangianDef}
103   \end{equation}
104 < then Eq.~\ref{introEquation:halmitonianPrinciple1} becomes
104 > Thus, Eq.~\ref{introEquation:halmitonianPrinciple1} becomes
105   \begin{equation}
106 < \delta \int_{t_1 }^{t_2 } {L dt = 0} ,
106 > \delta \int_{t_1 }^{t_2 } {L dt = 0} .
107   \label{introEquation:halmitonianPrinciple2}
108   \end{equation}
109  
# Line 138 | Line 135 | p_i  = \frac{{\partial L}}{{\partial q_i }}
135   p_i  = \frac{{\partial L}}{{\partial q_i }}
136   \label{introEquation:generalizedMomentaDot}
137   \end{equation}
141 –
138   With the help of the generalized momenta, we may now define a new
139   quantity $H$ by the equation
140   \begin{equation}
# Line 146 | Line 142 | where $ \dot q_1  \ldots \dot q_f $ are generalized ve
142   \label{introEquation:hamiltonianDefByLagrangian}
143   \end{equation}
144   where $ \dot q_1  \ldots \dot q_f $ are generalized velocities and
145 < $L$ is the Lagrangian function for the system.
146 <
151 < Differentiating Eq.~\ref{introEquation:hamiltonianDefByLagrangian},
152 < one can obtain
145 > $L$ is the Lagrangian function for the system. Differentiating
146 > Eq.~\ref{introEquation:hamiltonianDefByLagrangian}, one can obtain
147   \begin{equation}
148   dH = \sum\limits_k {\left( {p_k d\dot q_k  + \dot q_k dp_k  -
149   \frac{{\partial L}}{{\partial q_k }}dq_k  - \frac{{\partial
150   L}}{{\partial \dot q_k }}d\dot q_k } \right)}  - \frac{{\partial
151 < L}}{{\partial t}}dt \label{introEquation:diffHamiltonian1}
151 > L}}{{\partial t}}dt . \label{introEquation:diffHamiltonian1}
152   \end{equation}
153 < Making use of  Eq.~\ref{introEquation:generalizedMomenta}, the
154 < second and fourth terms in the parentheses cancel. Therefore,
153 > Making use of Eq.~\ref{introEquation:generalizedMomenta}, the second
154 > and fourth terms in the parentheses cancel. Therefore,
155   Eq.~\ref{introEquation:diffHamiltonian1} can be rewritten as
156   \begin{equation}
157   dH = \sum\limits_k {\left( {\dot q_k dp_k  - \dot p_k dq_k }
158 < \right)}  - \frac{{\partial L}}{{\partial t}}dt
158 > \right)}  - \frac{{\partial L}}{{\partial t}}dt .
159   \label{introEquation:diffHamiltonian2}
160   \end{equation}
161   By identifying the coefficients of $dq_k$, $dp_k$ and dt, we can
# Line 180 | Line 174 | t}}
174   t}}
175   \label{introEquation:motionHamiltonianTime}
176   \end{equation}
177 <
184 < Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
177 > where Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
178   Eq.~\ref{introEquation:motionHamiltonianMomentum} are Hamilton's
179   equation of motion. Due to their symmetrical formula, they are also
180   known as the canonical equations of motions \cite{Goldstein2001}.
# Line 195 | Line 188 | only works with 1st-order differential equations\cite{
188   statistical mechanics and quantum mechanics, since it treats the
189   coordinate and its time derivative as independent variables and it
190   only works with 1st-order differential equations\cite{Marion1990}.
198 –
191   In Newtonian Mechanics, a system described by conservative forces
192 < conserves the total energy \ref{introEquation:energyConservation}.
193 < It follows that Hamilton's equations of motion conserve the total
194 < Hamiltonian.
192 > conserves the total energy
193 > (Eq.~\ref{introEquation:energyConservation}). It follows that
194 > Hamilton's equations of motion conserve the total Hamiltonian
195   \begin{equation}
196   \frac{{dH}}{{dt}} = \sum\limits_i {\left( {\frac{{\partial
197   H}}{{\partial q_i }}\dot q_i  + \frac{{\partial H}}{{\partial p_i
198   }}\dot p_i } \right)}  = \sum\limits_i {\left( {\frac{{\partial
199   H}}{{\partial q_i }}\frac{{\partial H}}{{\partial p_i }} -
200   \frac{{\partial H}}{{\partial p_i }}\frac{{\partial H}}{{\partial
201 < q_i }}} \right) = 0} \label{introEquation:conserveHalmitonian}
201 > q_i }}} \right) = 0}. \label{introEquation:conserveHalmitonian}
202   \end{equation}
203  
204   \section{\label{introSection:statisticalMechanics}Statistical
# Line 227 | Line 219 | assigned to one of $6f$ mutually orthogonal axes, the
219   momentum variables. Consider a dynamic system of $f$ particles in a
220   cartesian space, where each of the $6f$ coordinates and momenta is
221   assigned to one of $6f$ mutually orthogonal axes, the phase space of
222 < this system is a $6f$ dimensional space. A point, $x = (q_1 , \ldots
223 < ,q_f ,p_1 , \ldots ,p_f )$, with a unique set of values of $6f$
224 < coordinates and momenta is a phase space vector.
225 <
222 > this system is a $6f$ dimensional space. A point, $x =
223 > (\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
224 > \over q} _1 , \ldots
225 > ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
226 > \over q} _f
227 > ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
228 > \over p} _1  \ldots
229 > ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
230 > \over p} _f )$ , with a unique set of values of $6f$ coordinates and
231 > momenta is a phase space vector.
232   %%%fix me
233 < A microscopic state or microstate of a classical system is
234 < specification of the complete phase space vector of a system at any
237 < instant in time. An ensemble is defined as a collection of systems
238 < sharing one or more macroscopic characteristics but each being in a
239 < unique microstate. The complete ensemble is specified by giving all
240 < systems or microstates consistent with the common macroscopic
241 < characteristics of the ensemble. Although the state of each
242 < individual system in the ensemble could be precisely described at
243 < any instance in time by a suitable phase space vector, when using
244 < ensembles for statistical purposes, there is no need to maintain
245 < distinctions between individual systems, since the numbers of
246 < systems at any time in the different states which correspond to
247 < different regions of the phase space are more interesting. Moreover,
248 < in the point of view of statistical mechanics, one would prefer to
249 < use ensembles containing a large enough population of separate
250 < members so that the numbers of systems in such different states can
251 < be regarded as changing continuously as we traverse different
252 < regions of the phase space. The condition of an ensemble at any time
233 >
234 > In statistical mechanics, the condition of an ensemble at any time
235   can be regarded as appropriately specified by the density $\rho$
236   with which representative points are distributed over the phase
237   space. The density distribution for an ensemble with $f$ degrees of
# Line 304 | Line 286 | can be used to describe the statistical properties of
286   statistical characteristics. As a function of macroscopic
287   parameters, such as temperature \textit{etc}, the partition function
288   can be used to describe the statistical properties of a system in
289 < thermodynamic equilibrium.
290 <
291 < As an ensemble of systems, each of which is known to be thermally
310 < isolated and conserve energy, the Microcanonical ensemble (NVE) has
311 < a partition function like,
289 > thermodynamic equilibrium. As an ensemble of systems, each of which
290 > is known to be thermally isolated and conserve energy, the
291 > Microcanonical ensemble (NVE) has a partition function like,
292   \begin{equation}
293 < \Omega (N,V,E) = e^{\beta TS} \label{introEquation:NVEPartition}.
293 > \Omega (N,V,E) = e^{\beta TS}. \label{introEquation:NVEPartition}
294   \end{equation}
295   A canonical ensemble (NVT)is an ensemble of systems, each of which
296   can share its energy with a large heat reservoir. The distribution
297   of the total energy amongst the possible dynamical states is given
298   by the partition function,
299   \begin{equation}
300 < \Omega (N,V,T) = e^{ - \beta A}
300 > \Omega (N,V,T) = e^{ - \beta A}.
301   \label{introEquation:NVTPartition}
302   \end{equation}
303   Here, $A$ is the Helmholtz free energy which is defined as $ A = U -
# Line 374 | Line 354 | simple form,
354   \frac{{\partial \rho }}{{\partial p_i }}\dot p_i } \right)}  = 0 .
355   \label{introEquation:liouvilleTheorem}
356   \end{equation}
377 –
357   Liouville's theorem states that the distribution function is
358   constant along any trajectory in phase space. In classical
359   statistical mechanics, since the number of members in an ensemble is
# Line 496 | Line 475 | developed to address this issue\cite{Dullweber1997, Mc
475   geometric integrators, which preserve various phase-flow invariants
476   such as symplectic structure, volume and time reversal symmetry, are
477   developed to address this issue\cite{Dullweber1997, McLachlan1998,
478 < Leimkuhler1999}. The velocity verlet method, which happens to be a
478 > Leimkuhler1999}. The velocity Verlet method, which happens to be a
479   simple example of symplectic integrator, continues to gain
480   popularity in the molecular dynamics community. This fact can be
481   partly explained by its geometric nature.
# Line 516 | Line 495 | $\omega(x, x) = 0$. The cross product operation in vec
495   $\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+
496   \lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and
497   $\omega(x, x) = 0$. The cross product operation in vector field is
498 < an example of symplectic form.
499 <
500 < One of the motivations to study \emph{symplectic manifolds} in
501 < Hamiltonian Mechanics is that a symplectic manifold can represent
502 < all possible configurations of the system and the phase space of the
503 < system can be described by it's cotangent bundle. Every symplectic
504 < manifold is even dimensional. For instance, in Hamilton equations,
526 < coordinate and momentum always appear in pairs.
498 > an example of symplectic form. One of the motivations to study
499 > \emph{symplectic manifolds} in Hamiltonian Mechanics is that a
500 > symplectic manifold can represent all possible configurations of the
501 > system and the phase space of the system can be described by it's
502 > cotangent bundle. Every symplectic manifold is even dimensional. For
503 > instance, in Hamilton equations, coordinate and momentum always
504 > appear in pairs.
505  
506   \subsection{\label{introSection:ODE}Ordinary Differential Equations}
507  
# Line 550 | Line 528 | system can be rewritten as,
528   \frac{d}{{dt}}x = J\nabla _x H(x)
529   \label{introEquation:compactHamiltonian}
530   \end{equation}In this case, $f$ is
531 < called a \emph{Hamiltonian vector field}.
532 <
555 < Another generalization of Hamiltonian dynamics is Poisson
556 < Dynamics\cite{Olver1986},
531 > called a \emph{Hamiltonian vector field}. Another generalization of
532 > Hamiltonian dynamics is Poisson Dynamics\cite{Olver1986},
533   \begin{equation}
534   \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
535   \end{equation}
# Line 591 | Line 567 | In most cases, it is not easy to find the exact flow $
567   \end{equation}
568  
569   In most cases, it is not easy to find the exact flow $\varphi_\tau$.
570 < Instead, we use a approximate map, $\psi_\tau$, which is usually
570 > Instead, we use an approximate map, $\psi_\tau$, which is usually
571   called integrator. The order of an integrator $\psi_\tau$ is $p$, if
572   the Taylor series of $\psi_\tau$ agree to order $p$,
573   \begin{equation}
574 < \psi_tau(x) = x + \tau f(x) + O(\tau^{p+1})
574 > \psi_\tau(x) = x + \tau f(x) + O(\tau^{p+1})
575   \end{equation}
576  
577   \subsection{\label{introSection:geometricProperties}Geometric Properties}
578  
579 < The hidden geometric properties\cite{Budd1999, Marsden1998} of ODE
580 < and its flow play important roles in numerical studies. Many of them
581 < can be found in systems which occur naturally in applications.
606 <
579 > The hidden geometric properties\cite{Budd1999, Marsden1998} of an
580 > ODE and its flow play important roles in numerical studies. Many of
581 > them can be found in systems which occur naturally in applications.
582   Let $\varphi$ be the flow of Hamiltonian vector field, $\varphi$ is
583   a \emph{symplectic} flow if it satisfies,
584   \begin{equation}
# Line 617 | Line 592 | symplectomorphism. As to the Poisson system,
592   \begin{equation}
593   {\varphi '}^T J \varphi ' = J \circ \varphi
594   \end{equation}
595 < is the property must be preserved by the integrator.
596 <
597 < It is possible to construct a \emph{volume-preserving} flow for a
598 < source free($ \nabla \cdot f = 0 $) ODE, if the flow satisfies $
599 < \det d\varphi  = 1$. One can show easily that a symplectic flow will
600 < be volume-preserving.
626 <
627 < Changing the variables $y = h(x)$ in a ODE\ref{introEquation:ODE}
628 < will result in a new system,
595 > is the property that must be preserved by the integrator. It is
596 > possible to construct a \emph{volume-preserving} flow for a source
597 > free ODE ($ \nabla \cdot f = 0 $), if the flow satisfies $ \det
598 > d\varphi  = 1$. One can show easily that a symplectic flow will be
599 > volume-preserving. Changing the variables $y = h(x)$ in an ODE
600 > (Eq.~\ref{introEquation:ODE}) will result in a new system,
601   \[
602   \dot y = \tilde f(y) = ((dh \cdot f)h^{ - 1} )(y).
603   \]
604   The vector filed $f$ has reversing symmetry $h$ if $f = - \tilde f$.
605   In other words, the flow of this vector field is reversible if and
606 < only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $.
607 <
636 < A \emph{first integral}, or conserved quantity of a general
606 > only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $. A
607 > \emph{first integral}, or conserved quantity of a general
608   differential function is a function $ G:R^{2d}  \to R^d $ which is
609   constant for all solutions of the ODE $\frac{{dx}}{{dt}} = f(x)$ ,
610   \[
# Line 646 | Line 617 | smooth function $G$ is given by,
617   which is the condition for conserving \emph{first integral}. For a
618   canonical Hamiltonian system, the time evolution of an arbitrary
619   smooth function $G$ is given by,
649 –
620   \begin{eqnarray}
621   \frac{{dG(x(t))}}{{dt}} & = & [\nabla _x G(x(t))]^T \dot x(t) \\
622                          & = & [\nabla _x G(x(t))]^T J\nabla _x H(x(t)). \\
623   \label{introEquation:firstIntegral1}
624   \end{eqnarray}
655 –
656 –
625   Using poisson bracket notion, Equation
626   \ref{introEquation:firstIntegral1} can be rewritten as
627   \[
# Line 666 | Line 634 | is a \emph{first integral}, which is due to the fact $
634   \]
635   As well known, the Hamiltonian (or energy) H of a Hamiltonian system
636   is a \emph{first integral}, which is due to the fact $\{ H,H\}  =
637 < 0$.
670 <
671 < When designing any numerical methods, one should always try to
637 > 0$. When designing any numerical methods, one should always try to
638   preserve the structural properties of the original ODE and its flow.
639  
640   \subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
641   A lot of well established and very effective numerical methods have
642   been successful precisely because of their symplecticities even
643   though this fact was not recognized when they were first
644 < constructed. The most famous example is the Verlet-leapfrog methods
644 > constructed. The most famous example is the Verlet-leapfrog method
645   in molecular dynamics. In general, symplectic integrators can be
646   constructed using one of four different methods.
647   \begin{enumerate}
# Line 709 | Line 675 | where each of the sub-flow is chosen such that each re
675   \label{introEquation:FlowDecomposition}
676   \end{equation}
677   where each of the sub-flow is chosen such that each represent a
678 < simpler integration of the system.
679 <
714 < Suppose that a Hamiltonian system takes the form,
678 > simpler integration of the system. Suppose that a Hamiltonian system
679 > takes the form,
680   \[
681   H = H_1 + H_2.
682   \]
# Line 752 | Line 717 | to its symmetric property,
717   \begin{equation}
718   \varphi _h^{ - 1} = \varphi _{ - h}.
719   \label{introEquation:timeReversible}
720 < \end{equation},appendixFig:architecture
720 > \end{equation}
721  
722 < \subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Example of Splitting Method}}
722 > \subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Examples of the Splitting Method}}
723   The classical equation for a system consisting of interacting
724   particles can be written in Hamiltonian form,
725   \[
726   H = T + V
727   \]
728   where $T$ is the kinetic energy and $V$ is the potential energy.
729 < Setting $H_1 = T, H_2 = V$ and applying Strang splitting, one
729 > Setting $H_1 = T, H_2 = V$ and applying the Strang splitting, one
730   obtains the following:
731   \begin{align}
732   q(\Delta t) &= q(0) + \dot{q}(0)\Delta t +
# Line 788 | Line 753 | q(\Delta t) &= q(0) + \Delta t\, \dot{q}\biggl (\frac{
753      \label{introEquation:Lp9b}\\%
754   %
755   \dot{q}(\Delta t) &= \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) +
756 <    \frac{\Delta t}{2m}\, F[q(0)]. \label{introEquation:Lp9c}
756 >    \frac{\Delta t}{2m}\, F[q(t)]. \label{introEquation:Lp9c}
757   \end{align}
758   From the preceding splitting, one can see that the integration of
759   the equations of motion would follow:
# Line 797 | Line 762 | the equations of motion would follow:
762  
763   \item Use the half step velocities to move positions one whole step, $\Delta t$.
764  
765 < \item Evaluate the forces at the new positions, $\mathbf{r}(\Delta t)$, and use the new forces to complete the velocity move.
765 > \item Evaluate the forces at the new positions, $\mathbf{q}(\Delta t)$, and use the new forces to complete the velocity move.
766  
767   \item Repeat from step 1 with the new position, velocities, and forces assuming the roles of the initial values.
768   \end{enumerate}
769 <
770 < Simply switching the order of splitting and composing, a new
771 < integrator, the \emph{position verlet} integrator, can be generated,
769 > By simply switching the order of the propagators in the splitting
770 > and composing a new integrator, the \emph{position verlet}
771 > integrator, can be generated,
772   \begin{align}
773   \dot q(\Delta t) &= \dot q(0) + \Delta tF(q(0))\left[ {q(0) +
774   \frac{{\Delta t}}{{2m}}\dot q(0)} \right], %
# Line 816 | Line 781 | q(\Delta t)} \right]. %
781  
782   \subsubsection{\label{introSection:errorAnalysis}\textbf{Error Analysis and Higher Order Methods}}
783  
784 < Baker-Campbell-Hausdorff formula can be used to determine the local
785 < error of splitting method in terms of commutator of the
784 > The Baker-Campbell-Hausdorff formula can be used to determine the
785 > local error of splitting method in terms of the commutator of the
786   operators(\ref{introEquation:exponentialOperator}) associated with
787 < the sub-flow. For operators $hX$ and $hY$ which are associate to
787 > the sub-flow. For operators $hX$ and $hY$ which are associated with
788   $\varphi_1(t)$ and $\varphi_2(t)$ respectively , we have
789   \begin{equation}
790   \exp (hX + hY) = \exp (hZ)
# Line 833 | Line 798 | Here, $[X,Y]$ is the commutators of operator $X$ and $
798   \[
799   [X,Y] = XY - YX .
800   \]
801 < Applying Baker-Campbell-Hausdorff formula\cite{Varadarajan1974} to
802 < Sprang splitting, we can obtain
801 > Applying the Baker-Campbell-Hausdorff formula\cite{Varadarajan1974}
802 > to the Sprang splitting, we can obtain
803   \begin{eqnarray*}
804   \exp (h X/2)\exp (h Y)\exp (h X/2) & = & \exp (h X + h Y + h^2 [X,Y]/4 + h^2 [Y,X]/4 \\
805                                     &   & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\
806                                     &   & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3[X,[X,Y]]/24 + \ldots )
807   \end{eqnarray*}
808 < Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0\], the dominant local
808 > Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0,\] the dominant local
809   error of Spring splitting is proportional to $h^3$. The same
810 < procedure can be applied to general splitting,  of the form
810 > procedure can be applied to a general splitting,  of the form
811   \begin{equation}
812   \varphi _{b_m h}^2  \circ \varphi _{a_m h}^1  \circ \varphi _{b_{m -
813   1} h}^2  \circ  \ldots  \circ \varphi _{a_1 h}^1 .
814   \end{equation}
815 < Careful choice of coefficient $a_1 \ldots b_m$ will lead to higher
816 < order method. Yoshida proposed an elegant way to compose higher
815 > A careful choice of coefficient $a_1 \ldots b_m$ will lead to higher
816 > order methods. Yoshida proposed an elegant way to compose higher
817   order methods based on symmetric splitting\cite{Yoshida1990}. Given
818   a symmetric second order base method $ \varphi _h^{(2)} $, a
819   fourth-order symmetric method can be constructed by composing,
# Line 861 | Line 826 | integrator $ \varphi _h^{(2n + 2)}$ can be composed by
826   integrator $ \varphi _h^{(2n + 2)}$ can be composed by
827   \begin{equation}
828   \varphi _h^{(2n + 2)}  = \varphi _{\alpha h}^{(2n)}  \circ \varphi
829 < _{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)}
829 > _{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)},
830   \end{equation}
831 < , if the weights are chosen as
831 > if the weights are chosen as
832   \[
833   \alpha  =  - \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }},\beta =
834   \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }} .
# Line 901 | Line 866 | initialization of a simulation. Sec.~\ref{introSection
866   These three individual steps will be covered in the following
867   sections. Sec.~\ref{introSec:initialSystemSettings} deals with the
868   initialization of a simulation. Sec.~\ref{introSection:production}
869 < will discusses issues in production run.
869 > will discusse issues in production run.
870   Sec.~\ref{introSection:Analysis} provides the theoretical tools for
871   trajectory analysis.
872  
# Line 914 | Line 879 | year, many more remain unknown due to the difficulties
879   databases, such as RCSB Protein Data Bank \textit{etc}. Although
880   thousands of crystal structures of molecules are discovered every
881   year, many more remain unknown due to the difficulties of
882 < purification and crystallization. Even for the molecule with known
883 < structure, some important information is missing. For example, the
882 > purification and crystallization. Even for molecules with known
883 > structure, some important information is missing. For example, a
884   missing hydrogen atom which acts as donor in hydrogen bonding must
885   be added. Moreover, in order to include electrostatic interaction,
886   one may need to specify the partial charges for individual atoms.
887   Under some circumstances, we may even need to prepare the system in
888 < a special setup. For instance, when studying transport phenomenon in
889 < membrane system, we may prepare the lipids in bilayer structure
890 < instead of placing lipids randomly in solvent, since we are not
891 < interested in self-aggregation and it takes a long time to happen.
888 > a special configuration. For instance, when studying transport
889 > phenomenon in membrane systems, we may prepare the lipids in a
890 > bilayer structure instead of placing lipids randomly in solvent,
891 > since we are not interested in the slow self-aggregation process.
892  
893   \subsubsection{\textbf{Minimization}}
894  
895   It is quite possible that some of molecules in the system from
896 < preliminary preparation may be overlapped with each other. This
897 < close proximity leads to high potential energy which consequently
898 < jeopardizes any molecular dynamics simulations. To remove these
899 < steric overlaps, one typically performs energy minimization to find
900 < a more reasonable conformation. Several energy minimization methods
901 < have been developed to exploit the energy surface and to locate the
902 < local minimum. While converging slowly near the minimum, steepest
903 < descent method is extremely robust when systems are far from
904 < harmonic. Thus, it is often used to refine structure from
905 < crystallographic data. Relied on the gradient or hessian, advanced
906 < methods like conjugate gradient and Newton-Raphson converge rapidly
907 < to a local minimum, while become unstable if the energy surface is
908 < far from quadratic. Another factor must be taken into account, when
896 > preliminary preparation may be overlapping with each other. This
897 > close proximity leads to high initial potential energy which
898 > consequently jeopardizes any molecular dynamics simulations. To
899 > remove these steric overlaps, one typically performs energy
900 > minimization to find a more reasonable conformation. Several energy
901 > minimization methods have been developed to exploit the energy
902 > surface and to locate the local minimum. While converging slowly
903 > near the minimum, steepest descent method is extremely robust when
904 > systems are strongly anharmonic. Thus, it is often used to refine
905 > structure from crystallographic data. Relied on the gradient or
906 > hessian, advanced methods like Newton-Raphson converge rapidly to a
907 > local minimum, but become unstable if the energy surface is far from
908 > quadratic. Another factor that must be taken into account, when
909   choosing energy minimization method, is the size of the system.
910   Steepest descent and conjugate gradient can deal with models of any
911 < size. Because of the limit of computation power to calculate hessian
912 < matrix and insufficient storage capacity to store them, most
913 < Newton-Raphson methods can not be used with very large models.
911 > size. Because of the limits on computer memory to store the hessian
912 > matrix and the computing power needed to diagonalized these
913 > matrices, most Newton-Raphson methods can not be used with very
914 > large systems.
915  
916   \subsubsection{\textbf{Heating}}
917  
918   Typically, Heating is performed by assigning random velocities
919 < according to a Gaussian distribution for a temperature. Beginning at
920 < a lower temperature and gradually increasing the temperature by
921 < assigning greater random velocities, we end up with setting the
922 < temperature of the system to a final temperature at which the
923 < simulation will be conducted. In heating phase, we should also keep
924 < the system from drifting or rotating as a whole. Equivalently, the
925 < net linear momentum and angular momentum of the system should be
926 < shifted to zero.
919 > according to a Maxwell-Boltzman distribution for a desired
920 > temperature. Beginning at a lower temperature and gradually
921 > increasing the temperature by assigning larger random velocities, we
922 > end up with setting the temperature of the system to a final
923 > temperature at which the simulation will be conducted. In heating
924 > phase, we should also keep the system from drifting or rotating as a
925 > whole. To do this, the net linear momentum and angular momentum of
926 > the system is shifted to zero after each resampling from the Maxwell
927 > -Boltzman distribution.
928  
929   \subsubsection{\textbf{Equilibration}}
930  
# Line 973 | Line 940 | way.
940  
941   \subsection{\label{introSection:production}Production}
942  
943 < Production run is the most important step of the simulation, in
943 > The production run is the most important step of the simulation, in
944   which the equilibrated structure is used as a starting point and the
945   motions of the molecules are collected for later analysis. In order
946   to capture the macroscopic properties of the system, the molecular
947 < dynamics simulation must be performed in correct and efficient way.
947 > dynamics simulation must be performed by sampling correctly and
948 > efficiently from the relevant thermodynamic ensemble.
949  
950   The most expensive part of a molecular dynamics simulation is the
951   calculation of non-bonded forces, such as van der Waals force and
952   Coulombic forces \textit{etc}. For a system of $N$ particles, the
953   complexity of the algorithm for pair-wise interactions is $O(N^2 )$,
954   which making large simulations prohibitive in the absence of any
955 < computation saving techniques.
955 > algorithmic tricks.
956  
957 < A natural approach to avoid system size issue is to represent the
957 > A natural approach to avoid system size issues is to represent the
958   bulk behavior by a finite number of the particles. However, this
959 < approach will suffer from the surface effect. To offset this,
960 < \textit{Periodic boundary condition} (see Fig.~\ref{introFig:pbc})
961 < is developed to simulate bulk properties with a relatively small
962 < number of particles. In this method, the simulation box is
963 < replicated throughout space to form an infinite lattice. During the
964 < simulation, when a particle moves in the primary cell, its image in
965 < other cells move in exactly the same direction with exactly the same
966 < orientation. Thus, as a particle leaves the primary cell, one of its
967 < images will enter through the opposite face.
959 > approach will suffer from the surface effect at the edges of the
960 > simulation. To offset this, \textit{Periodic boundary conditions}
961 > (see Fig.~\ref{introFig:pbc}) is developed to simulate bulk
962 > properties with a relatively small number of particles. In this
963 > method, the simulation box is replicated throughout space to form an
964 > infinite lattice. During the simulation, when a particle moves in
965 > the primary cell, its image in other cells move in exactly the same
966 > direction with exactly the same orientation. Thus, as a particle
967 > leaves the primary cell, one of its images will enter through the
968 > opposite face.
969   \begin{figure}
970   \centering
971   \includegraphics[width=\linewidth]{pbc.eps}
# Line 1008 | Line 977 | Another important technique to improve the efficiency
977  
978   %cutoff and minimum image convention
979   Another important technique to improve the efficiency of force
980 < evaluation is to apply cutoff where particles farther than a
981 < predetermined distance, are not included in the calculation
980 > evaluation is to apply spherical cutoff where particles farther than
981 > a predetermined distance are not included in the calculation
982   \cite{Frenkel1996}. The use of a cutoff radius will cause a
983   discontinuity in the potential energy curve. Fortunately, one can
984 < shift the potential to ensure the potential curve go smoothly to
985 < zero at the cutoff radius. Cutoff strategy works pretty well for
986 < Lennard-Jones interaction because of its short range nature.
987 < However, simply truncating the electrostatic interaction with the
988 < use of cutoff has been shown to lead to severe artifacts in
989 < simulations. Ewald summation, in which the slowly conditionally
990 < convergent Coulomb potential is transformed into direct and
991 < reciprocal sums with rapid and absolute convergence, has proved to
992 < minimize the periodicity artifacts in liquid simulations. Taking the
993 < advantages of the fast Fourier transform (FFT) for calculating
994 < discrete Fourier transforms, the particle mesh-based
984 > shift simple radial potential to ensure the potential curve go
985 > smoothly to zero at the cutoff radius. The cutoff strategy works
986 > well for Lennard-Jones interaction because of its short range
987 > nature. However, simply truncating the electrostatic interaction
988 > with the use of cutoffs has been shown to lead to severe artifacts
989 > in simulations. The Ewald summation, in which the slowly decaying
990 > Coulomb potential is transformed into direct and reciprocal sums
991 > with rapid and absolute convergence, has proved to minimize the
992 > periodicity artifacts in liquid simulations. Taking the advantages
993 > of the fast Fourier transform (FFT) for calculating discrete Fourier
994 > transforms, the particle mesh-based
995   methods\cite{Hockney1981,Shimada1993, Luty1994} are accelerated from
996 < $O(N^{3/2})$ to $O(N logN)$. An alternative approach is \emph{fast
997 < multipole method}\cite{Greengard1987, Greengard1994}, which treats
998 < Coulombic interaction exactly at short range, and approximate the
999 < potential at long range through multipolar expansion. In spite of
1000 < their wide acceptances at the molecular simulation community, these
1001 < two methods are hard to be implemented correctly and efficiently.
1002 < Instead, we use a damped and charge-neutralized Coulomb potential
1003 < method developed by Wolf and his coworkers\cite{Wolf1999}. The
1004 < shifted Coulomb potential for particle $i$ and particle $j$ at
1005 < distance $r_{rj}$ is given by:
996 > $O(N^{3/2})$ to $O(N logN)$. An alternative approach is the
997 > \emph{fast multipole method}\cite{Greengard1987, Greengard1994},
998 > which treats Coulombic interactions exactly at short range, and
999 > approximate the potential at long range through multipolar
1000 > expansion. In spite of their wide acceptance at the molecular
1001 > simulation community, these two methods are difficult to implement
1002 > correctly and efficiently. Instead, we use a damped and
1003 > charge-neutralized Coulomb potential method developed by Wolf and
1004 > his coworkers\cite{Wolf1999}. The shifted Coulomb potential for
1005 > particle $i$ and particle $j$ at distance $r_{rj}$ is given by:
1006   \begin{equation}
1007   V(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha
1008   r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow
# Line 1055 | Line 1024 | illustration of shifted Coulomb potential.}
1024  
1025   \subsection{\label{introSection:Analysis} Analysis}
1026  
1027 < Recently, advanced visualization technique are widely applied to
1027 > Recently, advanced visualization technique have become applied to
1028   monitor the motions of molecules. Although the dynamics of the
1029   system can be described qualitatively from animation, quantitative
1030 < trajectory analysis are more appreciable. According to the
1031 < principles of Statistical Mechanics,
1032 < Sec.~\ref{introSection:statisticalMechanics}, one can compute
1033 < thermodynamics properties, analyze fluctuations of structural
1034 < parameters, and investigate time-dependent processes of the molecule
1066 < from the trajectories.
1030 > trajectory analysis are more useful. According to the principles of
1031 > Statistical Mechanics, Sec.~\ref{introSection:statisticalMechanics},
1032 > one can compute thermodynamic properties, analyze fluctuations of
1033 > structural parameters, and investigate time-dependent processes of
1034 > the molecule from the trajectories.
1035  
1036 < \subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamics Properties}}
1036 > \subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamic Properties}}
1037  
1038 < Thermodynamics properties, which can be expressed in terms of some
1038 > Thermodynamic properties, which can be expressed in terms of some
1039   function of the coordinates and momenta of all particles in the
1040   system, can be directly computed from molecular dynamics. The usual
1041   way to measure the pressure is based on virial theorem of Clausius
# Line 1090 | Line 1058 | Structural Properties of a simple fluid can be describ
1058   \subsubsection{\label{introSection:structuralProperties}\textbf{Structural Properties}}
1059  
1060   Structural Properties of a simple fluid can be described by a set of
1061 < distribution functions. Among these functions,\emph{pair
1061 > distribution functions. Among these functions,the \emph{pair
1062   distribution function}, also known as \emph{radial distribution
1063 < function}, is of most fundamental importance to liquid-state theory.
1064 < Pair distribution function can be gathered by Fourier transforming
1065 < raw data from a series of neutron diffraction experiments and
1066 < integrating over the surface factor \cite{Powles1973}. The
1067 < experiment result can serve as a criterion to justify the
1068 < correctness of the theory. Moreover, various equilibrium
1069 < thermodynamic and structural properties can also be expressed in
1070 < terms of radial distribution function \cite{Allen1987}.
1071 <
1104 < A pair distribution functions $g(r)$ gives the probability that a
1063 > function}, is of most fundamental importance to liquid theory.
1064 > Experimentally, pair distribution function can be gathered by
1065 > Fourier transforming raw data from a series of neutron diffraction
1066 > experiments and integrating over the surface factor
1067 > \cite{Powles1973}. The experimental results can serve as a criterion
1068 > to justify the correctness of a liquid model. Moreover, various
1069 > equilibrium thermodynamic and structural properties can also be
1070 > expressed in terms of radial distribution function \cite{Allen1987}.
1071 > The pair distribution functions $g(r)$ gives the probability that a
1072   particle $i$ will be located at a distance $r$ from a another
1073   particle $j$ in the system
1074   \[
1075   g(r) = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j
1076 < \ne i} {\delta (r - r_{ij} )} } } \right\rangle.
1076 > \ne i} {\delta (r - r_{ij} )} } } \right\rangle = \frac{\rho
1077 > (r)}{\rho}.
1078   \]
1079   Note that the delta function can be replaced by a histogram in
1080 < computer simulation. Figure
1081 < \ref{introFigure:pairDistributionFunction} shows a typical pair
1082 < distribution function for the liquid argon system. The occurrence of
1115 < several peaks in the plot of $g(r)$ suggests that it is more likely
1116 < to find particles at certain radial values than at others. This is a
1117 < result of the attractive interaction at such distances. Because of
1118 < the strong repulsive forces at short distance, the probability of
1119 < locating particles at distances less than about 2.5{\AA} from each
1120 < other is essentially zero.
1080 > computer simulation. Peaks in $g(r)$ represent solvent shells, and
1081 > the height of these peaks gradually decreases to 1 as the liquid of
1082 > large distance approaches the bulk density.
1083  
1122 – %\begin{figure}
1123 – %\centering
1124 – %\includegraphics[width=\linewidth]{pdf.eps}
1125 – %\caption[Pair distribution function for the liquid argon
1126 – %]{Pair distribution function for the liquid argon}
1127 – %\label{introFigure:pairDistributionFunction}
1128 – %\end{figure}
1084  
1085   \subsubsection{\label{introSection:timeDependentProperties}\textbf{Time-dependent
1086   Properties}}
1087  
1088   Time-dependent properties are usually calculated using \emph{time
1089 < correlation function}, which correlates random variables $A$ and $B$
1090 < at two different time
1089 > correlation functions}, which correlate random variables $A$ and $B$
1090 > at two different times,
1091   \begin{equation}
1092   C_{AB} (t) = \left\langle {A(t)B(0)} \right\rangle.
1093   \label{introEquation:timeCorrelationFunction}
1094   \end{equation}
1095   If $A$ and $B$ refer to same variable, this kind of correlation
1096 < function is called \emph{auto correlation function}. One example of
1097 < auto correlation function is velocity auto-correlation function
1098 < which is directly related to transport properties of molecular
1099 < liquids:
1096 > function is called an \emph{autocorrelation function}. One example
1097 > of an auto correlation function is the velocity auto-correlation
1098 > function which is directly related to transport properties of
1099 > molecular liquids:
1100   \[
1101   D = \frac{1}{3}\int\limits_0^\infty  {\left\langle {v(t) \cdot v(0)}
1102   \right\rangle } dt
1103   \]
1104 < where $D$ is diffusion constant. Unlike velocity autocorrelation
1105 < function which is averaging over time origins and over all the
1106 < atoms, dipole autocorrelation are calculated for the entire system.
1107 < The dipole autocorrelation function is given by:
1104 > where $D$ is diffusion constant. Unlike the velocity autocorrelation
1105 > function, which is averaging over time origins and over all the
1106 > atoms, the dipole autocorrelation functions are calculated for the
1107 > entire system. The dipole autocorrelation function is given by:
1108   \[
1109   c_{dipole}  = \left\langle {u_{tot} (t) \cdot u_{tot} (t)}
1110   \right\rangle
# Line 1175 | Line 1130 | movement of the objects in 3D gaming engine or other p
1130   areas, from engineering, physics, to chemistry. For example,
1131   missiles and vehicle are usually modeled by rigid bodies.  The
1132   movement of the objects in 3D gaming engine or other physics
1133 < simulator is governed by the rigid body dynamics. In molecular
1134 < simulation, rigid body is used to simplify the model in
1135 < protein-protein docking study\cite{Gray2003}.
1133 > simulator is governed by rigid body dynamics. In molecular
1134 > simulations, rigid bodies are used to simplify protein-protein
1135 > docking studies\cite{Gray2003}.
1136  
1137   It is very important to develop stable and efficient methods to
1138 < integrate the equations of motion of orientational degrees of
1139 < freedom. Euler angles are the nature choice to describe the
1140 < rotational degrees of freedom. However, due to its singularity, the
1141 < numerical integration of corresponding equations of motion is very
1142 < inefficient and inaccurate. Although an alternative integrator using
1143 < different sets of Euler angles can overcome this
1144 < difficulty\cite{Barojas1973}, the computational penalty and the lost
1145 < of angular momentum conservation still remain. A singularity free
1146 < representation utilizing quaternions was developed by Evans in
1147 < 1977\cite{Evans1977}. Unfortunately, this approach suffer from the
1148 < nonseparable Hamiltonian resulted from quaternion representation,
1149 < which prevents the symplectic algorithm to be utilized. Another
1150 < different approach is to apply holonomic constraints to the atoms
1151 < belonging to the rigid body. Each atom moves independently under the
1152 < normal forces deriving from potential energy and constraint forces
1153 < which are used to guarantee the rigidness. However, due to their
1154 < iterative nature, SHAKE and Rattle algorithm converge very slowly
1155 < when the number of constraint increases\cite{Ryckaert1977,
1156 < Andersen1983}.
1138 > integrate the equations of motion for orientational degrees of
1139 > freedom. Euler angles are the natural choice to describe the
1140 > rotational degrees of freedom. However, due to $\frac {1}{sin
1141 > \theta}$ singularities, the numerical integration of corresponding
1142 > equations of motion is very inefficient and inaccurate. Although an
1143 > alternative integrator using multiple sets of Euler angles can
1144 > overcome this difficulty\cite{Barojas1973}, the computational
1145 > penalty and the loss of angular momentum conservation still remain.
1146 > A singularity-free representation utilizing quaternions was
1147 > developed by Evans in 1977\cite{Evans1977}. Unfortunately, this
1148 > approach uses a nonseparable Hamiltonian resulting from the
1149 > quaternion representation, which prevents the symplectic algorithm
1150 > to be utilized. Another different approach is to apply holonomic
1151 > constraints to the atoms belonging to the rigid body. Each atom
1152 > moves independently under the normal forces deriving from potential
1153 > energy and constraint forces which are used to guarantee the
1154 > rigidness. However, due to their iterative nature, the SHAKE and
1155 > Rattle algorithms also converge very slowly when the number of
1156 > constraints increases\cite{Ryckaert1977, Andersen1983}.
1157  
1158 < The break through in geometric literature suggests that, in order to
1158 > A break-through in geometric literature suggests that, in order to
1159   develop a long-term integration scheme, one should preserve the
1160 < symplectic structure of the flow. Introducing conjugate momentum to
1161 < rotation matrix $Q$ and re-formulating Hamiltonian's equation, a
1162 < symplectic integrator, RSHAKE\cite{Kol1997}, was proposed to evolve
1163 < the Hamiltonian system in a constraint manifold by iteratively
1164 < satisfying the orthogonality constraint $Q_T Q = 1$. An alternative
1165 < method using quaternion representation was developed by
1166 < Omelyan\cite{Omelyan1998}. However, both of these methods are
1167 < iterative and inefficient. In this section, we will present a
1160 > symplectic structure of the flow. By introducing a conjugate
1161 > momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's
1162 > equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was
1163 > proposed to evolve the Hamiltonian system in a constraint manifold
1164 > by iteratively satisfying the orthogonality constraint $Q^T Q = 1$.
1165 > An alternative method using the quaternion representation was
1166 > developed by Omelyan\cite{Omelyan1998}. However, both of these
1167 > methods are iterative and inefficient. In this section, we descibe a
1168   symplectic Lie-Poisson integrator for rigid body developed by
1169   Dullweber and his coworkers\cite{Dullweber1997} in depth.
1170  
1171 < \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
1172 < The motion of the rigid body is Hamiltonian with the Hamiltonian
1171 > \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Bodies}
1172 > The motion of a rigid body is Hamiltonian with the Hamiltonian
1173   function
1174   \begin{equation}
1175   H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) +
# Line 1228 | Line 1183 | where $I_{ii}$ is the diagonal element of the inertia
1183   I_{ii}^{ - 1}  = \frac{1}{2}\sum\limits_{i \ne j} {J_{jj}^{ - 1} }
1184   \]
1185   where $I_{ii}$ is the diagonal element of the inertia tensor. This
1186 < constrained Hamiltonian equation subjects to a holonomic constraint,
1186 > constrained Hamiltonian equation is subjected to a holonomic
1187 > constraint,
1188   \begin{equation}
1189   Q^T Q = 1, \label{introEquation:orthogonalConstraint}
1190   \end{equation}
1191 < which is used to ensure rotation matrix's orthogonality.
1192 < Differentiating \ref{introEquation:orthogonalConstraint} and using
1193 < Equation \ref{introEquation:RBMotionMomentum}, one may obtain,
1191 > which is used to ensure rotation matrix's unitarity. Differentiating
1192 > \ref{introEquation:orthogonalConstraint} and using Equation
1193 > \ref{introEquation:RBMotionMomentum}, one may obtain,
1194   \begin{equation}
1195   Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
1196   \label{introEquation:RBFirstOrderConstraint}
1197   \end{equation}
1242 –
1198   Using Equation (\ref{introEquation:motionHamiltonianCoordinate},
1199   \ref{introEquation:motionHamiltonianMomentum}), one can write down
1200   the equations of motion,
1246 –
1201   \begin{eqnarray}
1202   \frac{{dq}}{{dt}} & = & \frac{p}{m} \label{introEquation:RBMotionPosition}\\
1203   \frac{{dp}}{{dt}} & = & - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
1204   \frac{{dQ}}{{dt}} & = & PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
1205   \frac{{dP}}{{dt}} & = & - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}
1206   \end{eqnarray}
1253 –
1207   In general, there are two ways to satisfy the holonomic constraints.
1208 < We can use constraint force provided by lagrange multiplier on the
1209 < normal manifold to keep the motion on constraint space. Or we can
1210 < simply evolve the system in constraint manifold. These two methods
1211 < are proved to be equivalent. The holonomic constraint and equations
1212 < of motions define a constraint manifold for rigid body
1208 > We can use a constraint force provided by a Lagrange multiplier on
1209 > the normal manifold to keep the motion on constraint space. Or we
1210 > can simply evolve the system on the constraint manifold. These two
1211 > methods have been proved to be equivalent. The holonomic constraint
1212 > and equations of motions define a constraint manifold for rigid
1213 > bodies
1214   \[
1215   M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0}
1216   \right\}.
1217   \]
1264 –
1218   Unfortunately, this constraint manifold is not the cotangent bundle
1219 < $T_{\star}SO(3)$. However, it turns out that under symplectic
1219 > $T^* SO(3)$ which can be consider as a symplectic manifold on Lie
1220 > rotation group $SO(3)$. However, it turns out that under symplectic
1221   transformation, the cotangent space and the phase space are
1222 < diffeomorphic. Introducing
1222 > diffeomorphic. By introducing
1223   \[
1224   \tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
1225   \]
# Line 1275 | Line 1229 | T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \t
1229   T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \tilde Q =
1230   1,\tilde Q^T \tilde PJ^{ - 1}  + J^{ - 1} P^T \tilde Q = 0} \right\}
1231   \]
1278 –
1232   For a body fixed vector $X_i$ with respect to the center of mass of
1233   the rigid body, its corresponding lab fixed vector $X_0^{lab}$  is
1234   given as
# Line 1294 | Line 1247 | and
1247   \[
1248   \nabla _Q V(q,Q) = F(q,Q)X_i^t
1249   \]
1250 < respectively.
1251 <
1252 < As a common choice to describe the rotation dynamics of the rigid
1300 < body, angular momentum on body frame $\Pi  = Q^t P$ is introduced to
1301 < rewrite the equations of motion,
1250 > respectively. As a common choice to describe the rotation dynamics
1251 > of the rigid body, the angular momentum on the body fixed frame $\Pi
1252 > = Q^t P$ is introduced to rewrite the equations of motion,
1253   \begin{equation}
1254   \begin{array}{l}
1255 < \mathop \Pi \limits^ \bullet   = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda  \\
1256 < \mathop Q\limits^{{\rm{   }} \bullet }  = Q\Pi {\rm{ }}J^{ - 1}  \\
1255 > \dot \Pi  = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda,  \\
1256 > \dot Q  = Q\Pi {\rm{ }}J^{ - 1},  \\
1257   \end{array}
1258   \label{introEqaution:RBMotionPI}
1259   \end{equation}
1260 < , as well as holonomic constraints,
1260 > as well as holonomic constraints,
1261   \[
1262   \begin{array}{l}
1263 < \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
1264 < Q^T Q = 1 \\
1263 > \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0, \\
1264 > Q^T Q = 1 .\\
1265   \end{array}
1266   \]
1316 –
1267   For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a matrix $\hat v \in
1268   so(3)^ \star$, the hat-map isomorphism,
1269   \begin{equation}
# Line 1328 | Line 1278 | operations
1278   will let us associate the matrix products with traditional vector
1279   operations
1280   \[
1281 < \hat vu = v \times u
1281 > \hat vu = v \times u.
1282   \]
1283 < Using \ref{introEqaution:RBMotionPI}, one can construct a skew
1283 > Using Eq.~\ref{introEqaution:RBMotionPI}, one can construct a skew
1284   matrix,
1285 < \begin{equation}
1286 < (\mathop \Pi \limits^ \bullet   - \mathop \Pi \limits^ {\bullet  ^T}
1287 < ){\rm{ }} = {\rm{ }}(\Pi  - \Pi ^T ){\rm{ }}(J^{ - 1} \Pi  + \Pi J^{
1288 < - 1} ) + \sum\limits_i {[Q^T F_i (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]} -
1289 < (\Lambda  - \Lambda ^T ) . \label{introEquation:skewMatrixPI}
1290 < \end{equation}
1291 < Since $\Lambda$ is symmetric, the last term of Equation
1292 < \ref{introEquation:skewMatrixPI} is zero, which implies the Lagrange
1293 < multiplier $\Lambda$ is absent from the equations of motion. This
1294 < unique property eliminate the requirement of iterations which can
1295 < not be avoided in other methods\cite{Kol1997, Omelyan1998}.
1296 <
1297 < Applying hat-map isomorphism, we obtain the equation of motion for
1298 < angular momentum on body frame
1285 > \begin{eqnarray}
1286 > (\dot \Pi  - \dot \Pi ^T ){\rm{ }} &= &{\rm{ }}(\Pi  - \Pi ^T ){\rm{
1287 > }}(J^{ - 1} \Pi  + \Pi J^{ - 1} ) \notag \\
1288 > + \sum\limits_i {[Q^T F_i
1289 > (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]}  - (\Lambda  - \Lambda ^T ).
1290 > \label{introEquation:skewMatrixPI}
1291 > \end{eqnarray}
1292 > Since $\Lambda$ is symmetric, the last term of
1293 > Eq.~\ref{introEquation:skewMatrixPI} is zero, which implies the
1294 > Lagrange multiplier $\Lambda$ is absent from the equations of
1295 > motion. This unique property eliminates the requirement of
1296 > iterations which can not be avoided in other methods\cite{Kol1997,
1297 > Omelyan1998}. Applying the hat-map isomorphism, we obtain the
1298 > equation of motion for angular momentum on body frame
1299   \begin{equation}
1300   \dot \pi  = \pi  \times I^{ - 1} \pi  + \sum\limits_i {\left( {Q^T
1301   F_i (r,Q)} \right) \times X_i }.
# Line 1354 | Line 1304 | given by
1304   In the same manner, the equation of motion for rotation matrix is
1305   given by
1306   \[
1307 < \dot Q = Qskew(I^{ - 1} \pi )
1307 > \dot Q = Qskew(I^{ - 1} \pi ).
1308   \]
1309  
1310   \subsection{\label{introSection:SymplecticFreeRB}Symplectic
1311   Lie-Poisson Integrator for Free Rigid Body}
1312  
1313 < If there is not external forces exerted on the rigid body, the only
1314 < contribution to the rotational is from the kinetic potential (the
1315 < first term of \ref{introEquation:bodyAngularMotion}). The free rigid
1316 < body is an example of Lie-Poisson system with Hamiltonian function
1313 > If there are no external forces exerted on the rigid body, the only
1314 > contribution to the rotational motion is from the kinetic energy
1315 > (the first term of \ref{introEquation:bodyAngularMotion}). The free
1316 > rigid body is an example of a Lie-Poisson system with Hamiltonian
1317 > function
1318   \begin{equation}
1319   T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
1320   \label{introEquation:rotationalKineticRB}
# Line 1375 | Line 1326 | J(\pi ) = \left( {\begin{array}{*{20}c}
1326     0 & {\pi _3 } & { - \pi _2 }  \\
1327     { - \pi _3 } & 0 & {\pi _1 }  \\
1328     {\pi _2 } & { - \pi _1 } & 0  \\
1329 < \end{array}} \right)
1329 > \end{array}} \right).
1330   \end{equation}
1331   Thus, the dynamics of free rigid body is governed by
1332   \begin{equation}
1333 < \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
1333 > \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi ).
1334   \end{equation}
1384 –
1335   One may notice that each $T_i^r$ in Equation
1336   \ref{introEquation:rotationalKineticRB} can be solved exactly. For
1337   instance, the equations of motion due to $T_1^r$ are given by
# Line 1410 | Line 1360 | To reduce the cost of computing expensive functions in
1360   \end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
1361   \]
1362   To reduce the cost of computing expensive functions in $e^{\Delta
1363 < tR_1 }$, we can use Cayley transformation,
1363 > tR_1 }$, we can use Cayley transformation to obtain a single-aixs
1364 > propagator,
1365   \[
1366   e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
1367 < )
1367 > ).
1368   \]
1369   The flow maps for $T_2^r$ and $T_3^r$ can be found in the same
1370 < manner.
1371 <
1421 < In order to construct a second-order symplectic method, we split the
1422 < angular kinetic Hamiltonian function can into five terms
1370 > manner. In order to construct a second-order symplectic method, we
1371 > split the angular kinetic Hamiltonian function can into five terms
1372   \[
1373   T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
1374   ) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
1375 < (\pi _1 )
1376 < \].
1377 < Concatenating flows corresponding to these five terms, we can obtain
1378 < an symplectic integrator,
1375 > (\pi _1 ).
1376 > \]
1377 > By concatenating the propagators corresponding to these five terms,
1378 > we can obtain an symplectic integrator,
1379   \[
1380   \varphi _{\Delta t,T^r }  = \varphi _{\Delta t/2,\pi _1 }  \circ
1381   \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }
1382   \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
1383   _1 }.
1384   \]
1436 –
1385   The non-canonical Lie-Poisson bracket ${F, G}$ of two function
1386   $F(\pi )$ and $G(\pi )$ is defined by
1387   \[
1388   \{ F,G\} (\pi ) = [\nabla _\pi  F(\pi )]^T J(\pi )\nabla _\pi  G(\pi
1389 < )
1389 > ).
1390   \]
1391   If the Poisson bracket of a function $F$ with an arbitrary smooth
1392   function $G$ is zero, $F$ is a \emph{Casimir}, which is the
# Line 1449 | Line 1397 | then by the chain rule
1397   then by the chain rule
1398   \[
1399   \nabla _\pi  F(\pi ) = S'(\frac{{\parallel \pi \parallel ^2
1400 < }}{2})\pi
1400 > }}{2})\pi.
1401   \]
1402 < Thus $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel \pi
1402 > Thus, $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel
1403 > \pi
1404   \parallel ^2 }}{2})\pi  \times \pi  = 0 $. This explicit
1405 < Lie-Poisson integrator is found to be extremely efficient and stable
1406 < which can be explained by the fact the small angle approximation is
1407 < used and the norm of the angular momentum is conserved.
1405 > Lie-Poisson integrator is found to be both extremely efficient and
1406 > stable. These properties can be explained by the fact the small
1407 > angle approximation is used and the norm of the angular momentum is
1408 > conserved.
1409  
1410   \subsection{\label{introSection:RBHamiltonianSplitting} Hamiltonian
1411   Splitting for Rigid Body}
# Line 1463 | Line 1413 | energy and potential energy,
1413   The Hamiltonian of rigid body can be separated in terms of kinetic
1414   energy and potential energy,
1415   \[
1416 < H = T(p,\pi ) + V(q,Q)
1416 > H = T(p,\pi ) + V(q,Q).
1417   \]
1418   The equations of motion corresponding to potential energy and
1419   kinetic energy are listed in the below table,
1420   \begin{table}
1421 < \caption{Equations of motion due to Potential and Kinetic Energies}
1421 > \caption{EQUATIONS OF MOTION DUE TO POTENTIAL AND KINETIC ENERGIES}
1422   \begin{center}
1423   \begin{tabular}{|l|l|}
1424    \hline
# Line 1482 | Line 1432 | kinetic energy are listed in the below table,
1432   \end{tabular}
1433   \end{center}
1434   \end{table}
1435 < A second-order symplectic method is now obtained by the
1436 < composition of the flow maps,
1435 > A second-order symplectic method is now obtained by the composition
1436 > of the position and velocity propagators,
1437   \[
1438   \varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
1439   _{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
1440   \]
1441   Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
1442 < sub-flows which corresponding to force and torque respectively,
1442 > sub-propagators which corresponding to force and torque
1443 > respectively,
1444   \[
1445   \varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
1446   _{\Delta t/2,\tau }.
1447   \]
1448   Since the associated operators of $\varphi _{\Delta t/2,F} $ and
1449 < $\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
1450 < order inside $\varphi _{\Delta t/2,V}$ does not matter.
1451 <
1452 < Furthermore, kinetic potential can be separated to translational
1502 < kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
1449 > $\circ \varphi _{\Delta t/2,\tau }$ commute, the composition order
1450 > inside $\varphi _{\Delta t/2,V}$ does not matter. Furthermore, the
1451 > kinetic energy can be separated to translational kinetic term, $T^t
1452 > (p)$, and rotational kinetic term, $T^r (\pi )$,
1453   \begin{equation}
1454   T(p,\pi ) =T^t (p) + T^r (\pi ).
1455   \end{equation}
1456   where $ T^t (p) = \frac{1}{2}p^T m^{ - 1} p $ and $T^r (\pi )$ is
1457   defined by \ref{introEquation:rotationalKineticRB}. Therefore, the
1458 < corresponding flow maps are given by
1458 > corresponding propagators are given by
1459   \[
1460   \varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
1461   _{\Delta t,T^r }.
1462   \]
1463 < Finally, we obtain the overall symplectic flow maps for free moving
1464 < rigid body
1465 < \begin{equation}
1466 < \begin{array}{c}
1467 < \varphi _{\Delta t}  = \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \\
1468 <  \circ \varphi _{\Delta t,T^t }  \circ \varphi _{\Delta t/2,\pi _1 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi _1 }  \\
1519 <  \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1520 < \end{array}
1463 > Finally, we obtain the overall symplectic propagators for freely
1464 > moving rigid bodies
1465 > \begin{eqnarray*}
1466 > \varphi _{\Delta t}  &=& \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \\
1467 >  & & \circ \varphi _{\Delta t,T^t }  \circ \varphi _{\Delta t/2,\pi _1 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi _1 }  \\
1468 >  & & \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1469   \label{introEquation:overallRBFlowMaps}
1470 < \end{equation}
1470 > \end{eqnarray*}
1471  
1472   \section{\label{introSection:langevinDynamics}Langevin Dynamics}
1473   As an alternative to newtonian dynamics, Langevin dynamics, which
1474   mimics a simple heat bath with stochastic and dissipative forces,
1475   has been applied in a variety of studies. This section will review
1476 < the theory of Langevin dynamics simulation. A brief derivation of
1477 < generalized Langevin equation will be given first. Follow that, we
1478 < will discuss the physical meaning of the terms appearing in the
1479 < equation as well as the calculation of friction tensor from
1480 < hydrodynamics theory.
1476 > the theory of Langevin dynamics. A brief derivation of generalized
1477 > Langevin equation will be given first. Following that, we will
1478 > discuss the physical meaning of the terms appearing in the equation
1479 > as well as the calculation of friction tensor from hydrodynamics
1480 > theory.
1481  
1482   \subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
1483  
1484 < Harmonic bath model, in which an effective set of harmonic
1484 > A harmonic bath model, in which an effective set of harmonic
1485   oscillators are used to mimic the effect of a linearly responding
1486   environment, has been widely used in quantum chemistry and
1487   statistical mechanics. One of the successful applications of
1488 < Harmonic bath model is the derivation of Deriving Generalized
1489 < Langevin Dynamics. Lets consider a system, in which the degree of
1488 > Harmonic bath model is the derivation of the Generalized Langevin
1489 > Dynamics (GLE). Lets consider a system, in which the degree of
1490   freedom $x$ is assumed to couple to the bath linearly, giving a
1491   Hamiltonian of the form
1492   \begin{equation}
1493   H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
1494   \label{introEquation:bathGLE}.
1495   \end{equation}
1496 < Here $p$ is a momentum conjugate to $q$, $m$ is the mass associated
1497 < with this degree of freedom, $H_B$ is harmonic bath Hamiltonian,
1496 > Here $p$ is a momentum conjugate to $x$, $m$ is the mass associated
1497 > with this degree of freedom, $H_B$ is a harmonic bath Hamiltonian,
1498   \[
1499   H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
1500   }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  \omega _\alpha ^2 }
# Line 1554 | Line 1502 | $\omega _\alpha$ are the harmonic bath frequencies, $m
1502   \]
1503   where the index $\alpha$ runs over all the bath degrees of freedom,
1504   $\omega _\alpha$ are the harmonic bath frequencies, $m_\alpha$ are
1505 < the harmonic bath masses, and $\Delta U$ is bilinear system-bath
1505 > the harmonic bath masses, and $\Delta U$ is a bilinear system-bath
1506   coupling,
1507   \[
1508   \Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
1509   \]
1510 < where $g_\alpha$ are the coupling constants between the bath and the
1511 < coordinate $x$. Introducing
1510 > where $g_\alpha$ are the coupling constants between the bath
1511 > coordinates ($x_ \alpha$) and the system coordinate ($x$).
1512 > Introducing
1513   \[
1514   W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
1515   }}{{2m_\alpha  w_\alpha ^2 }}} x^2
1516 < \] and combining the last two terms in Equation
1517 < \ref{introEquation:bathGLE}, we may rewrite the Harmonic bath
1569 < Hamiltonian as
1516 > \]
1517 > and combining the last two terms in Eq.~\ref{introEquation:bathGLE}, we may rewrite the Harmonic bath Hamiltonian as
1518   \[
1519   H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
1520   {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
1521   w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
1522 < w_\alpha ^2 }}x} \right)^2 } \right\}}
1522 > w_\alpha ^2 }}x} \right)^2 } \right\}}.
1523   \]
1524   Since the first two terms of the new Hamiltonian depend only on the
1525   system coordinates, we can get the equations of motion for
1526 < Generalized Langevin Dynamics by Hamilton's equations
1579 < \ref{introEquation:motionHamiltonianCoordinate,
1580 < introEquation:motionHamiltonianMomentum},
1526 > Generalized Langevin Dynamics by Hamilton's equations,
1527   \begin{equation}
1528   m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
1529   \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   -
# Line 1590 | Line 1536 | m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x
1536   \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right).
1537   \label{introEquation:bathMotionGLE}
1538   \end{equation}
1593 –
1539   In order to derive an equation for $x$, the dynamics of the bath
1540   variables $x_\alpha$ must be solved exactly first. As an integral
1541   transform which is particularly useful in solving linear ordinary
1542 < differential equations, Laplace transform is the appropriate tool to
1543 < solve this problem. The basic idea is to transform the difficult
1542 > differential equations,the Laplace transform is the appropriate tool
1543 > to solve this problem. The basic idea is to transform the difficult
1544   differential equations into simple algebra problems which can be
1545 < solved easily. Then applying inverse Laplace transform, also known
1546 < as the Bromwich integral, we can retrieve the solutions of the
1547 < original problems.
1548 <
1604 < Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace
1605 < transform of f(t) is a new function defined as
1545 > solved easily. Then, by applying the inverse Laplace transform, also
1546 > known as the Bromwich integral, we can retrieve the solutions of the
1547 > original problems. Let $f(t)$ be a function defined on $ [0,\infty )
1548 > $. The Laplace transform of f(t) is a new function defined as
1549   \[
1550   L(f(t)) \equiv F(p) = \int_0^\infty  {f(t)e^{ - pt} dt}
1551   \]
1552   where  $p$ is real and  $L$ is called the Laplace Transform
1553   Operator. Below are some important properties of Laplace transform
1611 –
1554   \begin{eqnarray*}
1555   L(x + y)  & = & L(x) + L(y) \\
1556   L(ax)     & = & aL(x) \\
# Line 1616 | Line 1558 | Operator. Below are some important properties of Lapla
1558   L(\ddot x)& = & p^2 L(x) - px(0) - \dot x(0) \\
1559   L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right)& = & G(p)H(p) \\
1560   \end{eqnarray*}
1561 <
1620 <
1621 < Applying Laplace transform to the bath coordinates, we obtain
1561 > Applying the Laplace transform to the bath coordinates, we obtain
1562   \begin{eqnarray*}
1563   p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) & = & - \omega _\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) \\
1564   L(x_\alpha  ) & = & \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
1565   \end{eqnarray*}
1626 –
1566   By the same way, the system coordinates become
1567   \begin{eqnarray*}
1568 < mL(\ddot x) & = & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
1569 <  & & \mbox{} - \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x) - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0) - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}  \\
1568 > mL(\ddot x) & = &
1569 >  - \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x) - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0) - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}  \\
1570 >  & & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}}
1571   \end{eqnarray*}
1632 –
1572   With the help of some relatively important inverse Laplace
1573   transformations:
1574   \[
# Line 1639 | Line 1578 | transformations:
1578   L(1) = \frac{1}{p} \\
1579   \end{array}
1580   \]
1581 < , we obtain
1581 > we obtain
1582   \begin{eqnarray*}
1583   m\ddot x & =  & - \frac{{\partial W(x)}}{{\partial x}} -
1584   \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
# Line 1695 | Line 1634 | This property is what we expect from a truly random pr
1634   \end{array}
1635   \]
1636   This property is what we expect from a truly random process. As long
1637 < as the model, which is gaussian distribution in general, chosen for
1638 < $R(t)$ is a truly random process, the stochastic nature of the GLE
1700 < still remains.
1637 > as the model chosen for $R(t)$ was a gaussian distribution in
1638 > general, the stochastic nature of the GLE still remains.
1639  
1640   %dynamic friction kernel
1641   The convolution integral
# Line 1713 | Line 1651 | $\xi(t) = \Xi_0$. Hence, the convolution integral beco
1651   \[
1652   \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
1653   \]
1654 < and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1654 > and Eq.~\ref{introEuqation:GeneralizedLangevinDynamics} becomes
1655   \[
1656   m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
1657   \frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
1658   \]
1659 < which can be used to describe dynamic caging effect. The other
1660 < extreme is the bath that responds infinitely quickly to motions in
1661 < the system. Thus, $\xi (t)$ can be taken as a $delta$ function in
1662 < time:
1659 > which can be used to describe the effect of dynamic caging in
1660 > viscous solvents. The other extreme is the bath that responds
1661 > infinitely quickly to motions in the system. Thus, $\xi (t)$ can be
1662 > taken as a $delta$ function in time:
1663   \[
1664   \xi (t) = 2\xi _0 \delta (t)
1665   \]
# Line 1730 | Line 1668 | Hence, the convolution integral becomes
1668   \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
1669   {\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
1670   \]
1671 < and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1671 > and Eq.~\ref{introEuqation:GeneralizedLangevinDynamics} becomes
1672   \begin{equation}
1673   m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \xi _0 \dot
1674   x(t) + R(t) \label{introEquation:LangevinEquation}
# Line 1753 | Line 1691 | And since the $q$ coordinates are harmonic oscillators
1691   R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
1692   \]
1693   And since the $q$ coordinates are harmonic oscillators,
1756 –
1694   \begin{eqnarray*}
1695   \left\langle {q_\alpha ^2 } \right\rangle  & = & \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
1696   \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle & = & \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
# Line 1762 | Line 1699 | And since the $q$ coordinates are harmonic oscillators
1699    & = &\sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
1700    & = &kT\xi (t) \\
1701   \end{eqnarray*}
1765 –
1702   Thus, we recover the \emph{second fluctuation dissipation theorem}
1703   \begin{equation}
1704   \xi (t) = \left\langle {R(t)R(0)} \right\rangle

Diff Legend

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