--- trunk/tengDissertation/Introduction.tex	2006/06/26 13:42:53	2889
+++ trunk/tengDissertation/Introduction.tex	2006/07/18 05:47:33	2944
@@ -3,16 +3,17 @@ Closely related to Classical Mechanics, Molecular Dyna
 \section{\label{introSection:classicalMechanics}Classical
 Mechanics}
 
-Closely related to Classical Mechanics, Molecular Dynamics
-simulations are carried out by integrating the equations of motion
-for a given system of particles. There are three fundamental ideas
-behind classical mechanics. Firstly, one can determine the state of
-a mechanical system at any time of interest; Secondly, all the
-mechanical properties of the system at that time can be determined
-by combining the knowledge of the properties of the system with the
-specification of this state; Finally, the specification of the state
-when further combine with the laws of mechanics will also be
-sufficient to predict the future behavior of the system.
+Using equations of motion derived from Classical Mechanics,
+Molecular Dynamics simulations are carried out by integrating the
+equations of motion for a given system of particles. There are three
+fundamental ideas behind classical mechanics. Firstly, one can
+determine the state of a mechanical system at any time of interest;
+Secondly, all the mechanical properties of the system at that time
+can be determined by combining the knowledge of the properties of
+the system with the specification of this state; Finally, the
+specification of the state when further combined with the laws of
+mechanics will also be sufficient to predict the future behavior of
+the system.
 
 \subsection{\label{introSection:newtonian}Newtonian Mechanics}
 The discovery of Newton's three laws of mechanics which govern the
@@ -31,10 +32,9 @@ F_{ij} = -F_{ji}
 $F_{ji}$ be the force that particle $j$ exerts on particle $i$.
 Newton's third law states that
 \begin{equation}
-F_{ij} = -F_{ji}
+F_{ij} = -F_{ji}.
  \label{introEquation:newtonThirdLaw}
 \end{equation}
-
 Conservation laws of Newtonian Mechanics play very important roles
 in solving mechanics problems. The linear momentum of a particle is
 conserved if it is free or it experiences no force. The second
@@ -63,50 +63,48 @@ that if all forces are conservative, Energy
 \end{equation}
 If there are no external torques acting on a body, the angular
 momentum of it is conserved. The last conservation theorem state
-that if all forces are conservative, Energy
-\begin{equation}E = T + V \label{introEquation:energyConservation}
+that if all forces are conservative, energy is conserved,
+\begin{equation}E = T + V. \label{introEquation:energyConservation}
 \end{equation}
- is conserved. All of these conserved quantities are
-important factors to determine the quality of numerical integration
-schemes for rigid bodies \cite{Dullweber1997}.
+All of these conserved quantities are important factors to determine
+the quality of numerical integration schemes for rigid
+bodies.\cite{Dullweber1997}
 
 \subsection{\label{introSection:lagrangian}Lagrangian Mechanics}
 
-Newtonian Mechanics suffers from two important limitations: motions
-can only be described in cartesian coordinate systems. Moreover, It
-become impossible to predict analytically the properties of the
-system even if we know all of the details of the interaction. In
-order to overcome some of the practical difficulties which arise in
-attempts to apply Newton's equation to complex system, approximate
-numerical procedures may be developed.
+Newtonian Mechanics suffers from an important limitation: motion can
+only be described in cartesian coordinate systems which make it
+impossible to predict analytically the properties of the system even
+if we know all of the details of the interaction. In order to
+overcome some of the practical difficulties which arise in attempts
+to apply Newton's equation to complex systems, approximate numerical
+procedures may be developed.
 
 \subsubsection{\label{introSection:halmiltonPrinciple}\textbf{Hamilton's
 Principle}}
 
 Hamilton introduced the dynamical principle upon which it is
 possible to base all of mechanics and most of classical physics.
-Hamilton's Principle may be stated as follows,
-
-The actual trajectory, along which a dynamical system may move from
-one point to another within a specified time, is derived by finding
-the path which minimizes the time integral of the difference between
-the kinetic, $K$, and potential energies, $U$.
+Hamilton's Principle may be stated as follows: the trajectory, along
+which a dynamical system may move from one point to another within a
+specified time, is derived by finding the path which minimizes the
+time integral of the difference between the kinetic $K$, and
+potential energies $U$,
 \begin{equation}
-\delta \int_{t_1 }^{t_2 } {(K - U)dt = 0} ,
+\delta \int_{t_1 }^{t_2 } {(K - U)dt = 0}.
 \label{introEquation:halmitonianPrinciple1}
 \end{equation}
-
 For simple mechanical systems, where the forces acting on the
 different parts are derivable from a potential, the Lagrangian
 function $L$ can be defined as the difference between the kinetic
 energy of the system and its potential energy,
 \begin{equation}
-L \equiv K - U = L(q_i ,\dot q_i ) ,
+L \equiv K - U = L(q_i ,\dot q_i ).
 \label{introEquation:lagrangianDef}
 \end{equation}
-then Eq.~\ref{introEquation:halmitonianPrinciple1} becomes
+Thus, Eq.~\ref{introEquation:halmitonianPrinciple1} becomes
 \begin{equation}
-\delta \int_{t_1 }^{t_2 } {L dt = 0} ,
+\delta \int_{t_1 }^{t_2 } {L dt = 0} .
 \label{introEquation:halmitonianPrinciple2}
 \end{equation}
 
@@ -138,7 +136,6 @@ p_i  = \frac{{\partial L}}{{\partial q_i }}
 p_i  = \frac{{\partial L}}{{\partial q_i }}
 \label{introEquation:generalizedMomentaDot}
 \end{equation}
-
 With the help of the generalized momenta, we may now define a new
 quantity $H$ by the equation
 \begin{equation}
@@ -146,22 +143,20 @@ $L$ is the Lagrangian function for the system.
 \label{introEquation:hamiltonianDefByLagrangian}
 \end{equation}
 where $ \dot q_1  \ldots \dot q_f $ are generalized velocities and
-$L$ is the Lagrangian function for the system.
-
-Differentiating Eq.~\ref{introEquation:hamiltonianDefByLagrangian},
-one can obtain
+$L$ is the Lagrangian function for the system. Differentiating
+Eq.~\ref{introEquation:hamiltonianDefByLagrangian}, one can obtain
 \begin{equation}
 dH = \sum\limits_k {\left( {p_k d\dot q_k  + \dot q_k dp_k  -
 \frac{{\partial L}}{{\partial q_k }}dq_k  - \frac{{\partial
 L}}{{\partial \dot q_k }}d\dot q_k } \right)}  - \frac{{\partial
-L}}{{\partial t}}dt \label{introEquation:diffHamiltonian1}
+L}}{{\partial t}}dt . \label{introEquation:diffHamiltonian1}
 \end{equation}
-Making use of  Eq.~\ref{introEquation:generalizedMomenta}, the
-second and fourth terms in the parentheses cancel. Therefore,
+Making use of Eq.~\ref{introEquation:generalizedMomenta}, the second
+and fourth terms in the parentheses cancel. Therefore,
 Eq.~\ref{introEquation:diffHamiltonian1} can be rewritten as
 \begin{equation}
 dH = \sum\limits_k {\left( {\dot q_k dp_k  - \dot p_k dq_k }
-\right)}  - \frac{{\partial L}}{{\partial t}}dt
+\right)}  - \frac{{\partial L}}{{\partial t}}dt .
 \label{introEquation:diffHamiltonian2}
 \end{equation}
 By identifying the coefficients of $dq_k$, $dp_k$ and dt, we can
@@ -180,11 +175,10 @@ t}}
 t}}
 \label{introEquation:motionHamiltonianTime}
 \end{equation}
-
-Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
+where Eq.~\ref{introEquation:motionHamiltonianCoordinate} and
 Eq.~\ref{introEquation:motionHamiltonianMomentum} are Hamilton's
 equation of motion. Due to their symmetrical formula, they are also
-known as the canonical equations of motions \cite{Goldstein2001}.
+known as the canonical equations of motions.\cite{Goldstein2001}
 
 An important difference between Lagrangian approach and the
 Hamiltonian approach is that the Lagrangian is considered to be a
@@ -194,19 +188,18 @@ only works with 1st-order differential equations\cite{
 Hamiltonian Mechanics is more appropriate for application to
 statistical mechanics and quantum mechanics, since it treats the
 coordinate and its time derivative as independent variables and it
-only works with 1st-order differential equations\cite{Marion1990}.
-
+only works with 1st-order differential equations.\cite{Marion1990}
 In Newtonian Mechanics, a system described by conservative forces
-conserves the total energy \ref{introEquation:energyConservation}.
-It follows that Hamilton's equations of motion conserve the total
-Hamiltonian.
+conserves the total energy
+(Eq.~\ref{introEquation:energyConservation}). It follows that
+Hamilton's equations of motion conserve the total Hamiltonian
 \begin{equation}
 \frac{{dH}}{{dt}} = \sum\limits_i {\left( {\frac{{\partial
 H}}{{\partial q_i }}\dot q_i  + \frac{{\partial H}}{{\partial p_i
 }}\dot p_i } \right)}  = \sum\limits_i {\left( {\frac{{\partial
 H}}{{\partial q_i }}\frac{{\partial H}}{{\partial p_i }} -
 \frac{{\partial H}}{{\partial p_i }}\frac{{\partial H}}{{\partial
-q_i }}} \right) = 0} \label{introEquation:conserveHalmitonian}
+q_i }}} \right) = 0}. \label{introEquation:conserveHalmitonian}
 \end{equation}
 
 \section{\label{introSection:statisticalMechanics}Statistical
@@ -215,22 +208,29 @@ Statistical Mechanics concepts and theorem presented i
 The thermodynamic behaviors and properties of Molecular Dynamics
 simulation are governed by the principle of Statistical Mechanics.
 The following section will give a brief introduction to some of the
-Statistical Mechanics concepts and theorem presented in this
+Statistical Mechanics concepts and theorems presented in this
 dissertation.
 
 \subsection{\label{introSection:ensemble}Phase Space and Ensemble}
 
 Mathematically, phase space is the space which represents all
-possible states. Each possible state of the system corresponds to
-one unique point in the phase space. For mechanical systems, the
-phase space usually consists of all possible values of position and
-momentum variables. Consider a dynamic system of $f$ particles in a
-cartesian space, where each of the $6f$ coordinates and momenta is
-assigned to one of $6f$ mutually orthogonal axes, the phase space of
-this system is a $6f$ dimensional space. A point, $x = (\rightarrow
-q_1 , \ldots ,\rightarrow q_f ,\rightarrow p_1 , \ldots ,\rightarrow
-p_f )$, with a unique set of values of $6f$ coordinates and momenta
-is a phase space vector.
+possible states of a system. Each possible state of the system
+corresponds to one unique point in the phase space. For mechanical
+systems, the phase space usually consists of all possible values of
+position and momentum variables. Consider a dynamic system of $f$
+particles in a cartesian space, where each of the $6f$ coordinates
+and momenta is assigned to one of $6f$ mutually orthogonal axes, the
+phase space of this system is a $6f$ dimensional space. A point, $x
+=
+(\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
+\over q} _1 , \ldots
+,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
+\over q} _f
+,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
+\over p} _1  \ldots
+,\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}}
+\over p} _f )$ , with a unique set of values of $6f$ coordinates and
+momenta is a phase space vector.
 %%%fix me
 
 In statistical mechanics, the condition of an ensemble at any time
@@ -243,16 +243,15 @@ their locations which would change the density at any 
 \label{introEquation:densityDistribution}
 \end{equation}
 Governed by the principles of mechanics, the phase points change
-their locations which would change the density at any time at phase
+their locations which changes the density at any time at phase
 space. Hence, the density distribution is also to be taken as a
-function of the time.
-
-The number of systems $\delta N$ at time $t$ can be determined by,
+function of the time. The number of systems $\delta N$ at time $t$
+can be determined by,
 \begin{equation}
 \delta N = \rho (q,p,t)dq_1  \ldots dq_f dp_1  \ldots dp_f.
 \label{introEquation:deltaN}
 \end{equation}
-Assuming a large enough population of systems, we can sufficiently
+Assuming enough copies of the systems, we can sufficiently
 approximate $\delta N$ without introducing discontinuity when we go
 from one region in the phase space to another. By integrating over
 the whole phase space,
@@ -260,8 +259,8 @@ gives us an expression for the total number of the sys
 N = \int { \ldots \int {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f
 \label{introEquation:totalNumberSystem}
 \end{equation}
-gives us an expression for the total number of the systems. Hence,
-the probability per unit in the phase space can be obtained by,
+gives us an expression for the total number of copies. Hence, the
+probability per unit volume in the phase space can be obtained by,
 \begin{equation}
 \frac{{\rho (q,p,t)}}{N} = \frac{{\rho (q,p,t)}}{{\int { \ldots \int
 {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
@@ -270,53 +269,20 @@ momenta of the system. Even when the dynamics of the r
 With the help of Eq.~\ref{introEquation:unitProbability} and the
 knowledge of the system, it is possible to calculate the average
 value of any desired quantity which depends on the coordinates and
-momenta of the system. Even when the dynamics of the real system is
+momenta of the system. Even when the dynamics of the real system are
 complex, or stochastic, or even discontinuous, the average
-properties of the ensemble of possibilities as a whole remaining
-well defined. For a classical system in thermal equilibrium with its
+properties of the ensemble of possibilities as a whole remain well
+defined. For a classical system in thermal equilibrium with its
 environment, the ensemble average of a mechanical quantity, $\langle
 A(q , p) \rangle_t$, takes the form of an integral over the phase
 space of the system,
 \begin{equation}
 \langle  A(q , p) \rangle_t = \frac{{\int { \ldots \int {A(q,p)\rho
 (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}{{\int { \ldots \int {\rho
-(q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}
+(q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
 \label{introEquation:ensembelAverage}
 \end{equation}
 
-There are several different types of ensembles with different
-statistical characteristics. As a function of macroscopic
-parameters, such as temperature \textit{etc}, the partition function
-can be used to describe the statistical properties of a system in
-thermodynamic equilibrium.
-
-As an ensemble of systems, each of which is known to be thermally
-isolated and conserve energy, the Microcanonical ensemble (NVE) has
-a partition function like,
-\begin{equation}
-\Omega (N,V,E) = e^{\beta TS} \label{introEquation:NVEPartition}.
-\end{equation}
-A canonical ensemble (NVT)is an ensemble of systems, each of which
-can share its energy with a large heat reservoir. The distribution
-of the total energy amongst the possible dynamical states is given
-by the partition function,
-\begin{equation}
-\Omega (N,V,T) = e^{ - \beta A}
-\label{introEquation:NVTPartition}
-\end{equation}
-Here, $A$ is the Helmholtz free energy which is defined as $ A = U -
-TS$. Since most experiments are carried out under constant pressure
-condition, the isothermal-isobaric ensemble (NPT) plays a very
-important role in molecular simulations. The isothermal-isobaric
-ensemble allow the system to exchange energy with a heat bath of
-temperature $T$ and to change the volume as well. Its partition
-function is given as
-\begin{equation}
-\Delta (N,P,T) =  - e^{\beta G}.
- \label{introEquation:NPTPartition}
-\end{equation}
-Here, $G = U - TS + PV$, and $G$ is called Gibbs free energy.
-
 \subsection{\label{introSection:liouville}Liouville's theorem}
 
 Liouville's theorem is the foundation on which statistical mechanics
@@ -358,12 +324,11 @@ simple form,
 \frac{{\partial \rho }}{{\partial p_i }}\dot p_i } \right)}  = 0 .
 \label{introEquation:liouvilleTheorem}
 \end{equation}
-
 Liouville's theorem states that the distribution function is
 constant along any trajectory in phase space. In classical
-statistical mechanics, since the number of members in an ensemble is
-huge and constant, we can assume the local density has no reason
-(other than classical mechanics) to change,
+statistical mechanics, since the number of system copies in an
+ensemble is huge and constant, we can assume the local density has
+no reason (other than classical mechanics) to change,
 \begin{equation}
 \frac{{\partial \rho }}{{\partial t}} = 0.
 \label{introEquation:stationary}
@@ -393,8 +358,8 @@ With the help of stationary assumption
 \frac{{d(\delta N)}}{{dt}} = \frac{{d\rho }}{{dt}}\delta v + \rho
 \frac{d}{{dt}}(\delta v) = 0.
 \end{equation}
-With the help of stationary assumption
-(\ref{introEquation:stationary}), we obtain the principle of the
+With the help of the stationary assumption
+(Eq.~\ref{introEquation:stationary}), we obtain the principle of
 \emph{conservation of volume in phase space},
 \begin{equation}
 \frac{d}{{dt}}(\delta v) = \frac{d}{{dt}}\int { \ldots \int {dq_1 }
@@ -404,10 +369,10 @@ Liouville's theorem can be expresses in a variety of d
 
 \subsubsection{\label{introSection:liouvilleInOtherForms}\textbf{Liouville's Theorem in Other Forms}}
 
-Liouville's theorem can be expresses in a variety of different forms
+Liouville's theorem can be expressed in a variety of different forms
 which are convenient within different contexts. For any two function
 $F$ and $G$ of the coordinates and momenta of a system, the Poisson
-bracket ${F, G}$ is defined as
+bracket $\{F,G\}$ is defined as
 \begin{equation}
 \left\{ {F,G} \right\} = \sum\limits_i {\left( {\frac{{\partial
 F}}{{\partial q_i }}\frac{{\partial G}}{{\partial p_i }} -
@@ -415,9 +380,9 @@ Substituting equations of motion in Hamiltonian formal
 q_i }}} \right)}.
 \label{introEquation:poissonBracket}
 \end{equation}
-Substituting equations of motion in Hamiltonian formalism(
-Eq.~\ref{introEquation:motionHamiltonianCoordinate} ,
-Eq.~\ref{introEquation:motionHamiltonianMomentum} ) into
+Substituting equations of motion in Hamiltonian formalism
+(Eq.~\ref{introEquation:motionHamiltonianCoordinate} ,
+Eq.~\ref{introEquation:motionHamiltonianMomentum}) into
 (Eq.~\ref{introEquation:liouvilleTheorem}), we can rewrite
 Liouville's theorem using Poisson bracket notion,
 \begin{equation}
@@ -438,7 +403,7 @@ expressed as
 \left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - iL\rho
 \label{introEquation:liouvilleTheoremInOperator}
 \end{equation}
-
+which can help define a propagator $\rho (t) = e^{-iLt} \rho (0)$.
 \subsection{\label{introSection:ergodic}The Ergodic Hypothesis}
 
 Various thermodynamic properties can be calculated from Molecular
@@ -447,11 +412,11 @@ period of them which is different from the average beh
 simulation and the quality of the underlying model. However, both
 experiments and computer simulations are usually performed during a
 certain time interval and the measurements are averaged over a
-period of them which is different from the average behavior of
+period of time which is different from the average behavior of
 many-body system in Statistical Mechanics. Fortunately, the Ergodic
 Hypothesis makes a connection between time average and the ensemble
 average. It states that the time average and average over the
-statistical ensemble are identical \cite{Frenkel1996, Leach2001}.
+statistical ensemble are identical:\cite{Frenkel1996, Leach2001}
 \begin{equation}
 \langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty }
 \frac{1}{t}\int\limits_0^t {A(q(t),p(t))dt = \int\limits_\Gamma
@@ -460,68 +425,74 @@ sufficiently long time (longer than relaxation time), 
 where $\langle  A(q , p) \rangle_t$ is an equilibrium value of a
 physical quantity and $\rho (p(t), q(t))$ is the equilibrium
 distribution function. If an observation is averaged over a
-sufficiently long time (longer than relaxation time), all accessible
-microstates in phase space are assumed to be equally probed, giving
-a properly weighted statistical average. This allows the researcher
-freedom of choice when deciding how best to measure a given
-observable. In case an ensemble averaged approach sounds most
-reasonable, the Monte Carlo techniques\cite{Metropolis1949} can be
+sufficiently long time (longer than the relaxation time), all
+accessible microstates in phase space are assumed to be equally
+probed, giving a properly weighted statistical average. This allows
+the researcher freedom of choice when deciding how best to measure a
+given observable. In case an ensemble averaged approach sounds most
+reasonable, the Monte Carlo methods\cite{Metropolis1949} can be
 utilized. Or if the system lends itself to a time averaging
 approach, the Molecular Dynamics techniques in
 Sec.~\ref{introSection:molecularDynamics} will be the best
-choice\cite{Frenkel1996}.
+choice.\cite{Frenkel1996}
 
 \section{\label{introSection:geometricIntegratos}Geometric Integrators}
 A variety of numerical integrators have been proposed to simulate
 the motions of atoms in MD simulation. They usually begin with
-initial conditionals and move the objects in the direction governed
-by the differential equations. However, most of them ignore the
-hidden physical laws contained within the equations. Since 1990,
-geometric integrators, which preserve various phase-flow invariants
-such as symplectic structure, volume and time reversal symmetry, are
-developed to address this issue\cite{Dullweber1997, McLachlan1998,
-Leimkuhler1999}. The velocity Verlet method, which happens to be a
+initial conditions and move the objects in the direction governed by
+the differential equations. However, most of them ignore the hidden
+physical laws contained within the equations. Since 1990, geometric
+integrators, which preserve various phase-flow invariants such as
+symplectic structure, volume and time reversal symmetry, were
+developed to address this issue.\cite{Dullweber1997, McLachlan1998,
+Leimkuhler1999} The velocity Verlet method, which happens to be a
 simple example of symplectic integrator, continues to gain
 popularity in the molecular dynamics community. This fact can be
 partly explained by its geometric nature.
 
-\subsection{\label{introSection:symplecticManifold}Symplectic Manifolds}
+\subsection{\label{introSection:symplecticManifold}Manifolds and Bundles}
 A \emph{manifold} is an abstract mathematical space. It looks
 locally like Euclidean space, but when viewed globally, it may have
 more complicated structure. A good example of manifold is the
 surface of Earth. It seems to be flat locally, but it is round if
 viewed as a whole. A \emph{differentiable manifold} (also known as
 \emph{smooth manifold}) is a manifold on which it is possible to
-apply calculus on \emph{differentiable manifold}. A \emph{symplectic
-manifold} is defined as a pair $(M, \omega)$ which consists of a
-\emph{differentiable manifold} $M$ and a close, non-degenerated,
+apply calculus.\cite{Hirsch1997} A \emph{symplectic manifold} is
+defined as a pair $(M, \omega)$ which consists of a
+\emph{differentiable manifold} $M$ and a close, non-degenerate,
 bilinear symplectic form, $\omega$. A symplectic form on a vector
 space $V$ is a function $\omega(x, y)$ which satisfies
 $\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+
 \lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and
-$\omega(x, x) = 0$. The cross product operation in vector field is
-an example of symplectic form.
+$\omega(x, x) = 0$.\cite{McDuff1998} The cross product operation in
+vector field is an example of symplectic form.
+Given vector spaces $V$ and $W$ over same field $F$, $f: V \to W$ is a linear transformation if 
+\begin{eqnarray*}
+f(x+y) & = & f(x) + f(y) \\
+f(ax) & = & af(x)	
+\end{eqnarray*}
+are always satisfied for any two vectors $x$ and $y$ in $V$ and any scalar $a$ in $F$. One can define the dual vector space $V^*$ of $V$ if any two built-in linear transformations $\phi$ and $\psi$ in $V^*$ satisfy the following definition of addition and scalar multiplication:
+\begin{eqnarray*}
+(\phi+\psi)(x) & = & \phi(x)+\psi(x) \\
+(a\phi)(x) & = & a \phi(x)
+\end{eqnarray*}
+for all $a$ in $F$ and $x$ in $V$. For a manifold $M$, one can define a tangent vector of a tangent space $TM_q$ at every point $q$
+\begin{equation}
+\dot q = \mathop {\lim }\limits_{t \to 0} \frac{{\phi (t) - \phi (0)}}{t} 	
+\end{equation}
+where $\phi(0)=q$ and $\phi(t) \in M$. One may also define a cotangent space $T^*M_q$ as the dual space of the tangent space $TM_q$. The tangent space and the cotangent space are isomorphic to each other, since they are both real vector spaces with same dimension.
+The union of tangent spaces at every point of $M$ is called the tangent bundle of $M$ and is denoted by $TM$, while cotangent bundle $T^*M$ is defined as the union of the cotangent spaces to $M$.\cite{Jost2002} For a Hamiltonian system with configuration manifold $V$, the $(q,\dot q)$ phase space is the tangent bundle of the configuration manifold $V$, while the cotangent bundle is represented by $(q,p)$. 
 
-One of the motivations to study \emph{symplectic manifolds} in
-Hamiltonian Mechanics is that a symplectic manifold can represent
-all possible configurations of the system and the phase space of the
-system can be described by it's cotangent bundle. Every symplectic
-manifold is even dimensional. For instance, in Hamilton equations,
-coordinate and momentum always appear in pairs.
-
 \subsection{\label{introSection:ODE}Ordinary Differential Equations}
 
 For an ordinary differential system defined as
 \begin{equation}
 \dot x = f(x)
 \end{equation}
-where $x = x(q,p)^T$, this system is a canonical Hamiltonian, if
+where $x = x(q,p)$, this system is a canonical Hamiltonian, if
+$f(x) = J\nabla _x H(x)$. Here, $H = H (q, p)$ is Hamiltonian
+function and $J$ is the skew-symmetric matrix
 \begin{equation}
-f(r) = J\nabla _x H(r).
-\end{equation}
-$H = H (q, p)$ is Hamiltonian function and $J$ is the skew-symmetric
-matrix
-\begin{equation}
 J = \left( {\begin{array}{*{20}c}
    0 & I  \\
    { - I} & 0  \\
@@ -531,30 +502,30 @@ system can be rewritten as,
 where $I$ is an identity matrix. Using this notation, Hamiltonian
 system can be rewritten as,
 \begin{equation}
-\frac{d}{{dt}}x = J\nabla _x H(x)
+\frac{d}{{dt}}x = J\nabla _x H(x).
 \label{introEquation:compactHamiltonian}
 \end{equation}In this case, $f$ is
-called a \emph{Hamiltonian vector field}.
-
-Another generalization of Hamiltonian dynamics is Poisson
-Dynamics\cite{Olver1986},
+called a \emph{Hamiltonian vector field}. Another generalization of
+Hamiltonian dynamics is Poisson Dynamics,\cite{Olver1986}
 \begin{equation}
 \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
 \end{equation}
-The most obvious change being that matrix $J$ now depends on $x$.
+where the most obvious change being that matrix $J$ now depends on
+$x$.
 
-\subsection{\label{introSection:exactFlow}Exact Flow}
+\subsection{\label{introSection:exactFlow}Exact Propagator}
 
-Let $x(t)$ be the exact solution of the ODE system,
+Let $x(t)$ be the exact solution of the ODE
+system,
 \begin{equation}
-\frac{{dx}}{{dt}} = f(x) \label{introEquation:ODE}
-\end{equation}
-The exact flow(solution) $\varphi_\tau$ is defined by
-\[
-x(t+\tau) =\varphi_\tau(x(t))
+\frac{{dx}}{{dt}} = f(x), \label{introEquation:ODE}
+\end{equation} we can
+define its exact propagator $\varphi_\tau$:
+\[ x(t+\tau)
+=\varphi_\tau(x(t))
 \]
 where $\tau$ is a fixed time step and $\varphi$ is a map from phase
-space to itself. The flow has the continuous group property,
+space to itself. The propagator has the continuous group property,
 \begin{equation}
 \varphi _{\tau _1 }  \circ \varphi _{\tau _2 }  = \varphi _{\tau _1
 + \tau _2 } .
@@ -563,101 +534,84 @@ Therefore, the exact flow is self-adjoint,
 \begin{equation}
 \varphi _\tau   \circ \varphi _{ - \tau }  = I
 \end{equation}
-Therefore, the exact flow is self-adjoint,
+Therefore, the exact propagator is self-adjoint,
 \begin{equation}
 \varphi _\tau   = \varphi _{ - \tau }^{ - 1}.
 \end{equation}
-The exact flow can also be written in terms of the of an operator,
+In most cases, it is not easy to find the exact propagator
+$\varphi_\tau$. Instead, we use an approximate map, $\psi_\tau$,
+which is usually called an integrator. The order of an integrator
+$\psi_\tau$ is $p$, if the Taylor series of $\psi_\tau$ agree to
+order $p$,
 \begin{equation}
-\varphi _\tau  (x) = e^{\tau \sum\limits_i {f_i (x)\frac{\partial
-}{{\partial x_i }}} } (x) \equiv \exp (\tau f)(x).
-\label{introEquation:exponentialOperator}
-\end{equation}
-
-In most cases, it is not easy to find the exact flow $\varphi_\tau$.
-Instead, we use an approximate map, $\psi_\tau$, which is usually
-called integrator. The order of an integrator $\psi_\tau$ is $p$, if
-the Taylor series of $\psi_\tau$ agree to order $p$,
-\begin{equation}
 \psi_\tau(x) = x + \tau f(x) + O(\tau^{p+1})
 \end{equation}
 
 \subsection{\label{introSection:geometricProperties}Geometric Properties}
 
 The hidden geometric properties\cite{Budd1999, Marsden1998} of an
-ODE and its flow play important roles in numerical studies. Many of
-them can be found in systems which occur naturally in applications.
-
-Let $\varphi$ be the flow of Hamiltonian vector field, $\varphi$ is
-a \emph{symplectic} flow if it satisfies,
+ODE and its propagator play important roles in numerical studies.
+Many of them can be found in systems which occur naturally in
+applications. Let $\varphi$ be the propagator of Hamiltonian vector
+field, $\varphi$ is a \emph{symplectic} propagator if it satisfies,
 \begin{equation}
 {\varphi '}^T J \varphi ' = J.
 \end{equation}
 According to Liouville's theorem, the symplectic volume is invariant
-under a Hamiltonian flow, which is the basis for classical
-statistical mechanics. Furthermore, the flow of a Hamiltonian vector
-field on a symplectic manifold can be shown to be a
+under a Hamiltonian propagator, which is the basis for classical
+statistical mechanics. Furthermore, the propagator of a Hamiltonian
+vector field on a symplectic manifold can be shown to be a
 symplectomorphism. As to the Poisson system,
 \begin{equation}
 {\varphi '}^T J \varphi ' = J \circ \varphi
 \end{equation}
-is the property that must be preserved by the integrator.
-
-It is possible to construct a \emph{volume-preserving} flow for a
-source free ODE ($ \nabla \cdot f = 0 $), if the flow satisfies $
-\det d\varphi  = 1$. One can show easily that a symplectic flow will
-be volume-preserving.
-
-Changing the variables $y = h(x)$ in an ODE
-(Eq.~\ref{introEquation:ODE}) will result in a new system,
+is the property that must be preserved by the integrator. It is
+possible to construct a \emph{volume-preserving} propagator for a
+source free ODE ($ \nabla \cdot f = 0 $), if the propagator
+satisfies $ \det d\varphi  = 1$. One can show easily that a
+symplectic propagator will be volume-preserving. Changing the
+variables $y = h(x)$ in an ODE (Eq.~\ref{introEquation:ODE}) will
+result in a new system,
 \[
 \dot y = \tilde f(y) = ((dh \cdot f)h^{ - 1} )(y).
 \]
 The vector filed $f$ has reversing symmetry $h$ if $f = - \tilde f$.
-In other words, the flow of this vector field is reversible if and
-only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $.
-
-A \emph{first integral}, or conserved quantity of a general
-differential function is a function $ G:R^{2d}  \to R^d $ which is
-constant for all solutions of the ODE $\frac{{dx}}{{dt}} = f(x)$ ,
+In other words, the propagator of this vector field is reversible if
+and only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $. A
+conserved quantity of a general differential function is a function
+$ G:R^{2d}  \to R^d $ which is constant for all solutions of the ODE
+$\frac{{dx}}{{dt}} = f(x)$ ,
 \[
 \frac{{dG(x(t))}}{{dt}} = 0.
 \]
-Using chain rule, one may obtain,
+Using the chain rule, one may obtain,
 \[
-\sum\limits_i {\frac{{dG}}{{dx_i }}} f_i (x) = f \bullet \nabla G,
+\sum\limits_i {\frac{{dG}}{{dx_i }}} f_i (x) = f \cdot \nabla G,
 \]
-which is the condition for conserving \emph{first integral}. For a
-canonical Hamiltonian system, the time evolution of an arbitrary
-smooth function $G$ is given by,
-
+which is the condition for conserved quantities. For a canonical
+Hamiltonian system, the time evolution of an arbitrary smooth
+function $G$ is given by,
 \begin{eqnarray}
-\frac{{dG(x(t))}}{{dt}} & = & [\nabla _x G(x(t))]^T \dot x(t) \\
-                        & = & [\nabla _x G(x(t))]^T J\nabla _x H(x(t)). \\
+\frac{{dG(x(t))}}{{dt}} & = & [\nabla _x G(x(t))]^T \dot x(t) \notag\\
+                        & = & [\nabla _x G(x(t))]^T J\nabla _x H(x(t)).
 \label{introEquation:firstIntegral1}
 \end{eqnarray}
-
-
-Using poisson bracket notion, Equation
-\ref{introEquation:firstIntegral1} can be rewritten as
+Using poisson bracket notion, Eq.~\ref{introEquation:firstIntegral1}
+can be rewritten as
 \[
 \frac{d}{{dt}}G(x(t)) = \left\{ {G,H} \right\}(x(t)).
 \]
-Therefore, the sufficient condition for $G$ to be the \emph{first
-integral} of a Hamiltonian system is
-\[
-\left\{ {G,H} \right\} = 0.
-\]
-As well known, the Hamiltonian (or energy) H of a Hamiltonian system
-is a \emph{first integral}, which is due to the fact $\{ H,H\}  =
-0$.
-
+Therefore, the sufficient condition for $G$ to be a conserved
+quantity of a Hamiltonian system is $\left\{ {G,H} \right\} = 0.$ As
+is well known, the Hamiltonian (or energy) H of a Hamiltonian system
+is a conserved quantity, which is due to the fact $\{ H,H\}  = 0$.
 When designing any numerical methods, one should always try to
-preserve the structural properties of the original ODE and its flow.
+preserve the structural properties of the original ODE and its
+propagator.
 
 \subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
 A lot of well established and very effective numerical methods have
-been successful precisely because of their symplecticities even
+been successful precisely because of their symplectic nature even
 though this fact was not recognized when they were first
 constructed. The most famous example is the Verlet-leapfrog method
 in molecular dynamics. In general, symplectic integrators can be
@@ -668,43 +622,43 @@ constructed using one of four different methods.
 \item Runge-Kutta methods
 \item Splitting methods
 \end{enumerate}
-
-Generating function\cite{Channell1990} tends to lead to methods
+Generating functions\cite{Channell1990} tend to lead to methods
 which are cumbersome and difficult to use. In dissipative systems,
 variational methods can capture the decay of energy
-accurately\cite{Kane2000}. Since their geometrically unstable nature
+accurately.\cite{Kane2000} Since they are geometrically unstable
 against non-Hamiltonian perturbations, ordinary implicit Runge-Kutta
-methods are not suitable for Hamiltonian system. Recently, various
-high-order explicit Runge-Kutta methods
-\cite{Owren1992,Chen2003}have been developed to overcome this
-instability. However, due to computational penalty involved in
-implementing the Runge-Kutta methods, they have not attracted much
-attention from the Molecular Dynamics community. Instead, splitting
-methods have been widely accepted since they exploit natural
-decompositions of the system\cite{Tuckerman1992, McLachlan1998}.
+methods are not suitable for Hamiltonian
+system.\cite{Cartwright1992} Recently, various high-order explicit
+Runge-Kutta methods \cite{Owren1992,Chen2003} have been developed to
+overcome this instability. However, due to computational penalty
+involved in implementing the Runge-Kutta methods, they have not
+attracted much attention from the Molecular Dynamics community.
+Instead, splitting methods have been widely accepted since they
+exploit natural decompositions of the system.\cite{McLachlan1998,
+Tuckerman1992}
 
 \subsubsection{\label{introSection:splittingMethod}\textbf{Splitting Methods}}
 
 The main idea behind splitting methods is to decompose the discrete
-$\varphi_h$ as a composition of simpler flows,
+$\varphi_h$ as a composition of simpler propagators,
 \begin{equation}
 \varphi _h  = \varphi _{h_1 }  \circ \varphi _{h_2 }  \ldots  \circ
 \varphi _{h_n }
 \label{introEquation:FlowDecomposition}
 \end{equation}
-where each of the sub-flow is chosen such that each represent a
-simpler integration of the system.
-
-Suppose that a Hamiltonian system takes the form,
+where each of the sub-propagator is chosen such that each represent
+a simpler integration of the system. Suppose that a Hamiltonian
+system takes the form,
 \[
 H = H_1 + H_2.
 \]
 Here, $H_1$ and $H_2$ may represent different physical processes of
 the system. For instance, they may relate to kinetic and potential
 energy respectively, which is a natural decomposition of the
-problem. If $H_1$ and $H_2$ can be integrated using exact flows
-$\varphi_1(t)$ and $\varphi_2(t)$, respectively, a simple first
-order expression is then given by the Lie-Trotter formula
+problem. If $H_1$ and $H_2$ can be integrated using exact
+propagators $\varphi_1(t)$ and $\varphi_2(t)$, respectively, a
+simple first order expression is then given by the Lie-Trotter
+formula\cite{Trotter1959}
 \begin{equation}
 \varphi _h  = \varphi _{1,h}  \circ \varphi _{2,h},
 \label{introEquation:firstOrderSplitting}
@@ -713,8 +667,8 @@ It is easy to show that any composition of symplectic 
 continuous $\varphi _i$ over a time $h$. By definition, as
 $\varphi_i(t)$ is the exact solution of a Hamiltonian system, it
 must follow that each operator $\varphi_i(t)$ is a symplectic map.
-It is easy to show that any composition of symplectic flows yields a
-symplectic map,
+It is easy to show that any composition of symplectic propagators
+yields a symplectic map,
 \begin{equation}
 (\varphi '\phi ')^T J\varphi '\phi ' = \phi '^T \varphi '^T J\varphi
 '\phi ' = \phi '^T J\phi ' = J,
@@ -722,15 +676,15 @@ splitting in this context automatically generates a sy
 \end{equation}
 where $\phi$ and $\psi$ both are symplectic maps. Thus operator
 splitting in this context automatically generates a symplectic map.
-
-The Lie-Trotter splitting(\ref{introEquation:firstOrderSplitting})
-introduces local errors proportional to $h^2$, while Strang
-splitting gives a second-order decomposition,
+The Lie-Trotter
+splitting(Eq.~\ref{introEquation:firstOrderSplitting}) introduces
+local errors proportional to $h^2$, while the Strang splitting gives
+a second-order decomposition,\cite{Strang1968}
 \begin{equation}
 \varphi _h  = \varphi _{1,h/2}  \circ \varphi _{2,h}  \circ \varphi
 _{1,h/2} , \label{introEquation:secondOrderSplitting}
 \end{equation}
-which has a local error proportional to $h^3$. The Sprang
+which has a local error proportional to $h^3$. The Strang
 splitting's popularity in molecular simulation community attribute
 to its symmetric property,
 \begin{equation}
@@ -758,9 +712,9 @@ symplectic(\ref{introEquation:SymplecticFlowCompositio
 \end{align}
 where $F(t)$ is the force at time $t$. This integration scheme is
 known as \emph{velocity verlet} which is
-symplectic(\ref{introEquation:SymplecticFlowComposition}),
-time-reversible(\ref{introEquation:timeReversible}) and
-volume-preserving (\ref{introEquation:volumePreserving}). These
+symplectic(Eq.~\ref{introEquation:SymplecticFlowComposition}),
+time-reversible(Eq.~\ref{introEquation:timeReversible}) and
+volume-preserving (Eq.~\ref{introEquation:volumePreserving}). These
 geometric properties attribute to its long-time stability and its
 popularity in the community. However, the most commonly used
 velocity verlet integration scheme is written as below,
@@ -781,11 +735,10 @@ the equations of motion would follow:
 
 \item Use the half step velocities to move positions one whole step, $\Delta t$.
 
-\item Evaluate the forces at the new positions, $\mathbf{q}(\Delta t)$, and use the new forces to complete the velocity move.
+\item Evaluate the forces at the new positions, $q(\Delta t)$, and use the new forces to complete the velocity move.
 
 \item Repeat from step 1 with the new position, velocities, and forces assuming the roles of the initial values.
 \end{enumerate}
-
 By simply switching the order of the propagators in the splitting
 and composing a new integrator, the \emph{position verlet}
 integrator, can be generated,
@@ -801,11 +754,13 @@ The Baker-Campbell-Hausdorff formula can be used to de
 
 \subsubsection{\label{introSection:errorAnalysis}\textbf{Error Analysis and Higher Order Methods}}
 
-The Baker-Campbell-Hausdorff formula can be used to determine the
-local error of splitting method in terms of the commutator of the
-operators(\ref{introEquation:exponentialOperator}) associated with
-the sub-flow. For operators $hX$ and $hY$ which are associated with
-$\varphi_1(t)$ and $\varphi_2(t)$ respectively , we have
+The Baker-Campbell-Hausdorff formula\cite{Gilmore1974} can be used
+to determine the local error of a splitting method in terms of the
+commutator of the
+operators(Eq.~\ref{introEquation:exponentialOperator}) associated
+with the sub-propagator. For operators $hX$ and $hY$ which are
+associated with $\varphi_1(t)$ and $\varphi_2(t)$ respectively , we
+have
 \begin{equation}
 \exp (hX + hY) = \exp (hZ)
 \end{equation}
@@ -814,27 +769,28 @@ Here, $[X,Y]$ is the commutators of operator $X$ and $
 hZ = hX + hY + \frac{{h^2 }}{2}[X,Y] + \frac{{h^3 }}{2}\left(
 {[X,[X,Y]] + [Y,[Y,X]]} \right) +  \ldots .
 \end{equation}
-Here, $[X,Y]$ is the commutators of operator $X$ and $Y$ given by
+Here, $[X,Y]$ is the commutator of operator $X$ and $Y$ given by
 \[
 [X,Y] = XY - YX .
 \]
 Applying the Baker-Campbell-Hausdorff formula\cite{Varadarajan1974}
-to the Sprang splitting, we can obtain
+to the Strang splitting, we can obtain
 \begin{eqnarray*}
 \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 \\
                                    &   & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\
-                                   &   & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3[X,[X,Y]]/24 + \ldots )
+                                   &   & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3[X,[X,Y]]/24 + \ldots
+                                   ).
 \end{eqnarray*}
-Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0,\] the dominant local
-error of Spring splitting is proportional to $h^3$. The same
-procedure can be applied to a general splitting,  of the form
+Since $ [X,Y] + [Y,X] = 0$ and $ [X,X] = 0$, the dominant local
+error of Strang splitting is proportional to $h^3$. The same
+procedure can be applied to a general splitting of the form
 \begin{equation}
 \varphi _{b_m h}^2  \circ \varphi _{a_m h}^1  \circ \varphi _{b_{m -
 1} h}^2  \circ  \ldots  \circ \varphi _{a_1 h}^1 .
 \end{equation}
 A careful choice of coefficient $a_1 \ldots b_m$ will lead to higher
 order methods. Yoshida proposed an elegant way to compose higher
-order methods based on symmetric splitting\cite{Yoshida1990}. Given
+order methods based on symmetric splitting.\cite{Yoshida1990} Given
 a symmetric second order base method $ \varphi _h^{(2)} $, a
 fourth-order symmetric method can be constructed by composing,
 \[
@@ -862,14 +818,14 @@ simulations. For instance, instantaneous temperature o
 dynamical information. The basic idea of molecular dynamics is that
 macroscopic properties are related to microscopic behavior and
 microscopic behavior can be calculated from the trajectories in
-simulations. For instance, instantaneous temperature of an
-Hamiltonian system of $N$ particle can be measured by
+simulations. For instance, instantaneous temperature of a
+Hamiltonian system of $N$ particles can be measured by
 \[
 T = \sum\limits_{i = 1}^N {\frac{{m_i v_i^2 }}{{fk_B }}}
 \]
 where $m_i$ and $v_i$ are the mass and velocity of $i$th particle
 respectively, $f$ is the number of degrees of freedom, and $k_B$ is
-the boltzman constant.
+the Boltzman constant.
 
 A typical molecular dynamics run consists of three essential steps:
 \begin{enumerate}
@@ -886,9 +842,9 @@ will discusse issues in production run.
 These three individual steps will be covered in the following
 sections. Sec.~\ref{introSec:initialSystemSettings} deals with the
 initialization of a simulation. Sec.~\ref{introSection:production}
-will discusse issues in production run.
+ discusses issues of production runs.
 Sec.~\ref{introSection:Analysis} provides the theoretical tools for
-trajectory analysis.
+analysis of trajectories.
 
 \subsection{\label{introSec:initialSystemSettings}Initialization}
 
@@ -900,9 +856,9 @@ structure, some important information is missing. For 
 thousands of crystal structures of molecules are discovered every
 year, many more remain unknown due to the difficulties of
 purification and crystallization. Even for molecules with known
-structure, some important information is missing. For example, a
+structures, some important information is missing. For example, a
 missing hydrogen atom which acts as donor in hydrogen bonding must
-be added. Moreover, in order to include electrostatic interaction,
+be added. Moreover, in order to include electrostatic interactions,
 one may need to specify the partial charges for individual atoms.
 Under some circumstances, we may even need to prepare the system in
 a special configuration. For instance, when studying transport
@@ -920,31 +876,30 @@ near the minimum, steepest descent method is extremely
 minimization to find a more reasonable conformation. Several energy
 minimization methods have been developed to exploit the energy
 surface and to locate the local minimum. While converging slowly
-near the minimum, steepest descent method is extremely robust when
+near the minimum, the steepest descent method is extremely robust when
 systems are strongly anharmonic. Thus, it is often used to refine
-structure from crystallographic data. Relied on the gradient or
-hessian, advanced methods like Newton-Raphson converge rapidly to a
-local minimum, but become unstable if the energy surface is far from
+structures from crystallographic data. Relying on the Hessian,
+advanced methods like Newton-Raphson converge rapidly to a local
+minimum, but become unstable if the energy surface is far from
 quadratic. Another factor that must be taken into account, when
 choosing energy minimization method, is the size of the system.
 Steepest descent and conjugate gradient can deal with models of any
 size. Because of the limits on computer memory to store the hessian
-matrix and the computing power needed to diagonalized these
-matrices, most Newton-Raphson methods can not be used with very
-large systems.
+matrix and the computing power needed to diagonalize these matrices,
+most Newton-Raphson methods can not be used with very large systems.
 
 \subsubsection{\textbf{Heating}}
 
-Typically, Heating is performed by assigning random velocities
+Typically, heating is performed by assigning random velocities
 according to a Maxwell-Boltzman distribution for a desired
 temperature. Beginning at a lower temperature and gradually
 increasing the temperature by assigning larger random velocities, we
-end up with setting the temperature of the system to a final
-temperature at which the simulation will be conducted. In heating
-phase, we should also keep the system from drifting or rotating as a
-whole. To do this, the net linear momentum and angular momentum of
-the system is shifted to zero after each resampling from the Maxwell
--Boltzman distribution.
+end up setting the temperature of the system to a final temperature
+at which the simulation will be conducted. In the heating phase, we
+should also keep the system from drifting or rotating as a whole. To
+do this, the net linear momentum and angular momentum of the system
+is shifted to zero after each resampling from the Maxwell -Boltzman
+distribution.
 
 \subsubsection{\textbf{Equilibration}}
 
@@ -955,7 +910,7 @@ as a means to arrive at an equilibrated structure in a
 properties \textit{etc}, become independent of time. Strictly
 speaking, minimization and heating are not necessary, provided the
 equilibration process is long enough. However, these steps can serve
-as a means to arrive at an equilibrated structure in an effective
+as a mean to arrive at an equilibrated structure in an effective
 way.
 
 \subsection{\label{introSection:production}Production}
@@ -971,21 +926,19 @@ which making large simulations prohibitive in the abse
 calculation of non-bonded forces, such as van der Waals force and
 Coulombic forces \textit{etc}. For a system of $N$ particles, the
 complexity of the algorithm for pair-wise interactions is $O(N^2 )$,
-which making large simulations prohibitive in the absence of any
-algorithmic tricks.
-
-A natural approach to avoid system size issues is to represent the
-bulk behavior by a finite number of the particles. However, this
-approach will suffer from the surface effect at the edges of the
-simulation. To offset this, \textit{Periodic boundary conditions}
-(see Fig.~\ref{introFig:pbc}) is developed to simulate bulk
-properties with a relatively small number of particles. In this
-method, the simulation box is replicated throughout space to form an
-infinite lattice. During the simulation, when a particle moves in
-the primary cell, its image in other cells move in exactly the same
-direction with exactly the same orientation. Thus, as a particle
-leaves the primary cell, one of its images will enter through the
-opposite face.
+which makes large simulations prohibitive in the absence of any
+algorithmic tricks. A natural approach to avoid system size issues
+is to represent the bulk behavior by a finite number of the
+particles. However, this approach will suffer from surface effects
+at the edges of the simulation. To offset this, \textit{Periodic
+boundary conditions} (see Fig.~\ref{introFig:pbc}) were developed to
+simulate bulk properties with a relatively small number of
+particles. In this method, the simulation box is replicated
+throughout space to form an infinite lattice. During the simulation,
+when a particle moves in the primary cell, its image in other cells
+move in exactly the same direction with exactly the same
+orientation. Thus, as a particle leaves the primary cell, one of its
+images will enter through the opposite face.
 \begin{figure}
 \centering
 \includegraphics[width=\linewidth]{pbc.eps}
@@ -997,11 +950,11 @@ evaluation is to apply spherical cutoff where particle
 
 %cutoff and minimum image convention
 Another important technique to improve the efficiency of force
-evaluation is to apply spherical cutoff where particles farther than
-a predetermined distance are not included in the calculation
-\cite{Frenkel1996}. The use of a cutoff radius will cause a
-discontinuity in the potential energy curve. Fortunately, one can
-shift simple radial potential to ensure the potential curve go
+evaluation is to apply spherical cutoffs where particles farther
+than a predetermined distance are not included in the
+calculation.\cite{Frenkel1996} The use of a cutoff radius will cause
+a discontinuity in the potential energy curve. Fortunately, one can
+shift a simple radial potential to ensure the potential curve go
 smoothly to zero at the cutoff radius. The cutoff strategy works
 well for Lennard-Jones interaction because of its short range
 nature. However, simply truncating the electrostatic interaction
@@ -1009,9 +962,9 @@ periodicity artifacts in liquid simulations. Taking th
 in simulations. The Ewald summation, in which the slowly decaying
 Coulomb potential is transformed into direct and reciprocal sums
 with rapid and absolute convergence, has proved to minimize the
-periodicity artifacts in liquid simulations. Taking the advantages
-of the fast Fourier transform (FFT) for calculating discrete Fourier
-transforms, the particle mesh-based
+periodicity artifacts in liquid simulations. Taking advantage of
+fast Fourier transform (FFT) techniques for calculating discrete
+Fourier transforms, the particle mesh-based
 methods\cite{Hockney1981,Shimada1993, Luty1994} are accelerated from
 $O(N^{3/2})$ to $O(N logN)$. An alternative approach is the
 \emph{fast multipole method}\cite{Greengard1987, Greengard1994},
@@ -1021,13 +974,13 @@ his coworkers\cite{Wolf1999}. The shifted Coulomb pote
 simulation community, these two methods are difficult to implement
 correctly and efficiently. Instead, we use a damped and
 charge-neutralized Coulomb potential method developed by Wolf and
-his coworkers\cite{Wolf1999}. The shifted Coulomb potential for
+his coworkers.\cite{Wolf1999} The shifted Coulomb potential for
 particle $i$ and particle $j$ at distance $r_{rj}$ is given by:
 \begin{equation}
 V(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha
 r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow
 R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha
-r_{ij})}{r_{ij}}\right\}. \label{introEquation:shiftedCoulomb}
+r_{ij})}{r_{ij}}\right\}, \label{introEquation:shiftedCoulomb}
 \end{equation}
 where $\alpha$ is the convergence parameter. Due to the lack of
 inherent periodicity and rapid convergence,this method is extremely
@@ -1044,14 +997,15 @@ Recently, advanced visualization technique have become
 
 \subsection{\label{introSection:Analysis} Analysis}
 
-Recently, advanced visualization technique have become applied to
+Recently, advanced visualization techniques have been applied to
 monitor the motions of molecules. Although the dynamics of the
 system can be described qualitatively from animation, quantitative
-trajectory analysis are more useful. According to the principles of
-Statistical Mechanics, Sec.~\ref{introSection:statisticalMechanics},
-one can compute thermodynamic properties, analyze fluctuations of
-structural parameters, and investigate time-dependent processes of
-the molecule from the trajectories.
+trajectory analysis is more useful. According to the principles of
+Statistical Mechanics in
+Sec.~\ref{introSection:statisticalMechanics}, one can compute
+thermodynamic properties, analyze fluctuations of structural
+parameters, and investigate time-dependent processes of the molecule
+from the trajectories.
 
 \subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamic Properties}}
 
@@ -1081,22 +1035,21 @@ Experimentally, pair distribution function can be gath
 distribution functions. Among these functions,the \emph{pair
 distribution function}, also known as \emph{radial distribution
 function}, is of most fundamental importance to liquid theory.
-Experimentally, pair distribution function can be gathered by
+Experimentally, pair distribution functions can be gathered by
 Fourier transforming raw data from a series of neutron diffraction
-experiments and integrating over the surface factor
-\cite{Powles1973}. The experimental results can serve as a criterion
-to justify the correctness of a liquid model. Moreover, various
-equilibrium thermodynamic and structural properties can also be
-expressed in terms of radial distribution function \cite{Allen1987}.
-
-The pair distribution functions $g(r)$ gives the probability that a
-particle $i$ will be located at a distance $r$ from a another
-particle $j$ in the system
-\[
+experiments and integrating over the surface
+factor.\cite{Powles1973} The experimental results can serve as a
+criterion to justify the correctness of a liquid model. Moreover,
+various equilibrium thermodynamic and structural properties can also
+be expressed in terms of the radial distribution
+function.\cite{Allen1987} The pair distribution functions $g(r)$
+gives the probability that a particle $i$ will be located at a
+distance $r$ from a another particle $j$ in the system
+\begin{equation}
 g(r) = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j
 \ne i} {\delta (r - r_{ij} )} } } \right\rangle = \frac{\rho
 (r)}{\rho}.
-\]
+\end{equation}
 Note that the delta function can be replaced by a histogram in
 computer simulation. Peaks in $g(r)$ represent solvent shells, and
 the height of these peaks gradually decreases to 1 as the liquid of
@@ -1114,92 +1067,90 @@ function is called an \emph{autocorrelation function}.
 \label{introEquation:timeCorrelationFunction}
 \end{equation}
 If $A$ and $B$ refer to same variable, this kind of correlation
-function is called an \emph{autocorrelation function}. One example
-of an auto correlation function is the velocity auto-correlation
+functions are called \emph{autocorrelation functions}. One typical example is the velocity autocorrelation
 function which is directly related to transport properties of
 molecular liquids:
-\[
+\begin{equation}
 D = \frac{1}{3}\int\limits_0^\infty  {\left\langle {v(t) \cdot v(0)}
 \right\rangle } dt
-\]
+\end{equation}
 where $D$ is diffusion constant. Unlike the velocity autocorrelation
-function, which is averaging over time origins and over all the
-atoms, the dipole autocorrelation functions are calculated for the
+function, which is averaged over time origins and over all the
+atoms, the dipole autocorrelation functions is calculated for the
 entire system. The dipole autocorrelation function is given by:
-\[
+\begin{equation}
 c_{dipole}  = \left\langle {u_{tot} (t) \cdot u_{tot} (t)}
 \right\rangle
-\]
+\end{equation}
 Here $u_{tot}$ is the net dipole of the entire system and is given
 by
-\[
-u_{tot} (t) = \sum\limits_i {u_i (t)}
-\]
-In principle, many time correlation functions can be related with
+\begin{equation}
+u_{tot} (t) = \sum\limits_i {u_i (t)}.
+\end{equation}
+In principle, many time correlation functions can be related to
 Fourier transforms of the infrared, Raman, and inelastic neutron
 scattering spectra of molecular liquids. In practice, one can
-extract the IR spectrum from the intensity of dipole fluctuation at
-each frequency using the following relationship:
-\[
+extract the IR spectrum from the intensity of the molecular dipole
+fluctuation at each frequency using the following relationship:
+\begin{equation}
 \hat c_{dipole} (v) = \int_{ - \infty }^\infty  {c_{dipole} (t)e^{ -
-i2\pi vt} dt}
-\]
+i2\pi vt} dt}.
+\end{equation}
 
 \section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
 
 Rigid bodies are frequently involved in the modeling of different
-areas, from engineering, physics, to chemistry. For example,
-missiles and vehicle are usually modeled by rigid bodies.  The
-movement of the objects in 3D gaming engine or other physics
-simulator is governed by rigid body dynamics. In molecular
+areas, including engineering, physics and chemistry. For example,
+missiles and vehicles are usually modeled by rigid bodies.  The
+movement of the objects in 3D gaming engines or other physics
+simulators is governed by rigid body dynamics. In molecular
 simulations, rigid bodies are used to simplify protein-protein
-docking studies\cite{Gray2003}.
+docking studies.\cite{Gray2003}
 
 It is very important to develop stable and efficient methods to
 integrate the equations of motion for orientational degrees of
 freedom. Euler angles are the natural choice to describe the
 rotational degrees of freedom. However, due to $\frac {1}{sin
 \theta}$ singularities, the numerical integration of corresponding
-equations of motion is very inefficient and inaccurate. Although an
-alternative integrator using multiple sets of Euler angles can
-overcome this difficulty\cite{Barojas1973}, the computational
-penalty and the loss of angular momentum conservation still remain.
-A singularity-free representation utilizing quaternions was
-developed by Evans in 1977\cite{Evans1977}. Unfortunately, this
-approach uses a nonseparable Hamiltonian resulting from the
-quaternion representation, which prevents the symplectic algorithm
-to be utilized. Another different approach is to apply holonomic
-constraints to the atoms belonging to the rigid body. Each atom
-moves independently under the normal forces deriving from potential
-energy and constraint forces which are used to guarantee the
-rigidness. However, due to their iterative nature, the SHAKE and
-Rattle algorithms also converge very slowly when the number of
-constraints increases\cite{Ryckaert1977, Andersen1983}.
+equations of these motion is very inefficient and inaccurate.
+Although an alternative integrator using multiple sets of Euler
+angles can overcome this difficulty\cite{Barojas1973}, the
+computational penalty and the loss of angular momentum conservation
+still remain. A singularity-free representation utilizing
+quaternions was developed by Evans in 1977.\cite{Evans1977}
+Unfortunately, this approach used a nonseparable Hamiltonian
+resulting from the quaternion representation, which prevented the
+symplectic algorithm from being utilized. Another different approach
+is to apply holonomic constraints to the atoms belonging to the
+rigid body. Each atom moves independently under the normal forces
+deriving from potential energy and constraint forces which are used
+to guarantee the rigidness. However, due to their iterative nature,
+the SHAKE and Rattle algorithms also converge very slowly when the
+number of constraints increases.\cite{Ryckaert1977, Andersen1983}
 
 A break-through in geometric literature suggests that, in order to
 develop a long-term integration scheme, one should preserve the
-symplectic structure of the flow. By introducing a conjugate
+symplectic structure of the propagator. By introducing a conjugate
 momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's
 equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was
 proposed to evolve the Hamiltonian system in a constraint manifold
 by iteratively satisfying the orthogonality constraint $Q^T Q = 1$.
 An alternative method using the quaternion representation was
-developed by Omelyan\cite{Omelyan1998}. However, both of these
+developed by Omelyan.\cite{Omelyan1998} However, both of these
 methods are iterative and inefficient. In this section, we descibe a
-symplectic Lie-Poisson integrator for rigid body developed by
+symplectic Lie-Poisson integrator for rigid bodies developed by
 Dullweber and his coworkers\cite{Dullweber1997} in depth.
 
 \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Bodies}
-The motion of a rigid body is Hamiltonian with the Hamiltonian
-function
+The Hamiltonian of a rigid body is given by
 \begin{equation}
 H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) +
 V(q,Q) + \frac{1}{2}tr[(QQ^T  - 1)\Lambda ].
 \label{introEquation:RBHamiltonian}
 \end{equation}
-Here, $q$ and $Q$  are the position and rotation matrix for the
-rigid-body, $p$ and $P$  are conjugate momenta to $q$  and $Q$ , and
-$J$, a diagonal matrix, is defined by
+Here, $q$ and $Q$  are the position vector and rotation matrix for
+the rigid-body, $p$ and $P$  are conjugate momenta to $q$  and $Q$ ,
+and $J$, a diagonal matrix, is defined by
 \[
 I_{ii}^{ - 1}  = \frac{1}{2}\sum\limits_{i \ne j} {J_{jj}^{ - 1} }
 \]
@@ -1209,54 +1160,48 @@ which is used to ensure rotation matrix's unitarity. D
 \begin{equation}
 Q^T Q = 1, \label{introEquation:orthogonalConstraint}
 \end{equation}
-which is used to ensure rotation matrix's unitarity. Differentiating
-\ref{introEquation:orthogonalConstraint} and using Equation
-\ref{introEquation:RBMotionMomentum}, one may obtain,
-\begin{equation}
-Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
-\label{introEquation:RBFirstOrderConstraint}
-\end{equation}
-
-Using Equation (\ref{introEquation:motionHamiltonianCoordinate},
-\ref{introEquation:motionHamiltonianMomentum}), one can write down
+which is used to ensure the rotation matrix's unitarity. Using
+Eq.~\ref{introEquation:motionHamiltonianCoordinate} and Eq.~
+\ref{introEquation:motionHamiltonianMomentum}, one can write down
 the equations of motion,
-
 \begin{eqnarray}
- \frac{{dq}}{{dt}} & = & \frac{p}{m} \label{introEquation:RBMotionPosition}\\
- \frac{{dp}}{{dt}} & = & - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
- \frac{{dQ}}{{dt}} & = & PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
+ \frac{{dq}}{{dt}} & = & \frac{p}{m}, \label{introEquation:RBMotionPosition}\\
+ \frac{{dp}}{{dt}} & = & - \nabla _q V(q,Q), \label{introEquation:RBMotionMomentum}\\
+ \frac{{dQ}}{{dt}} & = & PJ^{ - 1},  \label{introEquation:RBMotionRotation}\\
  \frac{{dP}}{{dt}} & = & - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}
 \end{eqnarray}
-
+Differentiating Eq.~\ref{introEquation:orthogonalConstraint} and
+using Eq.~\ref{introEquation:RBMotionMomentum}, one may obtain,
+\begin{equation}
+Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
+\label{introEquation:RBFirstOrderConstraint}
+\end{equation}
 In general, there are two ways to satisfy the holonomic constraints.
 We can use a constraint force provided by a Lagrange multiplier on
-the normal manifold to keep the motion on constraint space. Or we
-can simply evolve the system on the constraint manifold. These two
-methods have been proved to be equivalent. The holonomic constraint
-and equations of motions define a constraint manifold for rigid
-bodies
+the normal manifold to keep the motion on the constraint space. Or
+we can simply evolve the system on the constraint manifold. These
+two methods have been proved to be equivalent. The holonomic
+constraint and equations of motions define a constraint manifold for
+rigid bodies
 \[
 M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0}
 \right\}.
 \]
-
-Unfortunately, this constraint manifold is not the cotangent bundle
-$T^* SO(3)$ which can be consider as a symplectic manifold on Lie
-rotation group $SO(3)$. However, it turns out that under symplectic
-transformation, the cotangent space and the phase space are
-diffeomorphic. By introducing
+Unfortunately, this constraint manifold is not $T^* SO(3)$ which is
+a symplectic manifold on Lie rotation group $SO(3)$. However, it
+turns out that under symplectic transformation, the cotangent space
+and the phase space are diffeomorphic. By introducing
 \[
 \tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
 \]
-the mechanical system subject to a holonomic constraint manifold $M$
+the mechanical system subjected to a holonomic constraint manifold $M$
 can be re-formulated as a Hamiltonian system on the cotangent space
 \[
 T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \tilde Q =
 1,\tilde Q^T \tilde PJ^{ - 1}  + J^{ - 1} P^T \tilde Q = 0} \right\}
 \]
