ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/tengDissertation/Methodology.tex
(Generate patch)

Comparing trunk/tengDissertation/Methodology.tex (file contents):
Revision 2804 by tim, Tue Jun 6 19:47:27 2006 UTC vs.
Revision 2853 by tim, Sun Jun 11 02:23:00 2006 UTC

# Line 16 | Line 16 | two decades. Matubayasi and Nakahara developed a time-
16  
17   Integration schemes for rotational motion of the rigid molecules in
18   microcanonical ensemble have been extensively studied in the last
19 < two decades. Matubayasi and Nakahara developed a time-reversible
20 < integrator for rigid bodies in quaternion representation. Although
21 < it is not symplectic, this integrator still demonstrates a better
22 < long-time energy conservation than traditional methods because of
23 < the time-reversible nature. Extending Trotter-Suzuki to general
24 < system with a flat phase space, Miller and his colleagues devised an
25 < novel symplectic, time-reversible and volume-preserving integrator
26 < in quaternion representation, which was shown to be superior to the
27 < time-reversible integrator of Matubayasi and Nakahara. However, all
28 < of the integrators in quaternion representation suffer from the
19 > two decades. Matubayasi developed a time-reversible integrator for
20 > rigid bodies in quaternion representation. Although it is not
21 > symplectic, this integrator still demonstrates a better long-time
22 > energy conservation than traditional methods because of the
23 > time-reversible nature. Extending Trotter-Suzuki to general system
24 > with a flat phase space, Miller and his colleagues devised an novel
25 > symplectic, time-reversible and volume-preserving integrator in
26 > quaternion representation, which was shown to be superior to the
27 > Matubayasi's time-reversible integrator. However, all of the
28 > integrators in quaternion representation suffer from the
29   computational penalty of constructing a rotation matrix from
30   quaternions to evolve coordinates and velocities at every time step.
31   An alternative integration scheme utilizing rotation matrix directly
# Line 139 | Line 139 | in Fig.~\ref{timestep}.
139   average 7\% increase in computation time using the DLM method in
140   place of quaternions. This cost is more than justified when
141   comparing the energy conservation of the two methods as illustrated
142 < in Fig.~\ref{timestep}.
142 > in Fig.~\ref{methodFig:timestep}.
143  
144   \begin{figure}
145   \centering
# Line 356 | Line 356 | relaxation of the pressure to the target value.  To se
356   \end{equation}
357  
358   In eq.(\ref{eq:melchionna1}), $\tau_B$ is the time constant for
359 < relaxation of the pressure to the target value.  To set values for
360 < $\tau_B$ or $P_{\mathrm{target}}$ in a simulation, one would use the
361 < {\tt tauBarostat} and {\tt targetPressure} keywords in the {\tt
362 < .bass} file.  The units for {\tt tauBarostat} are fs, and the units
363 < for the {\tt targetPressure} are atmospheres.  Like in the NVT
359 > relaxation of the pressure to the target value. Like in the NVT
360   integrator, the integration of the equations of motion is carried
361   out in a velocity-Verlet style 2 part algorithm:
362  
# Line 473 | Line 469 | Bond constraints are applied at the end of both the {\
469   P_{\mathrm{target}} \mathcal{V}(t).
470   \end{equation}
471  
476 – Bond constraints are applied at the end of both the {\tt moveA} and
477 – {\tt moveB} portions of the algorithm.  Details on the constraint
478 – algorithms are given in section \ref{oopseSec:rattle}.
479 –
472   \subsection{\label{methodSection:NPTf}Constant-pressure integration with a
473   flexible box (NPTf)}
474  
# Line 608 | Line 600 | assume non-orthorhombic geometries.
600  
601   \subsection{\label{methodSection:otherSpecialEnsembles}Other Special Ensembles}
602  
603 < \subsubsection{\label{methodSection:NPAT}NPAT Ensemble}
603 > \subsubsection{\label{methodSection:NPAT}\textbf{NPAT Ensemble}}
604  
605   A comprehensive understanding of structure¨Cfunction relations of
606   biological membrane system ultimately relies on structure and
# Line 631 | Line 623 | described for the NPTi integrator.
623   Note that the iterative schemes for NPAT are identical to those
624   described for the NPTi integrator.
625  
626 < \subsubsection{\label{methodSection:NPrT}NP$\gamma$T Ensemble}
626 > \subsubsection{\label{methodSection:NPrT}\textbf{NP$\gamma$T Ensemble}}
627  
628   Theoretically, the surface tension $\gamma$ of a stress free
629   membrane system should be zero since its surface free energy $G$ is
# Line 757 | Line 749 | F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\part
749   \begin{equation}
750   F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}=
751      -k_{\text{Harmonic}}(z(t)-z_{\text{cons}}).
760 – \end{equation}
761 –
762 –
763 – \section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies}
764 –
765 – %\subsection{\label{methodSection:temperature}Temperature Control}
766 –
767 – %\subsection{\label{methodSection:pressureControl}Pressure Control}
768 –
769 – %\section{\label{methodSection:hydrodynamics}Hydrodynamics}
770 –
771 – %applications of langevin dynamics
772 – As an excellent alternative to newtonian dynamics, Langevin
773 – dynamics, which mimics a simple heat bath with stochastic and
774 – dissipative forces, has been applied in a variety of studies. The
775 – stochastic treatment of the solvent enables us to carry out
776 – substantially longer time simulation. Implicit solvent Langevin
777 – dynamics simulation of met-enkephalin not only outperforms explicit
778 – solvent simulation on computation efficiency, but also agrees very
779 – well with explicit solvent simulation on dynamics
780 – properties\cite{Shen2002}. Recently, applying Langevin dynamics with
781 – UNRES model, Liow and his coworkers suggest that protein folding
782 – pathways can be possibly exploited within a reasonable amount of
783 – time\cite{Liwo2005}. The stochastic nature of the Langevin dynamics
784 – also enhances the sampling of the system and increases the
785 – probability of crossing energy barrier\cite{Banerjee2004, Cui2003}.
786 – Combining Langevin dynamics with Kramers's theory, Klimov and
787 – Thirumalai identified the free-energy barrier by studying the
788 – viscosity dependence of the protein folding rates\cite{Klimov1997}.
789 – In order to account for solvent induced interactions missing from
790 – implicit solvent model, Kaya incorporated desolvation free energy
791 – barrier into implicit coarse-grained solvent model in protein
792 – folding/unfolding study and discovered a higher free energy barrier
793 – between the native and denatured states. Because of its stability
794 – against noise, Langevin dynamics is very suitable for studying
795 – remagnetization processes in various
796 – systems\cite{Garcia-Palacios1998,Berkov2002,Denisov2003}. For
797 – instance, the oscillation power spectrum of nanoparticles from
798 – Langevin dynamics simulation has the same peak frequencies for
799 – different wave vectors,which recovers the property of magnetic
800 – excitations in small finite structures\cite{Berkov2005a}. In an
801 – attempt to reduce the computational cost of simulation, multiple
802 – time stepping (MTS) methods have been introduced and have been of
803 – great interest to macromolecule and protein
804 – community\cite{Tuckerman1992}. Relying on the observation that
805 – forces between distant atoms generally demonstrate slower
806 – fluctuations than forces between close atoms, MTS method are
807 – generally implemented by evaluating the slowly fluctuating forces
808 – less frequently than the fast ones. Unfortunately, nonlinear
809 – instability resulting from increasing timestep in MTS simulation
810 – have became a critical obstruction preventing the long time
811 – simulation. Due to the coupling to the heat bath, Langevin dynamics
812 – has been shown to be able to damp out the resonance artifact more
813 – efficiently\cite{Sandu1999}.
814 –
815 – %review rigid body dynamics
816 – Rigid bodies are frequently involved in the modeling of different
817 – areas, from engineering, physics, to chemistry. For example,
818 – missiles and vehicle are usually modeled by rigid bodies.  The
819 – movement of the objects in 3D gaming engine or other physics
820 – simulator is governed by the rigid body dynamics. In molecular
821 – simulation, rigid body is used to simplify the model in
822 – protein-protein docking study\cite{Gray2003}.
823 –
824 – It is very important to develop stable and efficient methods to
825 – integrate the equations of motion of orientational degrees of
826 – freedom. Euler angles are the nature choice to describe the
827 – rotational degrees of freedom. However, due to its singularity, the
828 – numerical integration of corresponding equations of motion is very
829 – inefficient and inaccurate. Although an alternative integrator using
830 – different sets of Euler angles can overcome this
831 – difficulty\cite{Ryckaert1977, Andersen1983}, the computational
832 – penalty and the lost of angular momentum conservation still remain.
833 – In 1977, a singularity free representation utilizing quaternions was
834 – developed by Evans\cite{Evans1977}. Unfortunately, this approach
835 – suffer from the nonseparable Hamiltonian resulted from quaternion
836 – representation, which prevents the symplectic algorithm to be
837 – utilized. Another different approach is to apply holonomic
838 – constraints to the atoms belonging to the rigid
839 – body\cite{Barojas1973}. Each atom moves independently under the
840 – normal forces deriving from potential energy and constraint forces
841 – which are used to guarantee the rigidness. However, due to their
842 – iterative nature, SHAKE and Rattle algorithm converge very slowly
843 – when the number of constraint increases.
844 –
845 – The break through in geometric literature suggests that, in order to
846 – develop a long-term integration scheme, one should preserve the
847 – geometric structure of the flow. Matubayasi and Nakahara developed a
848 – time-reversible integrator for rigid bodies in quaternion
849 – representation. Although it is not symplectic, this integrator still
850 – demonstrates a better long-time energy conservation than traditional
851 – methods because of the time-reversible nature. Extending
852 – Trotter-Suzuki to general system with a flat phase space, Miller and
853 – his colleagues devised an novel symplectic, time-reversible and
854 – volume-preserving integrator in quaternion representation. However,
855 – all of the integrators in quaternion representation suffer from the
856 – computational penalty of constructing a rotation matrix from
857 – quaternions to evolve coordinates and velocities at every time step.
858 – An alternative integration scheme utilizing rotation matrix directly
859 – is RSHAKE\cite{Kol1997}, in which a conjugate momentum to rotation
860 – matrix is introduced to re-formulate the Hamiltonian's equation and
861 – the Hamiltonian is evolved in a constraint manifold by iteratively
862 – satisfying the orthogonality constraint. However, RSHAKE is
863 – inefficient because of the iterative procedure. An extremely
864 – efficient integration scheme in rotation matrix representation,
865 – which also preserves the same structural properties of the
866 – Hamiltonian flow as Miller's integrator, is proposed by Dullweber,
867 – Leimkuhler and McLachlan (DLM)\cite{Dullweber1997}.
868 –
869 – %review langevin/browninan dynamics for arbitrarily shaped rigid body
870 – Combining Langevin or Brownian dynamics with rigid body dynamics,
871 – one can study the slow processes in biomolecular systems. Modeling
872 – the DNA as a chain of rigid spheres beads, which subject to harmonic
873 – potentials as well as excluded volume potentials, Mielke and his
874 – coworkers discover rapid superhelical stress generations from the
875 – stochastic simulation of twin supercoiling DNA with response to
876 – induced torques\cite{Mielke2004}. Membrane fusion is another key
877 – biological process which controls a variety of physiological
878 – functions, such as release of neurotransmitters \textit{etc}. A
879 – typical fusion event happens on the time scale of millisecond, which
880 – is impracticable to study using all atomistic model with newtonian
881 – mechanics. With the help of coarse-grained rigid body model and
882 – stochastic dynamics, the fusion pathways were exploited by many
883 – researchers\cite{Noguchi2001,Noguchi2002,Shillcock2005}. Due to the
884 – difficulty of numerical integration of anisotropy rotation, most of
885 – the rigid body models are simply modeled by sphere, cylinder,
886 – ellipsoid or other regular shapes in stochastic simulations. In an
887 – effort to account for the diffusion anisotropy of the arbitrary
888 – particles, Fernandes and de la Torre improved the original Brownian
889 – dynamics simulation algorithm\cite{Ermak1978,Allison1991} by
890 – incorporating a generalized $6\times6$ diffusion tensor and
891 – introducing a simple rotation evolution scheme consisting of three
892 – consecutive rotations\cite{Fernandes2002}. Unfortunately, unexpected
893 – error and bias are introduced into the system due to the arbitrary
894 – order of applying the noncommuting rotation
895 – operators\cite{Beard2003}. Based on the observation the momentum
896 – relaxation time is much less than the time step, one may ignore the
897 – inertia in Brownian dynamics. However, assumption of the zero
898 – average acceleration is not always true for cooperative motion which
899 – is common in protein motion. An inertial Brownian dynamics (IBD) was
900 – proposed to address this issue by adding an inertial correction
901 – term\cite{Beard2001}. As a complement to IBD which has a lower bound
902 – in time step because of the inertial relaxation time, long-time-step
903 – inertial dynamics (LTID) can be used to investigate the inertial
904 – behavior of the polymer segments in low friction
905 – regime\cite{Beard2001}. LTID can also deal with the rotational
906 – dynamics for nonskew bodies without translation-rotation coupling by
907 – separating the translation and rotation motion and taking advantage
908 – of the analytical solution of hydrodynamics properties. However,
909 – typical nonskew bodies like cylinder and ellipsoid are inadequate to
910 – represent most complex macromolecule assemblies. These intricate
911 – molecules have been represented by a set of beads and their
912 – hydrodynamics properties can be calculated using variant
913 – hydrodynamic interaction tensors.
914 –
915 – The goal of the present work is to develop a Langevin dynamics
916 – algorithm for arbitrary rigid particles by integrating the accurate
917 – estimation of friction tensor from hydrodynamics theory into the
918 – sophisticated rigid body dynamics.
919 –
920 –
921 – \subsection{Friction Tensor}
922 –
923 – For an arbitrary rigid body moves in a fluid, it may experience
924 – friction force $f_r$ or friction torque $\tau _r$ along the opposite
925 – direction of the velocity $v$ or angular velocity $\omega$ at
926 – arbitrary origin $P$,
927 – \begin{equation}
928 – \left( \begin{array}{l}
929 – f_r  \\
930 – \tau _r  \\
931 – \end{array} \right) =  - \left( {\begin{array}{*{20}c}
932 –   {\Xi _{P,t} } & {\Xi _{P,c}^T }  \\
933 –   {\Xi _{P,c} } & {\Xi _{P,r} }  \\
934 – \end{array}} \right)\left( \begin{array}{l}
935 – \nu  \\
936 – \omega  \\
937 – \end{array} \right)
938 – \end{equation}
939 – where $\Xi _{P,t}t$ is the translation friction tensor, $\Xi _{P,r}$
940 – is the rotational friction tensor and $\Xi _{P,c}$ is the
941 – translation-rotation coupling tensor. The procedure of calculating
942 – friction tensor using hydrodynamic tensor and comparison between
943 – bead model and shell model were elaborated by Carrasco \textit{et
944 – al}\cite{Carrasco1999}. An important property of the friction tensor
945 – is that the translational friction tensor is independent of origin
946 – while the rotational and coupling are sensitive to the choice of the
947 – origin \cite{Brenner1967}, which can be described by
948 – \begin{equation}
949 – \begin{array}{c}
950 – \Xi _{P,t}  = \Xi _{O,t}  = \Xi _t  \\
951 – \Xi _{P,c}  = \Xi _{O,c}  - r_{OP}  \times \Xi _t  \\
952 – \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  \\
953 – \end{array}
954 – \end{equation}
955 – Where $O$ is another origin and $r_{OP}$ is the vector joining $O$
956 – and $P$. It is also worthy of mention that both of translational and
957 – rotational frictional tensors are always symmetric. In contrast,
958 – coupling tensor is only symmetric at center of reaction:
959 – \begin{equation}
960 – \Xi _{R,c}  = \Xi _{R,c}^T
961 – \end{equation}
962 – The proper location for applying friction force is the center of
963 – reaction, at which the trace of rotational resistance tensor reaches
964 – minimum.
965 –
966 – \subsection{Rigid body dynamics}
967 –
968 – The Hamiltonian of rigid body can be separated in terms of potential
969 – energy $V(r,A)$ and kinetic energy $T(p,\pi)$,
970 – \[
971 – H = V(r,A) + T(v,\pi )
972 – \]
973 – A second-order symplectic method is now obtained by the composition
974 – of the flow maps,
975 – \[
976 – \varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
977 – _{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
978 – \]
979 – Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
980 – sub-flows which corresponding to force and torque respectively,
981 – \[
982 – \varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
983 – _{\Delta t/2,\tau }.
984 – \]
985 – Since the associated operators of $\varphi _{\Delta t/2,F} $ and
986 – $\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
987 – order inside $\varphi _{\Delta t/2,V}$ does not matter.
988 –
989 – Furthermore, kinetic potential can be separated to translational
990 – kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
991 – \begin{equation}
992 – T(p,\pi ) =T^t (p) + T^r (\pi ).
752   \end{equation}
994 – where $ T^t (p) = \frac{1}{2}v^T m v $ and $T^r (\pi )$ is defined
995 – by \ref{introEquation:rotationalKineticRB}. Therefore, the
996 – corresponding flow maps are given by
997 – \[
998 – \varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
999 – _{\Delta t,T^r }.
1000 – \]
1001 – The free rigid body is an example of Lie-Poisson system with
1002 – Hamiltonian function
1003 – \begin{equation}
1004 – T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
1005 – \label{introEquation:rotationalKineticRB}
1006 – \end{equation}
1007 – where $T_i^r (\pi _i ) = \frac{{\pi _i ^2 }}{{2I_i }}$ and
1008 – Lie-Poisson structure matrix,
1009 – \begin{equation}
1010 – J(\pi ) = \left( {\begin{array}{*{20}c}
1011 –   0 & {\pi _3 } & { - \pi _2 }  \\
1012 –   { - \pi _3 } & 0 & {\pi _1 }  \\
1013 –   {\pi _2 } & { - \pi _1 } & 0  \\
1014 – \end{array}} \right)
1015 – \end{equation}
1016 – Thus, the dynamics of free rigid body is governed by
1017 – \begin{equation}
1018 – \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
1019 – \end{equation}
1020 – One may notice that each $T_i^r$ in Equation
1021 – \ref{introEquation:rotationalKineticRB} can be solved exactly. For
1022 – instance, the equations of motion due to $T_1^r$ are given by
1023 – \begin{equation}
1024 – \frac{d}{{dt}}\pi  = R_1 \pi ,\frac{d}{{dt}}A = AR_1
1025 – \label{introEqaution:RBMotionSingleTerm}
1026 – \end{equation}
1027 – where
1028 – \[ R_1  = \left( {\begin{array}{*{20}c}
1029 –   0 & 0 & 0  \\
1030 –   0 & 0 & {\pi _1 }  \\
1031 –   0 & { - \pi _1 } & 0  \\
1032 – \end{array}} \right).
1033 – \]
1034 – The solutions of Equation \ref{introEqaution:RBMotionSingleTerm} is
1035 – \[
1036 – \pi (\Delta t) = e^{\Delta tR_1 } \pi (0),A(\Delta t) =
1037 – A(0)e^{\Delta tR_1 }
1038 – \]
1039 – with
1040 – \[
1041 – e^{\Delta tR_1 }  = \left( {\begin{array}{*{20}c}
1042 –   0 & 0 & 0  \\
1043 –   0 & {\cos \theta _1 } & {\sin \theta _1 }  \\
1044 –   0 & { - \sin \theta _1 } & {\cos \theta _1 }  \\
1045 – \end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
1046 – \]
1047 – To reduce the cost of computing expensive functions in $e^{\Delta
1048 – tR_1 }$, we can use Cayley transformation,
1049 – \[
1050 – e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
1051 – )
1052 – \]
1053 – The flow maps for $T_2^r$ and $T_3^r$ can be found in the same
1054 – manner.
1055 –
1056 – In order to construct a second-order symplectic method, we split the
1057 – angular kinetic Hamiltonian function into five terms
1058 – \[
1059 – T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
1060 – ) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
1061 – (\pi _1 )
1062 – \].
1063 – Concatenating flows corresponding to these five terms, we can obtain
1064 – the flow map for free rigid body,
1065 – \[
1066 – \varphi _{\Delta t,T^r }  = \varphi _{\Delta t/2,\pi _1 }  \circ
1067 – \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }
1068 – \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
1069 – _1 }.
1070 – \]
1071 –
1072 – The equations of motion corresponding to potential energy and
1073 – kinetic energy are listed in the below table,
1074 – \begin{center}
1075 – \begin{tabular}{|l|l|}
1076 –  \hline
1077 –  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
1078 –  Potential & Kinetic \\
1079 –  $\frac{{dq}}{{dt}} = \frac{p}{m}$ & $\frac{d}{{dt}}q = p$ \\
1080 –  $\frac{d}{{dt}}p =  - \frac{{\partial V}}{{\partial q}}$ & $ \frac{d}{{dt}}p = 0$ \\
1081 –  $\frac{d}{{dt}}Q = 0$ & $ \frac{d}{{dt}}Q = Qskew(I^{ - 1} j)$ \\
1082 –  $ \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$\\
1083 –  \hline
1084 – \end{tabular}
1085 – \end{center}
1086 –
1087 – Finally, we obtain the overall symplectic flow maps for free moving
1088 – rigid body
1089 – \begin{align*}
1090 – \varphi _{\Delta t}  = &\varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau } \circ  \\
1091 –  &\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  \\
1092 –  &\varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1093 – \label{introEquation:overallRBFlowMaps}
1094 – \end{align*}
1095 –
1096 – \subsection{Langevin dynamics for rigid particles of arbitrary shape}
1097 –
1098 – Consider a Langevin equation of motions in generalized coordinates
1099 – \begin{equation}
1100 – M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)}  + F_{r,i} (t)
1101 – \label{LDGeneralizedForm}
1102 – \end{equation}
1103 – where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
1104 – and moment of inertial) matrix and $V_i$ is a generalized velocity,
1105 – $V_i = V_i(v_i,\omega _i)$. The right side of Eq.
1106 – (\ref{LDGeneralizedForm}) consists of three generalized forces in
1107 – lab-fixed frame, systematic force $F_{s,i}$, dissipative force
1108 – $F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
1109 – system in Newtownian mechanics typically refers to lab-fixed frame,
1110 – it is also convenient to handle the rotation of rigid body in
1111 – body-fixed frame. Thus the friction and random forces are calculated
1112 – in body-fixed frame and converted back to lab-fixed frame by:
1113 – \[
1114 – \begin{array}{l}
1115 – F_{f,i}^l (t) = A^T F_{f,i}^b (t), \\
1116 – F_{r,i}^l (t) = A^T F_{r,i}^b (t) \\
1117 – \end{array}.
1118 – \]
1119 – Here, the body-fixed friction force $F_{r,i}^b$ is proportional to
1120 – the body-fixed velocity at center of resistance $v_{R,i}^b$ and
1121 – angular velocity $\omega _i$,
1122 – \begin{equation}
1123 – F_{r,i}^b (t) = \left( \begin{array}{l}
1124 – f_{r,i}^b (t) \\
1125 – \tau _{r,i}^b (t) \\
1126 – \end{array} \right) =  - \left( {\begin{array}{*{20}c}
1127 –   {\Xi _{R,t} } & {\Xi _{R,c}^T }  \\
1128 –   {\Xi _{R,c} } & {\Xi _{R,r} }  \\
1129 – \end{array}} \right)\left( \begin{array}{l}
1130 – v_{R,i}^b (t) \\
1131 – \omega _i (t) \\
1132 – \end{array} \right),
1133 – \end{equation}
1134 – while the random force $F_{r,i}^l$ is a Gaussian stochastic variable
1135 – with zero mean and variance
1136 – \begin{equation}
1137 – \left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
1138 – \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
1139 – 2k_B T\Xi _R \delta (t - t').
1140 – \end{equation}
1141 – The equation of motion for $v_i$ can be written as
1142 – \begin{equation}
1143 – m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
1144 – f_{r,i}^l (t)
1145 – \end{equation}
1146 – Since the frictional force is applied at the center of resistance
1147 – which generally does not coincide with the center of mass, an extra
1148 – torque is exerted at the center of mass. Thus, the net body-fixed
1149 – frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is
1150 – given by
1151 – \begin{equation}
1152 – \tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b
1153 – \end{equation}
1154 – where $r_{MR}$ is the vector from the center of mass to the center
1155 – of the resistance. Instead of integrating angular velocity in
1156 – lab-fixed frame, we consider the equation of motion of angular
1157 – momentum in body-fixed frame
1158 – \begin{equation}
1159 – \dot \pi _i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b
1160 – (t) + \tau _{r,i}^b(t)
1161 – \end{equation}
1162 –
1163 – Embedding the friction terms into force and torque, one can
1164 – integrate the langevin equations of motion for rigid body of
1165 – arbitrary shape in a velocity-Verlet style 2-part algorithm, where
1166 – $h= \delta t$:
1167 –
1168 – {\tt part one:}
1169 – \begin{align*}
1170 – v_i (t + h/2) &\leftarrow v_i (t) + \frac{{hf_{t,i}^l (t)}}{{2m_i }} \\
1171 – \pi _i (t + h/2) &\leftarrow \pi _i (t) + \frac{{h\tau _{t,i}^b (t)}}{2} \\
1172 – r_i (t + h) &\leftarrow r_i (t) + hv_i (t + h/2) \\
1173 – A_i (t + h) &\leftarrow rotate\left( {h\pi _i (t + h/2)I^{ - 1} } \right) \\
1174 – \end{align*}
1175 – In this context, the $\mathrm{rotate}$ function is the reversible
1176 – product of five consecutive body-fixed rotations,
1177 – \begin{equation}
1178 – \mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
1179 – \mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y
1180 – / 2) \cdot \mathsf{G}_x(a_x /2),
1181 – \end{equation}
1182 – where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
1183 – rotates both the rotation matrix ($\mathsf{A}$) and the body-fixed
1184 – angular momentum ($\pi$) by an angle $\theta$ around body-fixed axis
1185 – $\alpha$,
1186 – \begin{equation}
1187 – \mathsf{G}_\alpha( \theta ) = \left\{
1188 – \begin{array}{lcl}
1189 – \mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
1190 – {\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
1191 – j}(0).
1192 – \end{array}
1193 – \right.
1194 – \end{equation}
1195 – $\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis
1196 – rotation matrix.  For example, in the small-angle limit, the
1197 – rotation matrix around the body-fixed x-axis can be approximated as
1198 – \begin{equation}
1199 – \mathsf{R}_x(\theta) \approx \left(
1200 – \begin{array}{ccc}
1201 – 1 & 0 & 0 \\
1202 – 0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
1203 – \theta^2 / 4} \\
1204 – 0 & \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 +
1205 – \theta^2 / 4}
1206 – \end{array}
1207 – \right).
1208 – \end{equation}
1209 – All other rotations follow in a straightforward manner.
1210 –
1211 – After the first part of the propagation, the friction and random
1212 – forces are generated at the center of resistance in body-fixed frame
1213 – and converted back into lab-fixed frame
1214 – \[
1215 – f_{t,i}^l (t + h) =  - \left( {\frac{{\partial V}}{{\partial r_i }}}
1216 – \right)_{r_i (t + h)}  + A_i^T (t + h)[F_{f,i}^b (t + h) + F_{r,i}^b
1217 – (t + h)],
1218 – \]
1219 – while the system torque in lab-fixed frame is transformed into
1220 – body-fixed frame,
1221 – \[
1222 – \tau _{t,i}^b (t + h) = A\tau _{s,i}^l (t) + \tau _{n,i}^b (t) +
1223 – \tau _{r,i}^b (t).
1224 – \]
1225 – Once the forces and torques have been obtained at the new time step,
1226 – the velocities can be advanced to the same time value.
1227 –
1228 – {\tt part two:}
1229 – \begin{align*}
1230 – v_i (t) &\leftarrow v_i (t + h/2) + \frac{{hf_{t,i}^l (t + h)}}{{2m_i }} \\
1231 – \pi _i (t) &\leftarrow \pi _i (t + h/2) + \frac{{h\tau _{t,i}^b (t + h)}}{2} \\
1232 – \end{align*}
1233 –
1234 – \subsection{Results and discussion}

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines