--- trunk/tengDissertation/Introduction.tex	2006/06/05 21:24:52	2793
+++ trunk/tengDissertation/Introduction.tex	2006/06/28 17:36:32	2904
@@ -6,7 +6,7 @@ behind classical mechanics. Firstly, One can determine
 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
+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
@@ -17,24 +17,23 @@ Newton��s first law defines a class of inertial frames
 \subsection{\label{introSection:newtonian}Newtonian Mechanics}
 The discovery of Newton's three laws of mechanics which govern the
 motion of particles is the foundation of the classical mechanics.
-Newton��s first law defines a class of inertial frames. Inertial
+Newton's first law defines a class of inertial frames. Inertial
 frames are reference frames where a particle not interacting with
 other bodies will move with constant speed in the same direction.
-With respect to inertial frames Newton��s second law has the form
+With respect to inertial frames, Newton's second law has the form
 \begin{equation}
-F = \frac {dp}{dt} = \frac {mv}{dt}
+F = \frac {dp}{dt} = \frac {mdv}{dt}
 \label{introEquation:newtonSecondLaw}
 \end{equation}
 A point mass interacting with other bodies moves with the
 acceleration along the direction of the force acting on it. Let
 $F_{ij}$ be the force that particle $i$ exerts on particle $j$, and
 $F_{ji}$ be the force that particle $j$ exerts on particle $i$.
-Newton��s third law states that
+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
@@ -46,7 +45,7 @@ N \equiv r \times F \label{introEquation:torqueDefinit
 \end{equation}
 The torque $\tau$ with respect to the same origin is defined to be
 \begin{equation}
-N \equiv r \times F \label{introEquation:torqueDefinition}
+\tau \equiv r \times F \label{introEquation:torqueDefinition}
 \end{equation}
 Differentiating Eq.~\ref{introEquation:angularMomentumDefinition},
 \[
@@ -59,66 +58,60 @@ thus,
 \]
 thus,
 \begin{equation}
-\dot L = r \times \dot p = N
+\dot L = r \times \dot p = \tau
 \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
-scheme for rigid body \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: it
-describes their motion in special cartesian coordinate systems.
-Another limitation of Newtonian mechanics becomes obvious when we
-try to describe systems with large numbers of particles. It becomes
-very difficult to predict the properties of the system by carrying
-out calculations involving the each individual interaction between
-all the particles, 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, alternative procedures may be developed.
+Newtonian Mechanics suffers from a important limitation: motions 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 system, approximate numerical
+procedures may be developed.
 
-\subsubsection{\label{introSection:halmiltonPrinciple}Hamilton's
-Principle}
+\subsubsection{\label{introSection:halmiltonPrinciple}\textbf{Hamilton's
+Principle}}
 
 Hamilton introduced the dynamical principle upon which it is
-possible to base all of mechanics and, indeed, most of classical
-physics. Hamilton's Principle may be stated as follow,
-
-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$ \cite{Tolman1979}.
+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$,
 \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 part are derivable from a potential and the velocities are
-small compared with that of light, the Lagrangian function $L$ can
-be define as the difference between the kinetic energy of the system
-and its potential energy,
+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}
 
-\subsubsection{\label{introSection:equationOfMotionLagrangian}The
-Equations of Motion in Lagrangian Mechanics}
+\subsubsection{\label{introSection:equationOfMotionLagrangian}\textbf{The
+Equations of Motion in Lagrangian Mechanics}}
 
-For a holonomic system of $f$ degrees of freedom, the equations of
-motion in the Lagrangian form is
+For a system of $f$ degrees of freedom, the equations of motion in
+the Lagrangian form is
 \begin{equation}
 \frac{d}{{dt}}\frac{{\partial L}}{{\partial \dot q_i }} -
 \frac{{\partial L}}{{\partial q_i }} = 0,{\rm{ }}i = 1, \ldots,f
@@ -132,8 +125,7 @@ independent of generalized velocities, the generalized
 Arising from Lagrangian Mechanics, Hamiltonian Mechanics was
 introduced by William Rowan Hamilton in 1833 as a re-formulation of
 classical mechanics. If the potential energy of a system is
-independent of generalized velocities, the generalized momenta can
-be defined as
+independent of velocities, the momenta can be defined as
 \begin{equation}
 p_i = \frac{\partial L}{\partial \dot q_i}
 \label{introEquation:generalizedMomenta}
@@ -143,7 +135,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}
@@ -151,32 +142,30 @@ $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
 find
 \begin{equation}
-\frac{{\partial H}}{{\partial p_k }} = q_k
+\frac{{\partial H}}{{\partial p_k }} = \dot {q_k}
 \label{introEquation:motionHamiltonianCoordinate}
 \end{equation}
 \begin{equation}
-\frac{{\partial H}}{{\partial q_k }} =  - p_k
+\frac{{\partial H}}{{\partial q_k }} =  - \dot {p_k}
 \label{introEquation:motionHamiltonianMomentum}
 \end{equation}
 and
@@ -185,34 +174,31 @@ 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}.
 
 An important difference between Lagrangian approach and the
 Hamiltonian approach is that the Lagrangian is considered to be a
-function of the generalized velocities $\dot q_i$ and the
-generalized coordinates $q_i$, while the Hamiltonian is considered
-to be a function of the generalized momenta $p_i$ and the conjugate
-generalized coordinate $q_i$. 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}.
-
+function of the generalized velocities $\dot q_i$ and coordinates
+$q_i$, while the Hamiltonian is considered to be a function of the
+generalized momenta $p_i$ and the conjugate coordinates $q_i$.
+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}.
 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
@@ -230,42 +216,33 @@ momentum variables. Consider a dynamic system in a car
 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 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 = (q_1 , \ldots ,q_f ,p_1 ,
-\ldots ,p_f )$, with a unique set of values of $6f$ coordinates 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
 
-A microscopic state or microstate of a classical system is
-specification of the complete phase space vector of a system at any
-instant in time. An ensemble is defined as a collection of systems
-sharing one or more macroscopic characteristics but each being in a
-unique microstate. The complete ensemble is specified by giving all
-systems or microstates consistent with the common macroscopic
-characteristics of the ensemble. Although the state of each
-individual system in the ensemble could be precisely described at
-any instance in time by a suitable phase space vector, when using
-ensembles for statistical purposes, there is no need to maintain
-distinctions between individual systems, since the numbers of
-systems at any time in the different states which correspond to
-different regions of the phase space are more interesting. Moreover,
-in the point of view of statistical mechanics, one would prefer to
-use ensembles containing a large enough population of separate
-members so that the numbers of systems in such different states can
-be regarded as changing continuously as we traverse different
-regions of the phase space. The condition of an ensemble at any time
+In statistical mechanics, the condition of an ensemble at any time
 can be regarded as appropriately specified by the density $\rho$
 with which representative points are distributed over the phase
-space. The density of distribution for an ensemble with $f$ degrees
-of freedom is defined as,
+space. The density distribution for an ensemble with $f$ degrees of
+freedom is defined as,
 \begin{equation}
 \rho  = \rho (q_1 , \ldots ,q_f ,p_1 , \ldots ,p_f ,t).
 \label{introEquation:densityDistribution}
 \end{equation}
 Governed by the principles of mechanics, the phase points change
