--- trunk/langevin/langevin.tex	2008/01/11 22:08:57	3310
+++ trunk/langevin/langevin.tex	2008/01/18 20:43:53	3316
@@ -48,224 +48,255 @@ As alternative to Newtonian dynamics, Langevin dynamic
 \section{Introduction}
 
 %applications of langevin dynamics
-As 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. The stochastic treatment
-of the solvent enables us to carry out substantially longer time
-simulations. Implicit solvent Langevin dynamics simulations of
-met-enkephalin not only outperform explicit solvent simulations for
-computational efficiency, but also agrees very well with explicit
-solvent simulations for dynamical properties.\cite{Shen2002}
-Recently, applying Langevin dynamics with the UNRES model, Liow and
-his coworkers suggest that protein folding pathways can be possibly
-explored within a reasonable amount of time.\cite{Liwo2005} The
-stochastic nature of the Langevin dynamics also enhances the
-sampling of the system and increases the probability of crossing
-energy barriers.\cite{Banerjee2004, Cui2003} Combining Langevin
-dynamics with Kramers's theory, Klimov and Thirumalai identified
-free-energy barriers by studying the viscosity dependence of the
-protein folding rates.\cite{Klimov1997} In order to account for
-solvent induced interactions missing from implicit solvent model,
-Kaya incorporated desolvation free energy barrier into implicit
-coarse-grained solvent model in protein folding/unfolding studies
-and discovered a higher free energy barrier between the native and
-denatured states. Because of its stability against noise, Langevin
-dynamics is very suitable for studying remagnetization processes in
-various systems.\cite{Palacios1998,Berkov2002,Denisov2003} For
+Langevin dynamics, which mimics a simple heat bath with stochastic and
+dissipative forces, has been applied in a variety of situations as an
+alternative to molecular dynamics with explicit solvent molecules.
+The stochastic treatment of the solvent allows the use of simulations
+with substantially longer time and length scales.  In general, the
+dynamic and structural properties obtained from Langevin simulations
+agree quite well with similar properties obtained from explicit
+solvent simulations.
+
+Recent examples of the usefulness of Langevin simulations include a
+study of met-enkephalin in which Langevin simulations predicted
+dynamical properties that were largely in agreement with explicit
+solvent simulations.\cite{Shen2002} By applying Langevin dynamics with
+the UNRES model, Liow and his coworkers suggest that protein folding
+pathways can be explored within a reasonable amount of
+time.\cite{Liwo2005}
+
+The stochastic nature of Langevin dynamics also enhances the sampling
+of the system and increases the probability of crossing energy
+barriers.\cite{Cui2003,Banerjee2004} Combining Langevin dynamics with
+Kramers's theory, Klimov and Thirumalai identified free-energy
+barriers by studying the viscosity dependence of the protein folding
+rates.\cite{Klimov1997} In order to account for solvent induced
+interactions missing from the implicit solvent model, Kaya
+incorporated a desolvation free energy barrier into protein
+folding/unfolding studies and discovered a higher free energy barrier
+between the native and denatured states.\cite{XXX}
+
+Because of its stability against noise, Langevin dynamics has also
+proven useful for studying remagnetization processes in various
+systems.\cite{Palacios1998,Berkov2002,Denisov2003} [Check: For
 instance, the oscillation power spectrum of nanoparticles from
-Langevin dynamics simulation has the same peak frequencies for
-different wave vectors, which recovers the property of magnetic
-excitations in small finite structures.\cite{Berkov2005a}
+Langevin dynamics has the same peak frequencies for different wave
+vectors, which recovers the property of magnetic excitations in small
+finite structures.\cite{Berkov2005a}]
 
-%review rigid body dynamics
-Rigid bodies are frequently involved in the modeling of different
-areas, from engineering, physics, to chemistry. For example,
-missiles and vehicle are usually modeled by rigid bodies.  The
-movement of the objects in 3D gaming engine or other physics
-simulator is governed by the rigid body dynamics. In molecular
-simulation, rigid body is used to simplify the model in
-protein-protein docking study{\cite{Gray2003}}.
+In typical LD simulations, the friction and random forces on
+individual atoms are taken from the Stokes-Einstein hydrodynamic
+approximation,
+\begin{eqnarray}
+m \dot{v}(t) & = & -\nabla U(x) - \xi m v(t) + R(t) \\
+\langle R(t) \rangle & = & 0 \\
+\langle R(t) R(t') \rangle & = & 2 k_B T \xi m \delta(t - t')
+\end{eqnarray}
+where $\xi \approx 6 \pi \eta a$.  Here $\eta$ is the viscosity of the
+implicit solvent, and $a$ is the hydrodynamic radius of the atom.
 