-
 For a body fixed vector $X_i$ with respect to the center of mass of
-the rigid body, its corresponding lab fixed vector $X_0^{lab}$  is
+the rigid body, its corresponding lab fixed vector $X_i^{lab}$  is
 given as
 \begin{equation}
 X_i^{lab} = Q X_i + q.
@@ -1273,28 +1218,19 @@ respectively.
 \[
 \nabla _Q V(q,Q) = F(q,Q)X_i^t
 \]
-respectively.
-
-As a common choice to describe the rotation dynamics of the rigid
-body, the angular momentum on the body fixed frame $\Pi  = Q^t P$ is
-introduced to rewrite the equations of motion,
+respectively. As a common choice to describe the rotation dynamics
+of the rigid body, the angular momentum on the body fixed frame $\Pi
+= Q^t P$ is introduced to rewrite the equations of motion,
 \begin{equation}
  \begin{array}{l}
- \dot \Pi  = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda  \\
- \dot Q  = Q\Pi {\rm{ }}J^{ - 1}  \\
+ \dot \Pi  = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda,  \\
+ \dot Q  = Q\Pi {\rm{ }}J^{ - 1},  \\
  \end{array}
  \label{introEqaution:RBMotionPI}
 \end{equation}
-, as well as holonomic constraints,
-\[
-\begin{array}{l}
- \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
- Q^T Q = 1 \\
- \end{array}
-\]
-
-For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a matrix $\hat v \in
-so(3)^ \star$, the hat-map isomorphism,
+as well as holonomic constraints $\Pi J^{ - 1}  + J^{ - 1} \Pi ^t  =
+0$ and $Q^T Q = 1$. For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a
+matrix $\hat v \in so(3)^ \star$, the hat-map isomorphism,
 \begin{equation}
 v(v_1 ,v_2 ,v_3 ) \Leftrightarrow \hat v = \left(
 {\begin{array}{*{20}c}
@@ -1307,26 +1243,22 @@ operations
 will let us associate the matrix products with traditional vector
 operations
 \[
-\hat vu = v \times u
+\hat vu = v \times u.
 \]
-Using \ref{introEqaution:RBMotionPI}, one can construct a skew
+Using Eq.~\ref{introEqaution:RBMotionPI}, one can construct a skew
 matrix,
-
-\begin{eqnarry*}
-(\dot \Pi  - \dot \Pi ^T ){\rm{ }} = {\rm{ }}(\Pi  - \Pi ^T ){\rm{
-}}(J^{ - 1} \Pi  + \Pi J^{ - 1} ) + \sum\limits_i {[Q^T F_i
-(r,Q)X_i^T  - X_i F_i (r,Q)^T Q]}  - (\Lambda  - \Lambda ^T ).
-\label{introEquation:skewMatrixPI}
-\end{eqnarray*}
-
-Since $\Lambda$ is symmetric, the last term of Equation
-\ref{introEquation:skewMatrixPI} is zero, which implies the Lagrange
-multiplier $\Lambda$ is absent from the equations of motion. This
-unique property eliminates the requirement of iterations which can
-not be avoided in other methods\cite{Kol1997, Omelyan1998}.
-
-Applying the hat-map isomorphism, we obtain the equation of motion
-for angular momentum on body frame
+\begin{eqnarray}
+(\dot \Pi  - \dot \Pi ^T )&= &(\Pi  - \Pi ^T )(J^{ - 1} \Pi  + \Pi J^{ - 1} ) \notag \\
+& & + \sum\limits_i {[Q^T F_i (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]}  -
+(\Lambda  - \Lambda ^T ). \label{introEquation:skewMatrixPI}
+\end{eqnarray}
+Since $\Lambda$ is symmetric, the last term of
+Eq.~\ref{introEquation:skewMatrixPI} is zero, which implies the
+Lagrange multiplier $\Lambda$ is absent from the equations of
+motion. This unique property eliminates the requirement of
+iterations which can not be avoided in other methods.\cite{Kol1997,
+Omelyan1998} Applying the hat-map isomorphism, we obtain the
+equation of motion for angular momentum in the body frame
 \begin{equation}
 \dot \pi  = \pi  \times I^{ - 1} \pi  + \sum\limits_i {\left( {Q^T
 F_i (r,Q)} \right) \times X_i }.
@@ -1335,11 +1267,11 @@ given by
 In the same manner, the equation of motion for rotation matrix is
 given by
 \[
-\dot Q = Qskew(I^{ - 1} \pi )
+\dot Q = Qskew(I^{ - 1} \pi ).
 \]
 
 \subsection{\label{introSection:SymplecticFreeRB}Symplectic
-Lie-Poisson Integrator for Free Rigid Body}
+Lie-Poisson Integrator for Free Rigid Bodies}
 
 If there are no external forces exerted on the rigid body, the only
 contribution to the rotational motion is from the kinetic energy
@@ -1357,28 +1289,27 @@ J(\pi ) = \left( {\begin{array}{*{20}c}
    0 & {\pi _3 } & { - \pi _2 }  \\
    { - \pi _3 } & 0 & {\pi _1 }  \\
    {\pi _2 } & { - \pi _1 } & 0  \\
-\end{array}} \right)
+\end{array}} \right).
 \end{equation}
 Thus, the dynamics of free rigid body is governed by
 \begin{equation}
