--- trunk/tengDissertation/Langevin.tex	2006/06/11 02:06:01	2851
+++ trunk/tengDissertation/Langevin.tex	2006/06/28 19:09:25	2905
@@ -1,5 +1,5 @@
 
-\chapter{\label{chapt:methodology}Langevin Dynamics for Rigid Bodies of Arbitrary Shape}
+\chapter{\label{chapt:langevin}LANGEVIN DYNAMICS FOR RIGID BODIES OF ARBITRARY SHAPE}
 
 \section{Introduction}
 
@@ -28,78 +28,24 @@ systems\cite{Garcia-Palacios1998,Berkov2002,Denisov200
 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{Garcia-Palacios1998,Berkov2002,Denisov2003}. 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}. In an
-attempt to reduce the computational cost of simulation, multiple
-time stepping (MTS) methods have been introduced and have been of
-great interest to macromolecule and protein
-community\cite{Tuckerman1992}. Relying on the observation that
-forces between distant atoms generally demonstrate slower
-fluctuations than forces between close atoms, MTS method are
-generally implemented by evaluating the slowly fluctuating forces
-less frequently than the fast ones. Unfortunately, nonlinear
-instability resulting from increasing timestep in MTS simulation
-have became a critical obstruction preventing the long time
-simulation. Due to the coupling to the heat bath, Langevin dynamics
-has been shown to be able to damp out the resonance artifact more
-efficiently\cite{Sandu1999}.
+systems\cite{Palacios1998,Berkov2002,Denisov2003}. 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}. In an attempt to reduce the
+computational cost of simulation, multiple time stepping (MTS)
+methods have been introduced and have been of great interest to
+macromolecule and protein community\cite{Tuckerman1992}. Relying on
+the observation that forces between distant atoms generally
+demonstrate slower fluctuations than forces between close atoms, MTS
+method are generally implemented by evaluating the slowly
+fluctuating forces less frequently than the fast ones.
+Unfortunately, nonlinear instability resulting from increasing
+timestep in MTS simulation have became a critical obstruction
+preventing the long time simulation. Due to the coupling to the heat
+bath, Langevin dynamics has been shown to be able to damp out the
+resonance artifact more efficiently\cite{Sandu1999}.
 
-%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}}.
-
-It is very important to develop stable and efficient methods to
-integrate the equations of motion of orientational degrees of
-freedom. Euler angles are the nature choice to describe the
-rotational degrees of freedom. However, due to its singularity, the
-numerical integration of corresponding equations of motion is very
-inefficient and inaccurate. Although an alternative integrator using
-different sets of Euler angles can overcome this difficulty\cite{},
-the computational penalty and the lost of angular momentum
-conservation still remain. In 1977, a singularity free
-representation utilizing quaternions was developed by
-Evans\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\cite{}. 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.
-
-The break through in geometric literature suggests that, in order to
-develop a long-term integration scheme, one should preserve the
-geometric structure of the flow. Matubayasi and Nakahara developed a
-time-reversible integrator for rigid bodies in quaternion
-representation. Although it is not symplectic, this integrator still
-demonstrates a better long-time energy conservation than traditional
-methods because of the time-reversible nature. Extending
-Trotter-Suzuki to general system with a flat phase space, Miller and
-his colleagues devised an novel symplectic, time-reversible and
-volume-preserving integrator in quaternion representation. However,
-all of the integrators in quaternion representation suffer from the
-computational penalty of constructing a rotation matrix from
-quaternions to evolve coordinates and velocities at every time step.
-An alternative integration scheme utilizing rotation matrix directly
-is RSHAKE , in which a conjugate momentum to rotation matrix is
-introduced to re-formulate the Hamiltonian's equation and the
-Hamiltonian is evolved in a constraint manifold by iteratively
-satisfying the orthogonality constraint. However, RSHAKE is
-inefficient because of the iterative procedure. An extremely
-efficient integration scheme in rotation matrix representation,
-which also preserves the same structural properties of the
-Hamiltonian flow as Miller's integrator, is proposed by Dullweber,
-Leimkuhler and McLachlan (DLM).
-
 %review langevin/browninan dynamics for arbitrarily shaped rigid body
 Combining Langevin or Brownian dynamics with rigid body dynamics,
 one can study the slow processes in biomolecular systems. Modeling