-It is very important to develop stable and efficient methods to
-integrate the equations of motion for orientational degrees of
-freedom. Euler angles are the natural choice to describe the
-rotational degrees of freedom. However, due to $\frac {1}{sin
-\theta}$ singularities, the numerical integration of corresponding
-equations of these motion is very inefficient and inaccurate.
-Although an alternative integrator using multiple sets of Euler
-angles can overcome this difficulty\cite{Barojas1973}, the
-computational penalty and the loss of angular momentum conservation
-still remain. A singularity-free representation utilizing
-quaternions was developed by Evans in 1977.\cite{Evans1977}
-Unfortunately, this approach used a nonseparable Hamiltonian
-resulting from the quaternion representation, which prevented the
-symplectic algorithm from being utilized. Another different approach
-is to apply holonomic constraints to the atoms belonging to the
-rigid body. Each atom moves independently under the normal forces
-deriving from potential energy and constraint forces which are used
-to guarantee the rigidness. However, due to their iterative nature,
-the SHAKE and Rattle algorithms also converge very slowly when the
-number of constraints increases.\cite{Ryckaert1977, Andersen1983}
+The use of rigid substructures,\cite{???} 
+coarse-graining,\cite{Ayton,Sun,Zannoni} and ellipsoidal
+representations of protein side chains~\cite{Schulten} has made the
+use of the Stokes-Einstein approximation problematic.  A rigid
+substructure moves as a single unit with orientational as well as
+translational degrees of freedom.  This requires a more general
+treatment of the hydrodynamics than the spherical approximation
+provides.  The atoms involved in a rigid or coarse-grained structure
+should properly have solvent-mediated interactions with each
+other. The theory of interactions {\it between} bodies moving through
+a fluid has been developed over the past century and has been applied
+to simulations of Brownian
+motion.\cite{MarshallNewton,GarciaDeLaTorre} There a need to have a
+more thorough treatment of hydrodynamics included in the library of
+methods available for performing Langevin simulations.
 
-A break-through in geometric literature suggests that, in order to
+\subsection{Rigid Body Dynamics}
+Rigid bodies are frequently involved in the modeling of large
+collections of particles that move as a single unit.  In molecular
+simulations, rigid bodies have been used to simplify protein-protein
+docking,\cite{Gray2003} and lipid bilayer simulations.\cite{Sun2008}
+Many of the water models in common use are also rigid-body
+models,\cite{TIPs,SPC/E} although they are typically evolved using
+constraints rather than rigid body equations of motion.
+
+Euler angles are a natural choice to describe the rotational
+degrees of freedom.  However, due to $1 \over \sin \theta$
+singularities, the numerical integration of corresponding equations of
+these motion can become inaccurate (and inefficient).  Although an
+alternative integrator using multiple sets of Euler angles can
+overcome this problem,\cite{Barojas1973} the computational penalty and
+the loss of angular momentum conservation remain. A singularity-free
+representation utilizing quaternions was developed by Evans in
+1977.\cite{Evans1977} Unfortunately, this approach uses a nonseparable
+Hamiltonian resulting from the quaternion representation, which
+prevented symplectic algorithms from being utilized until very
+recently.\cite{Miller2002} Another approach is the application of
+holonomic constraints to the atoms belonging to the rigid body. Each
+atom moves independently under the normal forces deriving from
+potential energy and constraint forces which are used to guarantee the
+rigidness. However, due to their iterative nature, the SHAKE and
+Rattle algorithms also converge very slowly when the number of
+constraints increases.\cite{Ryckaert1977,Andersen1983}
+
+A breakthrough in geometric literature suggests that, in order to
 develop a long-term integration scheme, one should preserve the
 symplectic structure of the propagator. By introducing a conjugate
 momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's
-equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was
-proposed to evolve the Hamiltonian system in a constraint manifold
-by iteratively satisfying the orthogonality constraint $Q^T Q = 1$.
-An alternative method using the quaternion representation was
-developed by Omelyan.\cite{Omelyan1998} However, both of these
-methods are iterative and inefficient. In this section, we descibe a
-symplectic Lie-Poisson integrator for rigid bodies developed by
-Dullweber and his coworkers\cite{Dullweber1997} in depth.
+equation, a symplectic integrator, RSHAKE,\cite{Kol1997} was proposed
+to evolve the Hamiltonian system in a constraint manifold by
+iteratively satisfying the orthogonality constraint $Q^T Q = 1$.  An
+alternative method using the quaternion representation was developed
+by Omelyan.\cite{Omelyan1998} However, both of these methods are
+iterative and suffer from some related inefficiencies. A symplectic
+Lie-Poisson integrator for rigid bodies developed by Dullweber {\it et
+al.}\cite{Dullweber1997} gets around most of the limitations mentioned
+above and has become the basis for our Langevin integrator.
 
-%review langevin/browninan dynamics for arbitrarily shaped rigid body
-Combining Langevin or Brownian dynamics with rigid body dynamics,
-one can study slow processes in biomolecular systems. Modeling DNA
-as a chain of rigid beads, which are subject to harmonic potentials
-as well as excluded volume potentials, Mielke and his coworkers
-discovered rapid superhelical stress generations from the stochastic
-simulation of twin supercoiling DNA with response to induced
-torques.\cite{Mielke2004} Membrane fusion is another key biological
-process which controls a variety of physiological functions, such as
-release of neurotransmitters \textit{etc}. A typical fusion event
-happens on the time scale of a millisecond, which is impractical to
-study using atomistic models with newtonian mechanics. With the help
-of coarse-grained rigid body model and stochastic dynamics, the
-fusion pathways were explored by many
-researchers.\cite{Noguchi2001,Noguchi2002,Shillcock2005} Due to the
-difficulty of numerical integration of anisotropic rotation, most of
-the rigid body models are simply modeled using spheres, cylinders,
-ellipsoids or other regular shapes in stochastic simulations. In an
-effort to account for the diffusion anisotropy of arbitrary
-particles, Fernandes and de la Torre improved the original Brownian
-dynamics simulation algorithm\cite{Ermak1978,Allison1991} by
-incorporating a generalized $6\times6$ diffusion tensor and
-introducing a simple rotation evolution scheme consisting of three
-consecutive rotations.\cite{Fernandes2002} Unfortunately, unexpected
-errors and biases are introduced into the system due to the
+
+\subsection{The Hydrodynamic tensor and Brownian dynamics}
+Combining Brownian dynamics with rigid substructures, one can study
+slow processes in biomolecular systems.  Modeling DNA as a chain of
+beads which are subject to harmonic potentials as well as excluded
+volume potentials, Mielke and his coworkers discovered rapid
+superhelical stress generations from the stochastic simulation of twin
+supercoiling DNA with response to induced torques.\cite{Mielke2004}
+Membrane fusion is another key biological process which controls a
+variety of physiological functions, such as release of
+neurotransmitters \textit{etc}. A typical fusion event happens on the
+time scale of a millisecond, which is impractical to study using
+atomistic models with newtonian mechanics. With the help of
+coarse-grained rigid body model and stochastic dynamics, the fusion
+pathways were explored by Noguchi and others.\cite{Noguchi2001,Noguchi2002,Shillcock2005}
+
+Due to the difficulty of numerically integrating anisotropic
+rotational motion, most of the coarse-grained rigid body models are
+treated as spheres, cylinders, ellipsoids or other regular shapes in
+stochastic simulations.  In an effort to account for the diffusion
+anisotropy of arbitrarily-shaped particles, Fernandes and Garc\'{i}a
+de la Torre improved the original Brownian dynamics simulation
+algorithm~\cite{Ermak1978,Allison1991} by incorporating a generalized
+$6\times6$ diffusion tensor and introducing a rotational evolution
+scheme consisting of three consecutive rotations.\cite{Fernandes2002}
+Unfortunately, biases are introduced into the system due to the
 arbitrary order of applying the noncommuting rotation
 operators.\cite{Beard2003} Based on the observation the momentum
 relaxation time is much less than the time step, one may ignore the