-their value which would change the density at any time at phase
-space. Hence, the density of distribution is also to be taken as a
+their locations which would change 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,
@@ -273,10 +250,10 @@ Assuming a large enough population of systems are expl
 \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 are exploited, 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,
+Assuming a large enough population of 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,
 \begin{equation}
 N = \int { \ldots \int {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f
 \label{introEquation:totalNumberSystem}
@@ -288,16 +265,16 @@ With the help of Equation(\ref{introEquation:unitProba
 {\rho (q,p,t)dq_1 } ...dq_f dp_1 } ...dp_f }}.
 \label{introEquation:unitProbability}
 \end{equation}
-With the help of Equation(\ref{introEquation:unitProbability}) and
-the knowledge of the system, it is possible to calculate the average
+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
 complex, or stochastic, or even discontinuous, the average
-properties of the ensemble of possibilities as a whole may still
-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,
+properties of the ensemble of possibilities as a whole remaining
+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
@@ -307,30 +284,29 @@ parameters, such as temperature \textit{etc}, partitio
 
 There are several different types of ensembles with different
 statistical characteristics. As a function of macroscopic
-parameters, such as temperature \textit{etc}, 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, Microcanonical ensemble(NVE) has a
-partition function like,
+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}.
+\Omega (N,V,E) = e^{\beta TS}. \label{introEquation:NVEPartition}
 \end{equation}
-A canonical ensemble(NVT)is an ensemble of systems, each of which
+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}
+\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 experiment are carried out under constant pressure
-condition, isothermal-isobaric ensemble(NPT) play a very important
-role in molecular simulation. 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
+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}
@@ -339,13 +315,13 @@ The Liouville's theorem is the foundation on which sta
 
 \subsection{\label{introSection:liouville}Liouville's theorem}
 
-The Liouville's theorem is the foundation on which statistical
-mechanics rests. It describes the time evolution of phase space
+Liouville's theorem is the foundation on which statistical mechanics
+rests. It describes the time evolution of the phase space
 distribution function. In order to calculate the rate of change of
-$\rho$, we begin from Equation(\ref{introEquation:deltaN}). If we
-consider the two faces perpendicular to the $q_1$ axis, which are
-located at $q_1$ and $q_1 + \delta q_1$, the number of phase points
-leaving the opposite face is given by the expression,
+$\rho$, we begin from Eq.~\ref{introEquation:deltaN}. If we consider
+the two faces perpendicular to the $q_1$ axis, which are located at
+$q_1$ and $q_1 + \delta q_1$, the number of phase points leaving the
+opposite face is given by the expression,
 \begin{equation}
 \left( {\rho  + \frac{{\partial \rho }}{{\partial q_1 }}\delta q_1 }
 \right)\left( {\dot q_1  + \frac{{\partial \dot q_1 }}{{\partial q_1
@@ -369,7 +345,7 @@ divining $ \delta q_1  \ldots \delta q_f \delta p_1  \
 + \frac{{\partial \dot p_i }}{{\partial p_i }}} \right)}  = 0 ,
 \end{equation}
 which cancels the first terms of the right hand side. Furthermore,
-divining $ \delta q_1  \ldots \delta q_f \delta p_1  \ldots \delta
+dividing $ \delta q_1  \ldots \delta q_f \delta p_1  \ldots \delta
 p_f $ in both sides, we can write out Liouville's theorem in a
 simple form,
 \begin{equation}
@@ -378,11 +354,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 particles in the system
-is huge, we may be able to believe the system is stationary,
+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,
 \begin{equation}
 \frac{{\partial \rho }}{{\partial t}} = 0.
 \label{introEquation:stationary}
@@ -395,14 +371,14 @@ distribution,
 \label{introEquation:densityAndHamiltonian}
 \end{equation}
 
-\subsubsection{\label{introSection:phaseSpaceConservation}Conservation of Phase Space}
+\subsubsection{\label{introSection:phaseSpaceConservation}\textbf{Conservation of Phase Space}}
 Lets consider a region in the phase space,
 \begin{equation}
 \delta v = \int { \ldots \int {dq_1 } ...dq_f dp_1 } ..dp_f .
 \end{equation}
 If this region is small enough, the density $\rho$ can be regarded
-as uniform over the whole phase space. Thus, the number of phase
-points inside this region is given by,
+as uniform over the whole integral. Thus, the number of phase points
+inside this region is given by,
 \begin{equation}
 \delta N = \rho \delta v = \rho \int { \ldots \int {dq_1 } ...dq_f
 dp_1 } ..dp_f.
@@ -414,14 +390,14 @@ With the help of stationary assumption
 \end{equation}
 With the help of stationary assumption
 (\ref{introEquation:stationary}), we obtain the principle of the
-\emph{conservation of extension in phase space},
+\emph{conservation of volume in phase space},
 \begin{equation}
 \frac{d}{{dt}}(\delta v) = \frac{d}{{dt}}\int { \ldots \int {dq_1 }
 ...dq_f dp_1 } ..dp_f  = 0.
 \label{introEquation:volumePreserving}
 \end{equation}
 
-\subsubsection{\label{introSection:liouvilleInOtherForms}Liouville's Theorem in Other Forms}
+\subsubsection{\label{introSection:liouvilleInOtherForms}\textbf{Liouville's Theorem in Other Forms}}
 
 Liouville's theorem can be expresses in a variety of different forms
 which are convenient within different contexts. For any two function
@@ -435,10 +411,10 @@ Substituting equations of motion in Hamiltonian formal
 \label{introEquation:poissonBracket}
 \end{equation}
 Substituting equations of motion in Hamiltonian formalism(
-\ref{introEquation:motionHamiltonianCoordinate} ,
-\ref{introEquation:motionHamiltonianMomentum} ) into
-(\ref{introEquation:liouvilleTheorem}), we can rewrite Liouville's
-theorem using Poisson bracket notion,
+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}
 \left( {\frac{{\partial \rho }}{{\partial t}}} \right) =  - \left\{
 {\rho ,H} \right\}.
@@ -463,13 +439,13 @@ simulation and the quality of the underlying model. Ho
 Various thermodynamic properties can be calculated from Molecular
 Dynamics simulation. By comparing experimental values with the
 calculated properties, one can determine the accuracy of the
-simulation and the quality of the underlying model. However, both of
-experiment and computer simulation are usually performed during a
+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
-many-body system in Statistical Mechanics. Fortunately, Ergodic
-Hypothesis is proposed to make a connection between time average and
-ensemble average. It states that time average and average over the
+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}.
 \begin{equation}
 \langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty }
@@ -491,61 +467,49 @@ A variety of numerical integrators were proposed to si
 choice\cite{Frenkel1996}.
 
 \section{\label{introSection:geometricIntegratos}Geometric Integrators}
-A variety of numerical integrators were proposed to simulate the
-motions. They usually begin with an initial conditionals and move
-the objects in the direction governed by the differential equations.
-However, most of them ignore the hidden physical law 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 simple example of
-symplectic integrator, continues to gain its popularity in molecular
-dynamics community. This fact can be partly explained by its
-geometric nature.
+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
+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 Manifold}
-A \emph{manifold} is an abstract mathematical space. It locally
-looks like Euclidean space, but when viewed globally, it may have
-more complicate 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 with an open cover in which
-the covering neighborhoods are all smoothly isomorphic to one
-another. In other words,it is possible to apply calculus on
-\emph{differentiable manifold}. A \emph{symplectic manifold} is
-defined as a pair $(M, \omega)$ which consisting of a
+\subsection{\label{introSection:symplecticManifold}Symplectic Manifolds}
+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,
 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$. Cross product operation in vector field is an
-example of symplectic form.
+$\omega(x, x) = 0$. The cross product operation in vector field is
+an example of symplectic form. 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.
 
-One of the motivations to study \emph{symplectic manifold} 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.
-
-Let  $(M,\omega)$ and $(N, \eta)$ be symplectic manifolds. A map
-\[
-f : M \rightarrow N
-\]
-is a \emph{symplectomorphism} if it is a \emph{diffeomorphims} and
-the \emph{pullback} of $\eta$ under f is equal to $\omega$.
-Canonical transformation is an example of symplectomorphism in
-classical mechanics.
-
 \subsection{\label{introSection:ODE}Ordinary Differential Equations}
 
-For a ordinary differential system defined as
+For an ordinary differential system defined as
 \begin{equation}
 \dot x = f(x)
 \end{equation}
-where $x = x(q,p)^T$, this system is canonical Hamiltonian, if
+where $x = x(q,p)^T$, this system is a canonical Hamiltonian, if
 \begin{equation}
 f(r) = J\nabla _x H(r).
 \end{equation}
@@ -564,10 +528,8 @@ called a \emph{Hamiltonian vector field}.
 \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}
@@ -605,19 +567,18 @@ Instead, we use a approximate map, $\psi_\tau$, which 
 \end{equation}
 
 In most cases, it is not easy to find the exact flow $\varphi_\tau$.
-Instead, we use a approximate map, $\psi_\tau$, which is usually
+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})
+\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 ODE
-and its flow play important roles in numerical studies. Many of them
-can be found in systems which occur naturally in applications.
-
+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,
 \begin{equation}
@@ -631,23 +592,19 @@ is the property must be preserved by the integrator.
 \begin{equation}
 {\varphi '}^T J \varphi ' = J \circ \varphi
 \end{equation}
-is the property must be preserved by the integrator.
-
-It is possible to construct a \emph{volume-preserving} flow for a
-source free($ \nabla \cdot f = 0 $) ODE, 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 a ODE\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} 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,
 \[
 \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
+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)$ ,
 \[
@@ -660,14 +617,11 @@ smooth function $G$ is given by,
 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,
-
 \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)). \\
 \label{introEquation:firstIntegral1}
 \end{eqnarray}
-
-
 Using poisson bracket notion, Equation
 \ref{introEquation:firstIntegral1} can be rewritten as
 \[
@@ -680,17 +634,15 @@ is a \emph{first integral}, which is due to the fact $
 \]
 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$.
-
-When designing any numerical methods, one should always try to
+0$. When designing any numerical methods, one should always try to
 preserve the structural properties of the original ODE and its flow.
 
 \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
 though this fact was not recognized when they were first
-constructed. The most famous example is leapfrog methods in
-molecular dynamics. In general, symplectic integrators can be
+constructed. The most famous example is the Verlet-leapfrog method
+in molecular dynamics. In general, symplectic integrators can be
 constructed using one of four different methods.
 \begin{enumerate}
 \item Generating functions
@@ -708,12 +660,12 @@ implementing the Runge-Kutta methods, they do not attr
 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 do not attract too much
-attention from Molecular Dynamics community. Instead, splitting have
-been widely accepted since they exploit natural decompositions of
-the system\cite{Tuckerman1992, McLachlan1998}.
+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}.
 
-\subsubsection{\label{introSection:splittingMethod}Splitting Method}
+\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,
@@ -723,9 +675,8 @@ simpler integration of the system.
 \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,
+simpler integration of the system. Suppose that a Hamiltonian system
+takes the form,
 \[
 H = H_1 + H_2.
 \]
@@ -734,7 +685,7 @@ order is then given by the Lie-Trotter formula
 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 is then given by the Lie-Trotter formula
+order expression is then given by the Lie-Trotter formula
 \begin{equation}
 \varphi _h  = \varphi _{1,h}  \circ \varphi _{2,h},
 \label{introEquation:firstOrderSplitting}
@@ -760,22 +711,22 @@ which has a local error proportional to $h^3$. Sprang 
 \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$. Sprang splitting's
-popularity in molecular simulation community attribute to its
-symmetric property,
+which has a local error proportional to $h^3$. The Sprang
+splitting's popularity in molecular simulation community attribute
+to its symmetric property,
 \begin{equation}
 \varphi _h^{ - 1} = \varphi _{ - h}.
 \label{introEquation:timeReversible}
 \end{equation}
 
-\subsubsection{\label{introSection:exampleSplittingMethod}Example of Splitting Method}
+\subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Examples of the Splitting Method}}
 The classical equation for a system consisting of interacting
 particles can be written in Hamiltonian form,
 \[
 H = T + V
 \]
 where $T$ is the kinetic energy and $V$ is the potential energy.
-Setting $H_1 = T, H_2 = V$ and applying Strang splitting, one
+Setting $H_1 = T, H_2 = V$ and applying the Strang splitting, one
 obtains the following:
 \begin{align}
 q(\Delta t) &= q(0) + \dot{q}(0)\Delta t +
@@ -802,7 +753,7 @@ q(\Delta t) &= q(0) + \Delta t\, \dot{q}\biggl (\frac{
     \label{introEquation:Lp9b}\\%
 %
 \dot{q}(\Delta t) &= \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) +
-    \frac{\Delta t}{2m}\, F[q(0)]. \label{introEquation:Lp9c}
+    \frac{\Delta t}{2m}\, F[q(t)]. \label{introEquation:Lp9c}
 \end{align}
 From the preceding splitting, one can see that the integration of
 the equations of motion would follow:
@@ -811,13 +762,13 @@ 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{r}(\Delta t)$, and use the new forces to complete the velocity move.
+\item Evaluate the forces at the new positions, $\mathbf{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}
-
-Simply switching the order of splitting and composing, a new
-integrator, the \emph{position verlet} integrator, can be generated,
+By simply switching the order of the propagators in the splitting
+and composing a new integrator, the \emph{position verlet}
+integrator, can be generated,
 \begin{align}
 \dot q(\Delta t) &= \dot q(0) + \Delta tF(q(0))\left[ {q(0) +
 \frac{{\Delta t}}{{2m}}\dot q(0)} \right], %
@@ -828,12 +779,12 @@ q(\Delta t)} \right]. %
  \label{introEquation:positionVerlet2}
 \end{align}
 
-\subsubsection{\label{introSection:errorAnalysis}Error Analysis and Higher Order Methods}
+\subsubsection{\label{introSection:errorAnalysis}\textbf{Error Analysis and Higher Order Methods}}
 
