--- trunk/tengDissertation/Introduction.tex	2006/06/04 20:18:07	2786
+++ trunk/tengDissertation/Introduction.tex	2006/06/05 21:00:46	2789
@@ -498,10 +498,11 @@ issue\cite{}. The velocity verlet method, which happen
 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{}. 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,7 +566,8 @@ 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}
@@ -612,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,
@@ -658,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
 \[
@@ -679,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}
@@ -697,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{Tuckerman1992}.
+the system\cite{Tuckerman1992, McLachlan1998}.
 
 \subsubsection{\label{introSection:splittingMethod}Splitting Method}
 
@@ -844,8 +847,8 @@ 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 \\
@@ -860,9 +863,9 @@ order methods based on symmetric splitting. Given a sy
 \end{equation}
 Careful choice of coefficient $a_1 \ldot b_m$ will lead to higher
 order method. Yoshida proposed an elegant way to compose higher
-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)}
@@ -983,7 +986,7 @@ Production run is the most important steps of the simu
 
 \subsection{\label{introSection:production}Production}
 
-Production run is the most important steps of the simulation, in
+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
@@ -999,22 +1002,22 @@ approach will suffer from the surface effect. To offse
 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} 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]{pbcFig.eps}
-%\caption[An illustration of periodic boundary conditions]{A 2-D
-%illustration of periodic boundary conditions. As one particle leaves
-%the right of the simulation box, an image of it enters the left.}
-%\label{introFig:pbc}
-%\end{figure}
+\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
@@ -1032,16 +1035,18 @@ discrete Fourier transforms, the particle mesh-based m
 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 are
-accelerated from $O(N^{3/2})$ to $O(N logN)$. An alternative
-approach is \emph{fast multipole method}, 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. The shifted Coulomb potential
-for particle $i$ and particle $j$ at distance $r_{rj}$ is given by:
+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
@@ -1051,12 +1056,13 @@ efficient and easy to implement.
 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]{pbcFig.eps}
-%\caption[An illustration of shifted Coulomb potential]{An illustration of shifted Coulomb potential.}
-%\label{introFigure:shiftedCoulomb}
-%\end{figure}
+\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
 
@@ -1184,7 +1190,7 @@ protein-protein docking study{\cite{Gray2003}}.
 movement of the objects in 3D gaming engine or other physics
 simulator is governed by the rigid body dynamics. In molecular
 simulation, rigid body is used to simplify the model in
-protein-protein docking study{\cite{Gray2003}}.
+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
@@ -1192,32 +1198,33 @@ 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 $Q$ and re-formulating Hamiltonian's equation, a
-symplectic integrator, RSHAKE, was proposed to evolve the
-Hamiltonian system in a constraint manifold by iteratively
+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.
-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.
+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:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
 The motion of the rigid body is Hamiltonian with the Hamiltonian
@@ -1350,7 +1357,7 @@ not be avoided in other methods\cite{}.
 \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{}.
+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
@@ -1617,30 +1624,27 @@ Operator. Below are some important properties of Lapla
 \]
 where  $p$ is real and  $L$ is called the Laplace Transform
 Operator. Below are some important properties of Laplace transform
-\begin{equation}
-\begin{array}{c}
- 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{array}
-\end{equation}
+
+\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{array}{c}
- 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{array}
-\]
+\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{array}{c}
- mL(\ddot x) =  - \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{array}
-\]
+\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:
@@ -1652,8 +1656,8 @@ transformations:
  \end{array}
 \]
 , we obtain
-\begin{align}
-m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
+\begin{eqnarray*}
+m\ddot x & = & - \frac{{\partial W(x)}}{{\partial x}} -
 \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
 }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\int_0^t {\cos (\omega
 _\alpha  t)\dot x(t - \tau )d\tau  - \left[ {g_\alpha  x_\alpha  (0)
@@ -1661,7 +1665,7 @@ _\alpha  }}\sin (\omega _\alpha  t)} } \right\}}
 (\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
+& = & \mbox{} - \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\{
@@ -1669,7 +1673,7 @@ 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}
+\end{eqnarray*}
 
 Introducing a \emph{dynamic friction kernel}
 \begin{equation}
@@ -1764,16 +1768,16 @@ And since the $q$ coordinates are harmonic oscillators
 R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
 \]
 And since the $q$ coordinates are harmonic oscillators,
-\[
-\begin{array}{c}
- \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{array}
-\]
+
+\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
@@ -1845,7 +1849,7 @@ coordinates by
 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
+coordinates by\cite{Perrin1934, Perrin1936}
 \[
 \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
 }} = 1
@@ -1890,17 +1894,17 @@ translational and rotational motion of rigid body\cite
 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\cite{}, 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\cite{} where the molecule of interest was modeled by
-combinations of spheres(beads)\cite{} 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$,
+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 }
 \]
@@ -1921,8 +1925,8 @@ introduced by Rotne and Prager\cite{} and improved by 
 \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{} and improved by Garc\'{i}a de
-la Torre and Bloomfield,
+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
@@ -2038,5 +2042,26 @@ where $x_OR$, $y_OR$, $z_OR$ are the components of the
  (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
  \end{array} \right).
 \]
+
+
+\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$.
