trunk/tengDissertation/Langevin.tex


\chapter{\label{chapt:langevin}LANGEVIN DYNAMICS FOR RIGID BODIES OF ARBITRARY SHAPE}

\section{Introduction}

%applications of langevin dynamics
As an excellent 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 simulation. Implicit solvent Langevin
dynamics simulation of met-enkephalin not only outperforms explicit
solvent simulation on computation efficiency, but also agrees very
well with explicit solvent simulation on dynamics
properties\cite{Shen2002}. Recently, applying Langevin dynamics with
UNRES model, Liow and his coworkers suggest that protein folding
pathways can be possibly exploited 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 barrier\cite{Banerjee2004, Cui2003}.
Combining Langevin dynamics with Kramers's theory, Klimov and
Thirumalai identified the free-energy barrier 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 study 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 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 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
the DNA as a chain of rigid spheres beads, which subject to harmonic
potentials as well as excluded volume potentials, Mielke and his
coworkers discover 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 millisecond, which
is impracticable to study using all atomistic model with newtonian
mechanics. With the help of coarse-grained rigid body model and
stochastic dynamics, the fusion pathways were exploited by many
researchers\cite{Noguchi2001,Noguchi2002,Shillcock2005}. Due to the
difficulty of numerical integration of anisotropy rotation, most of
the rigid body models are simply modeled by sphere, cylinder,
ellipsoid or other regular shapes in stochastic simulations. In an
effort to account for the diffusion anisotropy of the 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
error and bias 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, assumption of the 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
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 cylinder and ellipsoid are inadequate to
represent most complex macromolecule assemblies. These intricate
molecules have been represented by a set of beads and their
hydrodynamics properties can be calculated using variant
hydrodynamic interaction tensors.

The goal of the present work is to develop a Langevin dynamics
algorithm for arbitrary rigid particles by integrating the accurate
estimation of friction tensor from hydrodynamics theory into the
sophisticated rigid body dynamics.

\section{Computational Methods{\label{methodSec}}}

\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} }  \\
   {\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.

\subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}}

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  \\
   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.

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
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 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}\textbf{The Resistance Tensor for Arbitrary Shapes}}

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 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
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
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 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)
\]
, 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 (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}
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 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}
 \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$.

\subsection{Langevin Dynamics for Rigid Particles of Arbitrary Shape\label{LDRB}}

Consider a Langevin equation of motions in generalized coordinates
\begin{equation}
M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)}  + F_{r,i} (t)
\label{LDGeneralizedForm}
\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
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,
it is also convenient to handle the rotation of rigid body in
body-fixed frame. Thus the friction and random forces are calculated
in body-fixed frame and converted back to lab-fixed frame by:
\[
\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}.
\]
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$,
\begin{equation}
F_{r,i}^b (t) = \left( \begin{array}{l}
 f_{r,i}^b (t) \\
 \tau _{r,i}^b (t) \\
 \end{array} \right) =  - \left( {\begin{array}{*{20}c}
   {\Xi _{R,t} } & {\Xi _{R,c}^T }  \\
   {\Xi _{R,c} } & {\Xi _{R,r} }  \\
\end{array}} \right)\left( \begin{array}{l}
 v_{R,i}^b (t) \\
 \omega _i (t) \\
 \end{array} \right),
\end{equation}
while the random force $F_{r,i}^l$ is a Gaussian stochastic variable
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'). \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) +
f_{r,i}^l (t)
\end{equation}
Since the frictional force is applied at the center of resistance
which generally does not coincide with the center of mass, an extra
torque is exerted at the center of mass. Thus, the net body-fixed
frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is
given by
\begin{equation}
\tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b
\end{equation}
where $r_{MR}$ is the vector from the center of mass to the center
of the resistance. Instead of integrating angular velocity in
lab-fixed frame, we consider the equation of motion of angular
momentum in body-fixed frame
\begin{equation}
\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
$h= \delta t$:

{\tt moveA:}
\begin{align*}
{\bf v}\left(t + h / 2\right)  &\leftarrow  {\bf v}(t)
    + \frac{h}{2} \left( {\bf f}(t) / m \right), \\
%
{\bf r}(t + h) &\leftarrow {\bf r}(t)
    + h  {\bf v}\left(t + h / 2 \right), \\
%
{\bf j}\left(t + h / 2 \right)  &\leftarrow {\bf j}(t)
    + \frac{h}{2} {\bf \tau}^b(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}
\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y
/ 2) \cdot \mathsf{G}_x(a_x /2),
\end{equation}
where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
rotates both the rotation matrix ($\mathsf{A}$) and the body-fixed
angular momentum (${\bf j}$) by an angle $\theta$ around body-fixed
axis $\alpha$,
\begin{equation}
\mathsf{G}_\alpha( \theta ) = \left\{
\begin{array}{lcl}
\mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
{\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
j}(0).
\end{array}
\right.
\end{equation}
$\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis
rotation matrix.  For example, in the small-angle limit, the
rotation matrix around the body-fixed x-axis can be approximated as
\begin{equation}
\mathsf{R}_x(\theta) \approx \left(
\begin{array}{ccc}
1 & 0 & 0 \\
0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
\theta^2 / 4} \\
0 & \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 +
\theta^2 / 4}
\end{array}
\right).
\end{equation}
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

{\tt doForces:}
\begin{align*}
{\bf f}(t + h) &\leftarrow
    - \left(\frac{\partial V}{\partial {\bf r}}\right)_{{\bf r}(t + h)}, \\
%
{\bf \tau}^{s}(t + h) &\leftarrow {\bf u}(t + h)
    \times \frac{\partial V}{\partial {\bf u}}, \\
%
{\bf \tau}^{b}(t + h) &\leftarrow \mathsf{A}(t + h)
    \cdot {\bf \tau}^s(t + h).
\end{align*}

{\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
\right)
    + \frac{h}{2} \left( {\bf f}(t + h) / m \right), \\
%
{\bf j}\left(t + h \right)  &\leftarrow {\bf j}\left(t + h / 2
\right)
    + \frac{h}{2} {\bf \tau}^b(t + h) .
\end{align*}

\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}{|l|l|l|}
  \hline
  $\eta$ & $\text{T}_{\text{avg}}$ & $\text{T}_{\text{rms}}$ \\
  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 Molecule}

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 & 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).
\]


\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 function of NVE (blue) and Langevin dynamics
(red).} \label{langevin:twoBanana}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=\linewidth]{uacf.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}

\section{Conclusions}
Revision:	2889
Committed:	Mon Jun 26 13:42:53 2006 UTC (18 years ago) by tim
Content type:	application/x-tex
File size:	27267 byte(s)
Log Message:	Change the case caption of tables