-Baker-Campbell-Hausdorff formula can be used to determine the local
-error of splitting method in terms of commutator of the
+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 associate to
+the sub-flow. 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)
@@ -847,22 +798,22 @@ Applying Baker-Campbell-Hausdorff formula\cite{Varadar
 \[
 [X,Y] = XY - YX .
 \]
-Applying Baker-Campbell-Hausdorff formula\cite{Varadarajan1974} to
-Sprang splitting, we can obtain
+Applying the Baker-Campbell-Hausdorff formula\cite{Varadarajan1974}
+to the Sprang 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 )
 \end{eqnarray*}
-Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0\], the dominant local
+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 general splitting,  of the form
+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}
-Careful choice of coefficient $a_1 \ldot b_m$ will lead to higher
-order method. Yoshida proposed an elegant way to compose higher
+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
 a symmetric second order base method $ \varphi _h^{(2)} $, a
 fourth-order symmetric method can be constructed by composing,
@@ -875,9 +826,9 @@ _{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)}
 integrator $ \varphi _h^{(2n + 2)}$ can be composed by
 \begin{equation}
 \varphi _h^{(2n + 2)}  = \varphi _{\alpha h}^{(2n)}  \circ \varphi
-_{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)}
+_{\beta h}^{(2n)}  \circ \varphi _{\alpha h}^{(2n)},
 \end{equation}
-, if the weights are chosen as
+if the weights are chosen as
 \[
 \alpha  =  - \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }},\beta =
 \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }} .
@@ -914,65 +865,68 @@ initialization of a simulation. Sec.~\ref{introSec:pro
 \end{enumerate}
 These three individual steps will be covered in the following
 sections. Sec.~\ref{introSec:initialSystemSettings} deals with the
-initialization of a simulation. Sec.~\ref{introSec:production} will
-discusses issues in production run. Sec.~\ref{introSection:Analysis}
-provides the theoretical tools for trajectory analysis.
+initialization of a simulation. Sec.~\ref{introSection:production}
+will discusse issues in production run.
+Sec.~\ref{introSection:Analysis} provides the theoretical tools for
+trajectory analysis.
 
 \subsection{\label{introSec:initialSystemSettings}Initialization}
 
-\subsubsection{Preliminary preparation}
+\subsubsection{\textbf{Preliminary preparation}}
 
 When selecting the starting structure of a molecule for molecular
 simulation, one may retrieve its Cartesian coordinates from public
 databases, such as RCSB Protein Data Bank \textit{etc}. Although
 thousands of crystal structures of molecules are discovered every
 year, many more remain unknown due to the difficulties of
-purification and crystallization. Even for the molecule with known
-structure, some important information is missing. For example, the
+purification and crystallization. Even for molecules with known
+structure, 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,
 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 setup. For instance, when studying transport phenomenon in
-membrane system, we may prepare the lipids in bilayer structure
-instead of placing lipids randomly in solvent, since we are not
-interested in self-aggregation and it takes a long time to happen.
+a special configuration. For instance, when studying transport
+phenomenon in membrane systems, we may prepare the lipids in a
+bilayer structure instead of placing lipids randomly in solvent,
+since we are not interested in the slow self-aggregation process.
 
-\subsubsection{Minimization}
+\subsubsection{\textbf{Minimization}}
 
 It is quite possible that some of molecules in the system from
-preliminary preparation may be overlapped with each other. This
-close proximity leads to high potential energy which consequently
-jeopardizes any molecular dynamics simulations. To remove these
-steric overlaps, one typically performs energy 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 systems are far from
-harmonic. Thus, it is often used to refine structure from
-crystallographic data. Relied on the gradient or hessian, advanced
-methods like conjugate gradient and Newton-Raphson converge rapidly
-to a local minimum, while become unstable if the energy surface is
-far from quadratic. Another factor must be taken into account, when
+preliminary preparation may be overlapping with each other. This
+close proximity leads to high initial potential energy which
+consequently jeopardizes any molecular dynamics simulations. To
+remove these steric overlaps, one typically performs energy
+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
+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
+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 limit of computation power to calculate hessian
-matrix and insufficient storage capacity to store them, most
-Newton-Raphson methods can not be used with very large models.
+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.
 
-\subsubsection{Heating}
+\subsubsection{\textbf{Heating}}
 
 Typically, Heating is performed by assigning random velocities
-according to a Gaussian distribution for a temperature. Beginning at
-a lower temperature and gradually increasing the temperature by
-assigning greater 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. Equivalently, the
-net linear momentum and angular momentum of the system should be
-shifted to zero.
+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.
 
-\subsubsection{Equilibration}
+\subsubsection{\textbf{Equilibration}}
 
 The purpose of equilibration is to allow the system to evolve
 spontaneously for a period of time and reach equilibrium. The
@@ -986,30 +940,32 @@ Production run is the most important step of the simul
 
 \subsection{\label{introSection:production}Production}
 
-Production run is the most important step of the simulation, in
+The production run is the most important step of the simulation, in
 which the equilibrated structure is used as a starting point and the
 motions of the molecules are collected for later analysis. In order
 to capture the macroscopic properties of the system, the molecular
-dynamics simulation must be performed in correct and efficient way.
+dynamics simulation must be performed by sampling correctly and
+efficiently from the relevant thermodynamic ensemble.
 
 The most expensive part of a molecular dynamics simulation is the
 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
-computation saving techniques.
+algorithmic tricks.
 
-A natural approach to avoid system size issue is to represent the
+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. To offset this,
-\textit{Periodic boundary condition} (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.
+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.
 \begin{figure}
 \centering
 \includegraphics[width=\linewidth]{pbc.eps}
@@ -1021,32 +977,32 @@ evaluation is to apply cutoff where particles farther 
 
 %cutoff and minimum image convention
 Another important technique to improve the efficiency of force
-evaluation is to apply cutoff where particles farther than a
-predetermined distance, are not included in the calculation
+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 the potential to ensure the potential curve go smoothly to
-zero at the cutoff radius. Cutoff strategy works pretty well for
-Lennard-Jones interaction because of its short range nature.
-However, simply truncating the electrostatic interaction with the
-use of cutoff has been shown to lead to severe artifacts in
-simulations. Ewald summation, in which the slowly conditionally
-convergent 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
+shift 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
+with the use of cutoffs has been shown to lead to severe artifacts
+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
 methods\cite{Hockney1981,Shimada1993, Luty1994} are accelerated from
-$O(N^{3/2})$ to $O(N logN)$. An alternative approach is \emph{fast
-multipole method}\cite{Greengard1987, Greengard1994}, which treats
-Coulombic interaction exactly at short range, and approximate the
-potential at long range through multipolar expansion. In spite of
-their wide acceptances at the molecular simulation community, these
-two methods are hard to be implemented 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 particle $i$ and particle $j$ at
-distance $r_{rj}$ is given by:
+$O(N^{3/2})$ to $O(N logN)$. An alternative approach is the
+\emph{fast multipole method}\cite{Greengard1987, Greengard1994},
+which treats Coulombic interactions exactly at short range, and
+approximate the potential at long range through multipolar
+expansion. In spite of their wide acceptance at the molecular
+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
+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
@@ -1068,19 +1024,18 @@ Recently, advanced visualization technique are widely 
 
 \subsection{\label{introSection:Analysis} Analysis}
 
-Recently, advanced visualization technique are widely applied to
+Recently, advanced visualization technique have become applied to
 monitor the motions of molecules. Although the dynamics of the
 system can be described qualitatively from animation, quantitative
-trajectory analysis are more appreciable. According to the
-principles of Statistical Mechanics,
-Sec.~\ref{introSection:statisticalMechanics}, one can compute
-thermodynamics properties, analyze fluctuations of structural
-parameters, and investigate time-dependent processes of the molecule
-from the trajectories.
+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.
 
-\subsubsection{\label{introSection:thermodynamicsProperties}Thermodynamics Properties}
+\subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamic Properties}}
 
