--- trunk/tengDissertation/Introduction.tex	2006/04/13 04:47:47	2706
+++ trunk/tengDissertation/Introduction.tex	2006/06/05 21:24:52	2793
@@ -93,7 +93,7 @@ the kinetic, $K$, and potential energies, $U$ \cite{to
 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{tolman79}.
+the kinetic, $K$, and potential energies, $U$ \cite{Tolman1979}.
 \begin{equation}
 \delta \int_{t_1 }^{t_2 } {(K - U)dt = 0} ,
 \label{introEquation:halmitonianPrinciple1}
@@ -189,7 +189,7 @@ known as the canonical equations of motions \cite{Gold
 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{Goldstein01}.
+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
@@ -200,7 +200,7 @@ equations\cite{Marion90}.
 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{Marion90}.
+equations\cite{Marion1990}.
 
 In Newtonian Mechanics, a system described by conservative forces
 conserves the total energy \ref{introEquation:energyConservation}.
@@ -470,7 +470,7 @@ statistical ensemble are identical \cite{Frenkel1996, 
 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
-statistical ensemble are identical \cite{Frenkel1996, leach01:mm}.
+statistical ensemble are identical \cite{Frenkel1996, Leach2001}.
 \begin{equation}
 \langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty }
 \frac{1}{t}\int\limits_0^t {A(q(t),p(t))dt = \int\limits_\Gamma
@@ -484,7 +484,7 @@ reasonable, the Monte Carlo techniques\cite{metropolis
 a properly weighted statistical average. This allows the researcher
 freedom of choice when deciding how best to measure a given
 observable. In case an ensemble averaged approach sounds most
-reasonable, the Monte Carlo techniques\cite{metropolis:1949} can be
+reasonable, the Monte Carlo techniques\cite{Metropolis1949} can be
 utilized. Or if the system lends itself to a time averaging
 approach, the Molecular Dynamics techniques in
 Sec.~\ref{introSection:molecularDynamics} will be the best
@@ -498,10 +498,11 @@ issue. The velocity verlet method, which happens to be
 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. 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.
+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.
 
 \subsection{\label{introSection:symplecticManifold}Symplectic Manifold}
 A \emph{manifold} is an abstract mathematical space. It locally
@@ -565,27 +566,13 @@ Another generalization of Hamiltonian dynamics is Pois
 \end{equation}In this case, $f$ is
 called a \emph{Hamiltonian vector field}.
 
-Another generalization of Hamiltonian dynamics is Poisson Dynamics,
+Another generalization of Hamiltonian dynamics is Poisson
+Dynamics\cite{Olver1986},
 \begin{equation}
 \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
 \end{equation}
 The most obvious change being that matrix $J$ now depends on $x$.
-The free rigid body is an example of Poisson system (actually a
-Lie-Poisson system) with Hamiltonian function of angular kinetic
-energy.
-\begin{equation}
-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{equation}
 
-\begin{equation}
-H = \frac{1}{2}\left( {\frac{{\pi _1^2 }}{{I_1 }} + \frac{{\pi _2^2
-}}{{I_2 }} + \frac{{\pi _3^2 }}{{I_3 }}} \right)
-\end{equation}
-
 \subsection{\label{introSection:exactFlow}Exact Flow}
 
 Let $x(t)$ be the exact solution of the ODE system,
@@ -627,9 +614,9 @@ The hidden geometric properties of ODE and its flow pl
 
 \subsection{\label{introSection:geometricProperties}Geometric Properties}
 
-The hidden geometric properties 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 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,
@@ -673,13 +660,14 @@ 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{equation}
-\begin{array}{c}
- \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)). \\
- \end{array}
+
+\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{equation}
+\end{eqnarray}
+
+
 Using poisson bracket notion, Equation
 \ref{introEquation:firstIntegral1} can be rewritten as
 \[
@@ -694,8 +682,7 @@ is a \emph{first integral}, which is due to the fact $
 is a \emph{first integral}, which is due to the fact $\{ H,H\}  =
 0$.
 
-
- When designing any numerical methods, one should always try to
+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}
@@ -712,18 +699,19 @@ Generating function tends to lead to methods which are
 \item Splitting methods
 \end{enumerate}
 
-Generating function tends to lead to methods which are cumbersome
-and difficult to use. In dissipative systems, variational methods
-can capture the decay of energy accurately. Since their
-geometrically unstable nature against non-Hamiltonian perturbations,
-ordinary implicit Runge-Kutta methods are not suitable for
-Hamiltonian system. Recently, various high-order explicit
-Runge--Kutta methods have been developed to overcome this
+Generating function\cite{Channell1990} tends to lead to methods
+which are cumbersome and difficult to use. In dissipative systems,
+variational methods can capture the decay of energy
+accurately\cite{Kane2000}. Since their geometrically unstable nature
+against non-Hamiltonian perturbations, ordinary implicit Runge-Kutta
+methods are not suitable for Hamiltonian system. Recently, various
+high-order explicit Runge-Kutta methods
+\cite{Owren1992,Chen2003}have been developed to overcome this
 instability. However, due to computational penalty involved in
 implementing the Runge-Kutta methods, they 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{Tuckerman92}.
+the system\cite{Tuckerman1992, McLachlan1998}.
 
 \subsubsection{\label{introSection:splittingMethod}Splitting Method}
 
@@ -837,7 +825,7 @@ q(\Delta t)} \right]. %
 %
 q(\Delta t) &= q(0) + \frac{{\Delta t}}{2}\left[ {\dot q(0) + \dot
 q(\Delta t)} \right]. %