-inertia in Brownian dynamics. However, the assumption of zero
-average acceleration is not always true for cooperative motion which
-is common in protein motion. An inertial Brownian dynamics (IBD) was
-proposed to address this issue by adding an inertial correction
+inertia in Brownian dynamics.  However, the assumption of zero average
+acceleration is not always true for cooperative motion which is common
+in proteins. An inertial Brownian dynamics (IBD) was proposed to
+address this issue by adding an inertial correction
 term.\cite{Beard2000} As a complement to IBD which has a lower bound
 in time step because of the inertial relaxation time, long-time-step
 inertial dynamics (LTID) can be used to investigate the inertial
 behavior of the polymer segments in low friction
 regime.\cite{Beard2000} LTID can also deal with the rotational
 dynamics for nonskew bodies without translation-rotation coupling by
-separating the translation and rotation motion and taking advantage
-of the analytical solution of hydrodynamics properties. However,
-typical nonskew bodies like cylinders and ellipsoids are inadequate
-to represent most complex macromolecule assemblies. These intricate
-molecules have been represented by a set of beads and their
-hydrodynamic properties can be calculated using variants on the
-standard hydrodynamic interaction tensors.
+separating the translation and rotation motion and taking advantage of
+the analytical solution of hydrodynamics properties. However, typical
+nonskew bodies like cylinders and ellipsoids are inadequate to
+represent most complex macromolecular assemblies.
 
 The goal of the present work is to develop a Langevin dynamics
 algorithm for arbitrary-shaped rigid particles by integrating the
-accurate estimation of friction tensor from hydrodynamics theory
-into the sophisticated rigid body dynamics algorithms.
+accurate estimation of friction tensor from hydrodynamics theory into
+a symplectic rigid body dynamics propagator.  In the sections below,
+we review some of the theory of hydrodynamic tensors developed for
+Brownian simulations of rigid multi-particle systems, we then present
+our integration method for a set of generalized Langevin equations of
+motion, and we compare the behavior of the new Langevin integrator to
+dynamical quantities obtained via explicit solvent molecular dynamics.
 
-\subsection{\label{introSection:frictionTensor}Friction Tensor}
-Theoretically, the friction kernel can be determined using the
+\subsection{\label{introSection:frictionTensor}The Friction Tensor}
+Theoretically, a complete friction kernel can be determined using the
 velocity autocorrelation function. However, this approach becomes