-\frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
+\frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi ).
 \end{equation}
-
-One may notice that each $T_i^r$ in Equation
-\ref{introEquation:rotationalKineticRB} can be solved exactly. For
-instance, the equations of motion due to $T_1^r$ are given by
+One may notice that each $T_i^r$ in
+Eq.~\ref{introEquation:rotationalKineticRB} can be solved exactly.
+For instance, the equations of motion due to $T_1^r$ are given by
 \begin{equation}
 \frac{d}{{dt}}\pi  = R_1 \pi ,\frac{d}{{dt}}Q = QR_1
 \label{introEqaution:RBMotionSingleTerm}
 \end{equation}
-where
+with
 \[ R_1  = \left( {\begin{array}{*{20}c}
    0 & 0 & 0  \\
    0 & 0 & {\pi _1 }  \\
    0 & { - \pi _1 } & 0  \\
 \end{array}} \right).
 \]
-The solutions of Equation \ref{introEqaution:RBMotionSingleTerm} is
+The solutions of Eq.~\ref{introEqaution:RBMotionSingleTerm} is
 \[
 \pi (\Delta t) = e^{\Delta tR_1 } \pi (0),Q(\Delta t) =
 Q(0)e^{\Delta tR_1 }
@@ -1392,15 +1323,24 @@ tR_1 }$, we can use Cayley transformation to obtain a 
 \end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
 \]
 To reduce the cost of computing expensive functions in $e^{\Delta
-tR_1 }$, we can use Cayley transformation to obtain a single-aixs
-propagator,
-\[
-e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
-)
-\]
-The flow maps for $T_2^r$ and $T_3^r$ can be found in the same
+tR_1 }$, we can use the Cayley transformation to obtain a
+single-aixs propagator,
+\begin{eqnarray*}
+e^{\Delta tR_1 }  & \approx & (1 - \Delta tR_1 )^{ - 1} (1 + \Delta
+tR_1 ) \\
+%
+& \approx & \left( \begin{array}{ccc}
+1 & 0 & 0 \\
+0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
+\theta^2 / 4} \\
+0 & \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 +
+\theta^2 / 4}
+\end{array}
+\right).
+\end{eqnarray*}
+The propagators for $T_2^r$ and $T_3^r$ can be found in the same
 manner. In order to construct a second-order symplectic method, we