-Thermodynamics properties, which can be expressed in terms of some
+Thermodynamic properties, which can be expressed in terms of some
 function of the coordinates and momenta of all particles in the
 system, can be directly computed from molecular dynamics. The usual
 way to measure the pressure is based on virial theorem of Clausius
@@ -1100,69 +1055,56 @@ P = \frac{{Nk_B T}}{V} - \frac{1}{{3V}}\left\langle {\
 < j} {r{}_{ij} \cdot f_{ij} } } \right\rangle
 \end{equation}
 
-\subsubsection{\label{introSection:structuralProperties}Structural Properties}
+\subsubsection{\label{introSection:structuralProperties}\textbf{Structural Properties}}
 
 Structural Properties of a simple fluid can be described by a set of
-distribution functions. Among these functions,\emph{pair
+distribution functions. Among these functions,the \emph{pair
 distribution function}, also known as \emph{radial distribution
-function}, is of most fundamental importance to liquid-state theory.
-Pair distribution function can be gathered by Fourier transforming
-raw data from a series of neutron diffraction experiments and
-integrating over the surface factor \cite{Powles1973}. The
-experiment result can serve as a criterion to justify the
-correctness of the theory. Moreover, various equilibrium
-thermodynamic and structural properties can also be expressed in
-terms of radial distribution function \cite{Allen1987}.
-
-A pair distribution functions $g(r)$ gives the probability that a
+function}, is of most fundamental importance to liquid theory.
+Experimentally, pair distribution function 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
 \[
 g(r) = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j
-\ne i} {\delta (r - r_{ij} )} } } \right\rangle.
+\ne i} {\delta (r - r_{ij} )} } } \right\rangle = \frac{\rho
+(r)}{\rho}.
 \]
 Note that the delta function can be replaced by a histogram in
-computer simulation. Figure
-\ref{introFigure:pairDistributionFunction} shows a typical pair
-distribution function for the liquid argon system. The occurrence of
-several peaks in the plot of $g(r)$ suggests that it is more likely
-to find particles at certain radial values than at others. This is a
-result of the attractive interaction at such distances. Because of
-the strong repulsive forces at short distance, the probability of
-locating particles at distances less than about 2.5{\AA} from each
-other is essentially zero.
+computer simulation. Peaks in $g(r)$ represent solvent shells, and
+the height of these peaks gradually decreases to 1 as the liquid of
+large distance approaches the bulk density.
 
-%\begin{figure}
-%\centering
-%\includegraphics[width=\linewidth]{pdf.eps}
-%\caption[Pair distribution function for the liquid argon
-%]{Pair distribution function for the liquid argon}
-%\label{introFigure:pairDistributionFunction}
-%\end{figure}
 
-\subsubsection{\label{introSection:timeDependentProperties}Time-dependent
-Properties}
+\subsubsection{\label{introSection:timeDependentProperties}\textbf{Time-dependent
+Properties}}
 
 Time-dependent properties are usually calculated using \emph{time
-correlation function}, which correlates random variables $A$ and $B$
-at two different time
+correlation functions}, which correlate random variables $A$ and $B$
+at two different times,
 \begin{equation}
 C_{AB} (t) = \left\langle {A(t)B(0)} \right\rangle.
 \label{introEquation:timeCorrelationFunction}
 \end{equation}
 If $A$ and $B$ refer to same variable, this kind of correlation
-function is called \emph{auto correlation function}. One example of
-auto correlation function is velocity auto-correlation function
-which is directly related to transport properties of molecular
-liquids:
+function is called an \emph{autocorrelation function}. One example
+of an auto correlation function is the velocity auto-correlation
+function which is directly related to transport properties of
+molecular liquids:
 \[
 D = \frac{1}{3}\int\limits_0^\infty  {\left\langle {v(t) \cdot v(0)}
 \right\rangle } dt
 \]
-where $D$ is diffusion constant. Unlike velocity autocorrelation
-function which is averaging over time origins and over all the
-atoms, dipole autocorrelation are calculated for the entire system.
-The dipole autocorrelation function is given by:
+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
+entire system. The dipole autocorrelation function is given by:
 \[
 c_{dipole}  = \left\langle {u_{tot} (t) \cdot u_{tot} (t)}
 \right\rangle
@@ -1188,46 +1130,46 @@ simulator is governed by the rigid body dynamics. In m
 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 the rigid body dynamics. In molecular
-simulation, rigid body is used to simplify the model in
-protein-protein docking study\cite{Gray2003}.
+simulator is governed by rigid body dynamics. In molecular
+simulations, rigid bodies are used to simplify protein-protein
+docking studies\cite{Gray2003}.
 
 It is very important to develop stable and efficient methods to
-integrate the equations of motion of orientational degrees of
-freedom. Euler angles are the nature choice to describe the
-rotational degrees of freedom. However, due to its singularity, the
-numerical integration of corresponding equations of motion is very
-inefficient and inaccurate. Although an alternative integrator using
-different sets of Euler angles can overcome this
-difficulty\cite{Barojas1973}, the computational penalty and the lost
-of angular momentum conservation still remain. A singularity free
-representation utilizing quaternions was developed by Evans in
-1977\cite{Evans1977}. Unfortunately, this approach suffer from the
-nonseparable Hamiltonian resulted from 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, SHAKE and Rattle algorithm converge very slowly
-when the number of constraint increases\cite{Ryckaert1977,
-Andersen1983}.
+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}.
 
-The break through in geometric literature suggests that, in order to
+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. Introducing conjugate momentum to
-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 quaternion representation was developed by
-Omelyan\cite{Omelyan1998}. However, both of these methods are
-iterative and inefficient. In this section, we will present a
+symplectic structure of the flow. 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
+methods are iterative and inefficient. In this section, we descibe a
 symplectic Lie-Poisson integrator for rigid body developed by
 Dullweber and his coworkers\cite{Dullweber1997} in depth.
 
-\subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
-The motion of the rigid body is Hamiltonian with the Hamiltonian
+\subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Bodies}
+The motion of a rigid body is Hamiltonian with the Hamiltonian
 function
 \begin{equation}
 H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) +
@@ -1241,45 +1183,43 @@ constrained Hamiltonian equation subjects to a holonom
 I_{ii}^{ - 1}  = \frac{1}{2}\sum\limits_{i \ne j} {J_{jj}^{ - 1} }
 \]
 where $I_{ii}$ is the diagonal element of the inertia tensor. This
-constrained Hamiltonian equation subjects to a holonomic constraint,
+constrained Hamiltonian equation is subjected to a holonomic
+constraint,
 \begin{equation}
 Q^T Q = 1, \label{introEquation:orthogonalConstraint}
 \end{equation}
-which is used to ensure rotation matrix's orthogonality.
-Differentiating \ref{introEquation:orthogonalConstraint} and using
-Equation \ref{introEquation:RBMotionMomentum}, one may obtain,
+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
 the equations of motion,
