--- trunk/tengDissertation/Introduction.tex	2006/06/06 14:56:36	2801
+++ trunk/tengDissertation/Introduction.tex	2006/06/07 21:03:46	2819
@@ -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,19 +17,19 @@ 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}
  \label{introEquation:newtonThirdLaw}
@@ -46,7 +46,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,7 +59,7 @@ 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
@@ -68,42 +68,38 @@ scheme for rigid body \cite{Dullweber1997}.
 \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}.
+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 two important limitations: motions
+can only be described in cartesian coordinate systems. Moreover, It
+become impossible to predict analytically the properties of the
+system even if we know all of the details of the interaction. In
+order to overcome some of the practical difficulties which arise in
+attempts to apply Newton's equation to complex system, approximate
+numerical procedures may be developed.
 
-\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,
+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$ \cite{Tolman1979}.
+the kinetic, $K$, and potential energies, $U$.
 \begin{equation}
 \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 ) ,
 \label{introEquation:lagrangianDef}
@@ -114,8 +110,8 @@ then Eq.~\ref{introEquation:halmitonianPrinciple1} bec
 \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
@@ -132,8 +128,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}
@@ -172,11 +167,11 @@ find
 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
@@ -193,14 +188,13 @@ function of the generalized velocities $\dot q_i$ and 
 
 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}.
@@ -230,12 +224,12 @@ 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
-momenta is a phase space vector.
+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 = (q_1 , \ldots
+,q_f ,p_1 , \ldots ,p_f )$, with a unique set of values of $6f$
+coordinates and momenta is a phase space vector.
 
 A microscopic state or microstate of a classical system is
 specification of the complete phase space vector of a system at any
@@ -257,15 +251,15 @@ space. The density of distribution for an ensemble wit
 regions of the phase space. 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 +267,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}
@@ -293,11 +287,11 @@ properties of the ensemble of possibilities as a whole
 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,12 +301,12 @@ 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
+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, Microcanonical ensemble(NVE) has a
+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}.
@@ -326,11 +320,12 @@ TS$. Since most experiment are carried out under const
 \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,8 +334,8 @@ 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
@@ -369,7 +364,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}
@@ -395,14 +390,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 +409,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
@@ -463,13 +458,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 +486,50 @@ 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 manifold} in
+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.
 
-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}
@@ -689,8 +673,8 @@ constructed. The most famous example is leapfrog metho
 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 methods
+in molecular dynamics. In general, symplectic integrators can be
 constructed using one of four different methods.
 \begin{enumerate}
 \item Generating functions
@@ -708,12 +692,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,
@@ -734,7 +718,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,15 +744,15 @@ 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{Example of Splitting Method}}
 The classical equation for a system consisting of interacting
 particles can be written in Hamiltonian form,
 \[
@@ -828,7 +812,7 @@ 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
@@ -921,7 +905,7 @@ 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
@@ -939,7 +923,7 @@ interested in self-aggregation and it takes a long tim
 instead of placing lipids randomly in solvent, since we are not
 interested in self-aggregation and it takes a long time to happen.
 
-\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
@@ -961,7 +945,7 @@ Newton-Raphson methods can not be used with very large
 matrix and insufficient storage capacity to store them, most
 Newton-Raphson methods can not be used with very large models.
 
-\subsubsection{Heating}
+\subsubsection{\textbf{Heating}}
 
 Typically, Heating is performed by assigning random velocities
 according to a Gaussian distribution for a temperature. Beginning at
@@ -973,7 +957,7 @@ shifted to zero.
 net linear momentum and angular momentum of the system should be
 shifted to zero.
 
-\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
@@ -1079,7 +1063,7 @@ from the trajectories.
 parameters, and investigate time-dependent processes of the molecule
 from the trajectories.
 
-\subsubsection{\label{introSection:thermodynamicsProperties}Thermodynamics Properties}
+\subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamics Properties}}
 
 Thermodynamics properties, which can be expressed in terms of some
 function of the coordinates and momenta of all particles in the
@@ -1101,7 +1085,7 @@ 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
@@ -1141,8 +1125,8 @@ other is essentially zero.
 %\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$
@@ -1696,7 +1680,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
@@ -1755,7 +1739,7 @@ particles is given in Sec.~\ref{introSection:frictionT
 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,
 \[
@@ -1821,7 +1805,7 @@ toque.
 where $F_r$ is the friction force and $\tau _R$ is the friction
 toque.
 
-\subsubsection{\label{introSection:resistanceTensorRegular}The Resistance Tensor for Regular Shape}
+\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}}
 
 For a spherical particle, the translational and rotational friction
 constant can be calculated from Stoke's law,
@@ -1883,7 +1867,7 @@ and
  \end{array}.
 \]
 
-\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}The Resistance Tensor for Arbitrary Shape}
+\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{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