-split the angular kinetic Hamiltonian function can into five terms
+split the angular kinetic Hamiltonian function into five terms
 \[
 T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
 ) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
@@ -1414,25 +1354,24 @@ _1 }.
 \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
 _1 }.
 \]
-
-The non-canonical Lie-Poisson bracket ${F, G}$ of two function
-$F(\pi )$ and $G(\pi )$ is defined by
+The non-canonical Lie-Poisson bracket $\{F, G\}$ of two functions $F(\pi )$ and $G(\pi )$ is defined by
 \[
 \{ F,G\} (\pi ) = [\nabla _\pi  F(\pi )]^T J(\pi )\nabla _\pi  G(\pi
-)
+).
 \]
 If the Poisson bracket of a function $F$ with an arbitrary smooth
 function $G$ is zero, $F$ is a \emph{Casimir}, which is the
 conserved quantity in Poisson system. We can easily verify that the
 norm of the angular momentum, $\parallel \pi
-\parallel$, is a \emph{Casimir}. Let$ F(\pi ) = S(\frac{{\parallel
+\parallel$, is a \emph{Casimir}.\cite{McLachlan1993} Let $F(\pi ) = S(\frac{{\parallel
 \pi \parallel ^2 }}{2})$ for an arbitrary function $ S:R \to R$ ,
 then by the chain rule
 \[
 \nabla _\pi  F(\pi ) = S'(\frac{{\parallel \pi \parallel ^2
-}}{2})\pi
+}}{2})\pi.
 \]