-\[
-\begin{array}{c}
- \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{array}
-\]
-
+\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{{dP}}{{dt}} & = & - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}
+\end{eqnarray}
 In general, there are two ways to satisfy the holonomic constraints.
-We can use constraint force provided by lagrange multiplier on the
-normal manifold to keep the motion on constraint space. Or we can
-simply evolve the system in constraint manifold. These two methods
-are proved to be equivalent. The holonomic constraint and equations
-of motions define a constraint manifold for rigid body
+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
 \[
 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_{\star}SO(3)$. However, it turns out that under symplectic
+$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. Introducing
+diffeomorphic. By introducing
 \[
 \tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
 \]
@@ -1289,7 +1229,6 @@ T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \t
 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
 given as
@@ -1308,26 +1247,23 @@ 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, angular momentum on body 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}
- \mathop \Pi \limits^ \bullet   = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda  \\
- \mathop Q\limits^{{\rm{   }} \bullet }  = 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,
+as well as holonomic constraints,
 \[
 \begin{array}{l}
- \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
- Q^T Q = 1 \\
+ \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,
 \begin{equation}
@@ -1342,25 +1278,24 @@ 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{equation}
-(\mathop \Pi \limits^ \bullet   - \mathop \Pi \limits^ \bullet  ^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{equation}
-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 eliminate the requirement of iterations which can
-not be avoided in other methods\cite{Kol1997, Omelyan1998}.
-
-Applying hat-map isomorphism, we obtain the equation of motion for
-angular momentum on body frame
+\begin{eqnarray}
+(\dot \Pi  - \dot \Pi ^T ){\rm{ }} &= &{\rm{ }}(\Pi  - \Pi ^T ){\rm{
+}}(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 on 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 }.
@@ -1369,16 +1304,16 @@ 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}
 
-If there is not external forces exerted on the rigid body, the only
-contribution to the rotational is from the kinetic potential (the
-first term of \ref{ introEquation:bodyAngularMotion}). The free
-rigid body is an example of Lie-Poisson system with Hamiltonian
+If there are no external forces exerted on the rigid body, the only
+contribution to the rotational motion is from the kinetic energy
+(the first term of \ref{introEquation:bodyAngularMotion}). The free
+rigid body is an example of a Lie-Poisson system with Hamiltonian
 function
 \begin{equation}
 T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
@@ -1391,13 +1326,12 @@ 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
@@ -1426,35 +1360,33 @@ tR_1 }$, we can use Cayley transformation,
 \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,
+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
-manner.
-
-In order to construct a second-order symplectic method, we split the
-angular kinetic Hamiltonian function can into five terms
+manner. In order to construct a second-order symplectic method, we
+split the angular kinetic Hamiltonian function can 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
-(\pi _1 )
-\].
-Concatenating flows corresponding to these five terms, we can obtain
-an symplectic integrator,
+(\pi _1 ).
+\]
+By concatenating the propagators corresponding to these five terms,
+we can obtain an symplectic integrator,
 \[
 \varphi _{\Delta t,T^r }  = \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 }.
 \]
-
 The non-canonical Lie-Poisson bracket ${F, G}$ of two function
 $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
@@ -1465,13 +1397,15 @@ then by the chain rule
 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 extremely efficient and stable
-which can be explained by the fact the small angle approximation is
-used and the norm of the angular momentum is conserved.
+Lie-Poisson integrator is found to be both extremely efficient and
+stable. These properties can be explained by the fact the small
+angle approximation is used and the norm of the angular momentum is
+conserved.
 
 \subsection{\label{introSection:RBHamiltonianSplitting} Hamiltonian
 Splitting for Rigid Body}
@@ -1479,12 +1413,12 @@ H = T(p,\pi ) + V(q,Q)
 The Hamiltonian of rigid body can be separated in terms of kinetic
 energy and potential energy,
 \[
-H = T(p,\pi ) + V(q,Q)
+H = T(p,\pi ) + V(q,Q).
 \]
 The equations of motion corresponding to potential energy and
 kinetic energy are listed in the below table,
 \begin{table}
-\caption{Equations of motion due to Potential and Kinetic Energies}
+\caption{EQUATIONS OF MOTION DUE TO POTENTIAL AND KINETIC ENERGIES}
 \begin{center}
 \begin{tabular}{|l|l|}
   \hline
@@ -1498,71 +1432,69 @@ A second-order symplectic method is now obtained by th
 \end{tabular}
 \end{center}
 \end{table}
-A second-order symplectic method is now obtained by the
-composition of the flow maps,
+A second-order symplectic method is now obtained by the composition
+of the position and velocity propagators,
 \[
 \varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
 _{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
 \]
 Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
-sub-flows which corresponding to force and torque respectively,
+sub-propagators which corresponding to force and torque
+respectively,
 \[
 \varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
 _{\Delta t/2,\tau }.
 \]
 Since the associated operators of $\varphi _{\Delta t/2,F} $ and
-$\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
-order inside $\varphi _{\Delta t/2,V}$ does not matter.
-
-Furthermore, kinetic potential can be separated to translational
-kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
+$\circ \varphi _{\Delta t/2,\tau }$ commute, the composition order
+inside $\varphi _{\Delta t/2,V}$ does not matter. Furthermore, the
+kinetic energy can be separated to translational kinetic term, $T^t
+(p)$, and rotational kinetic term, $T^r (\pi )$,
 \begin{equation}
 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 flow maps are given by
+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 flow maps for free moving
-rigid body
-\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}
+Finally, we obtain the overall symplectic propagators for freely
+moving rigid bodies
+\begin{eqnarray*}
+ \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}  .\\
 \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 simulation. A brief derivation of
-generalized Langevin equation will be given first. Follow 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.
+the theory of Langevin dynamics. A brief derivation of 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.
 
 \subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
 
-Harmonic bath model, in which an effective set of harmonic
+A harmonic bath model, in which an effective set of harmonic
 oscillators are used to mimic the effect of a linearly responding
 environment, has been widely used in quantum chemistry and
 statistical mechanics. One of the successful applications of
-Harmonic bath model is the derivation of Deriving Generalized
-Langevin Dynamics. Lets consider a system, in which the degree of
+Harmonic bath model is the derivation of the Generalized Langevin
+Dynamics (GLE). Lets 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}
 H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
 \label{introEquation:bathGLE}.
 \end{equation}
-Here $p$ is a momentum conjugate to $q$, $m$ is the mass associated
-with this degree of freedom, $H_B$ is harmonic bath Hamiltonian,
+Here $p$ is a momentum conjugate to $x$, $m$ is the mass associated
+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 }
@@ -1570,30 +1502,28 @@ the harmonic bath masses, and $\Delta U$ is bilinear s
 \]
 where the index $\alpha$ runs over all the bath degrees of freedom,
 $\omega _\alpha$ are the harmonic bath frequencies, $m_\alpha$ are
-the harmonic bath masses, and $\Delta U$ is bilinear system-bath
+the harmonic bath masses, and $\Delta U$ is a bilinear system-bath
 coupling,
 \[
 \Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
 \]