- \label{introEquation:positionVerlet1}
+ \label{introEquation:positionVerlet2}
 \end{align}
 
 \subsubsection{\label{introSection:errorAnalysis}Error Analysis and Higher Order Methods}
@@ -846,7 +834,7 @@ $\varphi_1(t)$ and $\varphi_2(t$ respectively , we hav
 error of splitting method in terms of commutator of the
 operators(\ref{introEquation:exponentialOperator}) associated with
 the sub-flow. For operators $hX$ and $hY$ which are associate to
-$\varphi_1(t)$ and $\varphi_2(t$ respectively , we have
+$\varphi_1(t)$ and $\varphi_2(t)$ respectively , we have
 \begin{equation}
 \exp (hX + hY) = \exp (hZ)
 \end{equation}
@@ -859,13 +847,12 @@ Applying Baker-Campbell-Hausdorff formula to Sprang sp
 \[
 [X,Y] = XY - YX .
 \]
-Applying Baker-Campbell-Hausdorff formula to Sprang splitting, we
-can obtain
+Applying Baker-Campbell-Hausdorff formula\cite{Varadarajan1974} to
+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 & & \mbox{} +
-\ldots )
+\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
 error of Spring splitting is proportional to $h^3$. The same
@@ -874,11 +861,11 @@ Careful choice of coefficient $a_1 ,\ldot , b_m$ will 
 \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
+Careful choice of coefficient $a_1 \ldot b_m$ will lead to higher
 order method. Yoshida proposed an elegant way to compose higher
-order methods based on symmetric splitting. Given a symmetric second
-order base method $ \varphi _h^{(2)} $, a fourth-order symmetric
-method can be constructed by composing,
+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,
 \[
 \varphi _h^{(4)}  = \varphi _{\alpha h}^{(2)}  \circ \varphi _{\beta
 h}^{(2)}  \circ \varphi _{\alpha h}^{(2)}
@@ -898,17 +885,303 @@ As a special discipline of molecular modeling, Molecul
 
 \section{\label{introSection:molecularDynamics}Molecular Dynamics}
 
-As a special discipline of molecular modeling, Molecular dynamics
-has proven to be a powerful tool for studying the functions of
-biological systems, providing structural, thermodynamic and
-dynamical information.
+As one of the principal tools of molecular modeling, Molecular
+dynamics has proven to be a powerful tool for studying the functions
+of biological systems, providing structural, thermodynamic and
+dynamical information. The basic idea of molecular dynamics is that
+macroscopic properties are related to microscopic behavior and
+microscopic behavior can be calculated from the trajectories in
+simulations. For instance, instantaneous temperature of an
+Hamiltonian system of $N$ particle can be measured by
+\[
+T = \sum\limits_{i = 1}^N {\frac{{m_i v_i^2 }}{{fk_B }}}
+\]
+where $m_i$ and $v_i$ are the mass and velocity of $i$th particle
+respectively, $f$ is the number of degrees of freedom, and $k_B$ is
+the boltzman constant.
 
-\subsection{\label{introSec:mdInit}Initialization}
+A typical molecular dynamics run consists of three essential steps:
+\begin{enumerate}
+  \item Initialization
+    \begin{enumerate}
+    \item Preliminary preparation
+    \item Minimization
+    \item Heating
+    \item Equilibration
+    \end{enumerate}
+  \item Production
+  \item Analysis
+\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.
 
-\subsection{\label{introSec:forceEvaluation}Force Evaluation}
+\subsection{\label{introSec:initialSystemSettings}Initialization}
 
-\subsection{\label{introSection:mdIntegration} Integration of the Equations of Motion}
+\subsubsection{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
+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.
+
+\subsubsection{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
+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.
+
+\subsubsection{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.
 
+\subsubsection{Equilibration}
+
+The purpose of equilibration is to allow the system to evolve
+spontaneously for a period of time and reach equilibrium. The
+procedure is continued until various statistical properties, such as
+temperature, pressure, energy, volume and other structural
+properties \textit{etc}, become independent of time. Strictly
+speaking, minimization and heating are not necessary, provided the
+equilibration process is long enough. However, these steps can serve
+as a means to arrive at an equilibrated structure in an effective
+way.
+
+\subsection{\label{introSection:production}Production}
+
+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.
+
+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.
+
+A natural approach to avoid system size issue 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.
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{pbc.eps}
+\caption[An illustration of periodic boundary conditions]{A 2-D
+illustration of periodic boundary conditions. As one particle leaves
+the left of the simulation box, an image of it enters the right.}
+\label{introFig:pbc}
+\end{figure}
+
+%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
+\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
+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:
+\begin{equation}
+V(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha
+r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow
+R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha
+r_{ij})}{r_{ij}}\right\}. \label{introEquation:shiftedCoulomb}
+\end{equation}
+where $\alpha$ is the convergence parameter. Due to the lack of
+inherent periodicity and rapid convergence,this method is extremely
+efficient and easy to implement.
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{shifted_coulomb.eps}
+\caption[An illustration of shifted Coulomb potential]{An
+illustration of shifted Coulomb potential.}
+\label{introFigure:shiftedCoulomb}
+\end{figure}
+
+%multiple time step
+
+\subsection{\label{introSection:Analysis} Analysis}
+
+Recently, advanced visualization technique are widely 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.
+
+\subsubsection{\label{introSection:thermodynamicsProperties}Thermodynamics Properties}
+
+Thermodynamics 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
+which states that the virial is equal to $-3Nk_BT$. For a system
+with forces between particles, the total virial, $W$, contains the
+contribution from external pressure and interaction between the
+particles:
+\[
+W =  - 3PV + \left\langle {\sum\limits_{i < j} {r{}_{ij} \cdot
+f_{ij} } } \right\rangle
+\]
+where $f_{ij}$ is the force between particle $i$ and $j$ at a
+distance $r_{ij}$. Thus, the expression for the pressure is given
+by:
+\begin{equation}
+P = \frac{{Nk_B T}}{V} - \frac{1}{{3V}}\left\langle {\sum\limits_{i
+< j} {r{}_{ij} \cdot f_{ij} } } \right\rangle
+\end{equation}
+
+\subsubsection{\label{introSection:structuralProperties}Structural Properties}
+
+Structural Properties of a simple fluid can be described by a set of
+distribution functions. Among these functions,\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
+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.
+\]
+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.
+
+%\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}
+
+Time-dependent properties are usually calculated using \emph{time
+correlation function}, which correlates random variables $A$ and $B$
+at two different time
+\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:
+\[
+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:
+\[
+c_{dipole}  = \left\langle {u_{tot} (t) \cdot u_{tot} (t)}
+\right\rangle
+\]
+Here $u_{tot}$ is the net dipole of the entire system and is given
+by
+\[
+u_{tot} (t) = \sum\limits_i {u_i (t)}
+\]
+In principle, many time correlation functions can be related with
+Fourier transforms of the infrared, Raman, and inelastic neutron
+scattering spectra of molecular liquids. In practice, one can
+extract the IR spectrum from the intensity of dipole fluctuation at
+each frequency using the following relationship:
+\[
+\hat c_{dipole} (v) = \int_{ - \infty }^\infty  {c_{dipole} (t)e^{ -
+i2\pi vt} dt}
+\]
+
 \section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
 
 Rigid bodies are frequently involved in the modeling of different
@@ -917,7 +1190,7 @@ protein-protein docking study{\cite{Gray03}}.
 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{Gray03}}.
+protein-protein docking study\cite{Gray2003}.
 
 It is very important to develop stable and efficient methods to
 integrate the equations of motion of orientational degrees of
@@ -925,36 +1198,37 @@ different sets of Euler angles can overcome this diffi
 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{},
-the computational penalty and the lost of angular momentum
-conservation still remain. A singularity free representation
-utilizing quaternions was developed by Evans in 1977. 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.
+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}.
 
 The 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 $A$ and re-formulating Hamiltonian's equation, a
