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 2844 by tim, Fri Jun 9 19:06:37 2006 UTC vs.
Revision 2899 by tim, Tue Jun 27 04:04:32 2006 UTC

# Line 31 | Line 31 | F_{ij} = -F_{ji}
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 | that if all forces are conservative, Energy
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
75 > can only be described in cartesian coordinate systems. Moreover, it
76 > becomes 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
# Line 85 | Line 84 | Hamilton's Principle may be stated as follows,
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
# Line 113 | Line 110 | For a holonomic system of $f$ degrees of freedom, the
110   \subsubsection{\label{introSection:equationOfMotionLagrangian}\textbf{The
111   Equations of Motion in Lagrangian Mechanics}}
112  
113 < For a holonomic system of $f$ degrees of freedom, the equations of
114 < motion in the Lagrangian form is
113 > For a system of $f$ degrees of freedom, the equations of motion in
114 > the Lagrangian form is
115   \begin{equation}
116   \frac{d}{{dt}}\frac{{\partial L}}{{\partial \dot q_i }} -
117   \frac{{\partial L}}{{\partial q_i }} = 0,{\rm{ }}i = 1, \ldots,f
# 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 | $L$ is the Lagrangian function for the system.
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}
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 }
# 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
# Line 227 | Line 219 | this system is a $6f$ dimensional space. A point, $x =
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.
222 > this system is a $6f$ dimensional space. A point, $x = (\rightarrow
223 > q_1 , \ldots ,\rightarrow q_f ,\rightarrow p_1 , \ldots ,\rightarrow
224 > p_f )$, with a unique set of values of $6f$ coordinates and momenta
225 > is a phase space vector.
226 > %%%fix me
227  
228 < A microscopic state or microstate of a classical system is
235 < specification of the complete phase space vector of a system at any
236 < instant in time. An ensemble is defined as a collection of systems
237 < sharing one or more macroscopic characteristics but each being in a
238 < unique microstate. The complete ensemble is specified by giving all
239 < systems or microstates consistent with the common macroscopic
240 < characteristics of the ensemble. Although the state of each
241 < individual system in the ensemble could be precisely described at
242 < any instance in time by a suitable phase space vector, when using
243 < ensembles for statistical purposes, there is no need to maintain
244 < distinctions between individual systems, since the numbers of
245 < systems at any time in the different states which correspond to
246 < different regions of the phase space are more interesting. Moreover,
247 < in the point of view of statistical mechanics, one would prefer to
248 < use ensembles containing a large enough population of separate
249 < members so that the numbers of systems in such different states can
250 < be regarded as changing continuously as we traverse different
251 < regions of the phase space. The condition of an ensemble at any time
228 > In statistical mechanics, the condition of an ensemble at any time
229   can be regarded as appropriately specified by the density $\rho$
230   with which representative points are distributed over the phase
231   space. The density distribution for an ensemble with $f$ degrees of
# Line 282 | Line 259 | With the help of Equation(\ref{introEquation:unitProba
259   {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
260   \label{introEquation:unitProbability}
261   \end{equation}
262 < With the help of Equation(\ref{introEquation:unitProbability}) and
263 < the knowledge of the system, it is possible to calculate the average
262 > With the help of Eq.~\ref{introEquation:unitProbability} and the
263 > knowledge of the system, it is possible to calculate the average
264   value of any desired quantity which depends on the coordinates and
265   momenta of the system. Even when the dynamics of the real system is
266   complex, or stochastic, or even discontinuous, the average
# Line 303 | Line 280 | thermodynamic equilibrium.
280   statistical characteristics. As a function of macroscopic
281   parameters, such as temperature \textit{etc}, the partition function
282   can be used to describe the statistical properties of a system in
283 < thermodynamic equilibrium.
284 <
285 < As an ensemble of systems, each of which is known to be thermally
309 < isolated and conserve energy, the Microcanonical ensemble(NVE) has a
310 < partition function like,
283 > thermodynamic equilibrium. As an ensemble of systems, each of which
284 > is known to be thermally isolated and conserve energy, the
285 > Microcanonical ensemble (NVE) has a partition function like,
286   \begin{equation}
287 < \Omega (N,V,E) = e^{\beta TS} \label{introEquation:NVEPartition}.
287 > \Omega (N,V,E) = e^{\beta TS}. \label{introEquation:NVEPartition}.
288   \end{equation}
289 < A canonical ensemble(NVT)is an ensemble of systems, each of which
289 > A canonical ensemble (NVT)is an ensemble of systems, each of which
290   can share its energy with a large heat reservoir. The distribution
291   of the total energy amongst the possible dynamical states is given
292   by the partition function,
293   \begin{equation}
294 < \Omega (N,V,T) = e^{ - \beta A}
294 > \Omega (N,V,T) = e^{ - \beta A}.
295   \label{introEquation:NVTPartition}
296   \end{equation}
297   Here, $A$ is the Helmholtz free energy which is defined as $ A = U -
298   TS$. Since most experiments are carried out under constant pressure
299 < condition, the isothermal-isobaric ensemble(NPT) plays a very
299 > condition, the isothermal-isobaric ensemble (NPT) plays a very
300   important role in molecular simulations. The isothermal-isobaric
301   ensemble allow the system to exchange energy with a heat bath of
302   temperature $T$ and to change the volume as well. Its partition
# Line 337 | Line 312 | $\rho$, we begin from Equation(\ref{introEquation:delt
312   Liouville's theorem is the foundation on which statistical mechanics
313   rests. It describes the time evolution of the phase space
314   distribution function. In order to calculate the rate of change of
315 < $\rho$, we begin from Equation(\ref{introEquation:deltaN}). If we
316 < consider the two faces perpendicular to the $q_1$ axis, which are
317 < located at $q_1$ and $q_1 + \delta q_1$, the number of phase points
318 < leaving the opposite face is given by the expression,
315 > $\rho$, we begin from Eq.~\ref{introEquation:deltaN}. If we consider
316 > the two faces perpendicular to the $q_1$ axis, which are located at
317 > $q_1$ and $q_1 + \delta q_1$, the number of phase points leaving the
318 > opposite face is given by the expression,
319   \begin{equation}
320   \left( {\rho  + \frac{{\partial \rho }}{{\partial q_1 }}\delta q_1 }
321   \right)\left( {\dot q_1  + \frac{{\partial \dot q_1 }}{{\partial q_1
# Line 373 | Line 348 | simple form,
348   \frac{{\partial \rho }}{{\partial p_i }}\dot p_i } \right)}  = 0 .
349   \label{introEquation:liouvilleTheorem}
350   \end{equation}
376 –
351   Liouville's theorem states that the distribution function is
352   constant along any trajectory in phase space. In classical
353 < statistical mechanics, since the number of particles in the system
354 < is huge, we may be able to believe the system is stationary,
353 > statistical mechanics, since the number of members in an ensemble is
354 > huge and constant, we can assume the local density has no reason
355 > (other than classical mechanics) to change,
356   \begin{equation}
357   \frac{{\partial \rho }}{{\partial t}} = 0.
358   \label{introEquation:stationary}
# Line 430 | Line 405 | Substituting equations of motion in Hamiltonian formal
405   \label{introEquation:poissonBracket}
406   \end{equation}
407   Substituting equations of motion in Hamiltonian formalism(
408 < \ref{introEquation:motionHamiltonianCoordinate} ,
409 < \ref{introEquation:motionHamiltonianMomentum} ) into
410 < (\ref{introEquation:liouvilleTheorem}), we can rewrite Liouville's
411 < theorem using Poisson bracket notion,
408 > Eq.~\ref{introEquation:motionHamiltonianCoordinate} ,
409 > Eq.~\ref{introEquation:motionHamiltonianMomentum} ) into
410 > (Eq.~\ref{introEquation:liouvilleTheorem}), we can rewrite
411 > Liouville's theorem using Poisson bracket notion,
412   \begin{equation}
413   \left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - \left\{
414   {\rho ,H} \right\}.
# Line 494 | Line 469 | Leimkuhler1999}. The velocity verlet method, which hap
469   geometric integrators, which preserve various phase-flow invariants
470   such as symplectic structure, volume and time reversal symmetry, are
471   developed to address this issue\cite{Dullweber1997, McLachlan1998,
472 < Leimkuhler1999}. The velocity verlet method, which happens to be a
472 > Leimkuhler1999}. The velocity Verlet method, which happens to be a
473   simple example of symplectic integrator, continues to gain
474   popularity in the molecular dynamics community. This fact can be
475   partly explained by its geometric nature.
# Line 514 | Line 489 | an example of symplectic form.
489   $\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+
490   \lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and
491   $\omega(x, x) = 0$. The cross product operation in vector field is
492 < an example of symplectic form.
493 <
494 < One of the motivations to study \emph{symplectic manifolds} in
495 < Hamiltonian Mechanics is that a symplectic manifold can represent
496 < all possible configurations of the system and the phase space of the
497 < system can be described by it's cotangent bundle. Every symplectic
498 < manifold is even dimensional. For instance, in Hamilton equations,
524 < coordinate and momentum always appear in pairs.
492 > an example of symplectic form. One of the motivations to study
493 > \emph{symplectic manifolds} in Hamiltonian Mechanics is that a
494 > symplectic manifold can represent all possible configurations of the
495 > system and the phase space of the system can be described by it's
496 > cotangent bundle. Every symplectic manifold is even dimensional. For
497 > instance, in Hamilton equations, coordinate and momentum always
498 > appear in pairs.
499  
500   \subsection{\label{introSection:ODE}Ordinary Differential Equations}
501  
# Line 548 | Line 522 | called a \emph{Hamiltonian vector field}.
522   \frac{d}{{dt}}x = J\nabla _x H(x)
523   \label{introEquation:compactHamiltonian}
524   \end{equation}In this case, $f$ is
525 < called a \emph{Hamiltonian vector field}.
526 <
553 < Another generalization of Hamiltonian dynamics is Poisson
554 < Dynamics\cite{Olver1986},
525 > called a \emph{Hamiltonian vector field}. Another generalization of
526 > Hamiltonian dynamics is Poisson Dynamics\cite{Olver1986},
527   \begin{equation}
528   \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
529   \end{equation}
# Line 589 | Line 561 | Instead, we use a approximate map, $\psi_\tau$, which
561   \end{equation}
562  
563   In most cases, it is not easy to find the exact flow $\varphi_\tau$.
564 < Instead, we use a approximate map, $\psi_\tau$, which is usually
564 > Instead, we use an approximate map, $\psi_\tau$, which is usually
565   called integrator. The order of an integrator $\psi_\tau$ is $p$, if
566   the Taylor series of $\psi_\tau$ agree to order $p$,
567   \begin{equation}
568 < \psi_tau(x) = x + \tau f(x) + O(\tau^{p+1})
568 > \psi_\tau(x) = x + \tau f(x) + O(\tau^{p+1})
569   \end{equation}
570  
571   \subsection{\label{introSection:geometricProperties}Geometric Properties}
572  
573 < The hidden geometric properties\cite{Budd1999, Marsden1998} of ODE
574 < and its flow play important roles in numerical studies. Many of them
575 < can be found in systems which occur naturally in applications.
604 <
573 > The hidden geometric properties\cite{Budd1999, Marsden1998} of an
574 > ODE and its flow play important roles in numerical studies. Many of
575 > them can be found in systems which occur naturally in applications.
576   Let $\varphi$ be the flow of Hamiltonian vector field, $\varphi$ is
577   a \emph{symplectic} flow if it satisfies,
578   \begin{equation}
# Line 615 | Line 586 | is the property must be preserved by the integrator.
586   \begin{equation}
587   {\varphi '}^T J \varphi ' = J \circ \varphi
588   \end{equation}
589 < is the property must be preserved by the integrator.
590 <
591 < It is possible to construct a \emph{volume-preserving} flow for a
592 < source free($ \nabla \cdot f = 0 $) ODE, if the flow satisfies $
593 < \det d\varphi  = 1$. One can show easily that a symplectic flow will
594 < be volume-preserving.
624 <
625 < Changing the variables $y = h(x)$ in a ODE\ref{introEquation:ODE}
626 < will result in a new system,
589 > is the property that must be preserved by the integrator. It is
590 > possible to construct a \emph{volume-preserving} flow for a source
591 > free ODE ($ \nabla \cdot f = 0 $), if the flow satisfies $ \det
592 > d\varphi  = 1$. One can show easily that a symplectic flow will be
593 > volume-preserving. Changing the variables $y = h(x)$ in an ODE
594 > (Eq.~\ref{introEquation:ODE}) will result in a new system,
595   \[
596   \dot y = \tilde f(y) = ((dh \cdot f)h^{ - 1} )(y).
597   \]
598   The vector filed $f$ has reversing symmetry $h$ if $f = - \tilde f$.
599   In other words, the flow of this vector field is reversible if and
600 < only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $.
601 <
634 < A \emph{first integral}, or conserved quantity of a general
600 > only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $. A
601 > \emph{first integral}, or conserved quantity of a general
602   differential function is a function $ G:R^{2d}  \to R^d $ which is
603   constant for all solutions of the ODE $\frac{{dx}}{{dt}} = f(x)$ ,
604   \[
# Line 644 | Line 611 | smooth function $G$ is given by,
611   which is the condition for conserving \emph{first integral}. For a
612   canonical Hamiltonian system, the time evolution of an arbitrary
613   smooth function $G$ is given by,
647 –
614   \begin{eqnarray}
615   \frac{{dG(x(t))}}{{dt}} & = & [\nabla _x G(x(t))]^T \dot x(t) \\
616                          & = & [\nabla _x G(x(t))]^T J\nabla _x H(x(t)). \\
617   \label{introEquation:firstIntegral1}
618   \end{eqnarray}
653 –
654 –
619   Using poisson bracket notion, Equation
620   \ref{introEquation:firstIntegral1} can be rewritten as
621   \[
# Line 664 | Line 628 | is a \emph{first integral}, which is due to the fact $
628   \]
629   As well known, the Hamiltonian (or energy) H of a Hamiltonian system
630   is a \emph{first integral}, which is due to the fact $\{ H,H\}  =
631 < 0$.
668 <
669 < When designing any numerical methods, one should always try to
631 > 0$. When designing any numerical methods, one should always try to
632   preserve the structural properties of the original ODE and its flow.
633  
634   \subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
635   A lot of well established and very effective numerical methods have
636   been successful precisely because of their symplecticities even
637   though this fact was not recognized when they were first
638 < constructed. The most famous example is the Verlet-leapfrog methods
638 > constructed. The most famous example is the Verlet-leapfrog method
639   in molecular dynamics. In general, symplectic integrators can be
640   constructed using one of four different methods.
641   \begin{enumerate}
# Line 707 | Line 669 | simpler integration of the system.
669   \label{introEquation:FlowDecomposition}
670   \end{equation}
671   where each of the sub-flow is chosen such that each represent a
672 < simpler integration of the system.
673 <
712 < Suppose that a Hamiltonian system takes the form,
672 > simpler integration of the system. Suppose that a Hamiltonian system
673 > takes the form,
674   \[
675   H = H_1 + H_2.
676   \]
# Line 750 | Line 711 | to its symmetric property,
711   \begin{equation}
712   \varphi _h^{ - 1} = \varphi _{ - h}.
713   \label{introEquation:timeReversible}
714 < \end{equation},appendixFig:architecture
714 > \end{equation}
715  
716 < \subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Example of Splitting Method}}
716 > \subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Examples of the Splitting Method}}
717   The classical equation for a system consisting of interacting
718   particles can be written in Hamiltonian form,
719   \[
720   H = T + V
721   \]
722   where $T$ is the kinetic energy and $V$ is the potential energy.
723 < Setting $H_1 = T, H_2 = V$ and applying Strang splitting, one
723 > Setting $H_1 = T, H_2 = V$ and applying the Strang splitting, one
724   obtains the following:
725   \begin{align}
726   q(\Delta t) &= q(0) + \dot{q}(0)\Delta t +
# Line 786 | Line 747 | q(\Delta t) &= q(0) + \Delta t\, \dot{q}\biggl (\frac{
747      \label{introEquation:Lp9b}\\%
748   %
749   \dot{q}(\Delta t) &= \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) +
750 <    \frac{\Delta t}{2m}\, F[q(0)]. \label{introEquation:Lp9c}
750 >    \frac{\Delta t}{2m}\, F[q(t)]. \label{introEquation:Lp9c}
751   \end{align}
752   From the preceding splitting, one can see that the integration of
753   the equations of motion would follow:
# Line 795 | Line 756 | the equations of motion would follow:
756  
757   \item Use the half step velocities to move positions one whole step, $\Delta t$.
758  
759 < \item Evaluate the forces at the new positions, $\mathbf{r}(\Delta t)$, and use the new forces to complete the velocity move.
759 > \item Evaluate the forces at the new positions, $\mathbf{q}(\Delta t)$, and use the new forces to complete the velocity move.
760  
761   \item Repeat from step 1 with the new position, velocities, and forces assuming the roles of the initial values.
762   \end{enumerate}
763 <
764 < Simply switching the order of splitting and composing, a new
765 < integrator, the \emph{position verlet} integrator, can be generated,
763 > By simply switching the order of the propagators in the splitting
764 > and composing a new integrator, the \emph{position verlet}
765 > integrator, can be generated,
766   \begin{align}
767   \dot q(\Delta t) &= \dot q(0) + \Delta tF(q(0))\left[ {q(0) +
768   \frac{{\Delta t}}{{2m}}\dot q(0)} \right], %
# Line 814 | Line 775 | Baker-Campbell-Hausdorff formula can be used to determ
775  
776   \subsubsection{\label{introSection:errorAnalysis}\textbf{Error Analysis and Higher Order Methods}}
777  
778 < Baker-Campbell-Hausdorff formula can be used to determine the local
779 < error of splitting method in terms of commutator of the
778 > The Baker-Campbell-Hausdorff formula can be used to determine the
779 > local error of splitting method in terms of the commutator of the
780   operators(\ref{introEquation:exponentialOperator}) associated with
781 < the sub-flow. For operators $hX$ and $hY$ which are associate to
781 > the sub-flow. For operators $hX$ and $hY$ which are associated with
782   $\varphi_1(t)$ and $\varphi_2(t)$ respectively , we have
783   \begin{equation}
784   \exp (hX + hY) = \exp (hZ)
# Line 831 | Line 792 | Applying Baker-Campbell-Hausdorff formula\cite{Varadar
792   \[
793   [X,Y] = XY - YX .
794   \]
795 < Applying Baker-Campbell-Hausdorff formula\cite{Varadarajan1974} to
796 < Sprang splitting, we can obtain
795 > Applying the Baker-Campbell-Hausdorff formula\cite{Varadarajan1974}
796 > to the Sprang splitting, we can obtain
797   \begin{eqnarray*}
798   \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 \\
799                                     &   & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\
800                                     &   & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3[X,[X,Y]]/24 + \ldots )
801   \end{eqnarray*}
802 < Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0\], the dominant local
802 > Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0,\] the dominant local
803   error of Spring splitting is proportional to $h^3$. The same
804 < procedure can be applied to general splitting,  of the form
804 > procedure can be applied to a general splitting,  of the form
805   \begin{equation}
806   \varphi _{b_m h}^2  \circ \varphi _{a_m h}^1  \circ \varphi _{b_{m -
807   1} h}^2  \circ  \ldots  \circ \varphi _{a_1 h}^1 .
808   \end{equation}
809 < Careful choice of coefficient $a_1 \ldots b_m$ will lead to higher
810 < order method. Yoshida proposed an elegant way to compose higher
809 > A careful choice of coefficient $a_1 \ldots b_m$ will lead to higher
810 > order methods. Yoshida proposed an elegant way to compose higher
811   order methods based on symmetric splitting\cite{Yoshida1990}. Given
812   a symmetric second order base method $ \varphi _h^{(2)} $, a
813   fourth-order symmetric method can be constructed by composing,
# Line 859 | Line 820 | _{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)}
820   integrator $ \varphi _h^{(2n + 2)}$ can be composed by
821   \begin{equation}
822   \varphi _h^{(2n + 2)}  = \varphi _{\alpha h}^{(2n)}  \circ \varphi
823 < _{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)}
823 > _{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)},
824   \end{equation}
825 < , if the weights are chosen as
825 > if the weights are chosen as
826   \[
827   \alpha  =  - \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }},\beta =
828   \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }} .
# Line 899 | Line 860 | will discusses issues in production run.
860   These three individual steps will be covered in the following
861   sections. Sec.~\ref{introSec:initialSystemSettings} deals with the
862   initialization of a simulation. Sec.~\ref{introSection:production}
863 < will discusses issues in production run.
863 > will discusse issues in production run.
864   Sec.~\ref{introSection:Analysis} provides the theoretical tools for
865   trajectory analysis.
866  
# Line 912 | Line 873 | purification and crystallization. Even for the molecul
873   databases, such as RCSB Protein Data Bank \textit{etc}. Although
874   thousands of crystal structures of molecules are discovered every
875   year, many more remain unknown due to the difficulties of
876 < purification and crystallization. Even for the molecule with known
877 < structure, some important information is missing. For example, the
876 > purification and crystallization. Even for molecules with known
877 > structure, some important information is missing. For example, a
878   missing hydrogen atom which acts as donor in hydrogen bonding must
879   be added. Moreover, in order to include electrostatic interaction,
880   one may need to specify the partial charges for individual atoms.
881   Under some circumstances, we may even need to prepare the system in
882 < a special setup. For instance, when studying transport phenomenon in
883 < membrane system, we may prepare the lipids in bilayer structure
884 < instead of placing lipids randomly in solvent, since we are not
885 < interested in self-aggregation and it takes a long time to happen.
882 > a special configuration. For instance, when studying transport
883 > phenomenon in membrane systems, we may prepare the lipids in a
884 > bilayer structure instead of placing lipids randomly in solvent,
885 > since we are not interested in the slow self-aggregation process.
886  
887   \subsubsection{\textbf{Minimization}}
888  
889   It is quite possible that some of molecules in the system from
890 < preliminary preparation may be overlapped with each other. This
891 < close proximity leads to high potential energy which consequently
892 < jeopardizes any molecular dynamics simulations. To remove these
893 < steric overlaps, one typically performs energy minimization to find
894 < a more reasonable conformation. Several energy minimization methods
895 < have been developed to exploit the energy surface and to locate the
896 < local minimum. While converging slowly near the minimum, steepest
897 < descent method is extremely robust when systems are far from
898 < harmonic. Thus, it is often used to refine structure from
899 < crystallographic data. Relied on the gradient or hessian, advanced
900 < methods like conjugate gradient and Newton-Raphson converge rapidly
901 < to a local minimum, while become unstable if the energy surface is
902 < far from quadratic. Another factor must be taken into account, when
890 > preliminary preparation may be overlapping with each other. This
891 > close proximity leads to high initial potential energy which
892 > consequently jeopardizes any molecular dynamics simulations. To
893 > remove these steric overlaps, one typically performs energy
894 > minimization to find a more reasonable conformation. Several energy
895 > minimization methods have been developed to exploit the energy
896 > surface and to locate the local minimum. While converging slowly
897 > near the minimum, steepest descent method is extremely robust when
898 > systems are strongly anharmonic. Thus, it is often used to refine
899 > structure from crystallographic data. Relied on the gradient or
900 > hessian, advanced methods like Newton-Raphson converge rapidly to a
901 > local minimum, but become unstable if the energy surface is far from
902 > quadratic. Another factor that must be taken into account, when
903   choosing energy minimization method, is the size of the system.
904   Steepest descent and conjugate gradient can deal with models of any
905 < size. Because of the limit of computation power to calculate hessian
906 < matrix and insufficient storage capacity to store them, most
907 < Newton-Raphson methods can not be used with very large models.
905 > size. Because of the limits on computer memory to store the hessian
906 > matrix and the computing power needed to diagonalized these
907 > matrices, most Newton-Raphson methods can not be used with very
908 > large systems.
909  
910   \subsubsection{\textbf{Heating}}
911  
912   Typically, Heating is performed by assigning random velocities
913 < according to a Gaussian distribution for a temperature. Beginning at
914 < a lower temperature and gradually increasing the temperature by
915 < assigning greater random velocities, we end up with setting the
916 < temperature of the system to a final temperature at which the
917 < simulation will be conducted. In heating phase, we should also keep
918 < the system from drifting or rotating as a whole. Equivalently, the
919 < net linear momentum and angular momentum of the system should be
920 < shifted to zero.
913 > according to a Maxwell-Boltzman distribution for a desired
914 > temperature. Beginning at a lower temperature and gradually
915 > increasing the temperature by assigning larger random velocities, we
916 > end up with setting the temperature of the system to a final
917 > temperature at which the simulation will be conducted. In heating
918 > phase, we should also keep the system from drifting or rotating as a
919 > whole. To do this, the net linear momentum and angular momentum of
920 > the system is shifted to zero after each resampling from the Maxwell
921 > -Boltzman distribution.
922  
923   \subsubsection{\textbf{Equilibration}}
924  
# Line 971 | Line 934 | Production run is the most important step of the simul
934  
935   \subsection{\label{introSection:production}Production}
936  
937 < Production run is the most important step of the simulation, in
937 > The production run is the most important step of the simulation, in
938   which the equilibrated structure is used as a starting point and the
939   motions of the molecules are collected for later analysis. In order
940   to capture the macroscopic properties of the system, the molecular
941 < dynamics simulation must be performed in correct and efficient way.
941 > dynamics simulation must be performed by sampling correctly and
942 > efficiently from the relevant thermodynamic ensemble.
943  
944   The most expensive part of a molecular dynamics simulation is the
945   calculation of non-bonded forces, such as van der Waals force and
946   Coulombic forces \textit{etc}. For a system of $N$ particles, the
947   complexity of the algorithm for pair-wise interactions is $O(N^2 )$,
948   which making large simulations prohibitive in the absence of any
949 < computation saving techniques.
949 > algorithmic tricks.
950  
951 < A natural approach to avoid system size issue is to represent the
951 > A natural approach to avoid system size issues is to represent the
952   bulk behavior by a finite number of the particles. However, this
953 < approach will suffer from the surface effect. To offset this,
954 < \textit{Periodic boundary condition} (see Fig.~\ref{introFig:pbc})
955 < is developed to simulate bulk properties with a relatively small
956 < number of particles. In this method, the simulation box is
957 < replicated throughout space to form an infinite lattice. During the
958 < simulation, when a particle moves in the primary cell, its image in
959 < other cells move in exactly the same direction with exactly the same
960 < orientation. Thus, as a particle leaves the primary cell, one of its
961 < images will enter through the opposite face.
953 > approach will suffer from the surface effect at the edges of the
954 > simulation. To offset this, \textit{Periodic boundary conditions}
955 > (see Fig.~\ref{introFig:pbc}) is developed to simulate bulk
956 > properties with a relatively small number of particles. In this
957 > method, the simulation box is replicated throughout space to form an
958 > infinite lattice. During the simulation, when a particle moves in
959 > the primary cell, its image in other cells move in exactly the same
960 > direction with exactly the same orientation. Thus, as a particle
961 > leaves the primary cell, one of its images will enter through the
962 > opposite face.
963   \begin{figure}
964   \centering
965   \includegraphics[width=\linewidth]{pbc.eps}
# Line 1006 | Line 971 | evaluation is to apply cutoff where particles farther
971  
972   %cutoff and minimum image convention
973   Another important technique to improve the efficiency of force
974 < evaluation is to apply cutoff where particles farther than a
975 < predetermined distance, are not included in the calculation
974 > evaluation is to apply spherical cutoff where particles farther than
975 > a predetermined distance are not included in the calculation
976   \cite{Frenkel1996}. The use of a cutoff radius will cause a
977   discontinuity in the potential energy curve. Fortunately, one can
978 < shift the potential to ensure the potential curve go smoothly to
979 < zero at the cutoff radius. Cutoff strategy works pretty well for
980 < Lennard-Jones interaction because of its short range nature.
981 < However, simply truncating the electrostatic interaction with the
982 < use of cutoff has been shown to lead to severe artifacts in
983 < simulations. Ewald summation, in which the slowly conditionally
984 < convergent Coulomb potential is transformed into direct and
985 < reciprocal sums with rapid and absolute convergence, has proved to
986 < minimize the periodicity artifacts in liquid simulations. Taking the
987 < advantages of the fast Fourier transform (FFT) for calculating
988 < discrete Fourier transforms, the particle mesh-based
978 > shift simple radial potential to ensure the potential curve go
979 > smoothly to zero at the cutoff radius. The cutoff strategy works
980 > well for Lennard-Jones interaction because of its short range
981 > nature. However, simply truncating the electrostatic interaction
982 > with the use of cutoffs has been shown to lead to severe artifacts
983 > in simulations. The Ewald summation, in which the slowly decaying
984 > Coulomb potential is transformed into direct and reciprocal sums
985 > with rapid and absolute convergence, has proved to minimize the
986 > periodicity artifacts in liquid simulations. Taking the advantages
987 > of the fast Fourier transform (FFT) for calculating discrete Fourier
988 > transforms, the particle mesh-based
989   methods\cite{Hockney1981,Shimada1993, Luty1994} are accelerated from
990 < $O(N^{3/2})$ to $O(N logN)$. An alternative approach is \emph{fast
991 < multipole method}\cite{Greengard1987, Greengard1994}, which treats
992 < Coulombic interaction exactly at short range, and approximate the
993 < potential at long range through multipolar expansion. In spite of
994 < their wide acceptances at the molecular simulation community, these
995 < two methods are hard to be implemented correctly and efficiently.
996 < Instead, we use a damped and charge-neutralized Coulomb potential
997 < method developed by Wolf and his coworkers\cite{Wolf1999}. The
998 < shifted Coulomb potential for particle $i$ and particle $j$ at
999 < distance $r_{rj}$ is given by:
990 > $O(N^{3/2})$ to $O(N logN)$. An alternative approach is the
991 > \emph{fast multipole method}\cite{Greengard1987, Greengard1994},
992 > which treats Coulombic interactions exactly at short range, and
993 > approximate the potential at long range through multipolar
994 > expansion. In spite of their wide acceptance at the molecular
995 > simulation community, these two methods are difficult to implement
996 > correctly and efficiently. Instead, we use a damped and
997 > charge-neutralized Coulomb potential method developed by Wolf and
998 > his coworkers\cite{Wolf1999}. The shifted Coulomb potential for
999 > particle $i$ and particle $j$ at distance $r_{rj}$ is given by:
1000   \begin{equation}
1001   V(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha
1002   r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow
# Line 1053 | Line 1018 | Recently, advanced visualization technique are widely
1018  
1019   \subsection{\label{introSection:Analysis} Analysis}
1020  
1021 < Recently, advanced visualization technique are widely applied to
1021 > Recently, advanced visualization technique have become applied to
1022   monitor the motions of molecules. Although the dynamics of the
1023   system can be described qualitatively from animation, quantitative
1024 < trajectory analysis are more appreciable. According to the
1025 < principles of Statistical Mechanics,
1026 < Sec.~\ref{introSection:statisticalMechanics}, one can compute
1027 < thermodynamics properties, analyze fluctuations of structural
1028 < parameters, and investigate time-dependent processes of the molecule
1064 < from the trajectories.
1024 > trajectory analysis are more useful. According to the principles of
1025 > Statistical Mechanics, Sec.~\ref{introSection:statisticalMechanics},
1026 > one can compute thermodynamic properties, analyze fluctuations of
1027 > structural parameters, and investigate time-dependent processes of
1028 > the molecule from the trajectories.
1029  
1030 < \subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamics Properties}}
1030 > \subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamic Properties}}
1031  
1032 < Thermodynamics properties, which can be expressed in terms of some
1032 > Thermodynamic properties, which can be expressed in terms of some
1033   function of the coordinates and momenta of all particles in the
1034   system, can be directly computed from molecular dynamics. The usual
1035   way to measure the pressure is based on virial theorem of Clausius
# Line 1088 | Line 1052 | distribution functions. Among these functions,\emph{pa
1052   \subsubsection{\label{introSection:structuralProperties}\textbf{Structural Properties}}
1053  
1054   Structural Properties of a simple fluid can be described by a set of
1055 < distribution functions. Among these functions,\emph{pair
1055 > distribution functions. Among these functions,the \emph{pair
1056   distribution function}, also known as \emph{radial distribution
1057 < function}, is of most fundamental importance to liquid-state theory.
1058 < Pair distribution function can be gathered by Fourier transforming
1059 < raw data from a series of neutron diffraction experiments and
1060 < integrating over the surface factor \cite{Powles1973}. The
1061 < experiment result can serve as a criterion to justify the
1062 < correctness of the theory. Moreover, various equilibrium
1063 < thermodynamic and structural properties can also be expressed in
1064 < terms of radial distribution function \cite{Allen1987}.
1065 <
1102 < A pair distribution functions $g(r)$ gives the probability that a
1057 > function}, is of most fundamental importance to liquid theory.
1058 > Experimentally, pair distribution function can be gathered by
1059 > Fourier transforming raw data from a series of neutron diffraction
1060 > experiments and integrating over the surface factor
1061 > \cite{Powles1973}. The experimental results can serve as a criterion
1062 > to justify the correctness of a liquid model. Moreover, various
1063 > equilibrium thermodynamic and structural properties can also be
1064 > expressed in terms of radial distribution function \cite{Allen1987}.
1065 > The pair distribution functions $g(r)$ gives the probability that a
1066   particle $i$ will be located at a distance $r$ from a another
1067   particle $j$ in the system
1068   \[
1069   g(r) = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j
1070 < \ne i} {\delta (r - r_{ij} )} } } \right\rangle.
1070 > \ne i} {\delta (r - r_{ij} )} } } \right\rangle = \frac{\rho
1071 > (r)}{\rho}.
1072   \]
1073   Note that the delta function can be replaced by a histogram in
1074 < computer simulation. Figure
1075 < \ref{introFigure:pairDistributionFunction} shows a typical pair
1076 < distribution function for the liquid argon system. The occurrence of
1113 < several peaks in the plot of $g(r)$ suggests that it is more likely
1114 < to find particles at certain radial values than at others. This is a
1115 < result of the attractive interaction at such distances. Because of
1116 < the strong repulsive forces at short distance, the probability of
1117 < locating particles at distances less than about 2.5{\AA} from each
1118 < other is essentially zero.
1074 > computer simulation. Peaks in $g(r)$ represent solvent shells, and
1075 > the height of these peaks gradually decreases to 1 as the liquid of
1076 > large distance approaches the bulk density.
1077  
1120 – %\begin{figure}
1121 – %\centering
1122 – %\includegraphics[width=\linewidth]{pdf.eps}
1123 – %\caption[Pair distribution function for the liquid argon
1124 – %]{Pair distribution function for the liquid argon}
1125 – %\label{introFigure:pairDistributionFunction}
1126 – %\end{figure}
1078  
1079   \subsubsection{\label{introSection:timeDependentProperties}\textbf{Time-dependent
1080   Properties}}
1081  
1082   Time-dependent properties are usually calculated using \emph{time
1083 < correlation function}, which correlates random variables $A$ and $B$
1084 < at two different time
1083 > correlation functions}, which correlate random variables $A$ and $B$
1084 > at two different times,
1085   \begin{equation}
1086   C_{AB} (t) = \left\langle {A(t)B(0)} \right\rangle.
1087   \label{introEquation:timeCorrelationFunction}
1088   \end{equation}
1089   If $A$ and $B$ refer to same variable, this kind of correlation
1090 < function is called \emph{auto correlation function}. One example of
1091 < auto correlation function is velocity auto-correlation function
1092 < which is directly related to transport properties of molecular
1093 < liquids:
1090 > function is called an \emph{autocorrelation function}. One example
1091 > of an auto correlation function is the velocity auto-correlation
1092 > function which is directly related to transport properties of
1093 > molecular liquids:
1094   \[
1095   D = \frac{1}{3}\int\limits_0^\infty  {\left\langle {v(t) \cdot v(0)}
1096   \right\rangle } dt
1097   \]
1098 < where $D$ is diffusion constant. Unlike velocity autocorrelation
1099 < function which is averaging over time origins and over all the
1100 < atoms, dipole autocorrelation are calculated for the entire system.
1101 < The dipole autocorrelation function is given by:
1098 > where $D$ is diffusion constant. Unlike the velocity autocorrelation
1099 > function, which is averaging over time origins and over all the
1100 > atoms, the dipole autocorrelation functions are calculated for the
1101 > entire system. The dipole autocorrelation function is given by:
1102   \[
1103   c_{dipole}  = \left\langle {u_{tot} (t) \cdot u_{tot} (t)}
1104   \right\rangle
# Line 1173 | Line 1124 | simulator is governed by the rigid body dynamics. In m
1124   areas, from engineering, physics, to chemistry. For example,
1125   missiles and vehicle are usually modeled by rigid bodies.  The
1126   movement of the objects in 3D gaming engine or other physics
1127 < simulator is governed by the rigid body dynamics. In molecular
1128 < simulation, rigid body is used to simplify the model in
1129 < protein-protein docking study\cite{Gray2003}.
1127 > simulator is governed by rigid body dynamics. In molecular
1128 > simulations, rigid bodies are used to simplify protein-protein
1129 > docking studies\cite{Gray2003}.
1130  
1131   It is very important to develop stable and efficient methods to
1132 < integrate the equations of motion of orientational degrees of
1133 < freedom. Euler angles are the nature choice to describe the
1134 < rotational degrees of freedom. However, due to its singularity, the
1135 < numerical integration of corresponding equations of motion is very
1136 < inefficient and inaccurate. Although an alternative integrator using
1137 < different sets of Euler angles can overcome this
1138 < difficulty\cite{Barojas1973}, the computational penalty and the lost
1139 < of angular momentum conservation still remain. A singularity free
1140 < representation utilizing quaternions was developed by Evans in
1141 < 1977\cite{Evans1977}. Unfortunately, this approach suffer from the
1142 < nonseparable Hamiltonian resulted from quaternion representation,
1143 < which prevents the symplectic algorithm to be utilized. Another
1144 < different approach is to apply holonomic constraints to the atoms
1145 < belonging to the rigid body. Each atom moves independently under the
1146 < normal forces deriving from potential energy and constraint forces
1147 < which are used to guarantee the rigidness. However, due to their
1148 < iterative nature, SHAKE and Rattle algorithm converge very slowly
1149 < when the number of constraint increases\cite{Ryckaert1977,
1150 < Andersen1983}.
1132 > integrate the equations of motion for orientational degrees of
1133 > freedom. Euler angles are the natural choice to describe the
1134 > rotational degrees of freedom. However, due to $\frac {1}{sin
1135 > \theta}$ singularities, the numerical integration of corresponding
1136 > equations of motion is very inefficient and inaccurate. Although an
1137 > alternative integrator using multiple sets of Euler angles can
1138 > overcome this difficulty\cite{Barojas1973}, the computational
1139 > penalty and the loss of angular momentum conservation still remain.
1140 > A singularity-free representation utilizing quaternions was
1141 > developed by Evans in 1977\cite{Evans1977}. Unfortunately, this
1142 > approach uses a nonseparable Hamiltonian resulting from the
1143 > quaternion representation, which prevents the symplectic algorithm
1144 > to be utilized. Another different approach is to apply holonomic
1145 > constraints to the atoms belonging to the rigid body. Each atom
1146 > moves independently under the normal forces deriving from potential
1147 > energy and constraint forces which are used to guarantee the
1148 > rigidness. However, due to their iterative nature, the SHAKE and
1149 > Rattle algorithms also converge very slowly when the number of
1150 > constraints increases\cite{Ryckaert1977, Andersen1983}.
1151  
1152 < The break through in geometric literature suggests that, in order to
1152 > A break-through in geometric literature suggests that, in order to
1153   develop a long-term integration scheme, one should preserve the
1154 < symplectic structure of the flow. Introducing conjugate momentum to
1155 < rotation matrix $Q$ and re-formulating Hamiltonian's equation, a
1156 < symplectic integrator, RSHAKE\cite{Kol1997}, was proposed to evolve
1157 < the Hamiltonian system in a constraint manifold by iteratively
1158 < satisfying the orthogonality constraint $Q_T Q = 1$. An alternative
1159 < method using quaternion representation was developed by
1160 < Omelyan\cite{Omelyan1998}. However, both of these methods are
1161 < iterative and inefficient. In this section, we will present a
1154 > symplectic structure of the flow. By introducing a conjugate
1155 > momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's
1156 > equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was
1157 > proposed to evolve the Hamiltonian system in a constraint manifold
1158 > by iteratively satisfying the orthogonality constraint $Q^T Q = 1$.
1159 > An alternative method using the quaternion representation was
1160 > developed by Omelyan\cite{Omelyan1998}. However, both of these
1161 > methods are iterative and inefficient. In this section, we descibe a
1162   symplectic Lie-Poisson integrator for rigid body developed by
1163   Dullweber and his coworkers\cite{Dullweber1997} in depth.
1164  
1165 < \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
1166 < The motion of the rigid body is Hamiltonian with the Hamiltonian
1165 > \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Bodies}
1166 > The motion of a rigid body is Hamiltonian with the Hamiltonian
1167   function
1168   \begin{equation}
1169   H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) +
# Line 1226 | Line 1177 | constrained Hamiltonian equation subjects to a holonom
1177   I_{ii}^{ - 1}  = \frac{1}{2}\sum\limits_{i \ne j} {J_{jj}^{ - 1} }
1178   \]
1179   where $I_{ii}$ is the diagonal element of the inertia tensor. This
1180 < constrained Hamiltonian equation subjects to a holonomic constraint,
1180 > constrained Hamiltonian equation is subjected to a holonomic
1181 > constraint,
1182   \begin{equation}
1183   Q^T Q = 1, \label{introEquation:orthogonalConstraint}
1184   \end{equation}
1185 < which is used to ensure rotation matrix's orthogonality.
1186 < Differentiating \ref{introEquation:orthogonalConstraint} and using
1187 < Equation \ref{introEquation:RBMotionMomentum}, one may obtain,
1185 > which is used to ensure rotation matrix's unitarity. Differentiating
1186 > \ref{introEquation:orthogonalConstraint} and using Equation
1187 > \ref{introEquation:RBMotionMomentum}, one may obtain,
1188   \begin{equation}
1189   Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
1190   \label{introEquation:RBFirstOrderConstraint}
1191   \end{equation}
1240 –
1192   Using Equation (\ref{introEquation:motionHamiltonianCoordinate},
1193   \ref{introEquation:motionHamiltonianMomentum}), one can write down
1194   the equations of motion,
1244 –
1195   \begin{eqnarray}
1196   \frac{{dq}}{{dt}} & = & \frac{p}{m} \label{introEquation:RBMotionPosition}\\
1197   \frac{{dp}}{{dt}} & = & - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
1198   \frac{{dQ}}{{dt}} & = & PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
1199   \frac{{dP}}{{dt}} & = & - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}
1200   \end{eqnarray}
1251 –
1201   In general, there are two ways to satisfy the holonomic constraints.
1202 < We can use constraint force provided by lagrange multiplier on the
1203 < normal manifold to keep the motion on constraint space. Or we can
1204 < simply evolve the system in constraint manifold. These two methods
1205 < are proved to be equivalent. The holonomic constraint and equations
1206 < of motions define a constraint manifold for rigid body
1202 > We can use a constraint force provided by a Lagrange multiplier on
1203 > the normal manifold to keep the motion on constraint space. Or we
1204 > can simply evolve the system on the constraint manifold. These two
1205 > methods have been proved to be equivalent. The holonomic constraint
1206 > and equations of motions define a constraint manifold for rigid
1207 > bodies
1208   \[
1209   M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0}
1210   \right\}.
1211   \]
1262 –
1212   Unfortunately, this constraint manifold is not the cotangent bundle
1213 < $T_{\star}SO(3)$. However, it turns out that under symplectic
1213 > $T^* SO(3)$ which can be consider as a symplectic manifold on Lie
1214 > rotation group $SO(3)$. However, it turns out that under symplectic
1215   transformation, the cotangent space and the phase space are
1216 < diffeomorphic. Introducing
1216 > diffeomorphic. By introducing
1217   \[
1218   \tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
1219   \]
# Line 1273 | Line 1223 | T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \t
1223   T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \tilde Q =
1224   1,\tilde Q^T \tilde PJ^{ - 1}  + J^{ - 1} P^T \tilde Q = 0} \right\}
1225   \]
1276 –
1226   For a body fixed vector $X_i$ with respect to the center of mass of
1227   the rigid body, its corresponding lab fixed vector $X_0^{lab}$  is
1228   given as
# Line 1292 | Line 1241 | respectively.
1241   \[
1242   \nabla _Q V(q,Q) = F(q,Q)X_i^t
1243   \]
1244 < respectively.
1245 <
1246 < As a common choice to describe the rotation dynamics of the rigid
1298 < body, angular momentum on body frame $\Pi  = Q^t P$ is introduced to
1299 < rewrite the equations of motion,
1244 > respectively. As a common choice to describe the rotation dynamics
1245 > of the rigid body, the angular momentum on the body fixed frame $\Pi
1246 > = Q^t P$ is introduced to rewrite the equations of motion,
1247   \begin{equation}
1248   \begin{array}{l}
1249 < \mathop \Pi \limits^ \bullet   = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda  \\
1250 < \mathop Q\limits^{{\rm{   }} \bullet }  = Q\Pi {\rm{ }}J^{ - 1}  \\
1249 > \dot \Pi  = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda,  \\
1250 > \dot Q  = Q\Pi {\rm{ }}J^{ - 1},  \\
1251   \end{array}
1252   \label{introEqaution:RBMotionPI}
1253   \end{equation}
1254 < , as well as holonomic constraints,
1254 > as well as holonomic constraints,
1255   \[
1256   \begin{array}{l}
1257 < \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
1258 < Q^T Q = 1 \\
1257 > \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0, \\
1258 > Q^T Q = 1 .\\
1259   \end{array}
1260   \]
1314 –
1261   For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a matrix $\hat v \in
1262   so(3)^ \star$, the hat-map isomorphism,
1263   \begin{equation}
# Line 1326 | Line 1272 | operations
1272   will let us associate the matrix products with traditional vector
1273   operations
1274   \[
1275 < \hat vu = v \times u
1275 > \hat vu = v \times u.
1276   \]
1277 < Using \ref{introEqaution:RBMotionPI}, one can construct a skew
1277 > Using Eq.~\ref{introEqaution:RBMotionPI}, one can construct a skew
1278   matrix,
1279 + \begin{eqnarray}
1280 + (\dot \Pi  - \dot \Pi ^T ){\rm{ }} &= &{\rm{ }}(\Pi  - \Pi ^T ){\rm{
1281 + }}(J^{ - 1} \Pi  + \Pi J^{ - 1} ) \notag \\
1282 + + \sum\limits_i {[Q^T F_i
1283 + (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]}  - (\Lambda  - \Lambda ^T ).
1284 + \label{introEquation:skewMatrixPI}
1285 + \end{eqnarray}
1286 + Since $\Lambda$ is symmetric, the last term of
1287 + Eq.~\ref{introEquation:skewMatrixPI} is zero, which implies the
1288 + Lagrange multiplier $\Lambda$ is absent from the equations of
1289 + motion. This unique property eliminates the requirement of
1290 + iterations which can not be avoided in other methods\cite{Kol1997,
1291 + Omelyan1998}. Applying the hat-map isomorphism, we obtain the
1292 + equation of motion for angular momentum on body frame
1293   \begin{equation}
1334 – (\mathop \Pi \limits^ \bullet   - \mathop \Pi \limits^ {\bullet  ^T}
1335 – ){\rm{ }} = {\rm{ }}(\Pi  - \Pi ^T ){\rm{ }}(J^{ - 1} \Pi  + \Pi J^{
1336 – - 1} ) + \sum\limits_i {[Q^T F_i (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]} -
1337 – (\Lambda  - \Lambda ^T ) . \label{introEquation:skewMatrixPI}
1338 – \end{equation}
1339 – Since $\Lambda$ is symmetric, the last term of Equation
1340 – \ref{introEquation:skewMatrixPI} is zero, which implies the Lagrange
1341 – multiplier $\Lambda$ is absent from the equations of motion. This
1342 – unique property eliminate the requirement of iterations which can
1343 – not be avoided in other methods\cite{Kol1997, Omelyan1998}.
1344 –
1345 – Applying hat-map isomorphism, we obtain the equation of motion for
1346 – angular momentum on body frame
1347 – \begin{equation}
1294   \dot \pi  = \pi  \times I^{ - 1} \pi  + \sum\limits_i {\left( {Q^T
1295   F_i (r,Q)} \right) \times X_i }.
1296   \label{introEquation:bodyAngularMotion}
# Line 1352 | Line 1298 | given by
1298   In the same manner, the equation of motion for rotation matrix is
1299   given by
1300   \[
1301 < \dot Q = Qskew(I^{ - 1} \pi )
1301 > \dot Q = Qskew(I^{ - 1} \pi ).
1302   \]
1303  
1304   \subsection{\label{introSection:SymplecticFreeRB}Symplectic
1305   Lie-Poisson Integrator for Free Rigid Body}
1306  
1307 < If there is not external forces exerted on the rigid body, the only
1308 < contribution to the rotational is from the kinetic potential (the
1309 < first term of \ref{introEquation:bodyAngularMotion}). The free rigid
1310 < body is an example of Lie-Poisson system with Hamiltonian function
1307 > If there are no external forces exerted on the rigid body, the only
1308 > contribution to the rotational motion is from the kinetic energy
1309 > (the first term of \ref{introEquation:bodyAngularMotion}). The free
1310 > rigid body is an example of a Lie-Poisson system with Hamiltonian
1311 > function
1312   \begin{equation}
1313   T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
1314   \label{introEquation:rotationalKineticRB}
# Line 1373 | Line 1320 | J(\pi ) = \left( {\begin{array}{*{20}c}
1320     0 & {\pi _3 } & { - \pi _2 }  \\
1321     { - \pi _3 } & 0 & {\pi _1 }  \\
1322     {\pi _2 } & { - \pi _1 } & 0  \\
1323 < \end{array}} \right)
1323 > \end{array}} \right).
1324   \end{equation}
1325   Thus, the dynamics of free rigid body is governed by
1326   \begin{equation}
1327 < \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
1327 > \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi ).
1328   \end{equation}
1382 –
1329   One may notice that each $T_i^r$ in Equation
1330   \ref{introEquation:rotationalKineticRB} can be solved exactly. For
1331   instance, the equations of motion due to $T_1^r$ are given by
# Line 1408 | Line 1354 | tR_1 }$, we can use Cayley transformation,
1354   \end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
1355   \]
1356   To reduce the cost of computing expensive functions in $e^{\Delta
1357 < tR_1 }$, we can use Cayley transformation,
1357 > tR_1 }$, we can use Cayley transformation to obtain a single-aixs
1358 > propagator,
1359   \[
1360   e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
1361 < )
1361 > ).
1362   \]
1363   The flow maps for $T_2^r$ and $T_3^r$ can be found in the same
1364 < manner.
1365 <
1419 < In order to construct a second-order symplectic method, we split the
1420 < angular kinetic Hamiltonian function can into five terms
1364 > manner. In order to construct a second-order symplectic method, we
1365 > split the angular kinetic Hamiltonian function can into five terms
1366   \[
1367   T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
1368   ) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
1369 < (\pi _1 )
1370 < \].
1371 < Concatenating flows corresponding to these five terms, we can obtain
1372 < an symplectic integrator,
1369 > (\pi _1 ).
1370 > \]
1371 > By concatenating the propagators corresponding to these five terms,
1372 > we can obtain an symplectic integrator,
1373   \[
1374   \varphi _{\Delta t,T^r }  = \varphi _{\Delta t/2,\pi _1 }  \circ
1375   \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }
1376   \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
1377   _1 }.
1378   \]
1434 –
1379   The non-canonical Lie-Poisson bracket ${F, G}$ of two function
1380   $F(\pi )$ and $G(\pi )$ is defined by
1381   \[
1382   \{ F,G\} (\pi ) = [\nabla _\pi  F(\pi )]^T J(\pi )\nabla _\pi  G(\pi
1383 < )
1383 > ).
1384   \]
1385   If the Poisson bracket of a function $F$ with an arbitrary smooth
1386   function $G$ is zero, $F$ is a \emph{Casimir}, which is the
# Line 1447 | Line 1391 | then by the chain rule
1391   then by the chain rule
1392   \[
1393   \nabla _\pi  F(\pi ) = S'(\frac{{\parallel \pi \parallel ^2
1394 < }}{2})\pi
1394 > }}{2})\pi.
1395   \]
1396 < Thus $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel \pi
1396 > Thus, $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel
1397 > \pi
1398   \parallel ^2 }}{2})\pi  \times \pi  = 0 $. This explicit
1399 < Lie-Poisson integrator is found to be extremely efficient and stable
1400 < which can be explained by the fact the small angle approximation is
1401 < used and the norm of the angular momentum is conserved.
1399 > Lie-Poisson integrator is found to be both extremely efficient and
1400 > stable. These properties can be explained by the fact the small
1401 > angle approximation is used and the norm of the angular momentum is
1402 > conserved.
1403  
1404   \subsection{\label{introSection:RBHamiltonianSplitting} Hamiltonian
1405   Splitting for Rigid Body}
# Line 1461 | Line 1407 | H = T(p,\pi ) + V(q,Q)
1407   The Hamiltonian of rigid body can be separated in terms of kinetic
1408   energy and potential energy,
1409   \[
1410 < H = T(p,\pi ) + V(q,Q)
1410 > H = T(p,\pi ) + V(q,Q).
1411   \]
1412   The equations of motion corresponding to potential energy and
1413   kinetic energy are listed in the below table,
1414   \begin{table}
1415 < \caption{Equations of motion due to Potential and Kinetic Energies}
1415 > \caption{EQUATIONS OF MOTION DUE TO POTENTIAL AND KINETIC ENERGIES}
1416   \begin{center}
1417   \begin{tabular}{|l|l|}
1418    \hline
# Line 1480 | Line 1426 | A second-order symplectic method is now obtained by th
1426   \end{tabular}
1427   \end{center}
1428   \end{table}
1429 < A second-order symplectic method is now obtained by the
1430 < composition of the flow maps,
1429 > A second-order symplectic method is now obtained by the composition
1430 > of the position and velocity propagators,
1431   \[
1432   \varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
1433   _{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
1434   \]
1435   Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
1436 < sub-flows which corresponding to force and torque respectively,
1436 > sub-propagators which corresponding to force and torque
1437 > respectively,
1438   \[
1439   \varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
1440   _{\Delta t/2,\tau }.
1441   \]
1442   Since the associated operators of $\varphi _{\Delta t/2,F} $ and
1443 < $\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
1444 < order inside $\varphi _{\Delta t/2,V}$ does not matter.
1445 <
1446 < Furthermore, kinetic potential can be separated to translational
1500 < kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
1443 > $\circ \varphi _{\Delta t/2,\tau }$ commute, the composition order
1444 > inside $\varphi _{\Delta t/2,V}$ does not matter. Furthermore, the
1445 > kinetic energy can be separated to translational kinetic term, $T^t
1446 > (p)$, and rotational kinetic term, $T^r (\pi )$,
1447   \begin{equation}
1448   T(p,\pi ) =T^t (p) + T^r (\pi ).
1449   \end{equation}
1450   where $ T^t (p) = \frac{1}{2}p^T m^{ - 1} p $ and $T^r (\pi )$ is
1451   defined by \ref{introEquation:rotationalKineticRB}. Therefore, the
1452 < corresponding flow maps are given by
1452 > corresponding propagators are given by
1453   \[
1454   \varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
1455   _{\Delta t,T^r }.
1456   \]
1457 < Finally, we obtain the overall symplectic flow maps for free moving
1458 < rigid body
1459 < \begin{equation}
1460 < \begin{array}{c}
1461 < \varphi _{\Delta t}  = \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \\
1462 <  \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 }  \\
1517 <  \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1518 < \end{array}
1457 > Finally, we obtain the overall symplectic propagators for freely
1458 > moving rigid bodies
1459 > \begin{eqnarray*}
1460 > \varphi _{\Delta t}  &=& \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \\
1461 >  & & \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 }  \\
1462 >  & & \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1463   \label{introEquation:overallRBFlowMaps}
1464 < \end{equation}
1464 > \end{eqnarray*}
1465  
1466   \section{\label{introSection:langevinDynamics}Langevin Dynamics}
1467   As an alternative to newtonian dynamics, Langevin dynamics, which
1468   mimics a simple heat bath with stochastic and dissipative forces,
1469   has been applied in a variety of studies. This section will review
1470 < the theory of Langevin dynamics simulation. A brief derivation of
1471 < generalized Langevin equation will be given first. Follow that, we
1472 < will discuss the physical meaning of the terms appearing in the
1473 < equation as well as the calculation of friction tensor from
1474 < hydrodynamics theory.
1470 > the theory of Langevin dynamics. A brief derivation of generalized
1471 > Langevin equation will be given first. Following that, we will
1472 > discuss the physical meaning of the terms appearing in the equation
1473 > as well as the calculation of friction tensor from hydrodynamics
1474 > theory.
1475  
1476   \subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
1477  
1478 < Harmonic bath model, in which an effective set of harmonic
1478 > A harmonic bath model, in which an effective set of harmonic
1479   oscillators are used to mimic the effect of a linearly responding
1480   environment, has been widely used in quantum chemistry and
1481   statistical mechanics. One of the successful applications of
1482 < Harmonic bath model is the derivation of Deriving Generalized
1483 < Langevin Dynamics. Lets consider a system, in which the degree of
1482 > Harmonic bath model is the derivation of the Generalized Langevin
1483 > Dynamics (GLE). Lets consider a system, in which the degree of
1484   freedom $x$ is assumed to couple to the bath linearly, giving a
1485   Hamiltonian of the form
1486   \begin{equation}
1487   H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
1488   \label{introEquation:bathGLE}.
1489   \end{equation}
1490 < Here $p$ is a momentum conjugate to $q$, $m$ is the mass associated
1491 < with this degree of freedom, $H_B$ is harmonic bath Hamiltonian,
1490 > Here $p$ is a momentum conjugate to $x$, $m$ is the mass associated
1491 > with this degree of freedom, $H_B$ is a harmonic bath Hamiltonian,
1492   \[
1493   H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
1494   }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  \omega _\alpha ^2 }
# Line 1552 | Line 1496 | the harmonic bath masses, and $\Delta U$ is bilinear s
1496   \]
1497   where the index $\alpha$ runs over all the bath degrees of freedom,
1498   $\omega _\alpha$ are the harmonic bath frequencies, $m_\alpha$ are
1499 < the harmonic bath masses, and $\Delta U$ is bilinear system-bath
1499 > the harmonic bath masses, and $\Delta U$ is a bilinear system-bath
1500   coupling,
1501   \[
1502   \Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
1503   \]
1504 < where $g_\alpha$ are the coupling constants between the bath and the
1505 < coordinate $x$. Introducing
1504 > where $g_\alpha$ are the coupling constants between the bath
1505 > coordinates ($x_ \alpha$) and the system coordinate ($x$).
1506 > Introducing
1507   \[
1508   W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
1509   }}{{2m_\alpha  w_\alpha ^2 }}} x^2
1510 < \] and combining the last two terms in Equation
1511 < \ref{introEquation:bathGLE}, we may rewrite the Harmonic bath
1567 < Hamiltonian as
1510 > \]
1511 > and combining the last two terms in Eq.~\ref{introEquation:bathGLE}, we may rewrite the Harmonic bath Hamiltonian as
1512   \[
1513   H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
1514   {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
1515   w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
1516 < w_\alpha ^2 }}x} \right)^2 } \right\}}
1516 > w_\alpha ^2 }}x} \right)^2 } \right\}}.
1517   \]
1518   Since the first two terms of the new Hamiltonian depend only on the
1519   system coordinates, we can get the equations of motion for
1520 < Generalized Langevin Dynamics by Hamilton's equations
1577 < \ref{introEquation:motionHamiltonianCoordinate,
1578 < introEquation:motionHamiltonianMomentum},
1520 > Generalized Langevin Dynamics by Hamilton's equations,
1521   \begin{equation}
1522   m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
1523   \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   -
# Line 1588 | Line 1530 | m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x
1530   \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right).
1531   \label{introEquation:bathMotionGLE}
1532   \end{equation}
1591 –
1533   In order to derive an equation for $x$, the dynamics of the bath
1534   variables $x_\alpha$ must be solved exactly first. As an integral
1535   transform which is particularly useful in solving linear ordinary
1536 < differential equations, Laplace transform is the appropriate tool to
1537 < solve this problem. The basic idea is to transform the difficult
1536 > differential equations,the Laplace transform is the appropriate tool
1537 > to solve this problem. The basic idea is to transform the difficult
1538   differential equations into simple algebra problems which can be
1539 < solved easily. Then applying inverse Laplace transform, also known
1540 < as the Bromwich integral, we can retrieve the solutions of the
1541 < original problems.
1542 <
1602 < Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace
1603 < transform of f(t) is a new function defined as
1539 > solved easily. Then, by applying the inverse Laplace transform, also
1540 > known as the Bromwich integral, we can retrieve the solutions of the
1541 > original problems. Let $f(t)$ be a function defined on $ [0,\infty )
1542 > $. The Laplace transform of f(t) is a new function defined as
1543   \[
1544   L(f(t)) \equiv F(p) = \int_0^\infty  {f(t)e^{ - pt} dt}
1545   \]
1546   where  $p$ is real and  $L$ is called the Laplace Transform
1547   Operator. Below are some important properties of Laplace transform
1609 –
1548   \begin{eqnarray*}
1549   L(x + y)  & = & L(x) + L(y) \\
1550   L(ax)     & = & aL(x) \\
# Line 1614 | Line 1552 | Operator. Below are some important properties of Lapla
1552   L(\ddot x)& = & p^2 L(x) - px(0) - \dot x(0) \\
1553   L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right)& = & G(p)H(p) \\
1554   \end{eqnarray*}
1555 <
1618 <
1619 < Applying Laplace transform to the bath coordinates, we obtain
1555 > Applying the Laplace transform to the bath coordinates, we obtain
1556   \begin{eqnarray*}
1557   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) \\
1558   L(x_\alpha  ) & = & \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
1559   \end{eqnarray*}
1624 –
1560   By the same way, the system coordinates become
1561   \begin{eqnarray*}
1562 < mL(\ddot x) & = & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
1563 <  & & \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\}}  \\
1562 > mL(\ddot x) & = &
1563 >  - \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\}}  \\
1564 >  & & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}}
1565   \end{eqnarray*}
1630 –
1566   With the help of some relatively important inverse Laplace
1567   transformations:
1568   \[
# Line 1637 | Line 1572 | transformations:
1572   L(1) = \frac{1}{p} \\
1573   \end{array}
1574   \]
1575 < , we obtain
1575 > we obtain
1576   \begin{eqnarray*}
1577   m\ddot x & =  & - \frac{{\partial W(x)}}{{\partial x}} -
1578   \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
# Line 1693 | Line 1628 | as the model, which is gaussian distribution in genera
1628   \end{array}
1629   \]
1630   This property is what we expect from a truly random process. As long
1631 < as the model, which is gaussian distribution in general, chosen for
1632 < $R(t)$ is a truly random process, the stochastic nature of the GLE
1698 < still remains.
1631 > as the model chosen for $R(t)$ was a gaussian distribution in
1632 > general, the stochastic nature of the GLE still remains.
1633  
1634   %dynamic friction kernel
1635   The convolution integral
# Line 1711 | Line 1645 | and Equation \ref{introEuqation:GeneralizedLangevinDyn
1645   \[
1646   \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
1647   \]
1648 < and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1648 > and Eq.~\ref{introEuqation:GeneralizedLangevinDynamics} becomes
1649   \[
1650   m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
1651   \frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
1652   \]
1653 < which can be used to describe dynamic caging effect. The other
1654 < extreme is the bath that responds infinitely quickly to motions in
1655 < the system. Thus, $\xi (t)$ can be taken as a $delta$ function in
1656 < time:
1653 > which can be used to describe the effect of dynamic caging in
1654 > viscous solvents. The other extreme is the bath that responds
1655 > infinitely quickly to motions in the system. Thus, $\xi (t)$ can be
1656 > taken as a $delta$ function in time:
1657   \[
1658   \xi (t) = 2\xi _0 \delta (t)
1659   \]
# Line 1728 | Line 1662 | and Equation \ref{introEuqation:GeneralizedLangevinDyn
1662   \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
1663   {\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
1664   \]
1665 < and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1665 > and Eq.~\ref{introEuqation:GeneralizedLangevinDynamics} becomes
1666   \begin{equation}
1667   m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \xi _0 \dot
1668   x(t) + R(t) \label{introEquation:LangevinEquation}
1669   \end{equation}
1670   which is known as the Langevin equation. The static friction
1671   coefficient $\xi _0$ can either be calculated from spectral density
1672 < or be determined by Stokes' law for regular shaped particles.A
1672 > or be determined by Stokes' law for regular shaped particles. A
1673   briefly review on calculating friction tensor for arbitrary shaped
1674   particles is given in Sec.~\ref{introSection:frictionTensor}.
1675  
# Line 1751 | Line 1685 | And since the $q$ coordinates are harmonic oscillators
1685   R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
1686   \]
1687   And since the $q$ coordinates are harmonic oscillators,
1754 –
1688   \begin{eqnarray*}
1689   \left\langle {q_\alpha ^2 } \right\rangle  & = & \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
1690   \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle & = & \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
# Line 1760 | Line 1693 | And since the $q$ coordinates are harmonic oscillators
1693    & = &\sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
1694    & = &kT\xi (t) \\
1695   \end{eqnarray*}
1763 –
1696   Thus, we recover the \emph{second fluctuation dissipation theorem}
1697   \begin{equation}
1698   \xi (t) = \left\langle {R(t)R(0)} \right\rangle
# Line 1768 | Line 1700 | can model the random force and friction kernel.
1700   \end{equation}
1701   In effect, it acts as a constraint on the possible ways in which one
1702   can model the random force and friction kernel.
1771 –
1772 – \subsection{\label{introSection:frictionTensor} Friction Tensor}
1773 – Theoretically, the friction kernel can be determined using velocity
1774 – autocorrelation function. However, this approach become impractical
1775 – when the system become more and more complicate. Instead, various
1776 – approaches based on hydrodynamics have been developed to calculate
1777 – the friction coefficients. The friction effect is isotropic in
1778 – Equation, $\zeta$ can be taken as a scalar. In general, friction
1779 – tensor $\Xi$ is a $6\times 6$ matrix given by
1780 – \[
1781 – \Xi  = \left( {\begin{array}{*{20}c}
1782 –   {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
1783 –   {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
1784 – \end{array}} \right).
1785 – \]
1786 – Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are translational friction
1787 – tensor and rotational resistance (friction) tensor respectively,
1788 – while ${\Xi^{tr} }$ is translation-rotation coupling tensor and $
1789 – {\Xi^{rt} }$ is rotation-translation coupling tensor. When a
1790 – particle moves in a fluid, it may experience friction force or
1791 – torque along the opposite direction of the velocity or angular
1792 – velocity,
1793 – \[
1794 – \left( \begin{array}{l}
1795 – F_R  \\
1796 – \tau _R  \\
1797 – \end{array} \right) =  - \left( {\begin{array}{*{20}c}
1798 –   {\Xi ^{tt} } & {\Xi ^{rt} }  \\
1799 –   {\Xi ^{tr} } & {\Xi ^{rr} }  \\
1800 – \end{array}} \right)\left( \begin{array}{l}
1801 – v \\
1802 – w \\
1803 – \end{array} \right)
1804 – \]
1805 – where $F_r$ is the friction force and $\tau _R$ is the friction
1806 – toque.
1807 –
1808 – \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}}
1809 –
1810 – For a spherical particle, the translational and rotational friction
1811 – constant can be calculated from Stoke's law,
1812 – \[
1813 – \Xi ^{tt}  = \left( {\begin{array}{*{20}c}
1814 –   {6\pi \eta R} & 0 & 0  \\
1815 –   0 & {6\pi \eta R} & 0  \\
1816 –   0 & 0 & {6\pi \eta R}  \\
1817 – \end{array}} \right)
1818 – \]
1819 – and
1820 – \[
1821 – \Xi ^{rr}  = \left( {\begin{array}{*{20}c}
1822 –   {8\pi \eta R^3 } & 0 & 0  \\
1823 –   0 & {8\pi \eta R^3 } & 0  \\
1824 –   0 & 0 & {8\pi \eta R^3 }  \\
1825 – \end{array}} \right)
1826 – \]
1827 – where $\eta$ is the viscosity of the solvent and $R$ is the
1828 – hydrodynamics radius.
1829 –
1830 – Other non-spherical shape, such as cylinder and ellipsoid
1831 – \textit{etc}, are widely used as reference for developing new
1832 – hydrodynamics theory, because their properties can be calculated
1833 – exactly. In 1936, Perrin extended Stokes's law to general ellipsoid,
1834 – also called a triaxial ellipsoid, which is given in Cartesian
1835 – coordinates by\cite{Perrin1934, Perrin1936}
1836 – \[
1837 – \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
1838 – }} = 1
1839 – \]
1840 – where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
1841 – due to the complexity of the elliptic integral, only the ellipsoid
1842 – with the restriction of two axes having to be equal, \textit{i.e.}
1843 – prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
1844 – exactly. Introducing an elliptic integral parameter $S$ for prolate,
1845 – \[
1846 – S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
1847 – } }}{b},
1848 – \]
1849 – and oblate,
1850 – \[
1851 – S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
1852 – }}{a}
1853 – \],
1854 – one can write down the translational and rotational resistance
1855 – tensors
1856 – \[
1857 – \begin{array}{l}
1858 – \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}} \\
1859 – \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S + 2a}} \\
1860 – \end{array},
1861 – \]
1862 – and
1863 – \[
1864 – \begin{array}{l}
1865 – \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}} \\
1866 – \Xi _b^{rr}  = \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}} \\
1867 – \end{array}.
1868 – \]
1869 –
1870 – \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shape}}
1871 –
1872 – Unlike spherical and other regular shaped molecules, there is not
1873 – analytical solution for friction tensor of any arbitrary shaped
1874 – rigid molecules. The ellipsoid of revolution model and general
1875 – triaxial ellipsoid model have been used to approximate the
1876 – hydrodynamic properties of rigid bodies. However, since the mapping
1877 – from all possible ellipsoidal space, $r$-space, to all possible
1878 – combination of rotational diffusion coefficients, $D$-space is not
1879 – unique\cite{Wegener1979} as well as the intrinsic coupling between
1880 – translational and rotational motion of rigid body, general ellipsoid
1881 – is not always suitable for modeling arbitrarily shaped rigid
1882 – molecule. A number of studies have been devoted to determine the
1883 – friction tensor for irregularly shaped rigid bodies using more
1884 – advanced method where the molecule of interest was modeled by
1885 – combinations of spheres(beads)\cite{Carrasco1999} and the
1886 – hydrodynamics properties of the molecule can be calculated using the
1887 – hydrodynamic interaction tensor. Let us consider a rigid assembly of
1888 – $N$ beads immersed in a continuous medium. Due to hydrodynamics
1889 – interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different
1890 – than its unperturbed velocity $v_i$,
1891 – \[
1892 – v'_i  = v_i  - \sum\limits_{j \ne i} {T_{ij} F_j }
1893 – \]
1894 – where $F_i$ is the frictional force, and $T_{ij}$ is the
1895 – hydrodynamic interaction tensor. The friction force of $i$th bead is
1896 – proportional to its ``net'' velocity
1897 – \begin{equation}
1898 – F_i  = \zeta _i v_i  - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
1899 – \label{introEquation:tensorExpression}
1900 – \end{equation}
1901 – This equation is the basis for deriving the hydrodynamic tensor. In
1902 – 1930, Oseen and Burgers gave a simple solution to Equation
1903 – \ref{introEquation:tensorExpression}
1904 – \begin{equation}
1905 – T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
1906 – R_{ij}^T }}{{R_{ij}^2 }}} \right).
1907 – \label{introEquation:oseenTensor}
1908 – \end{equation}
1909 – Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$.
1910 – A second order expression for element of different size was
1911 – introduced by Rotne and Prager\cite{Rotne1969} and improved by
1912 – Garc\'{i}a de la Torre and Bloomfield\cite{Torre1977},
1913 – \begin{equation}
1914 – T_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I +
1915 – \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma
1916 – _i^2  + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} -
1917 – \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
1918 – \label{introEquation:RPTensorNonOverlapped}
1919 – \end{equation}
1920 – Both of the Equation \ref{introEquation:oseenTensor} and Equation
1921 – \ref{introEquation:RPTensorNonOverlapped} have an assumption $R_{ij}
1922 – \ge \sigma _i  + \sigma _j$. An alternative expression for
1923 – overlapping beads with the same radius, $\sigma$, is given by
1924 – \begin{equation}
1925 – T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
1926 – \frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I +
1927 – \frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
1928 – \label{introEquation:RPTensorOverlapped}
1929 – \end{equation}
1930 –
1931 – To calculate the resistance tensor at an arbitrary origin $O$, we
1932 – construct a $3N \times 3N$ matrix consisting of $N \times N$
1933 – $B_{ij}$ blocks
1934 – \begin{equation}
1935 – B = \left( {\begin{array}{*{20}c}
1936 –   {B_{11} } &  \ldots  & {B_{1N} }  \\
1937 –    \vdots  &  \ddots  &  \vdots   \\
1938 –   {B_{N1} } &  \cdots  & {B_{NN} }  \\
1939 – \end{array}} \right),
1940 – \end{equation}
1941 – where $B_{ij}$ is given by
1942 – \[
1943 – B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
1944 – )T_{ij}
1945 – \]
1946 – where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
1947 – $B$, we obtain
1948 –
1949 – \[
1950 – C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
1951 –   {C_{11} } &  \ldots  & {C_{1N} }  \\
1952 –    \vdots  &  \ddots  &  \vdots   \\
1953 –   {C_{N1} } &  \cdots  & {C_{NN} }  \\
1954 – \end{array}} \right)
1955 – \]
1956 – , which can be partitioned into $N \times N$ $3 \times 3$ block
1957 – $C_{ij}$. With the help of $C_{ij}$ and skew matrix $U_i$
1958 – \[
1959 – U_i  = \left( {\begin{array}{*{20}c}
1960 –   0 & { - z_i } & {y_i }  \\
1961 –   {z_i } & 0 & { - x_i }  \\
1962 –   { - y_i } & {x_i } & 0  \\
1963 – \end{array}} \right)
1964 – \]
1965 – where $x_i$, $y_i$, $z_i$ are the components of the vector joining
1966 – bead $i$ and origin $O$. Hence, the elements of resistance tensor at
1967 – arbitrary origin $O$ can be written as
1968 – \begin{equation}
1969 – \begin{array}{l}
1970 – \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\
1971 – \Xi _{}^{tr}  = \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
1972 – \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j  \\
1973 – \end{array}
1974 – \label{introEquation:ResistanceTensorArbitraryOrigin}
1975 – \end{equation}
1976 –
1977 – The resistance tensor depends on the origin to which they refer. The
1978 – proper location for applying friction force is the center of
1979 – resistance (reaction), at which the trace of rotational resistance
1980 – tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
1981 – resistance is defined as an unique point of the rigid body at which
1982 – the translation-rotation coupling tensor are symmetric,
1983 – \begin{equation}
1984 – \Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
1985 – \label{introEquation:definitionCR}
1986 – \end{equation}
1987 – Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
1988 – we can easily find out that the translational resistance tensor is
1989 – origin independent, while the rotational resistance tensor and
1990 – translation-rotation coupling resistance tensor depend on the
1991 – origin. Given resistance tensor at an arbitrary origin $O$, and a
1992 – vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
1993 – obtain the resistance tensor at $P$ by
1994 – \begin{equation}
1995 – \begin{array}{l}
1996 – \Xi _P^{tt}  = \Xi _O^{tt}  \\
1997 – \Xi _P^{tr}  = \Xi _P^{rt}  = \Xi _O^{tr}  - U_{OP} \Xi _O^{tt}  \\
1998 – \Xi _P^{rr}  = \Xi _O^{rr}  - U_{OP} \Xi _O^{tt} U_{OP}  + \Xi _O^{tr} U_{OP}  - U_{OP} \Xi _O^{{tr} ^{^T }}  \\
1999 – \end{array}
2000 – \label{introEquation:resistanceTensorTransformation}
2001 – \end{equation}
2002 – where
2003 – \[
2004 – U_{OP}  = \left( {\begin{array}{*{20}c}
2005 –   0 & { - z_{OP} } & {y_{OP} }  \\
2006 –   {z_i } & 0 & { - x_{OP} }  \\
2007 –   { - y_{OP} } & {x_{OP} } & 0  \\
2008 – \end{array}} \right)
2009 – \]
2010 – Using Equations \ref{introEquation:definitionCR} and
2011 – \ref{introEquation:resistanceTensorTransformation}, one can locate
2012 – the position of center of resistance,
2013 – \begin{eqnarray*}
2014 – \left( \begin{array}{l}
2015 – x_{OR}  \\
2016 – y_{OR}  \\
2017 – z_{OR}  \\
2018 – \end{array} \right) & = &\left( {\begin{array}{*{20}c}
2019 –   {(\Xi _O^{rr} )_{yy}  + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} }  \\
2020 –   { - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz}  + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} }  \\
2021 –   { - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx}  + (\Xi _O^{rr} )_{yy} }  \\
2022 – \end{array}} \right)^{ - 1}  \\
2023 –  & & \left( \begin{array}{l}
2024 – (\Xi _O^{tr} )_{yz}  - (\Xi _O^{tr} )_{zy}  \\
2025 – (\Xi _O^{tr} )_{zx}  - (\Xi _O^{tr} )_{xz}  \\
2026 – (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
2027 – \end{array} \right) \\
2028 – \end{eqnarray*}
2029 –
2030 –
2031 –
2032 – where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
2033 – joining center of resistance $R$ and origin $O$.

Diff Legend

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