-Thus $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel \pi
+Thus, $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel
+\pi
 \parallel ^2 }}{2})\pi  \times \pi  = 0 $. This explicit
 Lie-Poisson integrator is found to be both extremely efficient and
 stable. These properties can be explained by the fact the small
@@ -1443,14 +1382,12 @@ energy and potential energy,
 Splitting for Rigid Body}
 
 The Hamiltonian of rigid body can be separated in terms of kinetic
-energy and potential energy,
-\[
-H = T(p,\pi ) + V(q,Q)
-\]
-The equations of motion corresponding to potential energy and
-kinetic energy are listed in the below table,
+energy and potential energy, $H = T(p,\pi ) + V(q,Q)$. The equations
+of motion corresponding to potential energy and kinetic energy are
+listed in Table~\ref{introTable:rbEquations}.
 \begin{table}
 \caption{EQUATIONS OF MOTION DUE TO POTENTIAL AND KINETIC ENERGIES}
+\label{introTable:rbEquations}
 \begin{center}
 \begin{tabular}{|l|l|}
   \hline
@@ -1486,32 +1423,28 @@ defined by \ref{introEquation:rotationalKineticRB}. Th
 T(p,\pi ) =T^t (p) + T^r (\pi ).
 \end{equation}
 where $ T^t (p) = \frac{1}{2}p^T m^{ - 1} p $ and $T^r (\pi )$ is