-symplectic integrator, RSHAKE, was proposed to evolve the
-Hamiltonian system in a constraint manifold by iteratively
-satisfying the orthogonality constraint $A_t A = 1$. An alternative
-method using quaternion representation was developed by Omelyan.
-However, both of these methods are iterative and inefficient. In
-this section, we will present a symplectic Lie-Poisson integrator
-for rigid body developed by Dullweber and his coworkers\cite{}.
+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 Lie-Poisson integrator for rigid body developed by
+Dullweber and his coworkers\cite{Dullweber1997} in depth.
 
-\subsection{\label{introSection:lieAlgebra}Lie Algebra}
-
 \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
-
+The motion of the 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) +
 V(q,Q) + \frac{1}{2}tr[(QQ^T  - 1)\Lambda ].
@@ -969,13 +1243,13 @@ Q^T Q = 1$, \label{introEquation:orthogonalConstraint}
 where $I_{ii}$ is the diagonal element of the inertia tensor. This
 constrained Hamiltonian equation subjects to a holonomic constraint,
 \begin{equation}
-Q^T Q = 1$, \label{introEquation:orthogonalConstraint}
+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,
 \begin{equation}
-Q^t PJ^{ - 1}  + J^{ - 1} P^t Q = 0 . \\
+Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
 \label{introEquation:RBFirstOrderConstraint}
 \end{equation}
 