@@ -132,11 +78,11 @@ term\cite{Beard2001}. As a complement to IBD which has
 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
-term\cite{Beard2001}. As a complement to IBD which has a lower bound
+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{Beard2001}. LTID can also deal with the rotational
+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,
@@ -151,16 +97,15 @@ sophisticated rigid body dynamics.
 estimation of friction tensor from hydrodynamics theory into the
 sophisticated rigid body dynamics.
 
-\section{Method{\label{methodSec}}}
+\section{Computational Methods{\label{methodSec}}}
 
-\subsection{\label{introSection:frictionTensor} Friction Tensor}
-Theoretically, the friction kernel can be determined using velocity
-autocorrelation function. However, this approach become impractical
-when the system become more and more complicate. Instead, various
-approaches based on hydrodynamics have been developed to calculate
-the friction coefficients. The friction effect is isotropic in
-Equation, $\zeta$ can be taken as a scalar. In general, friction
-tensor $\Xi$ is a $6\times 6$ matrix given by
+\subsection{\label{introSection:frictionTensor}Friction Tensor}
+Theoretically, the friction kernel can be determined using the
+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. In general,
+friction tensor $\Xi$ is a $6\times 6$ matrix given by
 \[
 \Xi  = \left( {\begin{array}{*{20}c}
    {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
@@ -191,8 +136,9 @@ For a spherical particle, the translational and rotati
 
 \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}}
 
-For a spherical particle, the translational and rotational friction
-constant can be calculated from Stoke's law,
+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}
    {6\pi \eta R} & 0 & 0  \\
@@ -214,9 +160,9 @@ exactly. In 1936, Perrin extended Stokes's law to gene
 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}
+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
@@ -225,46 +171,48 @@ exactly. Introducing an elliptic integral parameter $S
 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,
+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,
+and oblate :
 \[
 S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
-}}{a}
-\],
+}}{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},
+ \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 _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}\textbf{The Resistance Tensor for Arbitrary Shape}}
+\subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shapes}}
 
-Unlike spherical and other regular shaped molecules, there is not
-analytical solution for friction tensor of any arbitrary shaped
+Unlike spherical and other simply shaped molecules, there is no
+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 space, $r$-space, to all possible
-combination of rotational diffusion coefficients, $D$-space is not
+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 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
+translational and rotational motion of rigid body, general
+ellipsoids are 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
@@ -283,8 +231,8 @@ This equation is the basis for deriving the hydrodynam
 \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}
+1930, Oseen and Burgers gave a simple solution to
+Eq.~\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}
@@ -300,9 +248,9 @@ Both of the Equation \ref{introEquation:oseenTensor} a
 \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
+Both of the Eq.~\ref{introEquation:oseenTensor} and
+Eq.~\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 -
@@ -310,7 +258,6 @@ T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( 
 \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
@@ -326,17 +273,16 @@ where $\delta _{ij}$ is Kronecker delta function. Inve
 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
-
+where $\delta _{ij}$ is Kronecker delta function. Inverting the $B$
+matrix, 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)
+\end{array}} \right),
 \]
-, which can be partitioned into $N \times N$ $3 \times 3$ block
+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}
@@ -348,31 +294,30 @@ arbitrary origin $O$ can be written as
 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} } } , \\
+\begin{eqnarray}
+ \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } \notag , \\
  \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}
+ \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j. \notag \\
 \label{introEquation:ResistanceTensorArbitraryOrigin}
-\end{equation}
+\end{eqnarray}
 
 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,