-defined by \ref{introEquation:rotationalKineticRB}. Therefore, the
-corresponding propagators are given by
+defined by Eq.~\ref{introEquation:rotationalKineticRB}. Therefore,
+the corresponding propagators are given by
 \[
 \varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
 _{\Delta t,T^r }.
 \]
 Finally, we obtain the overall symplectic propagators for freely
 moving rigid bodies
-\begin{equation}
-\begin{array}{c}
- \varphi _{\Delta t}  = \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \\
-  \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 }  \\
-  \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
- \end{array}
+\begin{eqnarray}
+ \varphi _{\Delta t}  &=& \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \notag\\
+  & & \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\\
+  & & \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .
 \label{introEquation:overallRBFlowMaps}
-\end{equation}
+\end{eqnarray}
 
 \section{\label{introSection:langevinDynamics}Langevin Dynamics}
 As an alternative to newtonian dynamics, Langevin dynamics, which
 mimics a simple heat bath with stochastic and dissipative forces,
 has been applied in a variety of studies. This section will review
-the theory of Langevin dynamics. A brief derivation of generalized
+the theory of Langevin dynamics. A brief derivation of the generalized
 Langevin equation will be given first. Following that, we will
-discuss the physical meaning of the terms appearing in the equation
-as well as the calculation of friction tensor from hydrodynamics
-theory.
+discuss the physical meaning of the terms appearing in the equation.
 
 \subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
 