@@ -987,152 +1261,412 @@ the equations of motion,
  \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}\\
+ \frac{{dP}}{{dt}} =  - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}\\
  \end{array}
 \]
 
-
+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
 \[
-M = \left\{ {(Q,P):Q^T Q = 1,Q^t PJ^{ - 1}  + J^{ - 1} P^t Q = 0}
-\right\} .
+M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0}
+\right\}.
 \]
 
-\subsection{\label{introSection:DLMMotionEquation}The Euler Equations of Rigid Body Motion}
-
-\subsection{\label{introSection:symplecticDiscretizationRB}Symplectic Discretization of Euler Equations}
-
-
-\section{\label{introSection:langevinDynamics}Langevin Dynamics}
-
-\subsection{\label{introSection:LDIntroduction}Introduction and application of Langevin Dynamics}
-
-\subsection{\label{introSection:generalizedLangevinDynamics}Generalized Langevin Dynamics}
-
-\begin{equation}
-H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
-\label{introEquation:bathGLE}
-\end{equation}
-where $H_B$ is harmonic bath Hamiltonian,
+Unfortunately, this constraint manifold is not the cotangent bundle
+$T_{\star}SO(3)$. However, it turns out that under symplectic
+transformation, the cotangent space and the phase space are
+diffeomorphic. Introducing
 \[
-H_B  =\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
-}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  w_\alpha ^2 } \right\}}
+\tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
 \]
