--- trunk/tengDissertation/Methodology.tex	2006/06/06 14:56:36	2801
+++ trunk/tengDissertation/Methodology.tex	2006/06/06 19:47:27	2804
@@ -668,21 +668,567 @@ $\gamma$ is set to zero.
 the box axes does not change. It should be noted that the NPTxyz
 integrator is a special case of $NP\gamma T$ if the surface tension
 $\gamma$ is set to zero.
-
-%\section{\label{methodSection:constraintMethod}Constraint Method}
 
-%\subsection{\label{methodSection:bondConstraint}Bond Constraint for Rigid Body}
+\section{\label{methodSection:zcons}Z-Constraint Method}
 
-%\subsection{\label{methodSection:zcons}Z-constraint Method}
+Based on the fluctuation-dissipation theorem, a force
+auto-correlation method was developed by Roux and Karplus to
+investigate the dynamics of ions inside ion channels\cite{Roux1991}.
+The time-dependent friction coefficient can be calculated from the
+deviation of the instantaneous force from its mean force.
+\begin{equation}
+\xi(z,t)=\langle\delta F(z,t)\delta F(z,0)\rangle/k_{B}T,
+\end{equation}
+where%
+\begin{equation}
+\delta F(z,t)=F(z,t)-\langle F(z,t)\rangle.
+\end{equation}
 
-\section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies}
+If the time-dependent friction decays rapidly, the static friction
+coefficient can be approximated by
+\begin{equation}
+\xi_{\text{static}}(z)=\int_{0}^{\infty}\langle\delta F(z,t)\delta
+F(z,0)\rangle dt.
+\end{equation}
+Allowing diffusion constant to then be calculated through the
+Einstein relation:\cite{Marrink1994}
+\begin{equation}
+D(z)=\frac{k_{B}T}{\xi_{\text{static}}(z)}=\frac{(k_{B}T)^{2}}{\int_{0}^{\infty
+}\langle\delta F(z,t)\delta F(z,0)\rangle dt}.%
+\end{equation}
 
-\subsection{\label{methodSection:temperature}Temperature Control}
-
-\subsection{\label{methodSection:pressureControl}Pressure Control}
+The Z-Constraint method, which fixes the z coordinates of the
+molecules with respect to the center of the mass of the system, has
+been a method suggested to obtain the forces required for the force
+auto-correlation calculation.\cite{Marrink1994} However, simply
+resetting the coordinate will move the center of the mass of the
+whole system. To avoid this problem, we reset the forces of
+z-constrained molecules as well as subtract the total constraint
+forces from the rest of the system after the force calculation at
+each time step instead of resetting the coordinate.
 
-\section{\label{methodSection:hydrodynamics}Hydrodynamics}
+After the force calculation, define $G_\alpha$ as
+\begin{equation}
+G_{\alpha} = \sum_i F_{\alpha i}, \label{oopseEq:zc1}
+\end{equation}
+where $F_{\alpha i}$ is the force in the z direction of atom $i$ in
+z-constrained molecule $\alpha$. The forces of the z constrained
+molecule are then set to:
+\begin{equation}
+F_{\alpha i} = F_{\alpha i} -
+    \frac{m_{\alpha i} G_{\alpha}}{\sum_i m_{\alpha i}}.
+\end{equation}
+Here, $m_{\alpha i}$ is the mass of atom $i$ in the z-constrained
+molecule. Having rescaled the forces, the velocities must also be
+rescaled to subtract out any center of mass velocity in the z
+direction.
+\begin{equation}
+v_{\alpha i} = v_{\alpha i} -
+    \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}},
+\end{equation}
+where $v_{\alpha i}$ is the velocity of atom $i$ in the z direction.
+Lastly, all of the accumulated z constrained forces must be
+subtracted from the system to keep the system center of mass from
+drifting.
+\begin{equation}
+F_{\beta i} = F_{\beta i} - \frac{m_{\beta i} \sum_{\alpha}
+G_{\alpha}}
+    {\sum_{\beta}\sum_i m_{\beta i}},
+\end{equation}
+where $\beta$ are all of the unconstrained molecules in the system.
+Similarly, the velocities of the unconstrained molecules must also
+be scaled.
+\begin{equation}
+v_{\beta i} = v_{\beta i} + \sum_{\alpha}
+    \frac{\sum_i m_{\alpha i} v_{\alpha i}}{\sum_i m_{\alpha i}}.
+\end{equation}
 
-%\section{\label{methodSection:coarseGrained}Coarse-Grained Modeling}
+At the very beginning of the simulation, the molecules may not be at
+their constrained positions. To move a z-constrained molecule to its
+specified position, a simple harmonic potential is used
+\begin{equation}
+U(t)=\frac{1}{2}k_{\text{Harmonic}}(z(t)-z_{\text{cons}})^{2},%
+\end{equation}
+where $k_{\text{Harmonic}}$ is the harmonic force constant, $z(t)$
+is the current $z$ coordinate of the center of mass of the
+constrained molecule, and $z_{\text{cons}}$ is the constrained
+position. The harmonic force operating on the z-constrained molecule
+at time $t$ can be calculated by
+\begin{equation}
+F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}=
+    -k_{\text{Harmonic}}(z(t)-z_{\text{cons}}).
+\end{equation}
+
+
+\section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies}
+
+%\subsection{\label{methodSection:temperature}Temperature Control}
+
+%\subsection{\label{methodSection:pressureControl}Pressure Control}
+
+%\section{\label{methodSection:hydrodynamics}Hydrodynamics}
+
+%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}.
 
-%\section{\label{methodSection:moleculeScale}Molecular-Scale Modeling}
+%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{Ryckaert1977, Andersen1983}, 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{Barojas1973}. 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\cite{Kol1997}, 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.
+
+
+\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*}
+
+\subsection{Results and discussion}