@@ -1520,7 +1453,7 @@ Dynamics (GLE). Lets consider a system, in which the d
 environment, has been widely used in quantum chemistry and
 statistical mechanics. One of the successful applications of
 Harmonic bath model is the derivation of the Generalized Langevin
-Dynamics (GLE). Lets consider a system, in which the degree of
+Dynamics (GLE). Consider a system, in which the degree of
 freedom $x$ is assumed to couple to the bath linearly, giving a
 Hamiltonian of the form
 \begin{equation}
@@ -1531,7 +1464,7 @@ H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_
 with this degree of freedom, $H_B$ is a harmonic bath Hamiltonian,
 \[
 H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
-}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  \omega _\alpha ^2 }
+}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  x_\alpha ^2 }
 \right\}}
 \]
 where the index $\alpha$ runs over all the bath degrees of freedom,
@@ -1547,14 +1480,13 @@ W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\a
 \[
 W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
-\] and combining the last two terms in Equation
-\ref{introEquation:bathGLE}, we may rewrite the Harmonic bath
-Hamiltonian as
+\]
+and combining the last two terms in Eq.~\ref{introEquation:bathGLE}, we may rewrite the Harmonic bath Hamiltonian as
 \[
 H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
 {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
 w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
-w_\alpha ^2 }}x} \right)^2 } \right\}}
+w_\alpha ^2 }}x} \right)^2 } \right\}}.
 \]
 Since the first two terms of the new Hamiltonian depend only on the
 system coordinates, we can get the equations of motion for
@@ -1571,25 +1503,21 @@ m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x
 \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right).
 \label{introEquation:bathMotionGLE}
 \end{equation}
-
 In order to derive an equation for $x$, the dynamics of the bath
 variables $x_\alpha$ must be solved exactly first. As an integral
 transform which is particularly useful in solving linear ordinary
 differential equations,the Laplace transform is the appropriate tool
 to solve this problem. The basic idea is to transform the difficult
 differential equations into simple algebra problems which can be
-solved easily. Then, by applying the inverse Laplace transform, also
-known as the Bromwich integral, we can retrieve the solutions of the
-original problems.
-
-Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace
-transform of f(t) is a new function defined as
+solved easily. Then, by applying the inverse Laplace transform, we
+can retrieve the solutions of the original problems. Let $f(t)$ be a
+function defined on $ [0,\infty ) $, the Laplace transform of $f(t)$
+is a new function defined as
 \[
 L(f(t)) \equiv F(p) = \int_0^\infty  {f(t)e^{ - pt} dt}
 \]
 where  $p$ is real and  $L$ is called the Laplace Transform
-Operator. Below are some important properties of Laplace transform
-
+Operator. Below are some important properties of the Laplace transform
 \begin{eqnarray*}
  L(x + y)  & = & L(x) + L(y) \\
  L(ax)     & = & aL(x) \\
@@ -1597,20 +1525,17 @@ Operator. Below are some important properties of Lapla
  L(\ddot x)& = & p^2 L(x) - px(0) - \dot x(0) \\
  L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right)& = & G(p)H(p) \\
  \end{eqnarray*}
-
-
 Applying the Laplace transform to the bath coordinates, we obtain
 \begin{eqnarray*}
-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) \\
-L(x_\alpha  ) & = & \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
+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), \\
+L(x_\alpha  ) & = & \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}. \\
 \end{eqnarray*}
-
-By the same way, the system coordinates become
+In the same way, the system coordinates become
 \begin{eqnarray*}
- mL(\ddot x) & = & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
-  & & \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\}}  \\
+ mL(\ddot x) & = &
+  - \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\}}  \\
+  & & - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}}.
 \end{eqnarray*}
-
 With the help of some relatively important inverse Laplace
 transformations:
 \[
@@ -1620,7 +1545,7 @@ transformations:
  L(1) = \frac{1}{p} \\
  \end{array}
 \]
-, we obtain
+we obtain
 \begin{eqnarray*}
 m\ddot x & =  & - \frac{{\partial W(x)}}{{\partial x}} -
 \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
@@ -1629,12 +1554,12 @@ x_\alpha (0) - \frac{{g_\alpha  }}{{m_\alpha  \omega _
 & & + \sum\limits_{\alpha  = 1}^N {\left\{ {\left[ {g_\alpha
 x_\alpha (0) - \frac{{g_\alpha  }}{{m_\alpha  \omega _\alpha  }}}
 \right]\cos (\omega _\alpha  t) + \frac{{g_\alpha  \dot x_\alpha
-(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)} \right\}}
-\end{eqnarray*}
-\begin{eqnarray*}
-m\ddot x & = & - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
-{\sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
-}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha
+(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)} \right\}}\\
+%
+& = & -
+\frac{{\partial W(x)}}{{\partial x}} - \int_0^t {\sum\limits_{\alpha
+= 1}^N {\left( { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha
+^2 }}} \right)\cos (\omega _\alpha
 t)\dot x(t - \tau )d} \tau }  \\
 & & + \sum\limits_{\alpha  = 1}^N {\left\{ {\left[ {g_\alpha
 x_\alpha (0) - \frac{{g_\alpha }}{{m_\alpha \omega _\alpha  }}}
@@ -1661,24 +1586,18 @@ which is known as the \emph{generalized Langevin equat
 (t)\dot x(t - \tau )d\tau }  + R(t)
 \label{introEuqation:GeneralizedLangevinDynamics}
 \end{equation}
-which is known as the \emph{generalized Langevin equation}.
+which is known as the \emph{generalized Langevin equation} (GLE).
 
 \subsubsection{\label{introSection:randomForceDynamicFrictionKernel}\textbf{Random Force and Dynamic Friction Kernel}}
 
 One may notice that $R(t)$ depends only on initial conditions, which
 implies it is completely deterministic within the context of a
 harmonic bath. However, it is easy to verify that $R(t)$ is totally
-uncorrelated to $x$ and $\dot x$,
-\[
-\begin{array}{l}
- \left\langle {x(t)R(t)} \right\rangle  = 0, \\
- \left\langle {\dot x(t)R(t)} \right\rangle  = 0. \\
- \end{array}
-\]
-This property is what we expect from a truly random process. As long
-as the model chosen for $R(t)$ was a gaussian distribution in
+uncorrelated to $x$ and $\dot x$, $\left\langle {x(t)R(t)}
+\right\rangle  = 0, \left\langle {\dot x(t)R(t)} \right\rangle  =
+0.$ This property is what we expect from a truly random process. As
+long as the model chosen for $R(t)$ was a gaussian distribution in
 general, the stochastic nature of the GLE still remains.
-
 %dynamic friction kernel
 The convolution integral
 \[
@@ -1693,7 +1612,7 @@ and Equation \ref{introEuqation:GeneralizedLangevinDyn
 \[
 \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
 \]
-and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
+and Eq.~\ref{introEuqation:GeneralizedLangevinDynamics} becomes
 \[
 m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
 \frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
@@ -1703,14 +1622,14 @@ taken as a $delta$ function in time:
 infinitely quickly to motions in the system. Thus, $\xi (t)$ can be
 taken as a $delta$ function in time:
 \[
-\xi (t) = 2\xi _0 \delta (t)
+\xi (t) = 2\xi _0 \delta (t).
 \]
 Hence, the convolution integral becomes
 \[
 \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
 {\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
 \]
-and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
+and Eq.~\ref{introEuqation:GeneralizedLangevinDynamics} becomes
 \begin{equation}
 m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \xi _0 \dot
 x(t) + R(t) \label{introEquation:LangevinEquation}
@@ -1718,35 +1637,33 @@ briefly review on calculating friction tensor for arbi
 which is known as the Langevin equation. The static friction
 coefficient $\xi _0$ can either be calculated from spectral density
 or be determined by Stokes' law for regular shaped particles. A
-briefly review on calculating friction tensor for arbitrary shaped
+brief review on calculating friction tensors for arbitrary shaped
 particles is given in Sec.~\ref{introSection:frictionTensor}.
 
 \subsubsection{\label{introSection:secondFluctuationDissipation}\textbf{The Second Fluctuation Dissipation Theorem}}
 
-Defining a new set of coordinates,
+Defining a new set of coordinates
 \[
 q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
-^2 }}x(0)
-\],
-we can rewrite $R(T)$ as
+^2 }}x(0),
+\]
+we can rewrite $R(t)$ as
 \[
 R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
 \]
 And since the $q$ coordinates are harmonic oscillators,
-
 \begin{eqnarray*}
  \left\langle {q_\alpha ^2 } \right\rangle  & = & \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
  \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle & = & \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
  \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle & = &\delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
  \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 } }  \\
   & = &\sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
-  & = &kT\xi (t) \\
+  & = &kT\xi (t)
 \end{eqnarray*}
-
 Thus, we recover the \emph{second fluctuation dissipation theorem}
 \begin{equation}
 \xi (t) = \left\langle {R(t)R(0)} \right\rangle
-\label{introEquation:secondFluctuationDissipation}.
+\label{introEquation:secondFluctuationDissipation},
 \end{equation}
-In effect, it acts as a constraint on the possible ways in which one
-can model the random force and friction kernel.
+which acts as a constraint on the possible ways in which one can
+model the random force and friction kernel.