-and $\Delta U$ is bilinear system-bath coupling,
+the mechanical system subject to a holonomic constraint manifold $M$
+can be re-formulated as a Hamiltonian system on the cotangent space
 \[
-\Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
+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\}
 \]
-Completing the square,
+
+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
+\begin{equation}
+X_i^{lab} = Q X_i + q.
+\end{equation}
+Therefore, potential energy $V(q,Q)$ is defined by
 \[
-H_B  + \Delta U = \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\}}  - \sum\limits_{\alpha  =
-1}^N {\frac{{g_\alpha ^2 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
+V(q,Q) = V(Q X_0 + q).
 \]
-and putting it back into Eq.~\ref{introEquation:bathGLE},
+Hence, the force and torque are given by
 \[
-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\}}
+\nabla _q V(q,Q) = F(q,Q) = \sum\limits_i {F_i (q,Q)},
 \]
-where
+and
 \[
-W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
-}}{{2m_\alpha  w_\alpha ^2 }}} x^2
+\nabla _Q V(q,Q) = F(q,Q)X_i^t
 \]
-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},
-\begin{align}
-\dot p &=  - \frac{{\partial H}}{{\partial x}}
-       &= m\ddot x
-       &= - \frac{{\partial W(x)}}{{\partial x}} - \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)}
-\label{introEquation:Lp5}
-\end{align}
-, and
-\begin{align}
-\dot p_\alpha   &=  - \frac{{\partial H}}{{\partial x_\alpha  }}
-                &= m\ddot x_\alpha
-                &= \- m_\alpha  w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha}}{{m_\alpha  w_\alpha ^2 }}x} \right)
-\end{align}
+respectively.
 
-\subsection{\label{introSection:laplaceTransform}The Laplace Transform}
-
+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,
+\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}  \\
+ \end{array}
+ \label{introEqaution:RBMotionPI}
+\end{equation}
+, as well as holonomic constraints,
 \[
-L(x) = \int_0^\infty  {x(t)e^{ - pt} dt}
+\begin{array}{l}
+ \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
+ Q^T Q = 1 \\
+ \end{array}
 \]
 
+For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a matrix $\hat v \in
+so(3)^ \star$, the hat-map isomorphism,
+\begin{equation}
+v(v_1 ,v_2 ,v_3 ) \Leftrightarrow \hat v = \left(
+{\begin{array}{*{20}c}
+   0 & { - v_3 } & {v_2 }  \\
+   {v_3 } & 0 & { - v_1 }  \\
+   { - v_2 } & {v_1 } & 0  \\
+\end{array}} \right),
+\label{introEquation:hatmapIsomorphism}
+\end{equation}
+will let us associate the matrix products with traditional vector
+operations
 \[
-L(x + y) = L(x) + L(y)
+\hat vu = v \times u
 \]
 
-\[
-L(ax) = aL(x)
-\]
+Using \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{equation}
+\dot \pi  = \pi  \times I^{ - 1} \pi  + \sum\limits_i {\left( {Q^T
+F_i (r,Q)} \right) \times X_i }.
+\label{introEquation:bodyAngularMotion}
+\end{equation}
+In the same manner, the equation of motion for rotation matrix is
+given by
 \[
-L(\dot x) = pL(x) - px(0)
+\dot Q = Qskew(I^{ - 1} \pi )
 \]
 