-impractical when the system becomes more and more complicated.
-Instead, various approaches based on hydrodynamics have been
-developed to calculate the friction coefficients. In general, the
-friction tensor $\Xi$ is a $6\times 6$ matrix given by
-\[
+impractical when the solute becomes complex.  Instead, various
+approaches based on hydrodynamics have been developed to calculate the
+friction coefficients. In general, the friction tensor $\Xi$ is a
+$6\times 6$ matrix given by
+\begin{equation}
 \Xi  = \left( {\begin{array}{*{20}c}
    {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
    {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
 \end{array}} \right).
-\]
-Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are $3 \times 3$
-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,
-\[
+\end{equation}
+Here, $\Xi^{tt}$ and $\Xi^{rr}$ are $3 \times 3$ translational and
+rotational resistance (friction) tensors 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 ($\mathbf{F}_f$) and torque
+($\mathbf{\tau}_f$) in opposition to the directions of the velocity
+($\mathbf{v}$) and body-fixed angular velocity ($\mathbf{\omega}$),
+\begin{equation}
 \left( \begin{array}{l}
- F_R  \\
- \tau _R  \\
+ \mathbf{F}_f  \\
+ \mathbf{\tau}_f  \\
  \end{array} \right) =  - \left( {\begin{array}{*{20}c}
-   {\Xi ^{tt} } & {\Xi ^{rt} }  \\
-   {\Xi ^{tr} } & {\Xi ^{rr} }  \\
+   \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
-torque.
+ \mathbf{v} \\
+ \mathbf{\omega} \\
+ \end{array} \right).
+\end{equation}
 
 \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shapes}}
-
-For a spherical particle with slip boundary conditions, the
-translational and rotational friction constant can be calculated
-from Stoke's law,
-\[
-\Xi ^{tt}  = \left( {\begin{array}{*{20}c}
+For a spherical particle under ``stick'' boundary conditions, the
+translational and rotational friction tensors can be calculated from
+Stoke's law,
+\begin{equation}
+\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)
-\]
+\end{equation}
 and
-\[
+\begin{equation}
 \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)
-\]
+\end{equation}
 where $\eta$ is the viscosity of the solvent and $R$ is the
 hydrodynamic radius.
 
 Other non-spherical shapes, such as cylinders and ellipsoids, are
-widely used as references for developing new hydrodynamics theory,
+widely used as references for developing new hydrodynamics theories,
 because their properties can be calculated exactly. In 1936, Perrin
 extended Stokes's law to general ellipsoids, 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 being 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
-ellipsoids :
-\[
-S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
-} }}{b},
-\]
-and oblate ellipsoids:
-\[
-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{eqnarray*}
- \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{eqnarray*}
-and
-\begin{eqnarray*}
- \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}}.
+by\cite{Perrin1934,Perrin1936}
+\begin{equation}
+\frac{x^2 }{a^2} + \frac{y^2}{b^2} + \frac{z^2 }{c^2} = 1
+\end{equation}
+where the semi-axes are of lengths $a$, $b$, and $c$. Due to the
+complexity of the elliptic integral, only uniaxial ellipsoids,
+{\it i.e.} prolate ($ a \ge b = c$) and oblate ($ a < b = c $), can
+be solved exactly. Introducing an elliptic integral parameter $S$ for
+prolate ellipsoids :
+\begin{equation}
+S = \frac{2}{\sqrt{a^2  - b^2}} \ln \frac{a + \sqrt{a^2  - b^2}}{b},
+\end{equation}
+and oblate ellipsoids:
+\begin{equation}
+S = \frac{2}{\sqrt {b^2  - a^2 }} \arctan \frac{\sqrt {b^2  - a^2}}{a},
+\end{equation}
+one can write down the translational and rotational resistance
+tensors for oblate,
+\begin{eqnarray*}
+ \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{eqnarray*}
+and prolate,
+\begin{eqnarray*}
+ \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{eqnarray*}
+ellipsoids. For both spherical and ellipsoidal particles, the
+translation-rotation and rotation-translation coupling tensors are
+zero.
 
 \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}}
 
@@ -273,21 +304,21 @@ hydrodynamic properties of rigid bodies. However, sinc
 analytical solution for the friction tensor for arbitrarily 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 spaces, $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 bodies, general
-ellipsoids are not always suitable for modeling arbitrarily shaped
-rigid molecules. A number of studies have been devoted to
+hydrodynamic properties of rigid bodies. However, the mapping from all
+possible ellipsoidal spaces, $r$-space, to all possible combination of
+rotational diffusion coefficients, $D$-space, is not
+unique.\cite{Wegener1979} Additionally, because there is intrinsic
+coupling between translational and rotational motion of rigid bodies,
+general ellipsoids are not always suitable for modeling arbitrarily
+shaped rigid molecules.  A number of studies have been devoted to
 determining the friction tensor for irregularly shaped rigid bodies
-using more advanced methods where the molecule of interest was
-modeled by a combinations of spheres\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 hydrodynamic
-interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different
-than its unperturbed velocity $v_i$,
+using more advanced methods where the molecule of interest was modeled
+by a combinations of spheres\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 hydrodynamic 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 }
 \]