-where $g_\alpha$ are the coupling constants between the bath and the
-coordinate $x$. Introducing
+where $g_\alpha$ are the coupling constants between the bath
+coordinates ($x_ \alpha$) and the system coordinate ($x$).
+Introducing
 \[
 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
-Generalized Langevin Dynamics by Hamilton's equations
-\ref{introEquation:motionHamiltonianCoordinate,
-introEquation:motionHamiltonianMomentum},
+Generalized Langevin Dynamics by Hamilton's equations,
 \begin{equation}
 m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
 \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   -
@@ -1606,25 +1536,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, Laplace transform is the appropriate tool to
-solve this problem. The basic idea is to transform the difficult
+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 applying 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, 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
 \[
 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
-
 \begin{eqnarray*}
  L(x + y)  & = & L(x) + L(y) \\
  L(ax)     & = & aL(x) \\
@@ -1632,20 +1558,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 Laplace transform to the bath coordinates, we obtain
+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 }} \\
 \end{eqnarray*}
-
 By 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:
 \[
@@ -1655,27 +1578,27 @@ transformations:
  L(1) = \frac{1}{p} \\
  \end{array}
 \]
-, we obtain
-\[
-m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
+we obtain
+\begin{eqnarray*}
+m\ddot x & =  & - \frac{{\partial W(x)}}{{\partial x}} -
 \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
 }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\int_0^t {\cos (\omega
-_\alpha  t)\dot x(t - \tau )d\tau  - \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\}}
-\]
-\[
-m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
+_\alpha  t)\dot x(t - \tau )d\tau } } \right\}}  \\
+& & + \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
-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  }}} \right]\cos (\omega _\alpha  t) +
-\frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega _\alpha  }}\sin
-(\omega _\alpha  t)} \right\}}
-\]
-
+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  }}}
+\right]\cos (\omega _\alpha  t) + \frac{{g_\alpha  \dot x_\alpha
+(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)} \right\}}
+\end{eqnarray*}
 Introducing a \emph{dynamic friction kernel}
 \begin{equation}
 \xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
@@ -1698,7 +1621,7 @@ which is known as the \emph{generalized Langevin equat
 \end{equation}
 which is known as the \emph{generalized Langevin equation}.
 
-\subsubsection{\label{introSection:randomForceDynamicFrictionKernel}Random Force and Dynamic Friction Kernel}
+\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
@@ -1711,9 +1634,8 @@ as the model, which is gaussian distribution in genera
  \end{array}
 \]
 This property is what we expect from a truly random process. As long
-as the model, which is gaussian distribution in general, chosen for
-$R(t)$ is a truly random process, the stochastic nature of the GLE
-still remains.
+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
@@ -1729,15 +1651,15 @@ 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),
 \]
-which can be used to describe dynamic caging effect. The other
-extreme is the bath that responds infinitely quickly to motions in
-the system. Thus, $\xi (t)$ can be taken as a $delta$ function in
-time:
+which can be used to describe the effect of dynamic caging in
+viscous solvents. The other extreme is the bath that responds
+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)
 \]
@@ -1746,18 +1668,18 @@ and Equation \ref{introEuqation:GeneralizedLangevinDyn
 \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}
 \end{equation}
 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
+or be determined by Stokes' law for regular shaped particles. A
 briefly review on calculating friction tensor for arbitrary shaped
 particles is given in Sec.~\ref{introSection:frictionTensor}.
 
-\subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
+\subsubsection{\label{introSection:secondFluctuationDissipation}\textbf{The Second Fluctuation Dissipation Theorem}}
 
 Defining a new set of coordinates,
 \[
@@ -1769,7 +1691,6 @@ And since the $q$ coordinates are harmonic oscillators
 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) \\
@@ -1778,7 +1699,6 @@ And since the $q$ coordinates are harmonic oscillators
   & = &\sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  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
@@ -1786,266 +1706,3 @@ can model the random force and friction kernel.
 \end{equation}
 In effect, it acts as a constraint on the possible ways in which one
 can model the random force and friction kernel.