+resistance (or center of reaction), at which the trace of rotational
+resistance tensor, $ \Xi ^{rr}$ reaches a minimum value.
+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},
+From 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
+origin. Given the 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}
@@ -390,9 +335,9 @@ Using Equations \ref{introEquation:definitionCR} and
    { - 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,
+Using Eq.~\ref{introEquation:definitionCR} and
+Eq.~\ref{introEquation:resistanceTensorTransformation}, one can
+locate the position of center of resistance,
 \begin{eqnarray*}
  \left( \begin{array}{l}
  x_{OR}  \\
@@ -409,11 +354,10 @@ the position of center of resistance,
  (\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$.
 
-\subsection{Langevin dynamics for rigid particles of arbitrary shape}
+\subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}
 
 Consider a Langevin equation of motions in generalized coordinates
 \begin{equation}
@@ -422,8 +366,8 @@ $V_i = V_i(v_i,\omega _i)$. The right side of Eq.
 \end{equation}
 where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
 and moment of inertial) matrix and $V_i$ is a generalized velocity,
-$V_i = V_i(v_i,\omega _i)$. The right side of Eq.
-(\ref{LDGeneralizedForm}) consists of three generalized forces in
+$V_i = V_i(v_i,\omega _i)$. The right side of
+Eq.~\ref{LDGeneralizedForm} consists of three generalized forces in
 lab-fixed frame, systematic force $F_{s,i}$, dissipative force
 $F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
 system in Newtownian mechanics typically refers to lab-fixed frame,
@@ -433,12 +377,12 @@ in body-fixed frame and converted back to lab-fixed fr
 \[
 \begin{array}{l}
  F_{f,i}^l (t) = A^T F_{f,i}^b (t), \\
- F_{r,i}^l (t) = A^T F_{r,i}^b (t) \\
- \end{array}.
+ F_{r,i}^l (t) = A^T F_{r,i}^b (t). \\
+ \end{array}
 \]
 Here, the body-fixed friction force $F_{r,i}^b$ is proportional to
 the body-fixed velocity at center of resistance $v_{R,i}^b$ and
-angular velocity $\omega _i$,
+angular velocity $\omega _i$
 \begin{equation}
 F_{r,i}^b (t) = \left( \begin{array}{l}
  f_{r,i}^b (t) \\
@@ -456,9 +400,8 @@ with zero mean and variance
 \begin{equation}
 \left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
 \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
-2k_B T\Xi _R \delta (t - t').
+2k_B T\Xi _R \delta (t - t'). \label{randomForce}
 \end{equation}
-
 The equation of motion for $v_i$ can be written as
 \begin{equation}
 m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
@@ -480,7 +423,6 @@ + \tau _{r,i}^b(t)
 \dot j_i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b (t)
 + \tau _{r,i}^b(t)
 \end{equation}
-
 Embedding the friction terms into force and torque, one can
 integrate the langevin equations of motion for rigid body of
 arbitrary shape in a velocity-Verlet style 2-part algorithm, where
@@ -500,7 +442,6 @@ $h= \delta t$:
 \mathsf{A}(t + h) &\leftarrow \mathrm{rotate}\left( h {\bf j}
     (t + h / 2) \cdot \overleftrightarrow{\mathsf{I}}^{-1} \right).
 \end{align*}
-
 In this context, the $\mathrm{rotate}$ function is the reversible
 product of the three body-fixed rotations,
 \begin{equation}
@@ -535,11 +476,10 @@ All other rotations follow in a straightforward manner
 \end{array}
 \right).
 \end{equation}
-All other rotations follow in a straightforward manner.
+All other rotations follow in a straightforward manner. After the
+first part of the propagation, the forces and body-fixed torques are
+calculated at the new positions and orientations
 
-After the first part of the propagation, the forces and body-fixed
-torques are calculated at the new positions and orientations
-
 {\tt doForces:}
 \begin{align*}
 {\bf f}(t + h) &\leftarrow
@@ -551,12 +491,9 @@ torques are calculated at the new positions and orient
 {\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
     \cdot {\bf \tau}^s(t + h).
 \end{align*}
+Once the forces and torques have been obtained at the new time step,
+the velocities can be advanced to the same time value.
 
-{\sc oopse} automatically updates ${\bf u}$ when the rotation matrix
-$\mathsf{A}$ is calculated in {\tt moveA}.  Once the forces and
-torques have been obtained at the new time step, the velocities can
-be advanced to the same time value.
-
 {\tt moveB:}
 \begin{align*}
 {\bf v}\left(t + h \right)  &\leftarrow  {\bf v}\left(t + h / 2
@@ -568,6 +505,155 @@ be advanced to the same time value.
     + \frac{h}{2} {\bf \tau}^b(t + h) .
 \end{align*}
 
-\section{Results and discussion}
+\section{Results and Discussion}
 
+The Langevin algorithm described in previous section has been
+implemented in {\sc oopse}\cite{Meineke2005} and applied to the
+studies of kinetics and thermodynamic properties in several systems.
+
+\subsection{Temperature Control}
+
+As shown in Eq.~\ref{randomForce}, random collisions associated with
+the solvent's thermal motions is controlled by the external
+temperature. The capability to maintain the temperature of the whole
+system was usually used to measure the stability and efficiency of
+the algorithm. In order to verify the stability of this new
+algorithm, a series of simulations are performed on system
+consisiting of 256 SSD water molecules with different viscosities.
+Initial configuration for the simulations is taken from a 1ns NVT
+simulation with a cubic box of 19.7166~\AA. All simulation are
+carried out with cutoff radius of 9~\AA and 2 fs time step for 1 ns
+with reference temperature at 300~K. Average temperature as a
+function of $\eta$ is shown in Table \ref{langevin:viscosity} where
+the temperatures range from 303.04~K to 300.47~K for $\eta = 0.01 -
+1$ poise. The better temperature control at higher viscosity can be
+explained by the finite size effect and relative slow relaxation
+rate at lower viscosity regime.
+\begin{table}
+\caption{AVERAGE TEMPERATURES FROM LANGEVIN DYNAMICS SIMULATIONS OF
+SSD WATER MOLECULES WITH REFERENCE TEMPERATURE AT 300~K.}
+\label{langevin:viscosity}
+\begin{center}
+\begin{tabular}{lll}
+  \hline
+  $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\
+  \hline
+  1    & 300.47 & 10.99 \\
+  0.1  & 301.19 & 11.136 \\
+  0.01 & 303.04 & 11.796 \\
+  \hline
+\end{tabular}
+\end{center}
+\end{table}
+
+Another set of calculation were performed to study the efficiency of
+temperature control using different temperature coupling schemes.
+The starting configuration is cooled to 173~K and evolved using NVE,
+NVT, and Langevin dynamic with time step of 2 fs.
+Fig.~\ref{langevin:temperature} shows the heating curve obtained as
+the systems reach equilibrium. The orange curve in
+Fig.~\ref{langevin:temperature} represents the simulation using
+Nos\'e-Hoover temperature scaling scheme with thermostat of 5 ps
+which gives reasonable tight coupling, while the blue one from
+Langevin dynamics with viscosity of 0.1 poise demonstrates a faster
+scaling to the desire temperature. In extremely lower friction
+regime (when $ \eta  \approx 0$), Langevin dynamics becomes normal
+NVE (see green curve in Fig.~\ref{langevin:temperature}) which loses
+the temperature control ability.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{temperature.eps}
+\caption[Plot of Temperature Fluctuation Versus Time]{Plot of
+temperature fluctuation versus time.} \label{langevin:temperature}
+\end{figure}
+
+\subsection{Langevin Dynamics of Banana Shaped Molecules}
+
+In order to verify that Langevin dynamics can mimic the dynamics of
+the systems absent of explicit solvents, we carried out two sets of
+simulations and compare their dynamic properties.
+Fig.~\ref{langevin:twoBanana} shows a snapshot of the simulation
+made of 256 pentane molecules and two banana shaped molecules at
+273~K. It has an equivalent implicit solvent system containing only
+two banana shaped molecules with viscosity of 0.289 center poise. To
+calculate the hydrodynamic properties of the banana shaped molecule,
+we create a rough shell model (see Fig.~\ref{langevin:roughShell}),
+in which the banana shaped molecule is represented as a ``shell''
+made of 2266 small identical beads with size of 0.3 \AA on the
+surface. Applying the procedure described in
+Sec.~\ref{introEquation:ResistanceTensorArbitraryOrigin}, we
+identified the center of resistance at $(0, 0.7482, -0.1988)$, as
+well as the resistance tensor,
+\[
+\left( {\begin{array}{*{20}c}
+0.9261 & 0 & 0&0&0.08585&0.2057\\
+0& 0.9270&-0.007063& 0.08585&0&0\\
+0&0.007063&0.7494&0.2057&0&0\\
+0&0.0858&0.2057& 58.64& 0&-8.5736\\
+0.08585&0&0&0&48.30&3.219&\\
+0.2057&0&0&0&3.219&10.7373\\
+\end{array}} \right).
+\]
+%\[
+%\left( {\begin{array}{*{20}c}
+%0.9261 & 1.310e-14 & -7.292e-15&5.067e-14&0.08585&0.2057\\
+%3.968e-14& 0.9270&-0.007063& 0.08585&6.764e-14&4.846e-14\\
+%-6.561e-16&-0.007063&0.7494&0.2057&4.846e-14&1.5036e-14\\
+%5.067e-14&0.0858&0.2057& 58.64& 8.563e-13&-8.5736\\
+%0.08585&6.764e-14&4.846e-14&1.555e-12&48.30&3.219&\\
+%0.2057&4.846e-14&1.5036e-14&-3.904e-13&3.219&10.7373\\
+%\end{array}} \right).
+%\]
+
+Curves of velocity auto-correlation functions in
+Fig.~\ref{langevin:vacf} were shown to match each other very well.
+However, because of the stochastic nature, simulation using Langevin
+dynamics was shown to decay slightly fast. In order to study the
+rotational motion of the molecules, we also calculated the auto-
+correlation function of the principle axis of the second GB
+particle, $u$.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{roughShell.eps}
+\caption[Rough shell model for banana shaped molecule]{Rough shell
+model for banana shaped molecule.} \label{langevin:roughShell}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{twoBanana.eps}
+\caption[Snapshot from Simulation of Two Banana Shaped Molecules and
+256 Pentane Molecules]{Snapshot from simulation of two Banana shaped
+molecules and 256 pentane molecules.} \label{langevin:twoBanana}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{vacf.eps}
+\caption[Plots of Velocity Auto-correlation Functions]{Velocity
+auto-correlation functions in NVE (blue) and Langevin dynamics
+(red).} \label{langevin:vacf}
+\end{figure}
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{uacf.eps}
+\caption[Auto-correlation functions of the principle axis of the
+middle GB particle]{Auto-correlation functions of the principle axis
+of the middle GB particle in NVE (blue) and Langevin dynamics
+(red).} \label{langevin:twoBanana}
+\end{figure}
+
 \section{Conclusions}
+
+We have presented a new Langevin algorithm by incorporating the
+hydrodynamics properties of arbitrary shaped molecules into an
+advanced symplectic integration scheme. The temperature control
+ability of this algorithm was demonstrated by a set of simulations
+with different viscosities. It was also shown to have significant
+advantage of producing rapid thermal equilibration over
+Nos\'{e}-Hoover method. Further studies in systems involving banana
+shaped molecules illustrated that the dynamic properties could be
+preserved by using this new algorithm as an implicit solvent model.