-\[
-L(\ddot x) = p^2 L(x) - px(0) - \dot x(0)
-\]
+\subsection{\label{introSection:SymplecticFreeRB}Symplectic
+Lie-Poisson Integrator for Free Rigid Body}
 
-\[
-L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p)
-\]
+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
+function
+\begin{equation}
+T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
+\label{introEquation:rotationalKineticRB}
+\end{equation}
+where $T_i^r (\pi _i ) = \frac{{\pi _i ^2 }}{{2I_i }}$ and
+Lie-Poisson structure matrix,
+\begin{equation}
+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{equation}
+Thus, the dynamics of free rigid body is governed by
+\begin{equation}
+\frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
+\end{equation}
 
-Some relatively important transformation,
+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
+\begin{equation}
+\frac{d}{{dt}}\pi  = R_1 \pi ,\frac{d}{{dt}}Q = QR_1
+\label{introEqaution:RBMotionSingleTerm}
+\end{equation}
+where
+\[ R_1  = \left( {\begin{array}{*{20}c}
+   0 & 0 & 0  \\
+   0 & 0 & {\pi _1 }  \\
+   0 & { - \pi _1 } & 0  \\
+\end{array}} \right).
+\]
+The solutions of Equation \ref{introEqaution:RBMotionSingleTerm} is
 \[
-L(\cos at) = \frac{p}{{p^2  + a^2 }}
+\pi (\Delta t) = e^{\Delta tR_1 } \pi (0),Q(\Delta t) =
+Q(0)e^{\Delta tR_1 }
 \]
-
+with
 \[
-L(\sin at) = \frac{a}{{p^2  + a^2 }}
+e^{\Delta tR_1 }  = \left( {\begin{array}{*{20}c}
+   0 & 0 & 0  \\
+   0 & {\cos \theta _1 } & {\sin \theta _1 }  \\
+   0 & { - \sin \theta _1 } & {\cos \theta _1 }  \\
+\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,
 \[
-L(1) = \frac{1}{p}
+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.
 
-First, the bath coordinates,
+In order to construct a second-order symplectic method, we split the
+angular kinetic Hamiltonian function can into five terms
 \[
-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)
+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,
+\[
+\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
 \[
-L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) +
-px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}
+\{ F,G\} (\pi ) = [\nabla _\pi  F(\pi )]^T J(\pi )\nabla _\pi  G(\pi
+)
 \]
-Then, the system coordinates,
-\begin{align}
-mL(\ddot x) &=  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
-\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{\frac{{g_\alpha
-}}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha
-(0)}}{{p^2  + \omega _\alpha ^2 }} - \frac{{g_\alpha ^2 }}{{m_\alpha
-}}\omega _\alpha ^2 L(x)} \right\}}
-%
- &= - \frac{1}{p}\frac{{\partial W(x)}}{{\partial 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\}}
-\end{align}
-Then, the inverse transform,
+If the Poisson bracket of a function $F$ with an arbitrary smooth
+function $G$ is zero, $F$ is a \emph{Casimir}, which is the
+conserved quantity in Poisson system. We can easily verify that the
+norm of the angular momentum, $\parallel \pi
+\parallel$, is a \emph{Casimir}. Let$ F(\pi ) = S(\frac{{\parallel
+\pi \parallel ^2 }}{2})$ for an arbitrary function $ S:R \to R$ ,
+then by the chain rule
+\[
+\nabla _\pi  F(\pi ) = S'(\frac{{\parallel \pi \parallel ^2
+}}{2})\pi
+\]
+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.
 
