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 2776 by tim, Thu May 25 21:32:14 2006 UTC vs.
Revision 2819 by tim, Wed Jun 7 21:03:46 2006 UTC

# Line 6 | Line 6 | behind classical mechanics. Firstly, One can determine
6   Closely related to Classical Mechanics, Molecular Dynamics
7   simulations are carried out by integrating the equations of motion
8   for a given system of particles. There are three fundamental ideas
9 < behind classical mechanics. Firstly, One can determine the state of
9 > behind classical mechanics. Firstly, one can determine the state of
10   a mechanical system at any time of interest; Secondly, all the
11   mechanical properties of the system at that time can be determined
12   by combining the knowledge of the properties of the system with the
# Line 17 | Line 17 | Newton¡¯s first law defines a class of inertial frames
17   \subsection{\label{introSection:newtonian}Newtonian Mechanics}
18   The discovery of Newton's three laws of mechanics which govern the
19   motion of particles is the foundation of the classical mechanics.
20 < Newton¡¯s first law defines a class of inertial frames. Inertial
20 > Newton's first law defines a class of inertial frames. Inertial
21   frames are reference frames where a particle not interacting with
22   other bodies will move with constant speed in the same direction.
23 < With respect to inertial frames Newton¡¯s second law has the form
23 > With respect to inertial frames, Newton's second law has the form
24   \begin{equation}
25 < F = \frac {dp}{dt} = \frac {mv}{dt}
25 > F = \frac {dp}{dt} = \frac {mdv}{dt}
26   \label{introEquation:newtonSecondLaw}
27   \end{equation}
28   A point mass interacting with other bodies moves with the
29   acceleration along the direction of the force acting on it. Let
30   $F_{ij}$ be the force that particle $i$ exerts on particle $j$, and
31   $F_{ji}$ be the force that particle $j$ exerts on particle $i$.
32 < Newton¡¯s third law states that
32 > Newton's third law states that
33   \begin{equation}
34   F_{ij} = -F_{ji}
35   \label{introEquation:newtonThirdLaw}
# Line 46 | Line 46 | N \equiv r \times F \label{introEquation:torqueDefinit
46   \end{equation}
47   The torque $\tau$ with respect to the same origin is defined to be
48   \begin{equation}
49 < N \equiv r \times F \label{introEquation:torqueDefinition}
49 > \tau \equiv r \times F \label{introEquation:torqueDefinition}
50   \end{equation}
51   Differentiating Eq.~\ref{introEquation:angularMomentumDefinition},
52   \[
# Line 59 | Line 59 | thus,
59   \]
60   thus,
61   \begin{equation}
62 < \dot L = r \times \dot p = N
62 > \dot L = r \times \dot p = \tau
63   \end{equation}
64   If there are no external torques acting on a body, the angular
65   momentum of it is conserved. The last conservation theorem state
# Line 68 | Line 68 | scheme for rigid body \cite{Dullweber1997}.
68   \end{equation}
69   is conserved. All of these conserved quantities are
70   important factors to determine the quality of numerical integration
71 < scheme for rigid body \cite{Dullweber1997}.
71 > schemes for rigid bodies \cite{Dullweber1997}.
72  
73   \subsection{\label{introSection:lagrangian}Lagrangian Mechanics}
74  
75 < Newtonian Mechanics suffers from two important limitations: it
76 < describes their motion in special cartesian coordinate systems.
77 < Another limitation of Newtonian mechanics becomes obvious when we
78 < try to describe systems with large numbers of particles. It becomes
79 < very difficult to predict the properties of the system by carrying
80 < out calculations involving the each individual interaction between
81 < all the particles, even if we know all of the details of the
82 < interaction. In order to overcome some of the practical difficulties
83 < which arise in attempts to apply Newton's equation to complex
84 < system, alternative procedures may be developed.
75 > Newtonian Mechanics suffers from two important limitations: motions
76 > can only be described in cartesian coordinate systems. Moreover, It
77 > become impossible to predict analytically the properties of the
78 > system even if we know all of the details of the interaction. In
79 > order to overcome some of the practical difficulties which arise in
80 > attempts to apply Newton's equation to complex system, approximate
81 > numerical procedures may be developed.
82  
83 < \subsubsection{\label{introSection:halmiltonPrinciple}Hamilton's
84 < Principle}
83 > \subsubsection{\label{introSection:halmiltonPrinciple}\textbf{Hamilton's
84 > Principle}}
85  
86   Hamilton introduced the dynamical principle upon which it is
87 < possible to base all of mechanics and, indeed, most of classical
88 < physics. Hamilton's Principle may be stated as follow,
87 > possible to base all of mechanics and most of classical physics.
88 > Hamilton's Principle may be stated as follows,
89  
90   The actual trajectory, along which a dynamical system may move from
91   one point to another within a specified time, is derived by finding
92   the path which minimizes the time integral of the difference between
93 < the kinetic, $K$, and potential energies, $U$ \cite{tolman79}.
93 > the kinetic, $K$, and potential energies, $U$.
94   \begin{equation}
95   \delta \int_{t_1 }^{t_2 } {(K - U)dt = 0} ,
96   \label{introEquation:halmitonianPrinciple1}
97   \end{equation}
98  
99   For simple mechanical systems, where the forces acting on the
100 < different part are derivable from a potential and the velocities are
101 < small compared with that of light, the Lagrangian function $L$ can
102 < be define as the difference between the kinetic energy of the system
106 < and its potential energy,
100 > different parts are derivable from a potential, the Lagrangian
101 > function $L$ can be defined as the difference between the kinetic
102 > energy of the system and its potential energy,
103   \begin{equation}
104   L \equiv K - U = L(q_i ,\dot q_i ) ,
105   \label{introEquation:lagrangianDef}
# Line 114 | Line 110 | then Eq.~\ref{introEquation:halmitonianPrinciple1} bec
110   \label{introEquation:halmitonianPrinciple2}
111   \end{equation}
112  
113 < \subsubsection{\label{introSection:equationOfMotionLagrangian}The
114 < Equations of Motion in Lagrangian Mechanics}
113 > \subsubsection{\label{introSection:equationOfMotionLagrangian}\textbf{The
114 > Equations of Motion in Lagrangian Mechanics}}
115  
116   For a holonomic system of $f$ degrees of freedom, the equations of
117   motion in the Lagrangian form is
# Line 132 | Line 128 | independent of generalized velocities, the generalized
128   Arising from Lagrangian Mechanics, Hamiltonian Mechanics was
129   introduced by William Rowan Hamilton in 1833 as a re-formulation of
130   classical mechanics. If the potential energy of a system is
131 < independent of generalized velocities, the generalized momenta can
136 < be defined as
131 > independent of velocities, the momenta can be defined as
132   \begin{equation}
133   p_i = \frac{\partial L}{\partial \dot q_i}
134   \label{introEquation:generalizedMomenta}
# Line 172 | Line 167 | find
167   By identifying the coefficients of $dq_k$, $dp_k$ and dt, we can
168   find
169   \begin{equation}
170 < \frac{{\partial H}}{{\partial p_k }} = q_k
170 > \frac{{\partial H}}{{\partial p_k }} = \dot {q_k}
171   \label{introEquation:motionHamiltonianCoordinate}
172   \end{equation}
173   \begin{equation}
174 < \frac{{\partial H}}{{\partial q_k }} =  - p_k
174 > \frac{{\partial H}}{{\partial q_k }} =  - \dot {p_k}
175   \label{introEquation:motionHamiltonianMomentum}
176   \end{equation}
177   and
# Line 189 | Line 184 | known as the canonical equations of motions \cite{Gold
184   Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
185   Eq.~\ref{introEquation:motionHamiltonianMomentum} are Hamilton's
186   equation of motion. Due to their symmetrical formula, they are also
187 < known as the canonical equations of motions \cite{Goldstein01}.
187 > known as the canonical equations of motions \cite{Goldstein2001}.
188  
189   An important difference between Lagrangian approach and the
190   Hamiltonian approach is that the Lagrangian is considered to be a
191 < function of the generalized velocities $\dot q_i$ and the
192 < generalized coordinates $q_i$, while the Hamiltonian is considered
193 < to be a function of the generalized momenta $p_i$ and the conjugate
194 < generalized coordinate $q_i$. Hamiltonian Mechanics is more
195 < appropriate for application to statistical mechanics and quantum
196 < mechanics, since it treats the coordinate and its time derivative as
197 < independent variables and it only works with 1st-order differential
203 < equations\cite{Marion90}.
191 > function of the generalized velocities $\dot q_i$ and coordinates
192 > $q_i$, while the Hamiltonian is considered to be a function of the
193 > generalized momenta $p_i$ and the conjugate coordinates $q_i$.
194 > Hamiltonian Mechanics is more appropriate for application to
195 > statistical mechanics and quantum mechanics, since it treats the
196 > coordinate and its time derivative as independent variables and it
197 > only works with 1st-order differential equations\cite{Marion1990}.
198  
199   In Newtonian Mechanics, a system described by conservative forces
200   conserves the total energy \ref{introEquation:energyConservation}.
# Line 230 | Line 224 | momentum variables. Consider a dynamic system in a car
224   possible states. Each possible state of the system corresponds to
225   one unique point in the phase space. For mechanical systems, the
226   phase space usually consists of all possible values of position and
227 < momentum variables. Consider a dynamic system in a cartesian space,
228 < where each of the $6f$ coordinates and momenta is assigned to one of
229 < $6f$ mutually orthogonal axes, the phase space of this system is a
230 < $6f$ dimensional space. A point, $x = (q_1 , \ldots ,q_f ,p_1 ,
231 < \ldots ,p_f )$, with a unique set of values of $6f$ coordinates and
232 < momenta is a phase space vector.
227 > momentum variables. Consider a dynamic system of $f$ particles in a
228 > cartesian space, where each of the $6f$ coordinates and momenta is
229 > assigned to one of $6f$ mutually orthogonal axes, the phase space of
230 > this system is a $6f$ dimensional space. A point, $x = (q_1 , \ldots
231 > ,q_f ,p_1 , \ldots ,p_f )$, with a unique set of values of $6f$
232 > coordinates and momenta is a phase space vector.
233  
234   A microscopic state or microstate of a classical system is
235   specification of the complete phase space vector of a system at any
# Line 257 | Line 251 | space. The density of distribution for an ensemble wit
251   regions of the phase space. The condition of an ensemble at any time
252   can be regarded as appropriately specified by the density $\rho$
253   with which representative points are distributed over the phase
254 < space. The density of distribution for an ensemble with $f$ degrees
255 < of freedom is defined as,
254 > space. The density distribution for an ensemble with $f$ degrees of
255 > freedom is defined as,
256   \begin{equation}
257   \rho  = \rho (q_1 , \ldots ,q_f ,p_1 , \ldots ,p_f ,t).
258   \label{introEquation:densityDistribution}
259   \end{equation}
260   Governed by the principles of mechanics, the phase points change
261 < their value which would change the density at any time at phase
262 < space. Hence, the density of distribution is also to be taken as a
261 > their locations which would change the density at any time at phase
262 > space. Hence, the density distribution is also to be taken as a
263   function of the time.
264  
265   The number of systems $\delta N$ at time $t$ can be determined by,
# Line 273 | Line 267 | Assuming a large enough population of systems are expl
267   \delta N = \rho (q,p,t)dq_1  \ldots dq_f dp_1  \ldots dp_f.
268   \label{introEquation:deltaN}
269   \end{equation}
270 < Assuming a large enough population of systems are exploited, we can
271 < sufficiently approximate $\delta N$ without introducing
272 < discontinuity when we go from one region in the phase space to
273 < another. By integrating over the whole phase space,
270 > Assuming a large enough population of systems, we can sufficiently
271 > approximate $\delta N$ without introducing discontinuity when we go
272 > from one region in the phase space to another. By integrating over
273 > the whole phase space,
274   \begin{equation}
275   N = \int { \ldots \int {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f
276   \label{introEquation:totalNumberSystem}
# Line 293 | Line 287 | properties of the ensemble of possibilities as a whole
287   value of any desired quantity which depends on the coordinates and
288   momenta of the system. Even when the dynamics of the real system is
289   complex, or stochastic, or even discontinuous, the average
290 < properties of the ensemble of possibilities as a whole may still
291 < remain well defined. For a classical system in thermal equilibrium
292 < with its environment, the ensemble average of a mechanical quantity,
293 < $\langle A(q , p) \rangle_t$, takes the form of an integral over the
294 < phase space of the system,
290 > properties of the ensemble of possibilities as a whole remaining
291 > well defined. For a classical system in thermal equilibrium with its
292 > environment, the ensemble average of a mechanical quantity, $\langle
293 > A(q , p) \rangle_t$, takes the form of an integral over the phase
294 > space of the system,
295   \begin{equation}
296   \langle  A(q , p) \rangle_t = \frac{{\int { \ldots \int {A(q,p)\rho
297   (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}{{\int { \ldots \int {\rho
# Line 307 | Line 301 | parameters, such as temperature \textit{etc}, partitio
301  
302   There are several different types of ensembles with different
303   statistical characteristics. As a function of macroscopic
304 < parameters, such as temperature \textit{etc}, partition function can
305 < be used to describe the statistical properties of a system in
304 > parameters, such as temperature \textit{etc}, the partition function
305 > can be used to describe the statistical properties of a system in
306   thermodynamic equilibrium.
307  
308   As an ensemble of systems, each of which is known to be thermally
309 < isolated and conserve energy, Microcanonical ensemble(NVE) has a
309 > isolated and conserve energy, the Microcanonical ensemble(NVE) has a
310   partition function like,
311   \begin{equation}
312   \Omega (N,V,E) = e^{\beta TS} \label{introEquation:NVEPartition}.
# Line 326 | Line 320 | TS$. Since most experiment are carried out under const
320   \label{introEquation:NVTPartition}
321   \end{equation}
322   Here, $A$ is the Helmholtz free energy which is defined as $ A = U -
323 < TS$. Since most experiment are carried out under constant pressure
324 < condition, isothermal-isobaric ensemble(NPT) play a very important
325 < role in molecular simulation. The isothermal-isobaric ensemble allow
326 < the system to exchange energy with a heat bath of temperature $T$
327 < and to change the volume as well. Its partition function is given as
323 > TS$. Since most experiments are carried out under constant pressure
324 > condition, the isothermal-isobaric ensemble(NPT) plays a very
325 > important role in molecular simulations. The isothermal-isobaric
326 > ensemble allow the system to exchange energy with a heat bath of
327 > temperature $T$ and to change the volume as well. Its partition
328 > function is given as
329   \begin{equation}
330   \Delta (N,P,T) =  - e^{\beta G}.
331   \label{introEquation:NPTPartition}
# Line 339 | Line 334 | The Liouville's theorem is the foundation on which sta
334  
335   \subsection{\label{introSection:liouville}Liouville's theorem}
336  
337 < The Liouville's theorem is the foundation on which statistical
338 < mechanics rests. It describes the time evolution of phase space
337 > Liouville's theorem is the foundation on which statistical mechanics
338 > rests. It describes the time evolution of the phase space
339   distribution function. In order to calculate the rate of change of
340   $\rho$, we begin from Equation(\ref{introEquation:deltaN}). If we
341   consider the two faces perpendicular to the $q_1$ axis, which are
# Line 369 | Line 364 | divining $ \delta q_1  \ldots \delta q_f \delta p_1  \
364   + \frac{{\partial \dot p_i }}{{\partial p_i }}} \right)}  = 0 ,
365   \end{equation}
366   which cancels the first terms of the right hand side. Furthermore,
367 < divining $ \delta q_1  \ldots \delta q_f \delta p_1  \ldots \delta
367 > dividing $ \delta q_1  \ldots \delta q_f \delta p_1  \ldots \delta
368   p_f $ in both sides, we can write out Liouville's theorem in a
369   simple form,
370   \begin{equation}
# Line 395 | Line 390 | distribution,
390   \label{introEquation:densityAndHamiltonian}
391   \end{equation}
392  
393 < \subsubsection{\label{introSection:phaseSpaceConservation}Conservation of Phase Space}
393 > \subsubsection{\label{introSection:phaseSpaceConservation}\textbf{Conservation of Phase Space}}
394   Lets consider a region in the phase space,
395   \begin{equation}
396   \delta v = \int { \ldots \int {dq_1 } ...dq_f dp_1 } ..dp_f .
397   \end{equation}
398   If this region is small enough, the density $\rho$ can be regarded
399 < as uniform over the whole phase space. Thus, the number of phase
400 < points inside this region is given by,
399 > as uniform over the whole integral. Thus, the number of phase points
400 > inside this region is given by,
401   \begin{equation}
402   \delta N = \rho \delta v = \rho \int { \ldots \int {dq_1 } ...dq_f
403   dp_1 } ..dp_f.
# Line 414 | Line 409 | With the help of stationary assumption
409   \end{equation}
410   With the help of stationary assumption
411   (\ref{introEquation:stationary}), we obtain the principle of the
412 < \emph{conservation of extension in phase space},
412 > \emph{conservation of volume in phase space},
413   \begin{equation}
414   \frac{d}{{dt}}(\delta v) = \frac{d}{{dt}}\int { \ldots \int {dq_1 }
415   ...dq_f dp_1 } ..dp_f  = 0.
416   \label{introEquation:volumePreserving}
417   \end{equation}
418  
419 < \subsubsection{\label{introSection:liouvilleInOtherForms}Liouville's Theorem in Other Forms}
419 > \subsubsection{\label{introSection:liouvilleInOtherForms}\textbf{Liouville's Theorem in Other Forms}}
420  
421   Liouville's theorem can be expresses in a variety of different forms
422   which are convenient within different contexts. For any two function
# Line 463 | Line 458 | simulation and the quality of the underlying model. Ho
458   Various thermodynamic properties can be calculated from Molecular
459   Dynamics simulation. By comparing experimental values with the
460   calculated properties, one can determine the accuracy of the
461 < simulation and the quality of the underlying model. However, both of
462 < experiment and computer simulation are usually performed during a
461 > simulation and the quality of the underlying model. However, both
462 > experiments and computer simulations are usually performed during a
463   certain time interval and the measurements are averaged over a
464   period of them which is different from the average behavior of
465 < many-body system in Statistical Mechanics. Fortunately, Ergodic
466 < Hypothesis is proposed to make a connection between time average and
467 < ensemble average. It states that time average and average over the
468 < statistical ensemble are identical \cite{Frenkel1996, leach01:mm}.
465 > many-body system in Statistical Mechanics. Fortunately, the Ergodic
466 > Hypothesis makes a connection between time average and the ensemble
467 > average. It states that the time average and average over the
468 > statistical ensemble are identical \cite{Frenkel1996, Leach2001}.
469   \begin{equation}
470   \langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty }
471   \frac{1}{t}\int\limits_0^t {A(q(t),p(t))dt = \int\limits_\Gamma
# Line 484 | Line 479 | reasonable, the Monte Carlo techniques\cite{metropolis
479   a properly weighted statistical average. This allows the researcher
480   freedom of choice when deciding how best to measure a given
481   observable. In case an ensemble averaged approach sounds most
482 < reasonable, the Monte Carlo techniques\cite{metropolis:1949} can be
482 > reasonable, the Monte Carlo techniques\cite{Metropolis1949} can be
483   utilized. Or if the system lends itself to a time averaging
484   approach, the Molecular Dynamics techniques in
485   Sec.~\ref{introSection:molecularDynamics} will be the best
486   choice\cite{Frenkel1996}.
487  
488   \section{\label{introSection:geometricIntegratos}Geometric Integrators}
489 < A variety of numerical integrators were proposed to simulate the
490 < motions. They usually begin with an initial conditionals and move
491 < the objects in the direction governed by the differential equations.
492 < However, most of them ignore the hidden physical law contained
493 < within the equations. Since 1990, geometric integrators, which
494 < preserve various phase-flow invariants such as symplectic structure,
495 < volume and time reversal symmetry, are developed to address this
496 < issue. The velocity verlet method, which happens to be a simple
497 < example of symplectic integrator, continues to gain its popularity
498 < in molecular dynamics community. This fact can be partly explained
499 < by its geometric nature.
489 > A variety of numerical integrators have been proposed to simulate
490 > the motions of atoms in MD simulation. They usually begin with
491 > initial conditionals and move the objects in the direction governed
492 > by the differential equations. However, most of them ignore the
493 > hidden physical laws contained within the equations. Since 1990,
494 > geometric integrators, which preserve various phase-flow invariants
495 > such as symplectic structure, volume and time reversal symmetry, are
496 > developed to address this issue\cite{Dullweber1997, McLachlan1998,
497 > Leimkuhler1999}. The velocity verlet method, which happens to be a
498 > simple example of symplectic integrator, continues to gain
499 > popularity in the molecular dynamics community. This fact can be
500 > partly explained by its geometric nature.
501  
502 < \subsection{\label{introSection:symplecticManifold}Symplectic Manifold}
503 < A \emph{manifold} is an abstract mathematical space. It locally
504 < looks like Euclidean space, but when viewed globally, it may have
505 < more complicate structure. A good example of manifold is the surface
506 < of Earth. It seems to be flat locally, but it is round if viewed as
507 < a whole. A \emph{differentiable manifold} (also known as
508 < \emph{smooth manifold}) is a manifold with an open cover in which
509 < the covering neighborhoods are all smoothly isomorphic to one
510 < another. In other words,it is possible to apply calculus on
515 < \emph{differentiable manifold}. A \emph{symplectic manifold} is
516 < defined as a pair $(M, \omega)$ which consisting of a
502 > \subsection{\label{introSection:symplecticManifold}Symplectic Manifolds}
503 > A \emph{manifold} is an abstract mathematical space. It looks
504 > locally like Euclidean space, but when viewed globally, it may have
505 > more complicated structure. A good example of manifold is the
506 > surface of Earth. It seems to be flat locally, but it is round if
507 > viewed as a whole. A \emph{differentiable manifold} (also known as
508 > \emph{smooth manifold}) is a manifold on which it is possible to
509 > apply calculus on \emph{differentiable manifold}. A \emph{symplectic
510 > manifold} is defined as a pair $(M, \omega)$ which consists of a
511   \emph{differentiable manifold} $M$ and a close, non-degenerated,
512   bilinear symplectic form, $\omega$. A symplectic form on a vector
513   space $V$ is a function $\omega(x, y)$ which satisfies
514   $\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+
515   \lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and
516 < $\omega(x, x) = 0$. Cross product operation in vector field is an
517 < example of symplectic form.
516 > $\omega(x, x) = 0$. The cross product operation in vector field is
517 > an example of symplectic form.
518  
519 < One of the motivations to study \emph{symplectic manifold} in
519 > One of the motivations to study \emph{symplectic manifolds} in
520   Hamiltonian Mechanics is that a symplectic manifold can represent
521   all possible configurations of the system and the phase space of the
522   system can be described by it's cotangent bundle. Every symplectic
523   manifold is even dimensional. For instance, in Hamilton equations,
524   coordinate and momentum always appear in pairs.
525  
532 – Let  $(M,\omega)$ and $(N, \eta)$ be symplectic manifolds. A map
533 – \[
534 – f : M \rightarrow N
535 – \]
536 – is a \emph{symplectomorphism} if it is a \emph{diffeomorphims} and
537 – the \emph{pullback} of $\eta$ under f is equal to $\omega$.
538 – Canonical transformation is an example of symplectomorphism in
539 – classical mechanics.
540 –
526   \subsection{\label{introSection:ODE}Ordinary Differential Equations}
527  
528 < For a ordinary differential system defined as
528 > For an ordinary differential system defined as
529   \begin{equation}
530   \dot x = f(x)
531   \end{equation}
532 < where $x = x(q,p)^T$, this system is canonical Hamiltonian, if
532 > where $x = x(q,p)^T$, this system is a canonical Hamiltonian, if
533   \begin{equation}
534   f(r) = J\nabla _x H(r).
535   \end{equation}
# Line 565 | Line 550 | Another generalization of Hamiltonian dynamics is Pois
550   \end{equation}In this case, $f$ is
551   called a \emph{Hamiltonian vector field}.
552  
553 < Another generalization of Hamiltonian dynamics is Poisson Dynamics,
553 > Another generalization of Hamiltonian dynamics is Poisson
554 > Dynamics\cite{Olver1986},
555   \begin{equation}
556   \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
557   \end{equation}
# Line 612 | Line 598 | The hidden geometric properties of ODE and its flow pl
598  
599   \subsection{\label{introSection:geometricProperties}Geometric Properties}
600  
601 < The hidden geometric properties of ODE and its flow play important
602 < roles in numerical studies. Many of them can be found in systems
603 < which occur naturally in applications.
601 > The hidden geometric properties\cite{Budd1999, Marsden1998} of ODE
602 > and its flow play important roles in numerical studies. Many of them
603 > can be found in systems which occur naturally in applications.
604  
605   Let $\varphi$ be the flow of Hamiltonian vector field, $\varphi$ is
606   a \emph{symplectic} flow if it satisfies,
# Line 658 | Line 644 | smooth function $G$ is given by,
644   which is the condition for conserving \emph{first integral}. For a
645   canonical Hamiltonian system, the time evolution of an arbitrary
646   smooth function $G$ is given by,
647 < \begin{equation}
648 < \begin{array}{c}
649 < \frac{{dG(x(t))}}{{dt}} = [\nabla _x G(x(t))]^T \dot x(t) \\
650 <  = [\nabla _x G(x(t))]^T J\nabla _x H(x(t)). \\
665 < \end{array}
647 >
648 > \begin{eqnarray}
649 > \frac{{dG(x(t))}}{{dt}} & = & [\nabla _x G(x(t))]^T \dot x(t) \\
650 >                        & = & [\nabla _x G(x(t))]^T J\nabla _x H(x(t)). \\
651   \label{introEquation:firstIntegral1}
652 < \end{equation}
652 > \end{eqnarray}
653 >
654 >
655   Using poisson bracket notion, Equation
656   \ref{introEquation:firstIntegral1} can be rewritten as
657   \[
# Line 679 | Line 666 | is a \emph{first integral}, which is due to the fact $
666   is a \emph{first integral}, which is due to the fact $\{ H,H\}  =
667   0$.
668  
669 <
683 < When designing any numerical methods, one should always try to
669 > When designing any numerical methods, one should always try to
670   preserve the structural properties of the original ODE and its flow.
671  
672   \subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
673   A lot of well established and very effective numerical methods have
674   been successful precisely because of their symplecticities even
675   though this fact was not recognized when they were first
676 < constructed. The most famous example is leapfrog methods in
677 < molecular dynamics. In general, symplectic integrators can be
676 > constructed. The most famous example is the Verlet-leapfrog methods
677 > in molecular dynamics. In general, symplectic integrators can be
678   constructed using one of four different methods.
679   \begin{enumerate}
680   \item Generating functions
# Line 697 | Line 683 | Generating function tends to lead to methods which are
683   \item Splitting methods
684   \end{enumerate}
685  
686 < Generating function tends to lead to methods which are cumbersome
687 < and difficult to use. In dissipative systems, variational methods
688 < can capture the decay of energy accurately. Since their
689 < geometrically unstable nature against non-Hamiltonian perturbations,
690 < ordinary implicit Runge-Kutta methods are not suitable for
691 < Hamiltonian system. Recently, various high-order explicit
692 < Runge--Kutta methods have been developed to overcome this
686 > Generating function\cite{Channell1990} tends to lead to methods
687 > which are cumbersome and difficult to use. In dissipative systems,
688 > variational methods can capture the decay of energy
689 > accurately\cite{Kane2000}. Since their geometrically unstable nature
690 > against non-Hamiltonian perturbations, ordinary implicit Runge-Kutta
691 > methods are not suitable for Hamiltonian system. Recently, various
692 > high-order explicit Runge-Kutta methods
693 > \cite{Owren1992,Chen2003}have been developed to overcome this
694   instability. However, due to computational penalty involved in
695 < implementing the Runge-Kutta methods, they do not attract too much
696 < attention from Molecular Dynamics community. Instead, splitting have
697 < been widely accepted since they exploit natural decompositions of
698 < the system\cite{Tuckerman92}.
695 > implementing the Runge-Kutta methods, they have not attracted much
696 > attention from the Molecular Dynamics community. Instead, splitting
697 > methods have been widely accepted since they exploit natural
698 > decompositions of the system\cite{Tuckerman1992, McLachlan1998}.
699  
700 < \subsubsection{\label{introSection:splittingMethod}Splitting Method}
700 > \subsubsection{\label{introSection:splittingMethod}\textbf{Splitting Methods}}
701  
702   The main idea behind splitting methods is to decompose the discrete
703   $\varphi_h$ as a composition of simpler flows,
# Line 731 | Line 718 | order is then given by the Lie-Trotter formula
718   energy respectively, which is a natural decomposition of the
719   problem. If $H_1$ and $H_2$ can be integrated using exact flows
720   $\varphi_1(t)$ and $\varphi_2(t)$, respectively, a simple first
721 < order is then given by the Lie-Trotter formula
721 > order expression is then given by the Lie-Trotter formula
722   \begin{equation}
723   \varphi _h  = \varphi _{1,h}  \circ \varphi _{2,h},
724   \label{introEquation:firstOrderSplitting}
# Line 757 | Line 744 | which has a local error proportional to $h^3$. Sprang
744   \varphi _h  = \varphi _{1,h/2}  \circ \varphi _{2,h}  \circ \varphi
745   _{1,h/2} , \label{introEquation:secondOrderSplitting}
746   \end{equation}
747 < which has a local error proportional to $h^3$. Sprang splitting's
748 < popularity in molecular simulation community attribute to its
749 < symmetric property,
747 > which has a local error proportional to $h^3$. The Sprang
748 > splitting's popularity in molecular simulation community attribute
749 > to its symmetric property,
750   \begin{equation}
751   \varphi _h^{ - 1} = \varphi _{ - h}.
752   \label{introEquation:timeReversible}
753   \end{equation}
754  
755 < \subsubsection{\label{introSection:exampleSplittingMethod}Example of Splitting Method}
755 > \subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Example of Splitting Method}}
756   The classical equation for a system consisting of interacting
757   particles can be written in Hamiltonian form,
758   \[
# Line 825 | Line 812 | q(\Delta t)} \right]. %
812   \label{introEquation:positionVerlet2}
813   \end{align}
814  
815 < \subsubsection{\label{introSection:errorAnalysis}Error Analysis and Higher Order Methods}
815 > \subsubsection{\label{introSection:errorAnalysis}\textbf{Error Analysis and Higher Order Methods}}
816  
817   Baker-Campbell-Hausdorff formula can be used to determine the local
818   error of splitting method in terms of commutator of the
# Line 844 | Line 831 | Applying Baker-Campbell-Hausdorff formula to Sprang sp
831   \[
832   [X,Y] = XY - YX .
833   \]
834 < Applying Baker-Campbell-Hausdorff formula to Sprang splitting, we
835 < can obtain
836 < \begin{equation}
837 < \exp (h X/2)\exp (h Y)\exp (h X/2) & = & \exp (h X + h Y + h^2
838 < [X,Y]/4 + h^2 [Y,X]/4 \\ & & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\
839 < & & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3 [X,[X,Y]]/24 & & \mbox{} +
840 < \ldots )
854 < \end{equation}
834 > Applying Baker-Campbell-Hausdorff formula\cite{Varadarajan1974} to
835 > Sprang splitting, we can obtain
836 > \begin{eqnarray*}
837 > \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 \\
838 >                                   &   & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\
839 >                                   &   & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3[X,[X,Y]]/24 + \ldots )
840 > \end{eqnarray*}
841   Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0\], the dominant local
842   error of Spring splitting is proportional to $h^3$. The same
843   procedure can be applied to general splitting,  of the form
# Line 859 | Line 845 | Careful choice of coefficient $a_1 ,\ldot , b_m$ will
845   \varphi _{b_m h}^2  \circ \varphi _{a_m h}^1  \circ \varphi _{b_{m -
846   1} h}^2  \circ  \ldots  \circ \varphi _{a_1 h}^1 .
847   \end{equation}
848 < Careful choice of coefficient $a_1 ,\ldot , b_m$ will lead to higher
848 > Careful choice of coefficient $a_1 \ldots b_m$ will lead to higher
849   order method. Yoshida proposed an elegant way to compose higher
850 < order methods based on symmetric splitting. Given a symmetric second
851 < order base method $ \varphi _h^{(2)} $, a fourth-order symmetric
852 < method can be constructed by composing,
850 > order methods based on symmetric splitting\cite{Yoshida1990}. Given
851 > a symmetric second order base method $ \varphi _h^{(2)} $, a
852 > fourth-order symmetric method can be constructed by composing,
853   \[
854   \varphi _h^{(4)}  = \varphi _{\alpha h}^{(2)}  \circ \varphi _{\beta
855   h}^{(2)}  \circ \varphi _{\alpha h}^{(2)}
# Line 912 | Line 898 | initialization of a simulation. Sec.~\ref{introSec:pro
898   \end{enumerate}
899   These three individual steps will be covered in the following
900   sections. Sec.~\ref{introSec:initialSystemSettings} deals with the
901 < initialization of a simulation. Sec.~\ref{introSec:production} will
902 < discusses issues in production run. Sec.~\ref{introSection:Analysis}
903 < provides the theoretical tools for trajectory analysis.
901 > initialization of a simulation. Sec.~\ref{introSection:production}
902 > will discusses issues in production run.
903 > Sec.~\ref{introSection:Analysis} provides the theoretical tools for
904 > trajectory analysis.
905  
906   \subsection{\label{introSec:initialSystemSettings}Initialization}
907  
908 < \subsubsection{Preliminary preparation}
908 > \subsubsection{\textbf{Preliminary preparation}}
909  
910   When selecting the starting structure of a molecule for molecular
911   simulation, one may retrieve its Cartesian coordinates from public
# Line 936 | Line 923 | interested in self-aggregation and it takes a long tim
923   instead of placing lipids randomly in solvent, since we are not
924   interested in self-aggregation and it takes a long time to happen.
925  
926 < \subsubsection{Minimization}
926 > \subsubsection{\textbf{Minimization}}
927  
928   It is quite possible that some of molecules in the system from
929   preliminary preparation may be overlapped with each other. This
# Line 958 | Line 945 | Newton-Raphson methods can not be used with very large
945   matrix and insufficient storage capacity to store them, most
946   Newton-Raphson methods can not be used with very large models.
947  
948 < \subsubsection{Heating}
948 > \subsubsection{\textbf{Heating}}
949  
950   Typically, Heating is performed by assigning random velocities
951   according to a Gaussian distribution for a temperature. Beginning at
# Line 970 | Line 957 | shifted to zero.
957   net linear momentum and angular momentum of the system should be
958   shifted to zero.
959  
960 < \subsubsection{Equilibration}
960 > \subsubsection{\textbf{Equilibration}}
961  
962   The purpose of equilibration is to allow the system to evolve
963   spontaneously for a period of time and reach equilibrium. The
# Line 984 | Line 971 | Production run is the most important steps of the simu
971  
972   \subsection{\label{introSection:production}Production}
973  
974 < Production run is the most important steps of the simulation, in
974 > Production run is the most important step of the simulation, in
975   which the equilibrated structure is used as a starting point and the
976   motions of the molecules are collected for later analysis. In order
977   to capture the macroscopic properties of the system, the molecular
# Line 1000 | Line 987 | approach will suffer from the surface effect. To offse
987   A natural approach to avoid system size issue is to represent the
988   bulk behavior by a finite number of the particles. However, this
989   approach will suffer from the surface effect. To offset this,
990 < \textit{Periodic boundary condition} is developed to simulate bulk
991 < properties with a relatively small number of particles. In this
992 < method, the simulation box is replicated throughout space to form an
993 < infinite lattice. During the simulation, when a particle moves in
994 < the primary cell, its image in other cells move in exactly the same
995 < direction with exactly the same orientation. Thus, as a particle
996 < leaves the primary cell, one of its images will enter through the
997 < opposite face.
998 < %\begin{figure}
999 < %\centering
1000 < %\includegraphics[width=\linewidth]{pbcFig.eps}
1001 < %\caption[An illustration of periodic boundary conditions]{A 2-D
1002 < %illustration of periodic boundary conditions. As one particle leaves
1003 < %the right of the simulation box, an image of it enters the left.}
1004 < %\label{introFig:pbc}
1005 < %\end{figure}
990 > \textit{Periodic boundary condition} (see Fig.~\ref{introFig:pbc})
991 > is developed to simulate bulk properties with a relatively small
992 > number of particles. In this method, the simulation box is
993 > replicated throughout space to form an infinite lattice. During the
994 > simulation, when a particle moves in the primary cell, its image in
995 > other cells move in exactly the same direction with exactly the same
996 > orientation. Thus, as a particle leaves the primary cell, one of its
997 > images will enter through the opposite face.
998 > \begin{figure}
999 > \centering
1000 > \includegraphics[width=\linewidth]{pbc.eps}
1001 > \caption[An illustration of periodic boundary conditions]{A 2-D
1002 > illustration of periodic boundary conditions. As one particle leaves
1003 > the left of the simulation box, an image of it enters the right.}
1004 > \label{introFig:pbc}
1005 > \end{figure}
1006  
1007   %cutoff and minimum image convention
1008   Another important technique to improve the efficiency of force
# Line 1033 | Line 1020 | discrete Fourier transforms, the particle mesh-based m
1020   reciprocal sums with rapid and absolute convergence, has proved to
1021   minimize the periodicity artifacts in liquid simulations. Taking the
1022   advantages of the fast Fourier transform (FFT) for calculating
1023 < discrete Fourier transforms, the particle mesh-based methods are
1024 < accelerated from $O(N^{3/2})$ to $O(N logN)$. An alternative
1025 < approach is \emph{fast multipole method}, which treats Coulombic
1026 < interaction exactly at short range, and approximate the potential at
1027 < long range through multipolar expansion. In spite of their wide
1028 < acceptances at the molecular simulation community, these two methods
1029 < are hard to be implemented correctly and efficiently. Instead, we
1030 < use a damped and charge-neutralized Coulomb potential method
1031 < developed by Wolf and his coworkers. The shifted Coulomb potential
1032 < for particle $i$ and particle $j$ at distance $r_{rj}$ is given by:
1023 > discrete Fourier transforms, the particle mesh-based
1024 > methods\cite{Hockney1981,Shimada1993, Luty1994} are accelerated from
1025 > $O(N^{3/2})$ to $O(N logN)$. An alternative approach is \emph{fast
1026 > multipole method}\cite{Greengard1987, Greengard1994}, which treats
1027 > Coulombic interaction exactly at short range, and approximate the
1028 > potential at long range through multipolar expansion. In spite of
1029 > their wide acceptances at the molecular simulation community, these
1030 > two methods are hard to be implemented correctly and efficiently.
1031 > Instead, we use a damped and charge-neutralized Coulomb potential
1032 > method developed by Wolf and his coworkers\cite{Wolf1999}. The
1033 > shifted Coulomb potential for particle $i$ and particle $j$ at
1034 > distance $r_{rj}$ is given by:
1035   \begin{equation}
1036   V(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha
1037   r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow
# Line 1052 | Line 1041 | efficient and easy to implement.
1041   where $\alpha$ is the convergence parameter. Due to the lack of
1042   inherent periodicity and rapid convergence,this method is extremely
1043   efficient and easy to implement.
1044 < %\begin{figure}
1045 < %\centering
1046 < %\includegraphics[width=\linewidth]{pbcFig.eps}
1047 < %\caption[An illustration of shifted Coulomb potential]{An illustration of shifted Coulomb potential.}
1048 < %\label{introFigure:shiftedCoulomb}
1049 < %\end{figure}
1044 > \begin{figure}
1045 > \centering
1046 > \includegraphics[width=\linewidth]{shifted_coulomb.eps}
1047 > \caption[An illustration of shifted Coulomb potential]{An
1048 > illustration of shifted Coulomb potential.}
1049 > \label{introFigure:shiftedCoulomb}
1050 > \end{figure}
1051  
1052   %multiple time step
1053  
# Line 1073 | Line 1063 | from the trajectories.
1063   parameters, and investigate time-dependent processes of the molecule
1064   from the trajectories.
1065  
1066 < \subsubsection{\label{introSection:thermodynamicsProperties}Thermodynamics Properties}
1066 > \subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamics Properties}}
1067  
1068   Thermodynamics properties, which can be expressed in terms of some
1069   function of the coordinates and momenta of all particles in the
# Line 1095 | Line 1085 | P = \frac{{Nk_B T}}{V} - \frac{1}{{3V}}\left\langle {\
1085   < j} {r{}_{ij} \cdot f_{ij} } } \right\rangle
1086   \end{equation}
1087  
1088 < \subsubsection{\label{introSection:structuralProperties}Structural Properties}
1088 > \subsubsection{\label{introSection:structuralProperties}\textbf{Structural Properties}}
1089  
1090   Structural Properties of a simple fluid can be described by a set of
1091   distribution functions. Among these functions,\emph{pair
# Line 1103 | Line 1093 | integrating over the surface factor \cite{Powles73}. T
1093   function}, is of most fundamental importance to liquid-state theory.
1094   Pair distribution function can be gathered by Fourier transforming
1095   raw data from a series of neutron diffraction experiments and
1096 < integrating over the surface factor \cite{Powles73}. The experiment
1097 < result can serve as a criterion to justify the correctness of the
1098 < theory. Moreover, various equilibrium thermodynamic and structural
1099 < properties can also be expressed in terms of radial distribution
1100 < function \cite{allen87:csl}.
1096 > integrating over the surface factor \cite{Powles1973}. The
1097 > experiment result can serve as a criterion to justify the
1098 > correctness of the theory. Moreover, various equilibrium
1099 > thermodynamic and structural properties can also be expressed in
1100 > terms of radial distribution function \cite{Allen1987}.
1101  
1102   A pair distribution functions $g(r)$ gives the probability that a
1103   particle $i$ will be located at a distance $r$ from a another
# Line 1135 | Line 1125 | other is essentially zero.
1125   %\label{introFigure:pairDistributionFunction}
1126   %\end{figure}
1127  
1128 < \subsubsection{\label{introSection:timeDependentProperties}Time-dependent
1129 < Properties}
1128 > \subsubsection{\label{introSection:timeDependentProperties}\textbf{Time-dependent
1129 > Properties}}
1130  
1131   Time-dependent properties are usually calculated using \emph{time
1132   correlation function}, which correlates random variables $A$ and $B$
# Line 1185 | Line 1175 | protein-protein docking study{\cite{Gray03}}.
1175   movement of the objects in 3D gaming engine or other physics
1176   simulator is governed by the rigid body dynamics. In molecular
1177   simulation, rigid body is used to simplify the model in
1178 < protein-protein docking study{\cite{Gray03}}.
1178 > protein-protein docking study\cite{Gray2003}.
1179  
1180   It is very important to develop stable and efficient methods to
1181   integrate the equations of motion of orientational degrees of
# Line 1193 | Line 1183 | different sets of Euler angles can overcome this diffi
1183   rotational degrees of freedom. However, due to its singularity, the
1184   numerical integration of corresponding equations of motion is very
1185   inefficient and inaccurate. Although an alternative integrator using
1186 < different sets of Euler angles can overcome this difficulty\cite{},
1187 < the computational penalty and the lost of angular momentum
1188 < conservation still remain. A singularity free representation
1189 < utilizing quaternions was developed by Evans in 1977. Unfortunately,
1190 < this approach suffer from the nonseparable Hamiltonian resulted from
1191 < quaternion representation, which prevents the symplectic algorithm
1192 < to be utilized. Another different approach is to apply holonomic
1193 < constraints to the atoms belonging to the rigid body. Each atom
1194 < moves independently under the normal forces deriving from potential
1195 < energy and constraint forces which are used to guarantee the
1196 < rigidness. However, due to their iterative nature, SHAKE and Rattle
1197 < algorithm converge very slowly when the number of constraint
1198 < increases.
1186 > different sets of Euler angles can overcome this
1187 > difficulty\cite{Barojas1973}, the computational penalty and the lost
1188 > of angular momentum conservation still remain. A singularity free
1189 > representation utilizing quaternions was developed by Evans in
1190 > 1977\cite{Evans1977}. Unfortunately, this approach suffer from the
1191 > nonseparable Hamiltonian resulted from quaternion representation,
1192 > which prevents the symplectic algorithm to be utilized. Another
1193 > different approach is to apply holonomic constraints to the atoms
1194 > belonging to the rigid body. Each atom moves independently under the
1195 > normal forces deriving from potential energy and constraint forces
1196 > which are used to guarantee the rigidness. However, due to their
1197 > iterative nature, SHAKE and Rattle algorithm converge very slowly
1198 > when the number of constraint increases\cite{Ryckaert1977,
1199 > Andersen1983}.
1200  
1201   The break through in geometric literature suggests that, in order to
1202   develop a long-term integration scheme, one should preserve the
1203   symplectic structure of the flow. Introducing conjugate momentum to
1204   rotation matrix $Q$ and re-formulating Hamiltonian's equation, a
1205 < symplectic integrator, RSHAKE, was proposed to evolve the
1206 < Hamiltonian system in a constraint manifold by iteratively
1205 > symplectic integrator, RSHAKE\cite{Kol1997}, was proposed to evolve
1206 > the Hamiltonian system in a constraint manifold by iteratively
1207   satisfying the orthogonality constraint $Q_T Q = 1$. An alternative
1208 < method using quaternion representation was developed by Omelyan.
1209 < However, both of these methods are iterative and inefficient. In
1210 < this section, we will present a symplectic Lie-Poisson integrator
1211 < for rigid body developed by Dullweber and his
1212 < coworkers\cite{Dullweber1997} in depth.
1208 > method using quaternion representation was developed by
1209 > Omelyan\cite{Omelyan1998}. However, both of these methods are
1210 > iterative and inefficient. In this section, we will present a
1211 > symplectic Lie-Poisson integrator for rigid body developed by
1212 > Dullweber and his coworkers\cite{Dullweber1997} in depth.
1213  
1214   \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
1215   The motion of the rigid body is Hamiltonian with the Hamiltonian
# Line 1250 | Line 1241 | the equations of motion,
1241   Using Equation (\ref{introEquation:motionHamiltonianCoordinate},
1242   \ref{introEquation:motionHamiltonianMomentum}), one can write down
1243   the equations of motion,
1253 – \[
1254 – \begin{array}{c}
1255 – \frac{{dq}}{{dt}} = \frac{p}{m} \label{introEquation:RBMotionPosition}\\
1256 – \frac{{dp}}{{dt}} =  - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
1257 – \frac{{dQ}}{{dt}} = PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
1258 – \frac{{dP}}{{dt}} =  - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}\\
1259 – \end{array}
1260 – \]
1244  
1245 + \begin{eqnarray}
1246 + \frac{{dq}}{{dt}} & = & \frac{p}{m} \label{introEquation:RBMotionPosition}\\
1247 + \frac{{dp}}{{dt}} & = & - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
1248 + \frac{{dQ}}{{dt}} & = & PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
1249 + \frac{{dP}}{{dt}} & = & - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}
1250 + \end{eqnarray}
1251 +
1252   In general, there are two ways to satisfy the holonomic constraints.
1253   We can use constraint force provided by lagrange multiplier on the
1254   normal manifold to keep the motion on constraint space. Or we can
# Line 1338 | Line 1328 | operations
1328   \[
1329   \hat vu = v \times u
1330   \]
1341 –
1331   Using \ref{introEqaution:RBMotionPI}, one can construct a skew
1332   matrix,
1333   \begin{equation}
1334 < (\mathop \Pi \limits^ \bullet   - \mathop \Pi \limits^ \bullet  ^T
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}
# Line 1351 | Line 1340 | not be avoided in other methods\cite{}.
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{}.
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
# Line 1371 | Line 1360 | first term of \ref{ introEquation:bodyAngularMotion}).
1360  
1361   If there is not external forces exerted on the rigid body, the only
1362   contribution to the rotational is from the kinetic potential (the
1363 < first term of \ref{ introEquation:bodyAngularMotion}). The free
1364 < rigid body is an example of Lie-Poisson system with Hamiltonian
1376 < function
1363 > first term of \ref{introEquation:bodyAngularMotion}). The free rigid
1364 > body is an example of Lie-Poisson system with Hamiltonian function
1365   \begin{equation}
1366   T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
1367   \label{introEquation:rotationalKineticRB}
# Line 1618 | Line 1606 | Operator. Below are some important properties of Lapla
1606   \]
1607   where  $p$ is real and  $L$ is called the Laplace Transform
1608   Operator. Below are some important properties of Laplace transform
1609 < \begin{equation}
1610 < \begin{array}{c}
1611 < L(x + y) = L(x) + L(y) \\
1612 < L(ax) = aL(x) \\
1613 < L(\dot x) = pL(x) - px(0) \\
1614 < L(\ddot x) = p^2 L(x) - px(0) - \dot x(0) \\
1615 < L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p) \\
1616 < \end{array}
1617 < \end{equation}
1618 <
1609 >
1610 > \begin{eqnarray*}
1611 > L(x + y)  & = & L(x) + L(y) \\
1612 > L(ax)     & = & aL(x) \\
1613 > L(\dot x) & = & pL(x) - px(0) \\
1614 > L(\ddot x)& = & p^2 L(x) - px(0) - \dot x(0) \\
1615 > L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right)& = & G(p)H(p) \\
1616 > \end{eqnarray*}
1617 >
1618 >
1619   Applying Laplace transform to the bath coordinates, we obtain
1620 < \[
1621 < \begin{array}{c}
1622 < 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) \\
1623 < L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
1624 < \end{array}
1637 < \]
1620 > \begin{eqnarray*}
1621 > 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) \\
1622 > L(x_\alpha  ) & = & \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
1623 > \end{eqnarray*}
1624 >
1625   By the same way, the system coordinates become
1626 < \[
1627 < \begin{array}{c}
1628 < mL(\ddot x) =  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
1629 <  - \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\}}  \\
1643 < \end{array}
1644 < \]
1626 > \begin{eqnarray*}
1627 > mL(\ddot x) & = & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
1628 >  & & \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\}}  \\
1629 > \end{eqnarray*}
1630  
1631   With the help of some relatively important inverse Laplace
1632   transformations:
# Line 1653 | Line 1638 | transformations:
1638   \end{array}
1639   \]
1640   , we obtain
1641 < \begin{align}
1642 < m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
1641 > \begin{eqnarray*}
1642 > m\ddot x & =  & - \frac{{\partial W(x)}}{{\partial x}} -
1643   \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
1644   }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\int_0^t {\cos (\omega
1645 < _\alpha  t)\dot x(t - \tau )d\tau  - \left[ {g_\alpha  x_\alpha  (0)
1646 < - \frac{{g_\alpha  }}{{m_\alpha  \omega _\alpha  }}} \right]\cos
1647 < (\omega _\alpha  t) - \frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega
1648 < _\alpha  }}\sin (\omega _\alpha  t)} } \right\}}
1649 < %
1650 < &= - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
1645 > _\alpha  t)\dot x(t - \tau )d\tau } } \right\}}  \\
1646 > & & + \sum\limits_{\alpha  = 1}^N {\left\{ {\left[ {g_\alpha
1647 > x_\alpha (0) - \frac{{g_\alpha  }}{{m_\alpha  \omega _\alpha  }}}
1648 > \right]\cos (\omega _\alpha  t) + \frac{{g_\alpha  \dot x_\alpha
1649 > (0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)} \right\}}
1650 > \end{eqnarray*}
1651 > \begin{eqnarray*}
1652 > m\ddot x & = & - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
1653   {\sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1654   }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha
1655 < t)\dot x(t - \tau )d} \tau }  + \sum\limits_{\alpha  = 1}^N {\left\{
1656 < {\left[ {g_\alpha  x_\alpha  (0) - \frac{{g_\alpha  }}{{m_\alpha
1657 < \omega _\alpha  }}} \right]\cos (\omega _\alpha  t) +
1658 < \frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega _\alpha  }}\sin
1659 < (\omega _\alpha  t)} \right\}}
1660 < \end{align}
1674 <
1655 > t)\dot x(t - \tau )d} \tau }  \\
1656 > & & + \sum\limits_{\alpha  = 1}^N {\left\{ {\left[ {g_\alpha
1657 > x_\alpha (0) - \frac{{g_\alpha }}{{m_\alpha \omega _\alpha  }}}
1658 > \right]\cos (\omega _\alpha  t) + \frac{{g_\alpha  \dot x_\alpha
1659 > (0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)} \right\}}
1660 > \end{eqnarray*}
1661   Introducing a \emph{dynamic friction kernel}
1662   \begin{equation}
1663   \xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
# Line 1694 | Line 1680 | which is known as the \emph{generalized Langevin equat
1680   \end{equation}
1681   which is known as the \emph{generalized Langevin equation}.
1682  
1683 < \subsubsection{\label{introSection:randomForceDynamicFrictionKernel}Random Force and Dynamic Friction Kernel}
1683 > \subsubsection{\label{introSection:randomForceDynamicFrictionKernel}\textbf{Random Force and Dynamic Friction Kernel}}
1684  
1685   One may notice that $R(t)$ depends only on initial conditions, which
1686   implies it is completely deterministic within the context of a
# Line 1753 | Line 1739 | particles is given in Sec.~\ref{introSection:frictionT
1739   briefly review on calculating friction tensor for arbitrary shaped
1740   particles is given in Sec.~\ref{introSection:frictionTensor}.
1741  
1742 < \subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
1742 > \subsubsection{\label{introSection:secondFluctuationDissipation}\textbf{The Second Fluctuation Dissipation Theorem}}
1743  
1744   Defining a new set of coordinates,
1745   \[
# Line 1765 | Line 1751 | And since the $q$ coordinates are harmonic oscillators
1751   R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
1752   \]
1753   And since the $q$ coordinates are harmonic oscillators,
1754 < \[
1755 < \begin{array}{c}
1756 < \left\langle {q_\alpha ^2 } \right\rangle  = \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
1757 < \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  = \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
1758 < \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle  = \delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
1759 < \left\langle {R(t)R(0)} \right\rangle  = \sum\limits_\alpha  {\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle } }  \\
1760 <  = \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
1761 <  = kT\xi (t) \\
1762 < \end{array}
1763 < \]
1754 >
1755 > \begin{eqnarray*}
1756 > \left\langle {q_\alpha ^2 } \right\rangle  & = & \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
1757 > \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle & = & \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
1758 > \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle & = &\delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
1759 > \left\langle {R(t)R(0)} \right\rangle & = & \sum\limits_\alpha  {\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle } }  \\
1760 >  & = &\sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
1761 >  & = &kT\xi (t) \\
1762 > \end{eqnarray*}
1763 >
1764   Thus, we recover the \emph{second fluctuation dissipation theorem}
1765   \begin{equation}
1766   \xi (t) = \left\langle {R(t)R(0)} \right\rangle
# Line 1819 | Line 1805 | toque.
1805   where $F_r$ is the friction force and $\tau _R$ is the friction
1806   toque.
1807  
1808 < \subsubsection{\label{introSection:resistanceTensorRegular}The Resistance Tensor for Regular Shape}
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,
# Line 1846 | Line 1832 | coordinates by
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
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
# Line 1881 | Line 1867 | and
1867   \end{array}.
1868   \]
1869  
1870 < \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}The Resistance Tensor for Arbitrary Shape}
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
# Line 1890 | Line 1876 | unique\cite{Wegener79} as well as the intrinsic coupli
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{Wegener79} as well as the intrinsic coupling between
1880 < translational and rotational motion of rigid body\cite{}, general
1881 < ellipsoid is not always suitable for modeling arbitrarily shaped
1882 < rigid molecule. A number of studies have been devoted to determine
1883 < the friction tensor for irregularly shaped rigid bodies using more
1884 < advanced method\cite{} where the molecule of interest was modeled by
1885 < combinations of spheres(beads)\cite{} and the hydrodynamics
1886 < properties of the molecule can be calculated using the hydrodynamic
1887 < interaction tensor. Let us consider a rigid assembly of $N$ beads
1888 < immersed in a continuous medium. Due to hydrodynamics interaction,
1889 < the ``net'' velocity of $i$th bead, $v'_i$ is different than its
1890 < unperturbed velocity $v_i$,
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   \]
# Line 1922 | Line 1908 | introduced by Rotne and Prager\cite{} and improved by
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{} and improved by Garc\'{i}a de
1912 < la Torre and Bloomfield,
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
# Line 2009 | Line 1995 | obtain the resistance tensor at $P$ by
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 }  \\
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}
# Line 2024 | Line 2010 | the position of center of resistance,
2010   Using Equations \ref{introEquation:definitionCR} and
2011   \ref{introEquation:resistanceTensorTransformation}, one can locate
2012   the position of center of resistance,
2013 < \[
2014 < \left( \begin{array}{l}
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}
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} \left( \begin{array}{l}
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 < \]
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