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 2881 by tim, Fri Jun 23 20:21:54 2006 UTC vs.
Revision 2912 by tim, Fri Jun 30 02:45:29 2006 UTC

# Line 3 | Line 3 | Mechanics}
3   \section{\label{introSection:classicalMechanics}Classical
4   Mechanics}
5  
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
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
13 < specification of this state; Finally, the specification of the state
14 < when further combine with the laws of mechanics will also be
15 < sufficient to predict the future behavior of the system.
6 > Using equations of motion derived from Classical Mechanics,
7 > Molecular Dynamics simulations are carried out by integrating the
8 > equations of motion for a given system of particles. There are three
9 > fundamental ideas behind classical mechanics. Firstly, one can
10 > determine the state of a mechanical system at any time of interest;
11 > Secondly, all the mechanical properties of the system at that time
12 > can be determined by combining the knowledge of the properties of
13 > the system with the specification of this state; Finally, the
14 > specification of the state when further combined with the laws of
15 > mechanics will also be sufficient to predict the future behavior of
16 > the system.
17  
18   \subsection{\label{introSection:newtonian}Newtonian Mechanics}
19   The discovery of Newton's three laws of mechanics which govern the
# Line 31 | Line 32 | Newton's third law states that
32   $F_{ji}$ be the force that particle $j$ exerts on particle $i$.
33   Newton's third law states that
34   \begin{equation}
35 < F_{ij} = -F_{ji}
35 > F_{ij} = -F_{ji}.
36   \label{introEquation:newtonThirdLaw}
37   \end{equation}
37 –
38   Conservation laws of Newtonian Mechanics play very important roles
39   in solving mechanics problems. The linear momentum of a particle is
40   conserved if it is free or it experiences no force. The second
# Line 63 | Line 63 | momentum of it is conserved. The last conservation the
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
66 < that if all forces are conservative, Energy
67 < \begin{equation}E = T + V \label{introEquation:energyConservation}
66 > that if all forces are conservative, energy is conserved,
67 > \begin{equation}E = T + V. \label{introEquation:energyConservation}
68   \end{equation}
69 < is conserved. All of these conserved quantities are
70 < important factors to determine the quality of numerical integration
71 < schemes for rigid bodies \cite{Dullweber1997}.
69 > All of these conserved quantities are important factors to determine
70 > the quality of numerical integration schemes for rigid bodies
71 > \cite{Dullweber1997}.
72  
73   \subsection{\label{introSection:lagrangian}Lagrangian Mechanics}
74  
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.
75 > Newtonian Mechanics suffers from an important limitation: motion can
76 > only be described in cartesian coordinate systems which make it
77 > impossible to predict analytically the properties of the system even
78 > if we know all of the details of the interaction. In order to
79 > overcome some of the practical difficulties which arise in attempts
80 > to apply Newton's equation to complex systems, approximate numerical
81 > procedures may be developed.
82  
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 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$.
88 > Hamilton's Principle may be stated as follows: the trajectory, along
89 > which a dynamical system may move from one point to another within a
90 > specified time, is derived by finding the path which minimizes the
91 > time integral of the difference between the kinetic $K$, and
92 > potential energies $U$,
93   \begin{equation}
94 < \delta \int_{t_1 }^{t_2 } {(K - U)dt = 0} ,
94 > \delta \int_{t_1 }^{t_2 } {(K - U)dt = 0}.
95   \label{introEquation:halmitonianPrinciple1}
96   \end{equation}
98 –
97   For simple mechanical systems, where the forces acting on the
98   different parts are derivable from a potential, the Lagrangian
99   function $L$ can be defined as the difference between the kinetic
100   energy of the system and its potential energy,
101   \begin{equation}
102 < L \equiv K - U = L(q_i ,\dot q_i ) ,
102 > L \equiv K - U = L(q_i ,\dot q_i ).
103   \label{introEquation:lagrangianDef}
104   \end{equation}
105 < then Eq.~\ref{introEquation:halmitonianPrinciple1} becomes
105 > Thus, Eq.~\ref{introEquation:halmitonianPrinciple1} becomes
106   \begin{equation}
107 < \delta \int_{t_1 }^{t_2 } {L dt = 0} ,
107 > \delta \int_{t_1 }^{t_2 } {L dt = 0} .
108   \label{introEquation:halmitonianPrinciple2}
109   \end{equation}
110  
# Line 138 | Line 136 | p_i  = \frac{{\partial L}}{{\partial q_i }}
136   p_i  = \frac{{\partial L}}{{\partial q_i }}
137   \label{introEquation:generalizedMomentaDot}
138   \end{equation}
141 –
139   With the help of the generalized momenta, we may now define a new
140   quantity $H$ by the equation
141   \begin{equation}
# Line 146 | Line 143 | where $ \dot q_1  \ldots \dot q_f $ are generalized ve
143   \label{introEquation:hamiltonianDefByLagrangian}
144   \end{equation}
145   where $ \dot q_1  \ldots \dot q_f $ are generalized velocities and
146 < $L$ is the Lagrangian function for the system.
147 <
151 < Differentiating Eq.~\ref{introEquation:hamiltonianDefByLagrangian},
152 < one can obtain
146 > $L$ is the Lagrangian function for the system. Differentiating
147 > Eq.~\ref{introEquation:hamiltonianDefByLagrangian}, one can obtain
148   \begin{equation}
149   dH = \sum\limits_k {\left( {p_k d\dot q_k  + \dot q_k dp_k  -
150   \frac{{\partial L}}{{\partial q_k }}dq_k  - \frac{{\partial
151   L}}{{\partial \dot q_k }}d\dot q_k } \right)}  - \frac{{\partial
152 < L}}{{\partial t}}dt \label{introEquation:diffHamiltonian1}
152 > L}}{{\partial t}}dt . \label{introEquation:diffHamiltonian1}
153   \end{equation}
154 < Making use of  Eq.~\ref{introEquation:generalizedMomenta}, the
155 < second and fourth terms in the parentheses cancel. Therefore,
154 > Making use of Eq.~\ref{introEquation:generalizedMomenta}, the second
155 > and fourth terms in the parentheses cancel. Therefore,
156   Eq.~\ref{introEquation:diffHamiltonian1} can be rewritten as
157   \begin{equation}
158   dH = \sum\limits_k {\left( {\dot q_k dp_k  - \dot p_k dq_k }
159 < \right)}  - \frac{{\partial L}}{{\partial t}}dt
159 > \right)}  - \frac{{\partial L}}{{\partial t}}dt .
160   \label{introEquation:diffHamiltonian2}
161   \end{equation}
162   By identifying the coefficients of $dq_k$, $dp_k$ and dt, we can
# Line 180 | Line 175 | t}}
175   t}}
176   \label{introEquation:motionHamiltonianTime}
177   \end{equation}
178 <
184 < Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
178 > where Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
179   Eq.~\ref{introEquation:motionHamiltonianMomentum} are Hamilton's
180   equation of motion. Due to their symmetrical formula, they are also
181   known as the canonical equations of motions \cite{Goldstein2001}.
# Line 195 | Line 189 | only works with 1st-order differential equations\cite{
189   statistical mechanics and quantum mechanics, since it treats the
190   coordinate and its time derivative as independent variables and it
191   only works with 1st-order differential equations\cite{Marion1990}.
198 –
192   In Newtonian Mechanics, a system described by conservative forces
193 < conserves the total energy \ref{introEquation:energyConservation}.
194 < It follows that Hamilton's equations of motion conserve the total
195 < Hamiltonian.
193 > conserves the total energy
194 > (Eq.~\ref{introEquation:energyConservation}). It follows that
195 > Hamilton's equations of motion conserve the total Hamiltonian
196   \begin{equation}
197   \frac{{dH}}{{dt}} = \sum\limits_i {\left( {\frac{{\partial
198   H}}{{\partial q_i }}\dot q_i  + \frac{{\partial H}}{{\partial p_i
199   }}\dot p_i } \right)}  = \sum\limits_i {\left( {\frac{{\partial
200   H}}{{\partial q_i }}\frac{{\partial H}}{{\partial p_i }} -
201   \frac{{\partial H}}{{\partial p_i }}\frac{{\partial H}}{{\partial
202 < q_i }}} \right) = 0} \label{introEquation:conserveHalmitonian}
202 > q_i }}} \right) = 0}. \label{introEquation:conserveHalmitonian}
203   \end{equation}
204  
205   \section{\label{introSection:statisticalMechanics}Statistical
# Line 221 | Line 214 | Mathematically, phase space is the space which represe
214   \subsection{\label{introSection:ensemble}Phase Space and Ensemble}
215  
216   Mathematically, phase space is the space which represents all
217 < possible states. Each possible state of the system corresponds to
218 < one unique point in the phase space. For mechanical systems, the
219 < phase space usually consists of all possible values of position and
220 < momentum variables. Consider a dynamic system of $f$ particles in a
221 < cartesian space, where each of the $6f$ coordinates and momenta is
222 < assigned to one of $6f$ mutually orthogonal axes, the phase space of
223 < this system is a $6f$ dimensional space. A point, $x = (q_1 , \ldots
224 < ,q_f ,p_1 , \ldots ,p_f )$, with a unique set of values of $6f$
225 < coordinates and momenta is a phase space vector.
226 <
217 > possible states of a system. Each possible state of the system
218 > corresponds to one unique point in the phase space. For mechanical
219 > systems, the phase space usually consists of all possible values of
220 > position and momentum variables. Consider a dynamic system of $f$
221 > particles in a cartesian space, where each of the $6f$ coordinates
222 > and momenta is assigned to one of $6f$ mutually orthogonal axes, the
223 > phase space of this system is a $6f$ dimensional space. A point, $x
224 > =
225 > (\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
226 > \over q} _1 , \ldots
227 > ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
228 > \over q} _f
229 > ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
230 > \over p} _1  \ldots
231 > ,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
232 > \over p} _f )$ , with a unique set of values of $6f$ coordinates and
233 > momenta is a phase space vector.
234   %%%fix me
235 < A microscopic state or microstate of a classical system is
236 < specification of the complete phase space vector of a system at any
237 < instant in time. An ensemble is defined as a collection of systems
238 < sharing one or more macroscopic characteristics but each being in a
239 < unique microstate. The complete ensemble is specified by giving all
240 < systems or microstates consistent with the common macroscopic
241 < characteristics of the ensemble. Although the state of each
242 < individual system in the ensemble could be precisely described at
243 < any instance in time by a suitable phase space vector, when using
244 < ensembles for statistical purposes, there is no need to maintain
245 < distinctions between individual systems, since the numbers of
246 < systems at any time in the different states which correspond to
247 < different regions of the phase space are more interesting. Moreover,
248 < in the point of view of statistical mechanics, one would prefer to
249 < use ensembles containing a large enough population of separate
250 < members so that the numbers of systems in such different states can
251 < be regarded as changing continuously as we traverse different
252 < regions of the phase space. The condition of an ensemble at any time
235 >
236 > In statistical mechanics, the condition of an ensemble at any time
237   can be regarded as appropriately specified by the density $\rho$
238   with which representative points are distributed over the phase
239   space. The density distribution for an ensemble with $f$ degrees of
# Line 259 | Line 243 | Governed by the principles of mechanics, the phase poi
243   \label{introEquation:densityDistribution}
244   \end{equation}
245   Governed by the principles of mechanics, the phase points change
246 < their locations which would change the density at any time at phase
246 > their locations which changes the density at any time at phase
247   space. Hence, the density distribution is also to be taken as a
248 < function of the time.
249 <
266 < The number of systems $\delta N$ at time $t$ can be determined by,
248 > function of the time. The number of systems $\delta N$ at time $t$
249 > can be determined by,
250   \begin{equation}
251   \delta N = \rho (q,p,t)dq_1  \ldots dq_f dp_1  \ldots dp_f.
252   \label{introEquation:deltaN}
253   \end{equation}
254 < Assuming a large enough population of systems, we can sufficiently
254 > Assuming enough copies of the systems, we can sufficiently
255   approximate $\delta N$ without introducing discontinuity when we go
256   from one region in the phase space to another. By integrating over
257   the whole phase space,
# Line 276 | Line 259 | N = \int { \ldots \int {\rho (q,p,t)dq_1 } ...dq_f dp_
259   N = \int { \ldots \int {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f
260   \label{introEquation:totalNumberSystem}
261   \end{equation}
262 < gives us an expression for the total number of the systems. Hence,
263 < the probability per unit in the phase space can be obtained by,
262 > gives us an expression for the total number of copies. Hence, the
263 > probability per unit volume in the phase space can be obtained by,
264   \begin{equation}
265   \frac{{\rho (q,p,t)}}{N} = \frac{{\rho (q,p,t)}}{{\int { \ldots \int
266   {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
# Line 286 | Line 269 | value of any desired quantity which depends on the coo
269   With the help of Eq.~\ref{introEquation:unitProbability} and the
270   knowledge of the system, it is possible to calculate the average
271   value of any desired quantity which depends on the coordinates and
272 < momenta of the system. Even when the dynamics of the real system is
272 > momenta of the system. Even when the dynamics of the real system are
273   complex, or stochastic, or even discontinuous, the average
274 < properties of the ensemble of possibilities as a whole remaining
275 < well defined. For a classical system in thermal equilibrium with its
274 > properties of the ensemble of possibilities as a whole remain well
275 > defined. For a classical system in thermal equilibrium with its
276   environment, the ensemble average of a mechanical quantity, $\langle
277   A(q , p) \rangle_t$, takes the form of an integral over the phase
278   space of the system,
279   \begin{equation}
280   \langle  A(q , p) \rangle_t = \frac{{\int { \ldots \int {A(q,p)\rho
281   (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}{{\int { \ldots \int {\rho
282 < (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}
282 > (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
283   \label{introEquation:ensembelAverage}
284   \end{equation}
285  
303 – There are several different types of ensembles with different
304 – statistical characteristics. As a function of macroscopic
305 – parameters, such as temperature \textit{etc}, the partition function
306 – can be used to describe the statistical properties of a system in
307 – thermodynamic equilibrium.
308 –
309 – As an ensemble of systems, each of which is known to be thermally
310 – isolated and conserve energy, the Microcanonical ensemble (NVE) has
311 – a partition function like,
312 – \begin{equation}
313 – \Omega (N,V,E) = e^{\beta TS} \label{introEquation:NVEPartition}.
314 – \end{equation}
315 – A canonical ensemble (NVT)is an ensemble of systems, each of which
316 – can share its energy with a large heat reservoir. The distribution
317 – of the total energy amongst the possible dynamical states is given
318 – by the partition function,
319 – \begin{equation}
320 – \Omega (N,V,T) = e^{ - \beta A}
321 – \label{introEquation:NVTPartition}
322 – \end{equation}
323 – Here, $A$ is the Helmholtz free energy which is defined as $ A = U -
324 – TS$. Since most experiments are carried out under constant pressure
325 – condition, the isothermal-isobaric ensemble (NPT) plays a very
326 – important role in molecular simulations. The isothermal-isobaric
327 – ensemble allow the system to exchange energy with a heat bath of
328 – temperature $T$ and to change the volume as well. Its partition
329 – function is given as
330 – \begin{equation}
331 – \Delta (N,P,T) =  - e^{\beta G}.
332 – \label{introEquation:NPTPartition}
333 – \end{equation}
334 – Here, $G = U - TS + PV$, and $G$ is called Gibbs free energy.
335 –
286   \subsection{\label{introSection:liouville}Liouville's theorem}
287  
288   Liouville's theorem is the foundation on which statistical mechanics
# Line 374 | Line 324 | simple form,
324   \frac{{\partial \rho }}{{\partial p_i }}\dot p_i } \right)}  = 0 .
325   \label{introEquation:liouvilleTheorem}
326   \end{equation}
377 –
327   Liouville's theorem states that the distribution function is
328   constant along any trajectory in phase space. In classical
329 < statistical mechanics, since the number of members in an ensemble is
330 < huge and constant, we can assume the local density has no reason
331 < (other than classical mechanics) to change,
329 > statistical mechanics, since the number of system copies in an
330 > ensemble is huge and constant, we can assume the local density has
331 > no reason (other than classical mechanics) to change,
332   \begin{equation}
333   \frac{{\partial \rho }}{{\partial t}} = 0.
334   \label{introEquation:stationary}
# Line 409 | Line 358 | dp_1 } ..dp_f.
358   \frac{{d(\delta N)}}{{dt}} = \frac{{d\rho }}{{dt}}\delta v + \rho
359   \frac{d}{{dt}}(\delta v) = 0.
360   \end{equation}
361 < With the help of stationary assumption
362 < (\ref{introEquation:stationary}), we obtain the principle of the
361 > With the help of the stationary assumption
362 > (Eq.~\ref{introEquation:stationary}), we obtain the principle of
363   \emph{conservation of volume in phase space},
364   \begin{equation}
365   \frac{d}{{dt}}(\delta v) = \frac{d}{{dt}}\int { \ldots \int {dq_1 }
# Line 420 | Line 369 | With the help of stationary assumption
369  
370   \subsubsection{\label{introSection:liouvilleInOtherForms}\textbf{Liouville's Theorem in Other Forms}}
371  
372 < Liouville's theorem can be expresses in a variety of different forms
372 > Liouville's theorem can be expressed in a variety of different forms
373   which are convenient within different contexts. For any two function
374   $F$ and $G$ of the coordinates and momenta of a system, the Poisson
375   bracket ${F, G}$ is defined as
# Line 431 | Line 380 | q_i }}} \right)}.
380   q_i }}} \right)}.
381   \label{introEquation:poissonBracket}
382   \end{equation}
383 < Substituting equations of motion in Hamiltonian formalism(
384 < Eq.~\ref{introEquation:motionHamiltonianCoordinate} ,
385 < Eq.~\ref{introEquation:motionHamiltonianMomentum} ) into
383 > Substituting equations of motion in Hamiltonian formalism
384 > (Eq.~\ref{introEquation:motionHamiltonianCoordinate} ,
385 > Eq.~\ref{introEquation:motionHamiltonianMomentum}) into
386   (Eq.~\ref{introEquation:liouvilleTheorem}), we can rewrite
387   Liouville's theorem using Poisson bracket notion,
388   \begin{equation}
# Line 454 | Line 403 | expressed as
403   \left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - iL\rho
404   \label{introEquation:liouvilleTheoremInOperator}
405   \end{equation}
406 <
406 > which can help define a propagator $\rho (t) = e^{-iLt} \rho (0)$.
407   \subsection{\label{introSection:ergodic}The Ergodic Hypothesis}
408  
409   Various thermodynamic properties can be calculated from Molecular
# Line 463 | Line 412 | certain time interval and the measurements are average
412   simulation and the quality of the underlying model. However, both
413   experiments and computer simulations are usually performed during a
414   certain time interval and the measurements are averaged over a
415 < period of them which is different from the average behavior of
415 > period of time which is different from the average behavior of
416   many-body system in Statistical Mechanics. Fortunately, the Ergodic
417   Hypothesis makes a connection between time average and the ensemble
418   average. It states that the time average and average over the
419 < statistical ensemble are identical \cite{Frenkel1996, Leach2001}.
419 > statistical ensemble are identical \cite{Frenkel1996, Leach2001}:
420   \begin{equation}
421   \langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty }
422   \frac{1}{t}\int\limits_0^t {A(q(t),p(t))dt = \int\limits_\Gamma
# Line 476 | Line 425 | distribution function. If an observation is averaged o
425   where $\langle  A(q , p) \rangle_t$ is an equilibrium value of a
426   physical quantity and $\rho (p(t), q(t))$ is the equilibrium
427   distribution function. If an observation is averaged over a
428 < sufficiently long time (longer than relaxation time), all accessible
429 < microstates in phase space are assumed to be equally probed, giving
430 < a properly weighted statistical average. This allows the researcher
431 < freedom of choice when deciding how best to measure a given
432 < observable. In case an ensemble averaged approach sounds most
433 < reasonable, the Monte Carlo techniques\cite{Metropolis1949} can be
428 > sufficiently long time (longer than the relaxation time), all
429 > accessible microstates in phase space are assumed to be equally
430 > probed, giving a properly weighted statistical average. This allows
431 > the researcher freedom of choice when deciding how best to measure a
432 > given observable. In case an ensemble averaged approach sounds most
433 > reasonable, the Monte Carlo methods\cite{Metropolis1949} can be
434   utilized. Or if the system lends itself to a time averaging
435   approach, the Molecular Dynamics techniques in
436   Sec.~\ref{introSection:molecularDynamics} will be the best
# Line 494 | Line 443 | geometric integrators, which preserve various phase-fl
443   by the differential equations. However, most of them ignore the
444   hidden physical laws contained within the equations. Since 1990,
445   geometric integrators, which preserve various phase-flow invariants
446 < such as symplectic structure, volume and time reversal symmetry, are
447 < developed to address this issue\cite{Dullweber1997, McLachlan1998,
448 < Leimkuhler1999}. The velocity Verlet method, which happens to be a
449 < simple example of symplectic integrator, continues to gain
450 < popularity in the molecular dynamics community. This fact can be
451 < partly explained by its geometric nature.
446 > such as symplectic structure, volume and time reversal symmetry,
447 > were developed to address this issue\cite{Dullweber1997,
448 > McLachlan1998, Leimkuhler1999}. The velocity Verlet method, which
449 > happens to be a simple example of symplectic integrator, continues
450 > to gain popularity in the molecular dynamics community. This fact
451 > can be partly explained by its geometric nature.
452  
453   \subsection{\label{introSection:symplecticManifold}Symplectic Manifolds}
454   A \emph{manifold} is an abstract mathematical space. It looks
# Line 508 | Line 457 | viewed as a whole. A \emph{differentiable manifold} (a
457   surface of Earth. It seems to be flat locally, but it is round if
458   viewed as a whole. A \emph{differentiable manifold} (also known as
459   \emph{smooth manifold}) is a manifold on which it is possible to
460 < apply calculus on \emph{differentiable manifold}. A \emph{symplectic
461 < manifold} is defined as a pair $(M, \omega)$ which consists of a
460 > apply calculus\cite{Hirsch1997}. A \emph{symplectic manifold} is
461 > defined as a pair $(M, \omega)$ which consists of a
462   \emph{differentiable manifold} $M$ and a close, non-degenerated,
463   bilinear symplectic form, $\omega$. A symplectic form on a vector
464   space $V$ is a function $\omega(x, y)$ which satisfies
465   $\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+
466   \lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and
467 < $\omega(x, x) = 0$. The cross product operation in vector field is
468 < an example of symplectic form.
467 > $\omega(x, x) = 0$\cite{McDuff1998}. The cross product operation in
468 > vector field is an example of symplectic form. One of the
469 > motivations to study \emph{symplectic manifolds} in Hamiltonian
470 > Mechanics is that a symplectic manifold can represent all possible
471 > configurations of the system and the phase space of the system can
472 > be described by it's cotangent bundle\cite{Jost2002}. Every
473 > symplectic manifold is even dimensional. For instance, in Hamilton
474 > equations, coordinate and momentum always appear in pairs.
475  
521 – One of the motivations to study \emph{symplectic manifolds} in
522 – Hamiltonian Mechanics is that a symplectic manifold can represent
523 – all possible configurations of the system and the phase space of the
524 – system can be described by it's cotangent bundle. Every symplectic
525 – manifold is even dimensional. For instance, in Hamilton equations,
526 – coordinate and momentum always appear in pairs.
527 –
476   \subsection{\label{introSection:ODE}Ordinary Differential Equations}
477  
478   For an ordinary differential system defined as
# Line 532 | Line 480 | where $x = x(q,p)^T$, this system is a canonical Hamil
480   \dot x = f(x)
481   \end{equation}
482   where $x = x(q,p)^T$, this system is a canonical Hamiltonian, if
483 + $f(x) = J\nabla _x H(x)$. Here, $H = H (q, p)$ is Hamiltonian
484 + function and $J$ is the skew-symmetric matrix
485   \begin{equation}
536 – f(r) = J\nabla _x H(r).
537 – \end{equation}
538 – $H = H (q, p)$ is Hamiltonian function and $J$ is the skew-symmetric
539 – matrix
540 – \begin{equation}
486   J = \left( {\begin{array}{*{20}c}
487     0 & I  \\
488     { - I} & 0  \\
# Line 547 | Line 492 | system can be rewritten as,
492   where $I$ is an identity matrix. Using this notation, Hamiltonian
493   system can be rewritten as,
494   \begin{equation}
495 < \frac{d}{{dt}}x = J\nabla _x H(x)
495 > \frac{d}{{dt}}x = J\nabla _x H(x).
496   \label{introEquation:compactHamiltonian}
497   \end{equation}In this case, $f$ is
498 < called a \emph{Hamiltonian vector field}.
499 <
555 < Another generalization of Hamiltonian dynamics is Poisson
556 < Dynamics\cite{Olver1986},
498 > called a \emph{Hamiltonian vector field}. Another generalization of
499 > Hamiltonian dynamics is Poisson Dynamics\cite{Olver1986},
500   \begin{equation}
501   \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
502   \end{equation}
503   The most obvious change being that matrix $J$ now depends on $x$.
504  
505 < \subsection{\label{introSection:exactFlow}Exact Flow}
505 > \subsection{\label{introSection:exactFlow}Exact Propagator}
506  
507 < Let $x(t)$ be the exact solution of the ODE system,
507 > Let $x(t)$ be the exact solution of the ODE
508 > system,
509   \begin{equation}
510 < \frac{{dx}}{{dt}} = f(x) \label{introEquation:ODE}
511 < \end{equation}
512 < The exact flow(solution) $\varphi_\tau$ is defined by
513 < \[
514 < x(t+\tau) =\varphi_\tau(x(t))
510 > \frac{{dx}}{{dt}} = f(x), \label{introEquation:ODE}
511 > \end{equation} we can
512 > define its exact propagator $\varphi_\tau$:
513 > \[ x(t+\tau)
514 > =\varphi_\tau(x(t))
515   \]
516   where $\tau$ is a fixed time step and $\varphi$ is a map from phase
517 < space to itself. The flow has the continuous group property,
517 > space to itself. The propagator has the continuous group property,
518   \begin{equation}
519   \varphi _{\tau _1 }  \circ \varphi _{\tau _2 }  = \varphi _{\tau _1
520   + \tau _2 } .
# Line 579 | Line 523 | In particular,
523   \begin{equation}
524   \varphi _\tau   \circ \varphi _{ - \tau }  = I
525   \end{equation}
526 < Therefore, the exact flow is self-adjoint,
526 > Therefore, the exact propagator is self-adjoint,
527   \begin{equation}
528   \varphi _\tau   = \varphi _{ - \tau }^{ - 1}.
529   \end{equation}
530 < The exact flow can also be written in terms of the of an operator,
530 > The exact propagator can also be written in terms of operator,
531   \begin{equation}
532   \varphi _\tau  (x) = e^{\tau \sum\limits_i {f_i (x)\frac{\partial
533   }{{\partial x_i }}} } (x) \equiv \exp (\tau f)(x).
534   \label{introEquation:exponentialOperator}
535   \end{equation}
536 <
537 < In most cases, it is not easy to find the exact flow $\varphi_\tau$.
538 < Instead, we use an approximate map, $\psi_\tau$, which is usually
539 < called integrator. The order of an integrator $\psi_\tau$ is $p$, if
540 < the Taylor series of $\psi_\tau$ agree to order $p$,
536 > In most cases, it is not easy to find the exact propagator
537 > $\varphi_\tau$. Instead, we use an approximate map, $\psi_\tau$,
538 > which is usually called an integrator. The order of an integrator
539 > $\psi_\tau$ is $p$, if the Taylor series of $\psi_\tau$ agree to
540 > order $p$,
541   \begin{equation}
542   \psi_\tau(x) = x + \tau f(x) + O(\tau^{p+1})
543   \end{equation}
# Line 601 | Line 545 | The hidden geometric properties\cite{Budd1999, Marsden
545   \subsection{\label{introSection:geometricProperties}Geometric Properties}
546  
547   The hidden geometric properties\cite{Budd1999, Marsden1998} of an
548 < ODE and its flow play important roles in numerical studies. Many of
549 < them can be found in systems which occur naturally in applications.
550 <
551 < Let $\varphi$ be the flow of Hamiltonian vector field, $\varphi$ is
608 < a \emph{symplectic} flow if it satisfies,
548 > ODE and its propagator play important roles in numerical studies.
549 > Many of them can be found in systems which occur naturally in
550 > applications. Let $\varphi$ be the propagator of Hamiltonian vector
551 > field, $\varphi$ is a \emph{symplectic} propagator if it satisfies,
552   \begin{equation}
553   {\varphi '}^T J \varphi ' = J.
554   \end{equation}
555   According to Liouville's theorem, the symplectic volume is invariant
556 < under a Hamiltonian flow, which is the basis for classical
557 < statistical mechanics. Furthermore, the flow of a Hamiltonian vector
558 < field on a symplectic manifold can be shown to be a
556 > under a Hamiltonian propagator, which is the basis for classical
557 > statistical mechanics. Furthermore, the propagator of a Hamiltonian
558 > vector field on a symplectic manifold can be shown to be a
559   symplectomorphism. As to the Poisson system,
560   \begin{equation}
561   {\varphi '}^T J \varphi ' = J \circ \varphi
562   \end{equation}
563 < is the property that must be preserved by the integrator.
564 <
565 < It is possible to construct a \emph{volume-preserving} flow for a
566 < source free ODE ($ \nabla \cdot f = 0 $), if the flow satisfies $
567 < \det d\varphi  = 1$. One can show easily that a symplectic flow will
568 < be volume-preserving.
569 <
627 < Changing the variables $y = h(x)$ in an ODE
628 < (Eq.~\ref{introEquation:ODE}) will result in a new system,
563 > is the property that must be preserved by the integrator. It is
564 > possible to construct a \emph{volume-preserving} propagator for a
565 > source free ODE ($ \nabla \cdot f = 0 $), if the propagator
566 > satisfies $ \det d\varphi  = 1$. One can show easily that a
567 > symplectic propagator will be volume-preserving. Changing the
568 > variables $y = h(x)$ in an ODE (Eq.~\ref{introEquation:ODE}) will
569 > result in a new system,
570   \[
571   \dot y = \tilde f(y) = ((dh \cdot f)h^{ - 1} )(y).
572   \]
573   The vector filed $f$ has reversing symmetry $h$ if $f = - \tilde f$.
574 < In other words, the flow of this vector field is reversible if and
575 < only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $.
576 <
577 < A \emph{first integral}, or conserved quantity of a general
578 < differential function is a function $ G:R^{2d}  \to R^d $ which is
638 < constant for all solutions of the ODE $\frac{{dx}}{{dt}} = f(x)$ ,
574 > In other words, the propagator of this vector field is reversible if
575 > and only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $. A
576 > conserved quantity of a general differential function is a function
577 > $ G:R^{2d}  \to R^d $ which is constant for all solutions of the ODE
578 > $\frac{{dx}}{{dt}} = f(x)$ ,
579   \[
580   \frac{{dG(x(t))}}{{dt}} = 0.
581   \]
582 < Using chain rule, one may obtain,
582 > Using the chain rule, one may obtain,
583   \[
584 < \sum\limits_i {\frac{{dG}}{{dx_i }}} f_i (x) = f \bullet \nabla G,
584 > \sum\limits_i {\frac{{dG}}{{dx_i }}} f_i (x) = f \cdot \nabla G,
585   \]
586 < which is the condition for conserving \emph{first integral}. For a
587 < canonical Hamiltonian system, the time evolution of an arbitrary
588 < smooth function $G$ is given by,
649 <
586 > which is the condition for conserved quantities. For a canonical
587 > Hamiltonian system, the time evolution of an arbitrary smooth
588 > function $G$ is given by,
589   \begin{eqnarray}
590 < \frac{{dG(x(t))}}{{dt}} & = & [\nabla _x G(x(t))]^T \dot x(t) \\
591 <                        & = & [\nabla _x G(x(t))]^T J\nabla _x H(x(t)). \\
590 > \frac{{dG(x(t))}}{{dt}} & = & [\nabla _x G(x(t))]^T \dot x(t) \notag\\
591 >                        & = & [\nabla _x G(x(t))]^T J\nabla _x H(x(t)).
592   \label{introEquation:firstIntegral1}
593   \end{eqnarray}
594 <
595 <
657 < Using poisson bracket notion, Equation
658 < \ref{introEquation:firstIntegral1} can be rewritten as
594 > Using poisson bracket notion, Eq.~\ref{introEquation:firstIntegral1}
595 > can be rewritten as
596   \[
597   \frac{d}{{dt}}G(x(t)) = \left\{ {G,H} \right\}(x(t)).
598   \]
599 < Therefore, the sufficient condition for $G$ to be the \emph{first
600 < integral} of a Hamiltonian system is
601 < \[
602 < \left\{ {G,H} \right\} = 0.
666 < \]
667 < As well known, the Hamiltonian (or energy) H of a Hamiltonian system
668 < is a \emph{first integral}, which is due to the fact $\{ H,H\}  =
669 < 0$.
670 <
599 > Therefore, the sufficient condition for $G$ to be a conserved
600 > quantity of a Hamiltonian system is $\left\{ {G,H} \right\} = 0.$ As
601 > is well known, the Hamiltonian (or energy) H of a Hamiltonian system
602 > is a conserved quantity, which is due to the fact $\{ H,H\}  = 0$.
603   When designing any numerical methods, one should always try to
604 < preserve the structural properties of the original ODE and its flow.
604 > preserve the structural properties of the original ODE and its
605 > propagator.
606  
607   \subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
608   A lot of well established and very effective numerical methods have
609 < been successful precisely because of their symplecticities even
609 > been successful precisely because of their symplectic nature even
610   though this fact was not recognized when they were first
611   constructed. The most famous example is the Verlet-leapfrog method
612   in molecular dynamics. In general, symplectic integrators can be
# Line 684 | Line 617 | constructed using one of four different methods.
617   \item Runge-Kutta methods
618   \item Splitting methods
619   \end{enumerate}
620 <
688 < Generating function\cite{Channell1990} tends to lead to methods
620 > Generating functions\cite{Channell1990} tend to lead to methods
621   which are cumbersome and difficult to use. In dissipative systems,
622   variational methods can capture the decay of energy
623 < accurately\cite{Kane2000}. Since their geometrically unstable nature
623 > accurately\cite{Kane2000}. Since they are geometrically unstable
624   against non-Hamiltonian perturbations, ordinary implicit Runge-Kutta
625   methods are not suitable for Hamiltonian system. Recently, various
626 < high-order explicit Runge-Kutta methods
627 < \cite{Owren1992,Chen2003}have been developed to overcome this
628 < instability. However, due to computational penalty involved in
629 < implementing the Runge-Kutta methods, they have not attracted much
630 < attention from the Molecular Dynamics community. Instead, splitting
631 < methods have been widely accepted since they exploit natural
632 < decompositions of the system\cite{Tuckerman1992, McLachlan1998}.
626 > high-order explicit Runge-Kutta methods \cite{Owren1992,Chen2003}
627 > have been developed to overcome this instability. However, due to
628 > computational penalty involved in implementing the Runge-Kutta
629 > methods, they have not attracted much attention from the Molecular
630 > Dynamics community. Instead, splitting methods have been widely
631 > accepted since they exploit natural decompositions of the
632 > system\cite{Tuckerman1992, McLachlan1998}.
633  
634   \subsubsection{\label{introSection:splittingMethod}\textbf{Splitting Methods}}
635  
636   The main idea behind splitting methods is to decompose the discrete
637 < $\varphi_h$ as a composition of simpler flows,
637 > $\varphi_h$ as a composition of simpler propagators,
638   \begin{equation}
639   \varphi _h  = \varphi _{h_1 }  \circ \varphi _{h_2 }  \ldots  \circ
640   \varphi _{h_n }
641   \label{introEquation:FlowDecomposition}
642   \end{equation}
643 < where each of the sub-flow is chosen such that each represent a
644 < simpler integration of the system.
645 <
714 < Suppose that a Hamiltonian system takes the form,
643 > where each of the sub-propagator is chosen such that each represent
644 > a simpler integration of the system. Suppose that a Hamiltonian
645 > system takes the form,
646   \[
647   H = H_1 + H_2.
648   \]
649   Here, $H_1$ and $H_2$ may represent different physical processes of
650   the system. For instance, they may relate to kinetic and potential
651   energy respectively, which is a natural decomposition of the
652 < problem. If $H_1$ and $H_2$ can be integrated using exact flows
653 < $\varphi_1(t)$ and $\varphi_2(t)$, respectively, a simple first
654 < order expression is then given by the Lie-Trotter formula
652 > problem. If $H_1$ and $H_2$ can be integrated using exact
653 > propagators $\varphi_1(t)$ and $\varphi_2(t)$, respectively, a
654 > simple first order expression is then given by the Lie-Trotter
655 > formula
656   \begin{equation}
657   \varphi _h  = \varphi _{1,h}  \circ \varphi _{2,h},
658   \label{introEquation:firstOrderSplitting}
# Line 729 | Line 661 | must follow that each operator $\varphi_i(t)$ is a sym
661   continuous $\varphi _i$ over a time $h$. By definition, as
662   $\varphi_i(t)$ is the exact solution of a Hamiltonian system, it
663   must follow that each operator $\varphi_i(t)$ is a symplectic map.
664 < It is easy to show that any composition of symplectic flows yields a
665 < symplectic map,
664 > It is easy to show that any composition of symplectic propagators
665 > yields a symplectic map,
666   \begin{equation}
667   (\varphi '\phi ')^T J\varphi '\phi ' = \phi '^T \varphi '^T J\varphi
668   '\phi ' = \phi '^T J\phi ' = J,
# Line 738 | Line 670 | splitting in this context automatically generates a sy
670   \end{equation}
671   where $\phi$ and $\psi$ both are symplectic maps. Thus operator
672   splitting in this context automatically generates a symplectic map.
673 <
674 < The Lie-Trotter splitting(\ref{introEquation:firstOrderSplitting})
675 < introduces local errors proportional to $h^2$, while Strang
676 < splitting gives a second-order decomposition,
673 > The Lie-Trotter
674 > splitting(Eq.~\ref{introEquation:firstOrderSplitting}) introduces
675 > local errors proportional to $h^2$, while the Strang splitting gives
676 > a second-order decomposition,
677   \begin{equation}
678   \varphi _h  = \varphi _{1,h/2}  \circ \varphi _{2,h}  \circ \varphi
679   _{1,h/2} , \label{introEquation:secondOrderSplitting}
680   \end{equation}
681 < which has a local error proportional to $h^3$. The Sprang
681 > which has a local error proportional to $h^3$. The Strang
682   splitting's popularity in molecular simulation community attribute
683   to its symmetric property,
684   \begin{equation}
685   \varphi _h^{ - 1} = \varphi _{ - h}.
686   \label{introEquation:timeReversible}
687 < \end{equation},appendixFig:architecture
687 > \end{equation}
688  
689   \subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Examples of the Splitting Method}}
690   The classical equation for a system consisting of interacting
# Line 774 | Line 706 | known as \emph{velocity verlet} which is
706   \end{align}
707   where $F(t)$ is the force at time $t$. This integration scheme is
708   known as \emph{velocity verlet} which is
709 < symplectic(\ref{introEquation:SymplecticFlowComposition}),
710 < time-reversible(\ref{introEquation:timeReversible}) and
711 < volume-preserving (\ref{introEquation:volumePreserving}). These
709 > symplectic(Eq.~\ref{introEquation:SymplecticFlowComposition}),
710 > time-reversible(Eq.~\ref{introEquation:timeReversible}) and
711 > volume-preserving (Eq.~\ref{introEquation:volumePreserving}). These
712   geometric properties attribute to its long-time stability and its
713   popularity in the community. However, the most commonly used
714   velocity verlet integration scheme is written as below,
# Line 801 | Line 733 | the equations of motion would follow:
733  
734   \item Repeat from step 1 with the new position, velocities, and forces assuming the roles of the initial values.
735   \end{enumerate}
804 –
736   By simply switching the order of the propagators in the splitting
737   and composing a new integrator, the \emph{position verlet}
738   integrator, can be generated,
# Line 818 | Line 749 | The Baker-Campbell-Hausdorff formula can be used to de
749   \subsubsection{\label{introSection:errorAnalysis}\textbf{Error Analysis and Higher Order Methods}}
750  
751   The Baker-Campbell-Hausdorff formula can be used to determine the
752 < local error of splitting method in terms of the commutator of the
753 < operators(\ref{introEquation:exponentialOperator}) associated with
754 < the sub-flow. For operators $hX$ and $hY$ which are associated with
755 < $\varphi_1(t)$ and $\varphi_2(t)$ respectively , we have
752 > local error of a splitting method in terms of the commutator of the
753 > operators(Eq.~\ref{introEquation:exponentialOperator}) associated with
754 > the sub-propagator. For operators $hX$ and $hY$ which are associated
755 > with $\varphi_1(t)$ and $\varphi_2(t)$ respectively , we have
756   \begin{equation}
757   \exp (hX + hY) = \exp (hZ)
758   \end{equation}
# Line 830 | Line 761 | hZ = hX + hY + \frac{{h^2 }}{2}[X,Y] + \frac{{h^3 }}{2
761   hZ = hX + hY + \frac{{h^2 }}{2}[X,Y] + \frac{{h^3 }}{2}\left(
762   {[X,[X,Y]] + [Y,[Y,X]]} \right) +  \ldots .
763   \end{equation}
764 < Here, $[X,Y]$ is the commutators of operator $X$ and $Y$ given by
764 > Here, $[X,Y]$ is the commutator of operator $X$ and $Y$ given by
765   \[
766   [X,Y] = XY - YX .
767   \]
768   Applying the Baker-Campbell-Hausdorff formula\cite{Varadarajan1974}
769 < to the Sprang splitting, we can obtain
769 > to the Strang splitting, we can obtain
770   \begin{eqnarray*}
771   \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 \\
772                                     &   & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\
773 <                                   &   & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3[X,[X,Y]]/24 + \ldots )
773 >                                   &   & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3[X,[X,Y]]/24 + \ldots
774 >                                   ).
775   \end{eqnarray*}
776 < Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0,\] the dominant local
777 < error of Spring splitting is proportional to $h^3$. The same
778 < procedure can be applied to a general splitting,  of the form
776 > Since $ [X,Y] + [Y,X] = 0$ and $ [X,X] = 0$, the dominant local
777 > error of Strang splitting is proportional to $h^3$. The same
778 > procedure can be applied to a general splitting of the form
779   \begin{equation}
780   \varphi _{b_m h}^2  \circ \varphi _{a_m h}^1  \circ \varphi _{b_{m -
781   1} h}^2  \circ  \ldots  \circ \varphi _{a_1 h}^1 .
# Line 878 | Line 810 | microscopic behavior can be calculated from the trajec
810   dynamical information. The basic idea of molecular dynamics is that
811   macroscopic properties are related to microscopic behavior and
812   microscopic behavior can be calculated from the trajectories in
813 < simulations. For instance, instantaneous temperature of an
814 < Hamiltonian system of $N$ particle can be measured by
813 > simulations. For instance, instantaneous temperature of a
814 > Hamiltonian system of $N$ particles can be measured by
815   \[
816   T = \sum\limits_{i = 1}^N {\frac{{m_i v_i^2 }}{{fk_B }}}
817   \]
818   where $m_i$ and $v_i$ are the mass and velocity of $i$th particle
819   respectively, $f$ is the number of degrees of freedom, and $k_B$ is
820 < the boltzman constant.
820 > the Boltzman constant.
821  
822   A typical molecular dynamics run consists of three essential steps:
823   \begin{enumerate}
# Line 902 | Line 834 | initialization of a simulation. Sec.~\ref{introSection
834   These three individual steps will be covered in the following
835   sections. Sec.~\ref{introSec:initialSystemSettings} deals with the
836   initialization of a simulation. Sec.~\ref{introSection:production}
837 < will discusse issues in production run.
837 > discusses issues of production runs.
838   Sec.~\ref{introSection:Analysis} provides the theoretical tools for
839 < trajectory analysis.
839 > analysis of trajectories.
840  
841   \subsection{\label{introSec:initialSystemSettings}Initialization}
842  
# Line 916 | Line 848 | purification and crystallization. Even for molecules w
848   thousands of crystal structures of molecules are discovered every
849   year, many more remain unknown due to the difficulties of
850   purification and crystallization. Even for molecules with known
851 < structure, some important information is missing. For example, a
851 > structures, some important information is missing. For example, a
852   missing hydrogen atom which acts as donor in hydrogen bonding must
853 < be added. Moreover, in order to include electrostatic interaction,
853 > be added. Moreover, in order to include electrostatic interactions,
854   one may need to specify the partial charges for individual atoms.
855   Under some circumstances, we may even need to prepare the system in
856   a special configuration. For instance, when studying transport
# Line 938 | Line 870 | systems are strongly anharmonic. Thus, it is often use
870   surface and to locate the local minimum. While converging slowly
871   near the minimum, steepest descent method is extremely robust when
872   systems are strongly anharmonic. Thus, it is often used to refine
873 < structure from crystallographic data. Relied on the gradient or
874 < hessian, advanced methods like Newton-Raphson converge rapidly to a
875 < local minimum, but become unstable if the energy surface is far from
873 > structures from crystallographic data. Relying on the Hessian,
874 > advanced methods like Newton-Raphson converge rapidly to a local
875 > minimum, but become unstable if the energy surface is far from
876   quadratic. Another factor that must be taken into account, when
877   choosing energy minimization method, is the size of the system.
878   Steepest descent and conjugate gradient can deal with models of any
879   size. Because of the limits on computer memory to store the hessian
880 < matrix and the computing power needed to diagonalized these
881 < matrices, most Newton-Raphson methods can not be used with very
950 < large systems.
880 > matrix and the computing power needed to diagonalize these matrices,
881 > most Newton-Raphson methods can not be used with very large systems.
882  
883   \subsubsection{\textbf{Heating}}
884  
885 < Typically, Heating is performed by assigning random velocities
885 > Typically, heating is performed by assigning random velocities
886   according to a Maxwell-Boltzman distribution for a desired
887   temperature. Beginning at a lower temperature and gradually
888   increasing the temperature by assigning larger random velocities, we
889 < end up with setting the temperature of the system to a final
890 < temperature at which the simulation will be conducted. In heating
891 < phase, we should also keep the system from drifting or rotating as a
892 < whole. To do this, the net linear momentum and angular momentum of
893 < the system is shifted to zero after each resampling from the Maxwell
894 < -Boltzman distribution.
889 > end up setting the temperature of the system to a final temperature
890 > at which the simulation will be conducted. In heating phase, we
891 > should also keep the system from drifting or rotating as a whole. To
892 > do this, the net linear momentum and angular momentum of the system
893 > is shifted to zero after each resampling from the Maxwell -Boltzman
894 > distribution.
895  
896   \subsubsection{\textbf{Equilibration}}
897  
# Line 971 | Line 902 | equilibration process is long enough. However, these s
902   properties \textit{etc}, become independent of time. Strictly
903   speaking, minimization and heating are not necessary, provided the
904   equilibration process is long enough. However, these steps can serve
905 < as a means to arrive at an equilibrated structure in an effective
905 > as a mean to arrive at an equilibrated structure in an effective
906   way.
907  
908   \subsection{\label{introSection:production}Production}
# Line 987 | Line 918 | complexity of the algorithm for pair-wise interactions
918   calculation of non-bonded forces, such as van der Waals force and
919   Coulombic forces \textit{etc}. For a system of $N$ particles, the
920   complexity of the algorithm for pair-wise interactions is $O(N^2 )$,
921 < which making large simulations prohibitive in the absence of any
922 < algorithmic tricks.
923 <
924 < A natural approach to avoid system size issues is to represent the
925 < bulk behavior by a finite number of the particles. However, this
926 < approach will suffer from the surface effect at the edges of the
927 < simulation. To offset this, \textit{Periodic boundary conditions}
928 < (see Fig.~\ref{introFig:pbc}) is developed to simulate bulk
929 < properties with a relatively small number of particles. In this
930 < method, the simulation box is replicated throughout space to form an
931 < infinite lattice. During the simulation, when a particle moves in
932 < the primary cell, its image in other cells move in exactly the same
933 < direction with exactly the same orientation. Thus, as a particle
1003 < leaves the primary cell, one of its images will enter through the
1004 < opposite face.
921 > which makes large simulations prohibitive in the absence of any
922 > algorithmic tricks. A natural approach to avoid system size issues
923 > is to represent the bulk behavior by a finite number of the
924 > particles. However, this approach will suffer from surface effects
925 > at the edges of the simulation. To offset this, \textit{Periodic
926 > boundary conditions} (see Fig.~\ref{introFig:pbc}) were developed to
927 > simulate bulk properties with a relatively small number of
928 > particles. In this method, the simulation box is replicated
929 > throughout space to form an infinite lattice. During the simulation,
930 > when a particle moves in the primary cell, its image in other cells
931 > move in exactly the same direction with exactly the same
932 > orientation. Thus, as a particle leaves the primary cell, one of its
933 > images will enter through the opposite face.
934   \begin{figure}
935   \centering
936   \includegraphics[width=\linewidth]{pbc.eps}
# Line 1013 | Line 942 | Another important technique to improve the efficiency
942  
943   %cutoff and minimum image convention
944   Another important technique to improve the efficiency of force
945 < evaluation is to apply spherical cutoff where particles farther than
946 < a predetermined distance are not included in the calculation
945 > evaluation is to apply spherical cutoffs where particles farther
946 > than a predetermined distance are not included in the calculation
947   \cite{Frenkel1996}. The use of a cutoff radius will cause a
948   discontinuity in the potential energy curve. Fortunately, one can
949 < shift simple radial potential to ensure the potential curve go
949 > shift a simple radial potential to ensure the potential curve go
950   smoothly to zero at the cutoff radius. The cutoff strategy works
951   well for Lennard-Jones interaction because of its short range
952   nature. However, simply truncating the electrostatic interaction
# Line 1043 | Line 972 | R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha
972   V(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha
973   r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow
974   R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha
975 < r_{ij})}{r_{ij}}\right\}. \label{introEquation:shiftedCoulomb}
975 > r_{ij})}{r_{ij}}\right\}, \label{introEquation:shiftedCoulomb}
976   \end{equation}
977   where $\alpha$ is the convergence parameter. Due to the lack of
978   inherent periodicity and rapid convergence,this method is extremely
# Line 1060 | Line 989 | illustration of shifted Coulomb potential.}
989  
990   \subsection{\label{introSection:Analysis} Analysis}
991  
992 < Recently, advanced visualization technique have become applied to
992 > Recently, advanced visualization techniques have been applied to
993   monitor the motions of molecules. Although the dynamics of the
994   system can be described qualitatively from animation, quantitative
995 < trajectory analysis are more useful. According to the principles of
996 < Statistical Mechanics, Sec.~\ref{introSection:statisticalMechanics},
997 < one can compute thermodynamic properties, analyze fluctuations of
998 < structural parameters, and investigate time-dependent processes of
999 < the molecule from the trajectories.
995 > trajectory analysis is more useful. According to the principles of
996 > Statistical Mechanics in
997 > Sec.~\ref{introSection:statisticalMechanics}, one can compute
998 > thermodynamic properties, analyze fluctuations of structural
999 > parameters, and investigate time-dependent processes of the molecule
1000 > from the trajectories.
1001  
1002   \subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamic Properties}}
1003  
# Line 1097 | Line 1027 | function}, is of most fundamental importance to liquid
1027   distribution functions. Among these functions,the \emph{pair
1028   distribution function}, also known as \emph{radial distribution
1029   function}, is of most fundamental importance to liquid theory.
1030 < Experimentally, pair distribution function can be gathered by
1030 > Experimentally, pair distribution functions can be gathered by
1031   Fourier transforming raw data from a series of neutron diffraction
1032   experiments and integrating over the surface factor
1033   \cite{Powles1973}. The experimental results can serve as a criterion
1034   to justify the correctness of a liquid model. Moreover, various
1035   equilibrium thermodynamic and structural properties can also be
1036 < expressed in terms of radial distribution function \cite{Allen1987}.
1037 <
1038 < The pair distribution functions $g(r)$ gives the probability that a
1039 < particle $i$ will be located at a distance $r$ from a another
1040 < particle $j$ in the system
1111 < \[
1036 > expressed in terms of the radial distribution function
1037 > \cite{Allen1987}. The pair distribution functions $g(r)$ gives the
1038 > probability that a particle $i$ will be located at a distance $r$
1039 > from a another particle $j$ in the system
1040 > \begin{equation}
1041   g(r) = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j
1042   \ne i} {\delta (r - r_{ij} )} } } \right\rangle = \frac{\rho
1043   (r)}{\rho}.
1044 < \]
1044 > \end{equation}
1045   Note that the delta function can be replaced by a histogram in
1046   computer simulation. Peaks in $g(r)$ represent solvent shells, and
1047   the height of these peaks gradually decreases to 1 as the liquid of
# Line 1130 | Line 1059 | If $A$ and $B$ refer to same variable, this kind of co
1059   \label{introEquation:timeCorrelationFunction}
1060   \end{equation}
1061   If $A$ and $B$ refer to same variable, this kind of correlation
1062 < function is called an \emph{autocorrelation function}. One example
1063 < of an auto correlation function is the velocity auto-correlation
1062 > functions are called \emph{autocorrelation functions}. One example
1063 > of auto correlation function is the velocity auto-correlation
1064   function which is directly related to transport properties of
1065   molecular liquids:
1066   \[
# Line 1139 | Line 1068 | where $D$ is diffusion constant. Unlike the velocity a
1068   \right\rangle } dt
1069   \]
1070   where $D$ is diffusion constant. Unlike the velocity autocorrelation
1071 < function, which is averaging over time origins and over all the
1072 < atoms, the dipole autocorrelation functions are calculated for the
1071 > function, which is averaged over time origins and over all the
1072 > atoms, the dipole autocorrelation functions is calculated for the
1073   entire system. The dipole autocorrelation function is given by:
1074   \[
1075   c_{dipole}  = \left\langle {u_{tot} (t) \cdot u_{tot} (t)}
# Line 1149 | Line 1078 | by
1078   Here $u_{tot}$ is the net dipole of the entire system and is given
1079   by
1080   \[
1081 < u_{tot} (t) = \sum\limits_i {u_i (t)}
1081 > u_{tot} (t) = \sum\limits_i {u_i (t)}.
1082   \]
1083 < In principle, many time correlation functions can be related with
1083 > In principle, many time correlation functions can be related to
1084   Fourier transforms of the infrared, Raman, and inelastic neutron
1085   scattering spectra of molecular liquids. In practice, one can
1086 < extract the IR spectrum from the intensity of dipole fluctuation at
1087 < each frequency using the following relationship:
1086 > extract the IR spectrum from the intensity of the molecular dipole
1087 > fluctuation at each frequency using the following relationship:
1088   \[
1089   \hat c_{dipole} (v) = \int_{ - \infty }^\infty  {c_{dipole} (t)e^{ -
1090 < i2\pi vt} dt}
1090 > i2\pi vt} dt}.
1091   \]
1092  
1093   \section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
1094  
1095   Rigid bodies are frequently involved in the modeling of different
1096   areas, from engineering, physics, to chemistry. For example,
1097 < missiles and vehicle are usually modeled by rigid bodies.  The
1098 < movement of the objects in 3D gaming engine or other physics
1099 < simulator is governed by rigid body dynamics. In molecular
1097 > missiles and vehicles are usually modeled by rigid bodies.  The
1098 > movement of the objects in 3D gaming engines or other physics
1099 > simulators is governed by rigid body dynamics. In molecular
1100   simulations, rigid bodies are used to simplify protein-protein
1101   docking studies\cite{Gray2003}.
1102  
# Line 1176 | Line 1105 | rotational degrees of freedom. However, due to $\frac
1105   freedom. Euler angles are the natural choice to describe the
1106   rotational degrees of freedom. However, due to $\frac {1}{sin
1107   \theta}$ singularities, the numerical integration of corresponding
1108 < equations of motion is very inefficient and inaccurate. Although an
1109 < alternative integrator using multiple sets of Euler angles can
1110 < overcome this difficulty\cite{Barojas1973}, the computational
1111 < penalty and the loss of angular momentum conservation still remain.
1112 < A singularity-free representation utilizing quaternions was
1113 < developed by Evans in 1977\cite{Evans1977}. Unfortunately, this
1114 < approach uses a nonseparable Hamiltonian resulting from the
1115 < quaternion representation, which prevents the symplectic algorithm
1116 < to be utilized. Another different approach is to apply holonomic
1117 < constraints to the atoms belonging to the rigid body. Each atom
1118 < moves independently under the normal forces deriving from potential
1119 < energy and constraint forces which are used to guarantee the
1120 < rigidness. However, due to their iterative nature, the SHAKE and
1121 < Rattle algorithms also converge very slowly when the number of
1122 < constraints increases\cite{Ryckaert1977, Andersen1983}.
1108 > equations of these motion is very inefficient and inaccurate.
1109 > Although an alternative integrator using multiple sets of Euler
1110 > angles can overcome this difficulty\cite{Barojas1973}, the
1111 > computational penalty and the loss of angular momentum conservation
1112 > still remain. A singularity-free representation utilizing
1113 > quaternions was developed by Evans in 1977\cite{Evans1977}.
1114 > Unfortunately, this approach used a nonseparable Hamiltonian
1115 > resulting from the quaternion representation, which prevented the
1116 > symplectic algorithm from being utilized. Another different approach
1117 > is to apply holonomic constraints to the atoms belonging to the
1118 > rigid body. Each atom moves independently under the normal forces
1119 > deriving from potential energy and constraint forces which are used
1120 > to guarantee the rigidness. However, due to their iterative nature,
1121 > the SHAKE and Rattle algorithms also converge very slowly when the
1122 > number of constraints increases\cite{Ryckaert1977, Andersen1983}.
1123  
1124   A break-through in geometric literature suggests that, in order to
1125   develop a long-term integration scheme, one should preserve the
1126 < symplectic structure of the flow. By introducing a conjugate
1126 > symplectic structure of the propagator. By introducing a conjugate
1127   momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's
1128   equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was
1129   proposed to evolve the Hamiltonian system in a constraint manifold
# Line 1202 | Line 1131 | methods are iterative and inefficient. In this section
1131   An alternative method using the quaternion representation was
1132   developed by Omelyan\cite{Omelyan1998}. However, both of these
1133   methods are iterative and inefficient. In this section, we descibe a
1134 < symplectic Lie-Poisson integrator for rigid body developed by
1134 > symplectic Lie-Poisson integrator for rigid bodies developed by
1135   Dullweber and his coworkers\cite{Dullweber1997} in depth.
1136  
1137   \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Bodies}
# Line 1213 | Line 1142 | V(q,Q) + \frac{1}{2}tr[(QQ^T  - 1)\Lambda ].
1142   V(q,Q) + \frac{1}{2}tr[(QQ^T  - 1)\Lambda ].
1143   \label{introEquation:RBHamiltonian}
1144   \end{equation}
1145 < Here, $q$ and $Q$  are the position and rotation matrix for the
1146 < rigid-body, $p$ and $P$  are conjugate momenta to $q$  and $Q$ , and
1147 < $J$, a diagonal matrix, is defined by
1145 > Here, $q$ and $Q$  are the position vector and rotation matrix for
1146 > the rigid-body, $p$ and $P$  are conjugate momenta to $q$  and $Q$ ,
1147 > and $J$, a diagonal matrix, is defined by
1148   \[
1149   I_{ii}^{ - 1}  = \frac{1}{2}\sum\limits_{i \ne j} {J_{jj}^{ - 1} }
1150   \]
# Line 1225 | Line 1154 | Q^T Q = 1, \label{introEquation:orthogonalConstraint}
1154   \begin{equation}
1155   Q^T Q = 1, \label{introEquation:orthogonalConstraint}
1156   \end{equation}
1157 < which is used to ensure rotation matrix's unitarity. Differentiating
1158 < \ref{introEquation:orthogonalConstraint} and using Equation
1159 < \ref{introEquation:RBMotionMomentum}, one may obtain,
1231 < \begin{equation}
1232 < Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
1233 < \label{introEquation:RBFirstOrderConstraint}
1234 < \end{equation}
1235 <
1236 < Using Equation (\ref{introEquation:motionHamiltonianCoordinate},
1237 < \ref{introEquation:motionHamiltonianMomentum}), one can write down
1157 > which is used to ensure the rotation matrix's unitarity. Using
1158 > Eq.~\ref{introEquation:motionHamiltonianCoordinate} and Eq.~
1159 > \ref{introEquation:motionHamiltonianMomentum}, one can write down
1160   the equations of motion,
1239 –
1161   \begin{eqnarray}
1162 < \frac{{dq}}{{dt}} & = & \frac{p}{m} \label{introEquation:RBMotionPosition}\\
1163 < \frac{{dp}}{{dt}} & = & - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
1164 < \frac{{dQ}}{{dt}} & = & PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
1162 > \frac{{dq}}{{dt}} & = & \frac{p}{m}, \label{introEquation:RBMotionPosition}\\
1163 > \frac{{dp}}{{dt}} & = & - \nabla _q V(q,Q), \label{introEquation:RBMotionMomentum}\\
1164 > \frac{{dQ}}{{dt}} & = & PJ^{ - 1},  \label{introEquation:RBMotionRotation}\\
1165   \frac{{dP}}{{dt}} & = & - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}
1166   \end{eqnarray}
1167 <
1167 > Differentiating Eq.~\ref{introEquation:orthogonalConstraint} and
1168 > using Eq.~\ref{introEquation:RBMotionMomentum}, one may obtain,
1169 > \begin{equation}
1170 > Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
1171 > \label{introEquation:RBFirstOrderConstraint}
1172 > \end{equation}
1173   In general, there are two ways to satisfy the holonomic constraints.
1174   We can use a constraint force provided by a Lagrange multiplier on
1175 < the normal manifold to keep the motion on constraint space. Or we
1176 < can simply evolve the system on the constraint manifold. These two
1177 < methods have been proved to be equivalent. The holonomic constraint
1178 < and equations of motions define a constraint manifold for rigid
1179 < bodies
1175 > the normal manifold to keep the motion on the constraint space. Or
1176 > we can simply evolve the system on the constraint manifold. These
1177 > two methods have been proved to be equivalent. The holonomic
1178 > constraint and equations of motions define a constraint manifold for
1179 > rigid bodies
1180   \[
1181   M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0}
1182   \right\}.
1183   \]
1184 <
1185 < Unfortunately, this constraint manifold is not the cotangent bundle
1186 < $T_{\star}SO(3)$. However, it turns out that under symplectic
1187 < transformation, the cotangent space and the phase space are
1262 < diffeomorphic. By introducing
1184 > Unfortunately, this constraint manifold is not $T^* SO(3)$ which is
1185 > a symplectic manifold on Lie rotation group $SO(3)$. However, it
1186 > turns out that under symplectic transformation, the cotangent space
1187 > and the phase space are diffeomorphic. By introducing
1188   \[
1189   \tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
1190   \]
1191 < the mechanical system subject to a holonomic constraint manifold $M$
1191 > the mechanical system subjected to a holonomic constraint manifold $M$
1192   can be re-formulated as a Hamiltonian system on the cotangent space
1193   \[
1194   T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \tilde Q =
1195   1,\tilde Q^T \tilde PJ^{ - 1}  + J^{ - 1} P^T \tilde Q = 0} \right\}
1196   \]
1272 –
1197   For a body fixed vector $X_i$ with respect to the center of mass of
1198   the rigid body, its corresponding lab fixed vector $X_0^{lab}$  is
1199   given as
# Line 1288 | Line 1212 | and
1212   \[
1213   \nabla _Q V(q,Q) = F(q,Q)X_i^t
1214   \]
1215 < respectively.
1216 <
1217 < As a common choice to describe the rotation dynamics of the rigid
1294 < body, the angular momentum on the body fixed frame $\Pi  = Q^t P$ is
1295 < introduced to rewrite the equations of motion,
1215 > respectively. As a common choice to describe the rotation dynamics
1216 > of the rigid body, the angular momentum on the body fixed frame $\Pi
1217 > = Q^t P$ is introduced to rewrite the equations of motion,
1218   \begin{equation}
1219   \begin{array}{l}
1220 < \mathop \Pi \limits^ \bullet   = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda  \\
1221 < \mathop Q\limits^{{\rm{   }} \bullet }  = Q\Pi {\rm{ }}J^{ - 1}  \\
1220 > \dot \Pi  = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda,  \\
1221 > \dot Q  = Q\Pi {\rm{ }}J^{ - 1},  \\
1222   \end{array}
1223   \label{introEqaution:RBMotionPI}
1224   \end{equation}
1225 < , as well as holonomic constraints,
1226 < \[
1227 < \begin{array}{l}
1306 < \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
1307 < Q^T Q = 1 \\
1308 < \end{array}
1309 < \]
1310 <
1311 < For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a matrix $\hat v \in
1312 < so(3)^ \star$, the hat-map isomorphism,
1225 > as well as holonomic constraints $\Pi J^{ - 1}  + J^{ - 1} \Pi ^t  =
1226 > 0$ and $Q^T Q = 1$. For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a
1227 > matrix $\hat v \in so(3)^ \star$, the hat-map isomorphism,
1228   \begin{equation}
1229   v(v_1 ,v_2 ,v_3 ) \Leftrightarrow \hat v = \left(
1230   {\begin{array}{*{20}c}
# Line 1322 | Line 1237 | operations
1237   will let us associate the matrix products with traditional vector
1238   operations
1239   \[
1240 < \hat vu = v \times u
1240 > \hat vu = v \times u.
1241   \]
1242 < Using \ref{introEqaution:RBMotionPI}, one can construct a skew
1242 > Using Eq.~\ref{introEqaution:RBMotionPI}, one can construct a skew
1243   matrix,
1244 + \begin{eqnarray}
1245 + (\dot \Pi  - \dot \Pi ^T )&= &(\Pi  - \Pi ^T )(J^{ - 1} \Pi  + \Pi J^{ - 1} ) \notag \\
1246 + & & + \sum\limits_i {[Q^T F_i (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]}  -
1247 + (\Lambda  - \Lambda ^T ). \label{introEquation:skewMatrixPI}
1248 + \end{eqnarray}
1249 + Since $\Lambda$ is symmetric, the last term of
1250 + Eq.~\ref{introEquation:skewMatrixPI} is zero, which implies the
1251 + Lagrange multiplier $\Lambda$ is absent from the equations of
1252 + motion. This unique property eliminates the requirement of
1253 + iterations which can not be avoided in other methods\cite{Kol1997,
1254 + Omelyan1998}. Applying the hat-map isomorphism, we obtain the
1255 + equation of motion for angular momentum in the body frame
1256   \begin{equation}
1330 – (\mathop \Pi \limits^ \bullet   - \mathop \Pi \limits^ {\bullet  ^T}
1331 – ){\rm{ }} = {\rm{ }}(\Pi  - \Pi ^T ){\rm{ }}(J^{ - 1} \Pi  + \Pi J^{
1332 – - 1} ) + \sum\limits_i {[Q^T F_i (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]} -
1333 – (\Lambda  - \Lambda ^T ) . \label{introEquation:skewMatrixPI}
1334 – \end{equation}
1335 – Since $\Lambda$ is symmetric, the last term of Equation
1336 – \ref{introEquation:skewMatrixPI} is zero, which implies the Lagrange
1337 – multiplier $\Lambda$ is absent from the equations of motion. This
1338 – unique property eliminates the requirement of iterations which can
1339 – not be avoided in other methods\cite{Kol1997, Omelyan1998}.
1340 –
1341 – Applying the hat-map isomorphism, we obtain the equation of motion
1342 – for angular momentum on body frame
1343 – \begin{equation}
1257   \dot \pi  = \pi  \times I^{ - 1} \pi  + \sum\limits_i {\left( {Q^T
1258   F_i (r,Q)} \right) \times X_i }.
1259   \label{introEquation:bodyAngularMotion}
# Line 1348 | Line 1261 | given by
1261   In the same manner, the equation of motion for rotation matrix is
1262   given by
1263   \[
1264 < \dot Q = Qskew(I^{ - 1} \pi )
1264 > \dot Q = Qskew(I^{ - 1} \pi ).
1265   \]
1266  
1267   \subsection{\label{introSection:SymplecticFreeRB}Symplectic
1268 < Lie-Poisson Integrator for Free Rigid Body}
1268 > Lie-Poisson Integrator for Free Rigid Bodies}
1269  
1270   If there are no external forces exerted on the rigid body, the only
1271   contribution to the rotational motion is from the kinetic energy
# Line 1370 | Line 1283 | J(\pi ) = \left( {\begin{array}{*{20}c}
1283     0 & {\pi _3 } & { - \pi _2 }  \\
1284     { - \pi _3 } & 0 & {\pi _1 }  \\
1285     {\pi _2 } & { - \pi _1 } & 0  \\
1286 < \end{array}} \right)
1286 > \end{array}} \right).
1287   \end{equation}
1288   Thus, the dynamics of free rigid body is governed by
1289   \begin{equation}
1290 < \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
1290 > \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi ).
1291   \end{equation}
1292 <
1293 < One may notice that each $T_i^r$ in Equation
1294 < \ref{introEquation:rotationalKineticRB} can be solved exactly. For
1382 < instance, the equations of motion due to $T_1^r$ are given by
1292 > One may notice that each $T_i^r$ in
1293 > Eq.~\ref{introEquation:rotationalKineticRB} can be solved exactly.
1294 > For instance, the equations of motion due to $T_1^r$ are given by
1295   \begin{equation}
1296   \frac{d}{{dt}}\pi  = R_1 \pi ,\frac{d}{{dt}}Q = QR_1
1297   \label{introEqaution:RBMotionSingleTerm}
1298   \end{equation}
1299 < where
1299 > with
1300   \[ R_1  = \left( {\begin{array}{*{20}c}
1301     0 & 0 & 0  \\
1302     0 & 0 & {\pi _1 }  \\
1303     0 & { - \pi _1 } & 0  \\
1304   \end{array}} \right).
1305   \]
1306 < The solutions of Equation \ref{introEqaution:RBMotionSingleTerm} is
1306 > The solutions of Eq.~\ref{introEqaution:RBMotionSingleTerm} is
1307   \[
1308   \pi (\Delta t) = e^{\Delta tR_1 } \pi (0),Q(\Delta t) =
1309   Q(0)e^{\Delta tR_1 }
# Line 1405 | Line 1317 | To reduce the cost of computing expensive functions in
1317   \end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
1318   \]
1319   To reduce the cost of computing expensive functions in $e^{\Delta
1320 < tR_1 }$, we can use Cayley transformation to obtain a single-aixs
1321 < propagator,
1322 < \[
1323 < e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
1324 < )
1325 < \]
1326 < The flow maps for $T_2^r$ and $T_3^r$ can be found in the same
1320 > tR_1 }$, we can use the Cayley transformation to obtain a
1321 > single-aixs propagator,
1322 > \begin{eqnarray*}
1323 > e^{\Delta tR_1 }  & \approx & (1 - \Delta tR_1 )^{ - 1} (1 + \Delta
1324 > tR_1 ) \\
1325 > %
1326 > & \approx & \left( \begin{array}{ccc}
1327 > 1 & 0 & 0 \\
1328 > 0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
1329 > \theta^2 / 4} \\
1330 > 0 & \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 +
1331 > \theta^2 / 4}
1332 > \end{array}
1333 > \right).
1334 > \end{eqnarray*}
1335 > The propagators for $T_2^r$ and $T_3^r$ can be found in the same
1336   manner. In order to construct a second-order symplectic method, we
1337 < split the angular kinetic Hamiltonian function can into five terms
1337 > split the angular kinetic Hamiltonian function into five terms
1338   \[
1339   T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
1340   ) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
# Line 1427 | Line 1348 | _1 }.
1348   \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
1349   _1 }.
1350   \]
1430 –
1351   The non-canonical Lie-Poisson bracket ${F, G}$ of two function
1352   $F(\pi )$ and $G(\pi )$ is defined by
1353   \[
1354   \{ F,G\} (\pi ) = [\nabla _\pi  F(\pi )]^T J(\pi )\nabla _\pi  G(\pi
1355 < )
1355 > ).
1356   \]
1357   If the Poisson bracket of a function $F$ with an arbitrary smooth
1358   function $G$ is zero, $F$ is a \emph{Casimir}, which is the
1359   conserved quantity in Poisson system. We can easily verify that the
1360   norm of the angular momentum, $\parallel \pi
1361 < \parallel$, is a \emph{Casimir}. Let$ F(\pi ) = S(\frac{{\parallel
1361 > \parallel$, is a \emph{Casimir}\cite{McLachlan1993}. Let$ F(\pi ) = S(\frac{{\parallel
1362   \pi \parallel ^2 }}{2})$ for an arbitrary function $ S:R \to R$ ,
1363   then by the chain rule
1364   \[
1365   \nabla _\pi  F(\pi ) = S'(\frac{{\parallel \pi \parallel ^2
1366 < }}{2})\pi
1366 > }}{2})\pi.
1367   \]
1368 < Thus $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel \pi
1368 > Thus, $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel
1369 > \pi
1370   \parallel ^2 }}{2})\pi  \times \pi  = 0 $. This explicit
1371   Lie-Poisson integrator is found to be both extremely efficient and
1372   stable. These properties can be explained by the fact the small
# Line 1456 | Line 1377 | The Hamiltonian of rigid body can be separated in term
1377   Splitting for Rigid Body}
1378  
1379   The Hamiltonian of rigid body can be separated in terms of kinetic
1380 < energy and potential energy,
1381 < \[
1382 < H = T(p,\pi ) + V(q,Q)
1462 < \]
1463 < The equations of motion corresponding to potential energy and
1464 < kinetic energy are listed in the below table,
1380 > energy and potential energy, $H = T(p,\pi ) + V(q,Q)$. The equations
1381 > of motion corresponding to potential energy and kinetic energy are
1382 > listed in Table~\ref{introTable:rbEquations}
1383   \begin{table}
1384 < \caption{Equations of motion due to Potential and Kinetic Energies}
1384 > \caption{EQUATIONS OF MOTION DUE TO POTENTIAL AND KINETIC ENERGIES}
1385 > \label{introTable:rbEquations}
1386   \begin{center}
1387   \begin{tabular}{|l|l|}
1388    \hline
# Line 1499 | Line 1418 | where $ T^t (p) = \frac{1}{2}p^T m^{ - 1} p $ and $T^r
1418   T(p,\pi ) =T^t (p) + T^r (\pi ).
1419   \end{equation}
1420   where $ T^t (p) = \frac{1}{2}p^T m^{ - 1} p $ and $T^r (\pi )$ is
1421 < defined by \ref{introEquation:rotationalKineticRB}. Therefore, the
1422 < corresponding propagators are given by
1421 > defined by Eq.~\ref{introEquation:rotationalKineticRB}. Therefore,
1422 > the corresponding propagators are given by
1423   \[
1424   \varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
1425   _{\Delta t,T^r }.
1426   \]
1427   Finally, we obtain the overall symplectic propagators for freely
1428   moving rigid bodies
1429 < \begin{equation}
1430 < \begin{array}{c}
1431 < \varphi _{\Delta t}  = \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \\
1432 <  \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 }  \\
1514 <  \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1515 < \end{array}
1429 > \begin{eqnarray}
1430 > \varphi _{\Delta t}  &=& \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \notag\\
1431 >  & & \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 }  \notag\\
1432 >  & & \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .
1433   \label{introEquation:overallRBFlowMaps}
1434 < \end{equation}
1434 > \end{eqnarray}
1435  
1436   \section{\label{introSection:langevinDynamics}Langevin Dynamics}
1437   As an alternative to newtonian dynamics, Langevin dynamics, which
# Line 1560 | Line 1477 | W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\a
1477   \[
1478   W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
1479   }}{{2m_\alpha  w_\alpha ^2 }}} x^2
1480 < \] and combining the last two terms in Equation
1481 < \ref{introEquation:bathGLE}, we may rewrite the Harmonic bath
1565 < Hamiltonian as
1480 > \]
1481 > and combining the last two terms in Eq.~\ref{introEquation:bathGLE}, we may rewrite the Harmonic bath Hamiltonian as
1482   \[
1483   H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
1484   {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
1485   w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
1486 < w_\alpha ^2 }}x} \right)^2 } \right\}}
1486 > w_\alpha ^2 }}x} \right)^2 } \right\}}.
1487   \]
1488   Since the first two terms of the new Hamiltonian depend only on the
1489   system coordinates, we can get the equations of motion for
# Line 1584 | Line 1500 | m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x
1500   \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right).
1501   \label{introEquation:bathMotionGLE}
1502   \end{equation}
1587 –
1503   In order to derive an equation for $x$, the dynamics of the bath
1504   variables $x_\alpha$ must be solved exactly first. As an integral
1505   transform which is particularly useful in solving linear ordinary
1506   differential equations,the Laplace transform is the appropriate tool
1507   to solve this problem. The basic idea is to transform the difficult
1508   differential equations into simple algebra problems which can be
1509 < solved easily. Then, by applying the inverse Laplace transform, also
1510 < known as the Bromwich integral, we can retrieve the solutions of the
1511 < original problems.
1512 <
1598 < Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace
1599 < transform of f(t) is a new function defined as
1509 > solved easily. Then, by applying the inverse Laplace transform, we
1510 > can retrieve the solutions of the original problems. Let $f(t)$ be a
1511 > function defined on $ [0,\infty ) $, the Laplace transform of $f(t)$
1512 > is a new function defined as
1513   \[
1514   L(f(t)) \equiv F(p) = \int_0^\infty  {f(t)e^{ - pt} dt}
1515   \]
1516   where  $p$ is real and  $L$ is called the Laplace Transform
1517   Operator. Below are some important properties of Laplace transform
1605 –
1518   \begin{eqnarray*}
1519   L(x + y)  & = & L(x) + L(y) \\
1520   L(ax)     & = & aL(x) \\
# Line 1610 | Line 1522 | Operator. Below are some important properties of Lapla
1522   L(\ddot x)& = & p^2 L(x) - px(0) - \dot x(0) \\
1523   L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right)& = & G(p)H(p) \\
1524   \end{eqnarray*}
1613 –
1614 –
1525   Applying the Laplace transform to the bath coordinates, we obtain
1526   \begin{eqnarray*}
1527 < 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) \\
1528 < L(x_\alpha  ) & = & \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
1527 > 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), \\
1528 > L(x_\alpha  ) & = & \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}. \\
1529   \end{eqnarray*}
1530 <
1621 < By the same way, the system coordinates become
1530 > In the same way, the system coordinates become
1531   \begin{eqnarray*}
1532 < mL(\ddot x) & = & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
1533 <  & & \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\}}  \\
1532 > mL(\ddot x) & = &
1533 >  - \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\}}  \\
1534 >  & & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}}.
1535   \end{eqnarray*}
1626 –
1536   With the help of some relatively important inverse Laplace
1537   transformations:
1538   \[
# Line 1633 | Line 1542 | transformations:
1542   L(1) = \frac{1}{p} \\
1543   \end{array}
1544   \]
1545 < , we obtain
1545 > we obtain
1546   \begin{eqnarray*}
1547   m\ddot x & =  & - \frac{{\partial W(x)}}{{\partial x}} -
1548   \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
# Line 1642 | Line 1551 | x_\alpha (0) - \frac{{g_\alpha  }}{{m_\alpha  \omega _
1551   & & + \sum\limits_{\alpha  = 1}^N {\left\{ {\left[ {g_\alpha
1552   x_\alpha (0) - \frac{{g_\alpha  }}{{m_\alpha  \omega _\alpha  }}}
1553   \right]\cos (\omega _\alpha  t) + \frac{{g_\alpha  \dot x_\alpha
1554 < (0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)} \right\}}
1555 < \end{eqnarray*}
1556 < \begin{eqnarray*}
1557 < m\ddot x & = & - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
1558 < {\sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1559 < }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha
1554 > (0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)} \right\}}\\
1555 > %
1556 > & = & -
1557 > \frac{{\partial W(x)}}{{\partial x}} - \int_0^t {\sum\limits_{\alpha
1558 > = 1}^N {\left( { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha
1559 > ^2 }}} \right)\cos (\omega _\alpha
1560   t)\dot x(t - \tau )d} \tau }  \\
1561   & & + \sum\limits_{\alpha  = 1}^N {\left\{ {\left[ {g_\alpha
1562   x_\alpha (0) - \frac{{g_\alpha }}{{m_\alpha \omega _\alpha  }}}
# Line 1681 | Line 1590 | harmonic bath. However, it is easy to verify that $R(t
1590   One may notice that $R(t)$ depends only on initial conditions, which
1591   implies it is completely deterministic within the context of a
1592   harmonic bath. However, it is easy to verify that $R(t)$ is totally
1593 < uncorrelated to $x$ and $\dot x$,
1594 < \[
1595 < \begin{array}{l}
1596 < \left\langle {x(t)R(t)} \right\rangle  = 0, \\
1688 < \left\langle {\dot x(t)R(t)} \right\rangle  = 0. \\
1689 < \end{array}
1690 < \]
1691 < This property is what we expect from a truly random process. As long
1692 < as the model chosen for $R(t)$ was a gaussian distribution in
1593 > uncorrelated to $x$ and $\dot x$,$\left\langle {x(t)R(t)}
1594 > \right\rangle  = 0, \left\langle {\dot x(t)R(t)} \right\rangle  =
1595 > 0.$ This property is what we expect from a truly random process. As
1596 > long as the model chosen for $R(t)$ was a gaussian distribution in
1597   general, the stochastic nature of the GLE still remains.
1694 –
1598   %dynamic friction kernel
1599   The convolution integral
1600   \[
# Line 1706 | Line 1609 | $\xi(t) = \Xi_0$. Hence, the convolution integral beco
1609   \[
1610   \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
1611   \]
1612 < and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1612 > and Eq.~\ref{introEuqation:GeneralizedLangevinDynamics} becomes
1613   \[
1614   m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
1615   \frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
# Line 1723 | Line 1626 | Hence, the convolution integral becomes
1626   \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
1627   {\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
1628   \]
1629 < and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1629 > and Eq.~\ref{introEuqation:GeneralizedLangevinDynamics} becomes
1630   \begin{equation}
1631   m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \xi _0 \dot
1632   x(t) + R(t) \label{introEquation:LangevinEquation}
# Line 1731 | Line 1634 | or be determined by Stokes' law for regular shaped par
1634   which is known as the Langevin equation. The static friction
1635   coefficient $\xi _0$ can either be calculated from spectral density
1636   or be determined by Stokes' law for regular shaped particles. A
1637 < briefly review on calculating friction tensor for arbitrary shaped
1637 > brief review on calculating friction tensors for arbitrary shaped
1638   particles is given in Sec.~\ref{introSection:frictionTensor}.
1639  
1640   \subsubsection{\label{introSection:secondFluctuationDissipation}\textbf{The Second Fluctuation Dissipation Theorem}}
1641  
1642 < Defining a new set of coordinates,
1642 > Defining a new set of coordinates
1643   \[
1644   q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
1645 < ^2 }}x(0)
1646 < \],
1645 > ^2 }}x(0),
1646 > \]
1647   we can rewrite $R(T)$ as
1648   \[
1649   R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
1650   \]
1651   And since the $q$ coordinates are harmonic oscillators,
1749 –
1652   \begin{eqnarray*}
1653   \left\langle {q_\alpha ^2 } \right\rangle  & = & \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
1654   \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle & = & \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
1655   \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle & = &\delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
1656   \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 } }  \\
1657    & = &\sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
1658 <  & = &kT\xi (t) \\
1658 >  & = &kT\xi (t)
1659   \end{eqnarray*}
1758 –
1660   Thus, we recover the \emph{second fluctuation dissipation theorem}
1661   \begin{equation}
1662   \xi (t) = \left\langle {R(t)R(0)} \right\rangle
1663 < \label{introEquation:secondFluctuationDissipation}.
1663 > \label{introEquation:secondFluctuationDissipation},
1664   \end{equation}
1665 < In effect, it acts as a constraint on the possible ways in which one
1666 < can model the random force and friction kernel.
1665 > which acts as a constraint on the possible ways in which one can
1666 > model the random force and friction kernel.

Diff Legend

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