trunk/langevin/langevin.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
%\documentclass[aps,prb,preprint]{revtex4}
\documentclass[11pt]{article}
\usepackage{endfloat}
\usepackage{amsmath,bm}
\usepackage{amssymb}
\usepackage{epsf}
\usepackage{times}
\usepackage{mathptmx}
\usepackage{setspace}
\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{bibentry}
\usepackage{mathrsfs}
\usepackage[ref]{overcite}
\pagestyle{plain} \pagenumbering{arabic} \oddsidemargin 0.0cm
\evensidemargin 0.0cm \topmargin -21pt \headsep 10pt \textheight
9.0in \textwidth 6.5in \brokenpenalty=10000
\renewcommand{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note

\begin{document}

\title{Langevin Dynamics for Rigid Body of Arbitrary Shape }

\author{Teng Lin and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail:
gezelter@nd.edu} \\
Department of Chemistry and Biochemistry\\
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle \doublespacing

\begin{abstract}

\end{abstract}

\newpage

%\narrowtext

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                          BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%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)\cite{Dullweber1997}.

%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{Beard2001}. 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
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{Method{\label{methodSec}}}

\subsection{Friction Tensor}

For an arbitrary rigid body moves in a fluid, it may experience
friction force $f_r$ or friction torque $\tau _r$ along the opposite
direction of the velocity $v$ or angular velocity $\omega$ at
arbitrary origin $P$,
\begin{equation}
\left( \begin{array}{l}
 f_r  \\
 \tau _r  \\
 \end{array} \right) =  - \left( {\begin{array}{*{20}c}
   {\Xi _{P,t} } & {\Xi _{P,c}^T }  \\
   {\Xi _{P,c} } & {\Xi _{P,r} }  \\
\end{array}} \right)\left( \begin{array}{l}
 \nu  \\
 \omega  \\
 \end{array} \right)
\end{equation}
where $\Xi _{P,t}t$ is the translation friction tensor, $\Xi _{P,r}$
is the rotational friction tensor and $\Xi _{P,c}$ is the
translation-rotation coupling tensor. The procedure of calculating
friction tensor using hydrodynamic tensor and comparison between
bead model and shell model were elaborated by Carrasco \textit{et
al}\cite{Carrasco1999}. An important property of the friction tensor
is that the translational friction tensor is independent of origin
while the rotational and coupling are sensitive to the choice of the
origin \cite{Brenner1967}, which can be described by
\begin{equation}
\begin{array}{c}
 \Xi _{P,t}  = \Xi _{O,t}  = \Xi _t  \\
 \Xi _{P,c}  = \Xi _{O,c}  - r_{OP}  \times \Xi _t  \\
 \Xi _{P,r}  = \Xi _{O,r}  - r_{OP}  \times \Xi _t  \times r_{OP}  + \Xi _{O,c}  \times r_{OP}  - r_{OP}  \times \Xi _{O,c}^T  \\
 \end{array}
\end{equation}
Where $O$ is another origin and $r_{OP}$ is the vector joining $O$
and $P$. It is also worthy of mention that both of translational and
rotational frictional tensors are always symmetric. In contrast,
coupling tensor is only symmetric at center of reaction:
\begin{equation}
\Xi _{R,c}  = \Xi _{R,c}^T
\end{equation}
The proper location for applying friction force is the center of
reaction, at which the trace of rotational resistance tensor reaches
minimum.

\subsection{Rigid body dynamics}

The Hamiltonian of rigid body can be separated in terms of potential
energy $V(r,A)$ and kinetic energy $T(p,\pi)$,
\[
H = V(r,A) + T(v,\pi )
\]
A second-order symplectic method is now obtained by the composition
of the flow maps,
\[
\varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
_{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
\]
Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
sub-flows which corresponding to force and torque respectively,
\[
\varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
_{\Delta t/2,\tau }.
\]
Since the associated operators of $\varphi _{\Delta t/2,F} $ and
$\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
order inside $\varphi _{\Delta t/2,V}$ does not matter.

Furthermore, kinetic potential can be separated to translational
kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
\begin{equation}
T(p,\pi ) =T^t (p) + T^r (\pi ).
\end{equation}
where $ T^t (p) = \frac{1}{2}v^T m v $ and $T^r (\pi )$ is defined
by \ref{introEquation:rotationalKineticRB}. Therefore, the
corresponding flow maps are given by
\[
\varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
_{\Delta t,T^r }.
\]
The free rigid body is an example of Lie-Poisson system with
Hamiltonian function
\begin{equation}
T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
\label{introEquation:rotationalKineticRB}
\end{equation}
where $T_i^r (\pi _i ) = \frac{{\pi _i ^2 }}{{2I_i }}$ and
Lie-Poisson structure matrix,
\begin{equation}
J(\pi ) = \left( {\begin{array}{*{20}c}
   0 & {\pi _3 } & { - \pi _2 }  \\
   { - \pi _3 } & 0 & {\pi _1 }  \\
   {\pi _2 } & { - \pi _1 } & 0  \\
\end{array}} \right)
\end{equation}
Thus, the dynamics of free rigid body is governed by
\begin{equation}
\frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
\end{equation}
One may notice that each $T_i^r$ in Equation
\ref{introEquation:rotationalKineticRB} can be solved exactly. For
instance, the equations of motion due to $T_1^r$ are given by
\begin{equation}
\frac{d}{{dt}}\pi  = R_1 \pi ,\frac{d}{{dt}}A = AR_1
\label{introEqaution:RBMotionSingleTerm}
\end{equation}
where
\[ R_1  = \left( {\begin{array}{*{20}c}
   0 & 0 & 0  \\
   0 & 0 & {\pi _1 }  \\
   0 & { - \pi _1 } & 0  \\
\end{array}} \right).
\]
The solutions of Equation \ref{introEqaution:RBMotionSingleTerm} is
\[
\pi (\Delta t) = e^{\Delta tR_1 } \pi (0),A(\Delta t) =
A(0)e^{\Delta tR_1 }
\]
with
\[
e^{\Delta tR_1 }  = \left( {\begin{array}{*{20}c}
   0 & 0 & 0  \\
   0 & {\cos \theta _1 } & {\sin \theta _1 }  \\
   0 & { - \sin \theta _1 } & {\cos \theta _1 }  \\
\end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
\]
To reduce the cost of computing expensive functions in $e^{\Delta
tR_1 }$, we can use Cayley transformation,
\[
e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
)
\]
The flow maps for $T_2^r$ and $T_3^r$ can be found in the same
manner.

In order to construct a second-order symplectic method, we split the
angular kinetic Hamiltonian function into five terms
\[
T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
(\pi _1 )
\].
Concatenating flows corresponding to these five terms, we can obtain
the flow map for free rigid body,
\[
\varphi _{\Delta t,T^r }  = \varphi _{\Delta t/2,\pi _1 }  \circ
\varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }
\circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
_1 }.
\]