@@ -1061,16 +1092,47 @@ used recently as models for lipid molecules.\cite{SunG
 
 \subsection{Composite sphero-ellipsoids}
 Spherical heads perched on the ends of Gay-Berne ellipsoids have been
-used recently as models for lipid molecules.\cite{SunGezelter08,Ayton01}
-
+used recently as models for lipid
+molecules.\cite{SunGezelter08,Ayton01}
 MORE DETAILS
 
+A reference system composed of a single lipid rigid body embedded in a
+sea of 1929 solvent particles was created and run under standard
+(microcanonical) molecular dynamics.  The resulting viscosity of this
+mixture was 0.349 centipoise (as estimated using
+Eq. (\ref{eq:shear})).  To calculate the hydrodynamic properties of
+the lipid rigid body model, we created a rough shell (see
+Fig.~\ref{fig:roughShell}), in which the lipid is represented as a
+``shell'' made of 3550 identical beads (0.25 \AA\ in diameter)
+distributed on the surface.  Applying the procedure described in
+Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
+identified the center of resistance, ${\bf r} = $(0 \AA, 0 \AA, 1.46
+\AA).
 
+
 \subsection{Summary}
 According to our simulations, the langevin dynamics is a reliable
 theory to apply to replace the explicit solvents, especially for the
 translation properties. For large molecules, the rotation properties
 are also mimiced reasonablly well.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{graph}
+\caption[Mean squared displacements and orientational
+correlation functions for each of the model rigid bodies.]{The
+mean-squared displacements ($\langle r^2(t) \rangle$) and
+orientational correlation functions ($C_2(t)$) for each of the model
+rigid bodies studied.  The circles are the results for microcanonical
+simulations with explicit solvent molecules, while the other data sets
+are results for Langevin dynamics using the different hydrodynamic
+tensor approximations.  The Perrin model for the ellipsoids is
+considered the ``exact'' hydrodynamic behavior (this can also be said
+for the translational motion of the dumbbell operating under the bead
+model). In most cases, the various hydrodynamics models reproduce
+each other quantitatively.}
+\label{fig:results} 
+\end{figure}
 
 \begin{table*}
 \begin{minipage}{\linewidth}
@@ -1090,11 +1152,11 @@ sphere    & 0.261  & ?    & & 2.59 & exact       & 2.5
 \cline{2-3} \cline{5-7}
 model & $\eta$ (centipoise)  & D & & Analytical & method & Hydrodynamics & simulation \\ 
 \hline
-sphere    & 0.261  & ?    & & 2.59 & exact       & 2.59 & 2.56 \\
+sphere    & 0.279  & 3.06 & & 2.42 & exact       & 2.42 & 2.33 \\
 ellipsoid & 0.255  & 2.44 & & 2.34 & exact       & 2.34 & 2.37 \\
           & 0.255  & 2.44 & & 2.34 & rough shell & 2.36 & 2.28 \\
-dumbbell  & 0.322  & ?    & & 1.57 & bead model  & 1.57 & 1.57 \\
-          & 0.322  & ?    & & 1.57 & rough shell & ?    & ?    \\
+dumbbell  & 0.308  & 2.06 & & 1.64 & bead model  & 1.65 & 1.62 \\
+          & 0.308  & 2.06 & & 1.64 & rough shell & 1.59 & 1.62 \\
 banana    & 0.298  & 1.53 & &      & rough shell & 1.56 & 1.55 \\
 lipid     & 0.349  & 0.96 & &      & rough shell & 1.33 & 1.32 \\
 \end{tabular}
@@ -1119,11 +1181,11 @@ sphere    & 0.261  &      & & 9.06 & exact       & 9.0
 \cline{2-3} \cline{5-7}
 model & $\eta$ (centipoise) & $\tau$ & & Perrin & method & Hydrodynamic  & simulation \\ 
 \hline
-sphere    & 0.261  &      & & 9.06 & exact       & 9.06 & 9.11 \\
+sphere    & 0.279  &      & & 9.69 & exact       & 9.69 & 9.64 \\
 ellipsoid & 0.255  & 46.7 & & 22.0 & exact       & 22.0 & 22.2 \\
           & 0.255  & 46.7 & & 22.0 & rough shell & 22.6 & 22.2 \\
-dumbbell  & 0.322  & 14.0 & &      & bead model  & 52.3 & 52.8 \\
-          & 0.322  & 14.0 & &      & rough shell & ?    & ?    \\
+dumbbell  & 0.308  & 14.1 & &      & bead model  & 50.0 & 50.1 \\
+          & 0.308  & 14.1 & &      & rough shell & 41.5 & 41.3 \\
 banana    & 0.298  & 63.8 & &      & rough shell & 70.9 & 70.9 \\
 lipid     & 0.349  & 78.0 & &      & rough shell & 76.9 & 77.9 \\
 \hline