-
-\subsection{\label{introSection:frictionTensor} Friction Tensor}
-Theoretically, the friction kernel can be determined using velocity
-autocorrelation function. However, this approach become impractical
-when the system become more and more complicate. Instead, various
-approaches based on hydrodynamics have been developed to calculate
-the friction coefficients. The friction effect is isotropic in
-Equation, $\zeta$ can be taken as a scalar. In general, friction
-tensor $\Xi$ is a $6\times 6$ matrix given by
-\[
-\Xi  = \left( {\begin{array}{*{20}c}
-   {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
-   {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
-\end{array}} \right).
-\]
-Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are translational friction
-tensor and rotational resistance (friction) tensor respectively,
-while ${\Xi^{tr} }$ is translation-rotation coupling tensor and $
-{\Xi^{rt} }$ is rotation-translation coupling tensor. When a
-particle moves in a fluid, it may experience friction force or
-torque along the opposite direction of the velocity or angular
-velocity,
-\[
-\left( \begin{array}{l}
- F_R  \\
- \tau _R  \\
- \end{array} \right) =  - \left( {\begin{array}{*{20}c}
-   {\Xi ^{tt} } & {\Xi ^{rt} }  \\
-   {\Xi ^{tr} } & {\Xi ^{rr} }  \\
-\end{array}} \right)\left( \begin{array}{l}
- v \\
- w \\
- \end{array} \right)
-\]
-where $F_r$ is the friction force and $\tau _R$ is the friction
-toque.
-
-\subsubsection{\label{introSection:resistanceTensorRegular}The Resistance Tensor for Regular Shape}
-
-For a spherical particle, the translational and rotational friction
-constant can be calculated from Stoke's law,
-\[
-\Xi ^{tt}  = \left( {\begin{array}{*{20}c}
-   {6\pi \eta R} & 0 & 0  \\
-   0 & {6\pi \eta R} & 0  \\
-   0 & 0 & {6\pi \eta R}  \\
-\end{array}} \right)
-\]
-and
-\[
-\Xi ^{rr}  = \left( {\begin{array}{*{20}c}
-   {8\pi \eta R^3 } & 0 & 0  \\
-   0 & {8\pi \eta R^3 } & 0  \\
-   0 & 0 & {8\pi \eta R^3 }  \\
-\end{array}} \right)
-\]
-where $\eta$ is the viscosity of the solvent and $R$ is the
-hydrodynamics radius.
-
-Other non-spherical shape, such as cylinder and ellipsoid
-\textit{etc}, are widely used as reference for developing new
-hydrodynamics theory, because their properties can be calculated
-exactly. In 1936, Perrin extended Stokes's law to general ellipsoid,
-also called a triaxial ellipsoid, which is given in Cartesian
-coordinates by\cite{Perrin1934, Perrin1936}
-\[
-\frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
-}} = 1
-\]
-where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
-due to the complexity of the elliptic integral, only the ellipsoid
-with the restriction of two axes having to be equal, \textit{i.e.}
-prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
-exactly. Introducing an elliptic integral parameter $S$ for prolate,
-\[
-S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
-} }}{b},
-\]
-and oblate,
-\[
-S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
-}}{a}
-\],
-one can write down the translational and rotational resistance
-tensors
-\[
-\begin{array}{l}
- \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}} \\
- \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S + 2a}} \\
- \end{array},
-\]
-and
-\[
-\begin{array}{l}
- \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}} \\
- \Xi _b^{rr}  = \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}} \\
- \end{array}.
-\]
-
-\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}The Resistance Tensor for Arbitrary Shape}
-
-Unlike spherical and other regular shaped molecules, there is not
-analytical solution for friction tensor of any arbitrary shaped
-rigid molecules. The ellipsoid of revolution model and general
-triaxial ellipsoid model have been used to approximate the
-hydrodynamic properties of rigid bodies. However, since the mapping
-from all possible ellipsoidal space, $r$-space, to all possible
-combination of rotational diffusion coefficients, $D$-space is not
-unique\cite{Wegener1979} as well as the intrinsic coupling between
-translational and rotational motion of rigid body, general ellipsoid
-is not always suitable for modeling arbitrarily shaped rigid
-molecule. A number of studies have been devoted to determine the
-friction tensor for irregularly shaped rigid bodies using more
-advanced method where the molecule of interest was modeled by
-combinations of spheres(beads)\cite{Carrasco1999} and the
-hydrodynamics properties of the molecule can be calculated using the
-hydrodynamic interaction tensor. Let us consider a rigid assembly of
-$N$ beads immersed in a continuous medium. Due to hydrodynamics
-interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different
-than its unperturbed velocity $v_i$,
-\[
-v'_i  = v_i  - \sum\limits_{j \ne i} {T_{ij} F_j }
-\]
-where $F_i$ is the frictional force, and $T_{ij}$ is the
-hydrodynamic interaction tensor. The friction force of $i$th bead is
-proportional to its ``net'' velocity
-\begin{equation}
-F_i  = \zeta _i v_i  - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
-\label{introEquation:tensorExpression}
-\end{equation}
-This equation is the basis for deriving the hydrodynamic tensor. In
-1930, Oseen and Burgers gave a simple solution to Equation
-\ref{introEquation:tensorExpression}
-\begin{equation}
-T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
-R_{ij}^T }}{{R_{ij}^2 }}} \right).
-\label{introEquation:oseenTensor}
-\end{equation}
-Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$.
-A second order expression for element of different size was
-introduced by Rotne and Prager\cite{Rotne1969} and improved by
-Garc\'{i}a de la Torre and Bloomfield\cite{Torre1977},
-\begin{equation}
-T_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I +
-\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma
-_i^2  + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} -
-\frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
-\label{introEquation:RPTensorNonOverlapped}
-\end{equation}
-Both of the Equation \ref{introEquation:oseenTensor} and Equation
-\ref{introEquation:RPTensorNonOverlapped} have an assumption $R_{ij}
-\ge \sigma _i  + \sigma _j$. An alternative expression for
-overlapping beads with the same radius, $\sigma$, is given by
-\begin{equation}
-T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
-\frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I +
-\frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
-\label{introEquation:RPTensorOverlapped}
-\end{equation}
-
-To calculate the resistance tensor at an arbitrary origin $O$, we
-construct a $3N \times 3N$ matrix consisting of $N \times N$
-$B_{ij}$ blocks
-\begin{equation}
-B = \left( {\begin{array}{*{20}c}
-   {B_{11} } &  \ldots  & {B_{1N} }  \\
-    \vdots  &  \ddots  &  \vdots   \\
-   {B_{N1} } &  \cdots  & {B_{NN} }  \\
-\end{array}} \right),
-\end{equation}
-where $B_{ij}$ is given by
-\[
-B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
-)T_{ij}
-\]
-where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
-$B$, we obtain
-
-\[
-C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
-   {C_{11} } &  \ldots  & {C_{1N} }  \\
-    \vdots  &  \ddots  &  \vdots   \\
-   {C_{N1} } &  \cdots  & {C_{NN} }  \\
-\end{array}} \right)
-\]
-, which can be partitioned into $N \times N$ $3 \times 3$ block
-$C_{ij}$. With the help of $C_{ij}$ and skew matrix $U_i$
-\[
-U_i  = \left( {\begin{array}{*{20}c}
-   0 & { - z_i } & {y_i }  \\
-   {z_i } & 0 & { - x_i }  \\
-   { - y_i } & {x_i } & 0  \\
-\end{array}} \right)
-\]
-where $x_i$, $y_i$, $z_i$ are the components of the vector joining
-bead $i$ and origin $O$. Hence, the elements of resistance tensor at
-arbitrary origin $O$ can be written as
-\begin{equation}
-\begin{array}{l}
- \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\
- \Xi _{}^{tr}  = \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
- \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j  \\
- \end{array}
-\label{introEquation:ResistanceTensorArbitraryOrigin}
-\end{equation}
-
-The resistance tensor depends on the origin to which they refer. The
-proper location for applying friction force is the center of
-resistance (reaction), at which the trace of rotational resistance
-tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
-resistance is defined as an unique point of the rigid body at which
-the translation-rotation coupling tensor are symmetric,
-\begin{equation}
-\Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
-\label{introEquation:definitionCR}
-\end{equation}
-Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
-we can easily find out that the translational resistance tensor is
-origin independent, while the rotational resistance tensor and
-translation-rotation coupling resistance tensor depend on the
-origin. Given resistance tensor at an arbitrary origin $O$, and a
-vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
-obtain the resistance tensor at $P$ by
-\begin{equation}
-\begin{array}{l}
- \Xi _P^{tt}  = \Xi _O^{tt}  \\
- \Xi _P^{tr}  = \Xi _P^{rt}  = \Xi _O^{tr}  - U_{OP} \Xi _O^{tt}  \\
- \Xi _P^{rr}  = \Xi _O^{rr}  - U_{OP} \Xi _O^{tt} U_{OP}  + \Xi _O^{tr} U_{OP}  - U_{OP} \Xi _O^{tr} ^{^T }  \\
- \end{array}
- \label{introEquation:resistanceTensorTransformation}
-\end{equation}
-where
-\[
-U_{OP}  = \left( {\begin{array}{*{20}c}
-   0 & { - z_{OP} } & {y_{OP} }  \\
-   {z_i } & 0 & { - x_{OP} }  \\
-   { - y_{OP} } & {x_{OP} } & 0  \\
-\end{array}} \right)
-\]
-Using Equations \ref{introEquation:definitionCR} and
-\ref{introEquation:resistanceTensorTransformation}, one can locate
-the position of center of resistance,
-\begin{eqnarray*}
- \left( \begin{array}{l}
- x_{OR}  \\
- y_{OR}  \\
- z_{OR}  \\
- \end{array} \right) & = &\left( {\begin{array}{*{20}c}
-   {(\Xi _O^{rr} )_{yy}  + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} }  \\
-   { - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz}  + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} }  \\
-   { - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx}  + (\Xi _O^{rr} )_{yy} }  \\
-\end{array}} \right)^{ - 1}  \\
-  & & \left( \begin{array}{l}
- (\Xi _O^{tr} )_{yz}  - (\Xi _O^{tr} )_{zy}  \\
- (\Xi _O^{tr} )_{zx}  - (\Xi _O^{tr} )_{xz}  \\
- (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
- \end{array} \right) \\
-\end{eqnarray*}
-
-
-
-where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
-joining center of resistance $R$ and origin $O$.