The equations of motion corresponding to potential energy and
kinetic energy are listed in the below table,
\begin{center}
\begin{tabular}{|l|l|}
  \hline
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  Potential & Kinetic \\
  $\frac{{dq}}{{dt}} = \frac{p}{m}$ & $\frac{d}{{dt}}q = p$ \\
  $\frac{d}{{dt}}p =  - \frac{{\partial V}}{{\partial q}}$ & $ \frac{d}{{dt}}p = 0$ \\
  $\frac{d}{{dt}}Q = 0$ & $ \frac{d}{{dt}}Q = Qskew(I^{ - 1} j)$ \\
  $ \frac{d}{{dt}}\pi  = \sum\limits_i {\left( {Q^T F_i (r,Q)} \right) \times X_i }$ & $\frac{d}{{dt}}\pi  = \pi  \times I^{ - 1} \pi$\\
  \hline
\end{tabular}
\end{center}

Finally, we obtain the overall symplectic flow maps for free moving
rigid body
\begin{align*}
 \varphi _{\Delta t}  = &\varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau } \circ  \\
  &\varphi _{\Delta t,T^t }  \circ \varphi _{\Delta t/2,\pi _1 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }  \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi _1 } \circ  \\
  &\varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
\label{introEquation:overallRBFlowMaps}
\end{align*}

\subsection{Langevin dynamics for rigid particles of arbitrary shape}

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').
\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 \pi _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 part one:}
\begin{align*}
 v_i (t + h/2) &\leftarrow v_i (t) + \frac{{hf_{t,i}^l (t)}}{{2m_i }} \\
 \pi _i (t + h/2) &\leftarrow \pi _i (t) + \frac{{h\tau _{t,i}^b (t)}}{2} \\
 r_i (t + h) &\leftarrow r_i (t) + hv_i (t + h/2) \\
 A_i (t + h) &\leftarrow rotate\left( {h\pi _i (t + h/2)I^{ - 1} } \right) \\
\end{align*}
In this context, the $\mathrm{rotate}$ function is the reversible
product of five consecutive 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 ($\pi$) 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 friction and random
forces are generated at the center of resistance in body-fixed frame
and converted back into lab-fixed frame
\[
f_{t,i}^l (t + h) =  - \left( {\frac{{\partial V}}{{\partial r_i }}}
\right)_{r_i (t + h)}  + A_i^T (t + h)[F_{f,i}^b (t + h) + F_{r,i}^b
(t + h)],
\]
while the system torque in lab-fixed frame is transformed into
body-fixed frame,
\[
\tau _{t,i}^b (t + h) = A\tau _{s,i}^l (t) + \tau _{n,i}^b (t) +
\tau _{r,i}^b (t).
\]
Once the forces and torques have been obtained at the new time step,
the velocities can be advanced to the same time value.

{\tt part two:}
\begin{align*}
 v_i (t) &\leftarrow v_i (t + h/2) + \frac{{hf_{t,i}^l (t + h)}}{{2m_i }} \\
 \pi _i (t) &\leftarrow \pi _i (t + h/2) + \frac{{h\tau _{t,i}^b (t + h)}}{2} \\
\end{align*}

\section{Results and discussion}

\subsection{}

\subsection{Lipid bilayer}

\subsection{Liquid crystal}

\section{Conclusions}

\section{Acknowledgments}
Support for this project was provided by the National Science
Foundation under grant CHE-0134881. T.L. also acknowledges the
financial support from center of applied mathematics at University
of Notre Dame.
\newpage

\bibliographystyle{jcp2}
\bibliography{langevin}

\end{document}
Revision:	2746
Committed:	Wed May 3 16:44:46 2006 UTC (18 years, 4 months ago) by tim
Content type:	application/x-tex
File size:	22599 byte(s)
Log Message:	initial draft