-\begin{align}
-m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
+\subsection{\label{introSection:RBHamiltonianSplitting} Hamiltonian
+Splitting for Rigid Body}
+
+The Hamiltonian of rigid body can be separated in terms of kinetic
+energy and potential energy,
+\[
+H = T(p,\pi ) + V(q,Q)
+\]
+The equations of motion corresponding to potential energy and
+kinetic energy are listed in the below table,
+\begin{table}
+\caption{Equations of motion due to Potential and Kinetic Energies}
+\begin{center}
+\begin{tabular}{|l|l|}
+  \hline
+  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
+  Potential & Kinetic \\
+  $\frac{{dq}}{{dt}} = \frac{p}{m}$ & $\frac{d}{{dt}}q = p$ \\
+  $\frac{d}{{dt}}p =  - \frac{{\partial V}}{{\partial q}}$ & $ \frac{d}{{dt}}p = 0$ \\
+  $\frac{d}{{dt}}Q = 0$ & $ \frac{d}{{dt}}Q = Qskew(I^{ - 1} j)$ \\
+  $ \frac{d}{{dt}}\pi  = \sum\limits_i {\left( {Q^T F_i (r,Q)} \right) \times X_i }$ & $\frac{d}{{dt}}\pi  = \pi  \times I^{ - 1} \pi$\\
+  \hline
+\end{tabular}
+\end{center}
+\end{table}
+A second-order symplectic method is now obtained by the
+composition of the flow maps,
+\[
+\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,
+\[
+\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 )$,
+\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
+\[
+\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}
+\label{introEquation:overallRBFlowMaps}
+\end{equation}
+
+\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.
+
+\subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
+
+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
+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,
+\[
+H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
+}}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  \omega _\alpha ^2 }
+\right\}}
+\]
+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
+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
+\[
+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
+\[
+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\}}
+\]
+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},
+\begin{equation}
+m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
+\sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   -
+\frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)},
+\label{introEquation:coorMotionGLE}
+\end{equation}
+and
+\begin{equation}
+m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x_\alpha   -
+\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 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
+\[
+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) \\
+ L(\dot x) & = & pL(x) - px(0) \\
+ 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
+\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\}}  \\
+\end{eqnarray*}
+
+With the help of some relatively important inverse Laplace
+transformations:
+\[
+\begin{array}{c}
+ L(\cos at) = \frac{p}{{p^2  + a^2 }} \\
+ L(\sin at) = \frac{a}{{p^2  + a^2 }} \\
+ L(1) = \frac{1}{p} \\
+ \end{array}
+\]
+, we obtain
+\[
+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\}}
-%
-&= - \frac{{\partial W(x)}}{{\partial x}} - \int_0^t
+\]
+\[
+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\{
@@ -1140,74 +1674,378 @@ t)\dot x(t - \tau )d} \tau }  + \sum\limits_{\alpha  =
 \omega _\alpha  }}} \right]\cos (\omega _\alpha  t) +
 \frac{{g_\alpha  \dot x_\alpha  (0)}}{{\omega _\alpha  }}\sin
 (\omega _\alpha  t)} \right\}}
-\end{align}
+\]
 
+Introducing a \emph{dynamic friction kernel}
 \begin{equation}
+\xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
+}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
+\label{introEquation:dynamicFrictionKernelDefinition}
+\end{equation}
+and \emph{a random force}
+\begin{equation}
+R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
+- \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
+\right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
+(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t),
+\label{introEquation:randomForceDefinition}
+\end{equation}
+the equation of motion can be rewritten as
+\begin{equation}
 m\ddot x =  - \frac{{\partial W}}{{\partial x}} - \int_0^t {\xi
 (t)\dot x(t - \tau )d\tau }  + R(t)
 \label{introEuqation:GeneralizedLangevinDynamics}
 \end{equation}
-%where $ {\xi (t)}$ is friction kernel, $R(t)$ is random force and
-%$W$ is the potential of mean force. $W(x) =  - kT\ln p(x)$
+which is known as the \emph{generalized Langevin equation}.
+
+\subsubsection{\label{introSection:randomForceDynamicFrictionKernel}Random Force and Dynamic Friction Kernel}
+
+One may notice that $R(t)$ depends only on initial conditions, which
+implies it is completely deterministic within the context of a
+harmonic bath. However, it is easy to verify that $R(t)$ is totally
+uncorrelated to $x$ and $\dot x$,
 \[
-\xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
-}}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
+\begin{array}{l}
+ \left\langle {x(t)R(t)} \right\rangle  = 0, \\
+ \left\langle {\dot x(t)R(t)} \right\rangle  = 0. \\
+ \end{array}
 \]
-For an infinite harmonic bath, we can use the spectral density and
-an integral over frequencies.
+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.
 
+%dynamic friction kernel
+The convolution integral
 \[
-R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
-- \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
-\right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
-(0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)
+\int_0^t {\xi (t)\dot x(t - \tau )d\tau }
 \]
-The random forces depend only on initial conditions.
-
-\subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
-So we can define a new set of coordinates,
+depends on the entire history of the evolution of $x$, which implies
+that the bath retains memory of previous motions. In other words,
+the bath requires a finite time to respond to change in the motion
+of the system. For a sluggish bath which responds slowly to changes
+in the system coordinate, we may regard $\xi(t)$ as a constant
+$\xi(t) = \Xi_0$. Hence, the convolution integral becomes
 \[
-q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
-^2 }}x(0)
+\int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
 \]
-This makes
+and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
 \[
-R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}
+m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
+\frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
 \]
-And since the $q$ coordinates are harmonic oscillators,
+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:
 \[
-\begin{array}{l}
- \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  = \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
- \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle  = \delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
- \end{array}
+\xi (t) = 2\xi _0 \delta (t)
+\]
+Hence, the convolution integral becomes
+\[
+\int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
+{\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
+\]
+and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
+\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
+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}
+
+Defining a new set of coordinates,
+\[
+q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
+^2 }}x(0)
+\],
+we can rewrite $R(T)$ as
+\[
+R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
 \]
+And since the $q$ coordinates are harmonic oscillators,
 
-\begin{align}
-\left\langle {R(t)R(0)} \right\rangle  &= \sum\limits_\alpha
-{\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha
-(t)q_\beta  (0)} \right\rangle } }
-%
-&= \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)}
-\right\rangle \cos (\omega _\alpha  t)}
-%
-&= kT\xi (t)
-\end{align}
+\begin{eqnarray*}
+ \left\langle {q_\alpha ^2 } \right\rangle  & = & \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
+ \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle & = & \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
+ \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle & = &\delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
+ \left\langle {R(t)R(0)} \right\rangle & = & \sum\limits_\alpha  {\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle } }  \\
+  & = &\sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
+  & = &kT\xi (t) \\
+\end{eqnarray*}
 
+Thus, we recover the \emph{second fluctuation dissipation theorem}
 \begin{equation}
 \xi (t) = \left\langle {R(t)R(0)} \right\rangle
-\label{introEquation:secondFluctuationDissipation}
+\label{introEquation:secondFluctuationDissipation}.
 \end{equation}
+In effect, it acts as a constraint on the possible ways in which one
+can model the random force and friction kernel.
 
-\section{\label{introSection:hydroynamics}Hydrodynamics}
-
 \subsection{\label{introSection:frictionTensor} Friction Tensor}
-\subsection{\label{introSection:analyticalApproach}Analytical
-Approach}
+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.
 
-\subsection{\label{introSection:approximationApproach}Approximation
-Approach}
+\subsubsection{\label{introSection:resistanceTensorRegular}The Resistance Tensor for Regular Shape}
 
-\subsection{\label{introSection:centersRigidBody}Centers of Rigid
-Body}
+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.
 
-\section{\label{introSection:correlationFunctions}Correlation Functions}
+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$.