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 2850 by tim, Sun Jun 11 01:55:48 2006 UTC vs.
Revision 2851 by tim, Sun Jun 11 02:06:01 2006 UTC

# Line 756 | Line 756 | F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\part
756   \begin{equation}
757   F_{z_{\text{Harmonic}}}(t)=-\frac{\partial U(t)}{\partial z(t)}=
758      -k_{\text{Harmonic}}(z(t)-z_{\text{cons}}).
759 – \end{equation}
760 –
761 –
762 – \section{\label{methodSection:langevin}Integrators for Langevin Dynamics of Rigid Bodies}
763 –
764 – %\subsection{\label{methodSection:temperature}Temperature Control}
765 –
766 – %\subsection{\label{methodSection:pressureControl}Pressure Control}
767 –
768 – %\section{\label{methodSection:hydrodynamics}Hydrodynamics}
769 –
770 – %applications of langevin dynamics
771 – As an excellent alternative to newtonian dynamics, Langevin
772 – dynamics, which mimics a simple heat bath with stochastic and
773 – dissipative forces, has been applied in a variety of studies. The
774 – stochastic treatment of the solvent enables us to carry out
775 – substantially longer time simulation. Implicit solvent Langevin
776 – dynamics simulation of met-enkephalin not only outperforms explicit
777 – solvent simulation on computation efficiency, but also agrees very
778 – well with explicit solvent simulation on dynamics
779 – properties\cite{Shen2002}. Recently, applying Langevin dynamics with
780 – UNRES model, Liow and his coworkers suggest that protein folding
781 – pathways can be possibly exploited within a reasonable amount of
782 – time\cite{Liwo2005}. The stochastic nature of the Langevin dynamics
783 – also enhances the sampling of the system and increases the
784 – probability of crossing energy barrier\cite{Banerjee2004, Cui2003}.
785 – Combining Langevin dynamics with Kramers's theory, Klimov and
786 – Thirumalai identified the free-energy barrier by studying the
787 – viscosity dependence of the protein folding rates\cite{Klimov1997}.
788 – In order to account for solvent induced interactions missing from
789 – implicit solvent model, Kaya incorporated desolvation free energy
790 – barrier into implicit coarse-grained solvent model in protein
791 – folding/unfolding study and discovered a higher free energy barrier
792 – between the native and denatured states. Because of its stability
793 – against noise, Langevin dynamics is very suitable for studying
794 – remagnetization processes in various
795 – systems\cite{Palacios1998,Berkov2002,Denisov2003}. For instance, the
796 – oscillation power spectrum of nanoparticles from Langevin dynamics
797 – simulation has the same peak frequencies for different wave
798 – vectors,which recovers the property of magnetic excitations in small
799 – finite structures\cite{Berkov2005a}. In an attempt to reduce the
800 – computational cost of simulation, multiple time stepping (MTS)
801 – methods have been introduced and have been of great interest to
802 – macromolecule and protein community\cite{Tuckerman1992}. Relying on
803 – the observation that forces between distant atoms generally
804 – demonstrate slower fluctuations than forces between close atoms, MTS
805 – method are generally implemented by evaluating the slowly
806 – fluctuating forces less frequently than the fast ones.
807 – Unfortunately, nonlinear instability resulting from increasing
808 – timestep in MTS simulation have became a critical obstruction
809 – preventing the long time simulation. Due to the coupling to the heat
810 – bath, Langevin dynamics has been shown to be able to damp out the
811 – resonance artifact more efficiently\cite{Sandu1999}.
812 –
813 – It is very important to develop stable and efficient methods to
814 – integrate the equations of motion of orientational degrees of
815 – freedom. Euler angles are the nature choice to describe the
816 – rotational degrees of freedom. However, due to its singularity, the
817 – numerical integration of corresponding equations of motion is very
818 – inefficient and inaccurate. Although an alternative integrator using
819 – different sets of Euler angles can overcome this
820 – difficulty\cite{Ryckaert1977, Andersen1983}, the computational
821 – penalty and the lost of angular momentum conservation still remain.
822 – In 1977, a singularity free representation utilizing quaternions was
823 – developed by Evans\cite{Evans1977}. Unfortunately, this approach
824 – suffer from the nonseparable Hamiltonian resulted from quaternion
825 – representation, which prevents the symplectic algorithm to be
826 – utilized. Another different approach is to apply holonomic
827 – constraints to the atoms belonging to the rigid
828 – body\cite{Barojas1973}. Each atom moves independently under the
829 – normal forces deriving from potential energy and constraint forces
830 – which are used to guarantee the rigidness. However, due to their
831 – iterative nature, SHAKE and Rattle algorithm converge very slowly
832 – when the number of constraint increases.
833 –
834 – The break through in geometric literature suggests that, in order to
835 – develop a long-term integration scheme, one should preserve the
836 – geometric structure of the flow. Matubayasi developed a
837 – time-reversible integrator for rigid bodies in quaternion
838 – representation. Although it is not symplectic, this integrator still
839 – demonstrates a better long-time energy conservation than traditional
840 – methods because of the time-reversible nature. Extending
841 – Trotter-Suzuki to general system with a flat phase space, Miller and
842 – his colleagues devised an novel symplectic, time-reversible and
843 – volume-preserving integrator in quaternion representation. However,
844 – all of the integrators in quaternion representation suffer from the
845 – computational penalty of constructing a rotation matrix from
846 – quaternions to evolve coordinates and velocities at every time step.
847 – An alternative integration scheme utilizing rotation matrix directly
848 – is RSHAKE\cite{Kol1997}, in which a conjugate momentum to rotation
849 – matrix is introduced to re-formulate the Hamiltonian's equation and
850 – the Hamiltonian is evolved in a constraint manifold by iteratively
851 – satisfying the orthogonality constraint. However, RSHAKE is
852 – inefficient because of the iterative procedure. An extremely
853 – efficient integration scheme in rotation matrix representation,
854 – which also preserves the same structural properties of the
855 – Hamiltonian flow as Miller's integrator, is proposed by Dullweber,
856 – Leimkuhler and McLachlan (DLM)\cite{Dullweber1997}.
857 –
858 – %review langevin/browninan dynamics for arbitrarily shaped rigid body
859 – Combining Langevin or Brownian dynamics with rigid body dynamics,
860 – one can study the slow processes in biomolecular systems. Modeling
861 – the DNA as a chain of rigid spheres beads, which subject to harmonic
862 – potentials as well as excluded volume potentials, Mielke and his
863 – coworkers discover rapid superhelical stress generations from the
864 – stochastic simulation of twin supercoiling DNA with response to
865 – induced torques\cite{Mielke2004}. Membrane fusion is another key
866 – biological process which controls a variety of physiological
867 – functions, such as release of neurotransmitters \textit{etc}. A
868 – typical fusion event happens on the time scale of millisecond, which
869 – is impracticable to study using all atomistic model with newtonian
870 – mechanics. With the help of coarse-grained rigid body model and
871 – stochastic dynamics, the fusion pathways were exploited by many
872 – researchers\cite{Noguchi2001,Noguchi2002,Shillcock2005}. Due to the
873 – difficulty of numerical integration of anisotropy rotation, most of
874 – the rigid body models are simply modeled by sphere, cylinder,
875 – ellipsoid or other regular shapes in stochastic simulations. In an
876 – effort to account for the diffusion anisotropy of the arbitrary
877 – particles, Fernandes and de la Torre improved the original Brownian
878 – dynamics simulation algorithm\cite{Ermak1978,Allison1991} by
879 – incorporating a generalized $6\times6$ diffusion tensor and
880 – introducing a simple rotation evolution scheme consisting of three
881 – consecutive rotations\cite{Fernandes2002}. Unfortunately, unexpected
882 – error and bias are introduced into the system due to the arbitrary
883 – order of applying the noncommuting rotation
884 – operators\cite{Beard2003}. Based on the observation the momentum
885 – relaxation time is much less than the time step, one may ignore the
886 – inertia in Brownian dynamics. However, assumption of the zero
887 – average acceleration is not always true for cooperative motion which
888 – is common in protein motion. An inertial Brownian dynamics (IBD) was
889 – proposed to address this issue by adding an inertial correction
890 – term\cite{Beard2003}. As a complement to IBD which has a lower bound
891 – in time step because of the inertial relaxation time, long-time-step
892 – inertial dynamics (LTID) can be used to investigate the inertial
893 – behavior of the polymer segments in low friction
894 – regime\cite{Beard2003}. LTID can also deal with the rotational
895 – dynamics for nonskew bodies without translation-rotation coupling by
896 – separating the translation and rotation motion and taking advantage
897 – of the analytical solution of hydrodynamics properties. However,
898 – typical nonskew bodies like cylinder and ellipsoid are inadequate to
899 – represent most complex macromolecule assemblies. These intricate
900 – molecules have been represented by a set of beads and their
901 – hydrodynamics properties can be calculated using variant
902 – hydrodynamic interaction tensors.
903 –
904 – The goal of the present work is to develop a Langevin dynamics
905 – algorithm for arbitrary rigid particles by integrating the accurate
906 – estimation of friction tensor from hydrodynamics theory into the
907 – sophisticated rigid body dynamics.
908 –
909 – \subsection{\label{introSection:frictionTensor} Friction Tensor}
910 – Theoretically, the friction kernel can be determined using velocity
911 – autocorrelation function. However, this approach become impractical
912 – when the system become more and more complicate. Instead, various
913 – approaches based on hydrodynamics have been developed to calculate
914 – the friction coefficients. The friction effect is isotropic in
915 – Equation, $\zeta$ can be taken as a scalar. In general, friction
916 – tensor $\Xi$ is a $6\times 6$ matrix given by
917 – \[
918 – \Xi  = \left( {\begin{array}{*{20}c}
919 –   {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
920 –   {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
921 – \end{array}} \right).
922 – \]
923 – Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are translational friction
924 – tensor and rotational resistance (friction) tensor respectively,
925 – while ${\Xi^{tr} }$ is translation-rotation coupling tensor and $
926 – {\Xi^{rt} }$ is rotation-translation coupling tensor. When a
927 – particle moves in a fluid, it may experience friction force or
928 – torque along the opposite direction of the velocity or angular
929 – velocity,
930 – \[
931 – \left( \begin{array}{l}
932 – F_R  \\
933 – \tau _R  \\
934 – \end{array} \right) =  - \left( {\begin{array}{*{20}c}
935 –   {\Xi ^{tt} } & {\Xi ^{rt} }  \\
936 –   {\Xi ^{tr} } & {\Xi ^{rr} }  \\
937 – \end{array}} \right)\left( \begin{array}{l}
938 – v \\
939 – w \\
940 – \end{array} \right)
941 – \]
942 – where $F_r$ is the friction force and $\tau _R$ is the friction
943 – toque.
944 –
945 – \subsubsection{\label{introSection:resistanceTensorRegular}\textbf{The Resistance Tensor for Regular Shape}}
946 –
947 – For a spherical particle, the translational and rotational friction
948 – constant can be calculated from Stoke's law,
949 – \[
950 – \Xi ^{tt}  = \left( {\begin{array}{*{20}c}
951 –   {6\pi \eta R} & 0 & 0  \\
952 –   0 & {6\pi \eta R} & 0  \\
953 –   0 & 0 & {6\pi \eta R}  \\
954 – \end{array}} \right)
955 – \]
956 – and
957 – \[
958 – \Xi ^{rr}  = \left( {\begin{array}{*{20}c}
959 –   {8\pi \eta R^3 } & 0 & 0  \\
960 –   0 & {8\pi \eta R^3 } & 0  \\
961 –   0 & 0 & {8\pi \eta R^3 }  \\
962 – \end{array}} \right)
963 – \]
964 – where $\eta$ is the viscosity of the solvent and $R$ is the
965 – hydrodynamics radius.
966 –
967 – Other non-spherical shape, such as cylinder and ellipsoid
968 – \textit{etc}, are widely used as reference for developing new
969 – hydrodynamics theory, because their properties can be calculated
970 – exactly. In 1936, Perrin extended Stokes's law to general ellipsoid,
971 – also called a triaxial ellipsoid, which is given in Cartesian
972 – coordinates by\cite{Perrin1934, Perrin1936}
973 – \[
974 – \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
975 – }} = 1
976 – \]
977 – where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
978 – due to the complexity of the elliptic integral, only the ellipsoid
979 – with the restriction of two axes having to be equal, \textit{i.e.}
980 – prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
981 – exactly. Introducing an elliptic integral parameter $S$ for prolate,
982 – \[
983 – S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
984 – } }}{b},
985 – \]
986 – and oblate,
987 – \[
988 – S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
989 – }}{a}
990 – \],
991 – one can write down the translational and rotational resistance
992 – tensors
993 – \[
994 – \begin{array}{l}
995 – \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}} \\
996 – \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S + 2a}} \\
997 – \end{array},
998 – \]
999 – and
1000 – \[
1001 – \begin{array}{l}
1002 – \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}} \\
1003 – \Xi _b^{rr}  = \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}} \\
1004 – \end{array}.
1005 – \]
1006 –
1007 – \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}\textbf{The Resistance Tensor for Arbitrary Shape}}
1008 –
1009 – Unlike spherical and other regular shaped molecules, there is not
1010 – analytical solution for friction tensor of any arbitrary shaped
1011 – rigid molecules. The ellipsoid of revolution model and general
1012 – triaxial ellipsoid model have been used to approximate the
1013 – hydrodynamic properties of rigid bodies. However, since the mapping
1014 – from all possible ellipsoidal space, $r$-space, to all possible
1015 – combination of rotational diffusion coefficients, $D$-space is not
1016 – unique\cite{Wegener1979} as well as the intrinsic coupling between
1017 – translational and rotational motion of rigid body, general ellipsoid
1018 – is not always suitable for modeling arbitrarily shaped rigid
1019 – molecule. A number of studies have been devoted to determine the
1020 – friction tensor for irregularly shaped rigid bodies using more
1021 – advanced method where the molecule of interest was modeled by
1022 – combinations of spheres(beads)\cite{Carrasco1999} and the
1023 – hydrodynamics properties of the molecule can be calculated using the
1024 – hydrodynamic interaction tensor. Let us consider a rigid assembly of
1025 – $N$ beads immersed in a continuous medium. Due to hydrodynamics
1026 – interaction, the ``net'' velocity of $i$th bead, $v'_i$ is different
1027 – than its unperturbed velocity $v_i$,
1028 – \[
1029 – v'_i  = v_i  - \sum\limits_{j \ne i} {T_{ij} F_j }
1030 – \]
1031 – where $F_i$ is the frictional force, and $T_{ij}$ is the
1032 – hydrodynamic interaction tensor. The friction force of $i$th bead is
1033 – proportional to its ``net'' velocity
1034 – \begin{equation}
1035 – F_i  = \zeta _i v_i  - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
1036 – \label{introEquation:tensorExpression}
759   \end{equation}
1038 – This equation is the basis for deriving the hydrodynamic tensor. In
1039 – 1930, Oseen and Burgers gave a simple solution to Equation
1040 – \ref{introEquation:tensorExpression}
1041 – \begin{equation}
1042 – T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
1043 – R_{ij}^T }}{{R_{ij}^2 }}} \right). \label{introEquation:oseenTensor}
1044 – \end{equation}
1045 – Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$.
1046 – A second order expression for element of different size was
1047 – introduced by Rotne and Prager\cite{Rotne1969} and improved by
1048 – Garc\'{i}a de la Torre and Bloomfield\cite{Torre1977},
1049 – \begin{equation}
1050 – T_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I +
1051 – \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma
1052 – _i^2  + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} -
1053 – \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
1054 – \label{introEquation:RPTensorNonOverlapped}
1055 – \end{equation}
1056 – Both of the Equation \ref{introEquation:oseenTensor} and Equation
1057 – \ref{introEquation:RPTensorNonOverlapped} have an assumption $R_{ij}
1058 – \ge \sigma _i  + \sigma _j$. An alternative expression for
1059 – overlapping beads with the same radius, $\sigma$, is given by
1060 – \begin{equation}
1061 – T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
1062 – \frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I +
1063 – \frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
1064 – \label{introEquation:RPTensorOverlapped}
1065 – \end{equation}
1066 –
1067 – To calculate the resistance tensor at an arbitrary origin $O$, we
1068 – construct a $3N \times 3N$ matrix consisting of $N \times N$
1069 – $B_{ij}$ blocks
1070 – \begin{equation}
1071 – B = \left( {\begin{array}{*{20}c}
1072 –   {B_{11} } &  \ldots  & {B_{1N} }  \\
1073 –    \vdots  &  \ddots  &  \vdots   \\
1074 –   {B_{N1} } &  \cdots  & {B_{NN} }  \\
1075 – \end{array}} \right),
1076 – \end{equation}
1077 – where $B_{ij}$ is given by
1078 – \[
1079 – B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
1080 – )T_{ij}
1081 – \]
1082 – where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
1083 – $B$, we obtain
1084 –
1085 – \[
1086 – C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
1087 –   {C_{11} } &  \ldots  & {C_{1N} }  \\
1088 –    \vdots  &  \ddots  &  \vdots   \\
1089 –   {C_{N1} } &  \cdots  & {C_{NN} }  \\
1090 – \end{array}} \right)
1091 – \]
1092 – , which can be partitioned into $N \times N$ $3 \times 3$ block
1093 – $C_{ij}$. With the help of $C_{ij}$ and skew matrix $U_i$
1094 – \[
1095 – U_i  = \left( {\begin{array}{*{20}c}
1096 –   0 & { - z_i } & {y_i }  \\
1097 –   {z_i } & 0 & { - x_i }  \\
1098 –   { - y_i } & {x_i } & 0  \\
1099 – \end{array}} \right)
1100 – \]
1101 – where $x_i$, $y_i$, $z_i$ are the components of the vector joining
1102 – bead $i$ and origin $O$. Hence, the elements of resistance tensor at
1103 – arbitrary origin $O$ can be written as
1104 – \begin{equation}
1105 – \begin{array}{l}
1106 – \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\
1107 – \Xi _{}^{tr}  = \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
1108 – \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j  \\
1109 – \end{array}
1110 – \label{introEquation:ResistanceTensorArbitraryOrigin}
1111 – \end{equation}
1112 –
1113 – The resistance tensor depends on the origin to which they refer. The
1114 – proper location for applying friction force is the center of
1115 – resistance (reaction), at which the trace of rotational resistance
1116 – tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
1117 – resistance is defined as an unique point of the rigid body at which
1118 – the translation-rotation coupling tensor are symmetric,
1119 – \begin{equation}
1120 – \Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
1121 – \label{introEquation:definitionCR}
1122 – \end{equation}
1123 – Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
1124 – we can easily find out that the translational resistance tensor is
1125 – origin independent, while the rotational resistance tensor and
1126 – translation-rotation coupling resistance tensor depend on the
1127 – origin. Given resistance tensor at an arbitrary origin $O$, and a
1128 – vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
1129 – obtain the resistance tensor at $P$ by
1130 – \begin{equation}
1131 – \begin{array}{l}
1132 – \Xi _P^{tt}  = \Xi _O^{tt}  \\
1133 – \Xi _P^{tr}  = \Xi _P^{rt}  = \Xi _O^{tr}  - U_{OP} \Xi _O^{tt}  \\
1134 – \Xi _P^{rr}  = \Xi _O^{rr}  - U_{OP} \Xi _O^{tt} U_{OP}  + \Xi _O^{tr} U_{OP}  - U_{OP} \Xi _O^{{tr} ^{^T }}  \\
1135 – \end{array}
1136 – \label{introEquation:resistanceTensorTransformation}
1137 – \end{equation}
1138 – where
1139 – \[
1140 – U_{OP}  = \left( {\begin{array}{*{20}c}
1141 –   0 & { - z_{OP} } & {y_{OP} }  \\
1142 –   {z_i } & 0 & { - x_{OP} }  \\
1143 –   { - y_{OP} } & {x_{OP} } & 0  \\
1144 – \end{array}} \right)
1145 – \]
1146 – Using Equations \ref{introEquation:definitionCR} and
1147 – \ref{introEquation:resistanceTensorTransformation}, one can locate
1148 – the position of center of resistance,
1149 – \begin{eqnarray*}
1150 – \left( \begin{array}{l}
1151 – x_{OR}  \\
1152 – y_{OR}  \\
1153 – z_{OR}  \\
1154 – \end{array} \right) & = &\left( {\begin{array}{*{20}c}
1155 –   {(\Xi _O^{rr} )_{yy}  + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} }  \\
1156 –   { - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz}  + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} }  \\
1157 –   { - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx}  + (\Xi _O^{rr} )_{yy} }  \\
1158 – \end{array}} \right)^{ - 1}  \\
1159 –  & & \left( \begin{array}{l}
1160 – (\Xi _O^{tr} )_{yz}  - (\Xi _O^{tr} )_{zy}  \\
1161 – (\Xi _O^{tr} )_{zx}  - (\Xi _O^{tr} )_{xz}  \\
1162 – (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
1163 – \end{array} \right) \\
1164 – \end{eqnarray*}
1165 –
1166 – where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
1167 – joining center of resistance $R$ and origin $O$.
1168 –
1169 – \subsection{Langevin dynamics for rigid particles of arbitrary shape}
1170 –
1171 – Consider a Langevin equation of motions in generalized coordinates
1172 – \begin{equation}
1173 – M_i \dot V_i (t) = F_{s,i} (t) + F_{f,i(t)}  + F_{r,i} (t)
1174 – \label{LDGeneralizedForm}
1175 – \end{equation}
1176 – where $M_i$ is a $6\times6$ generalized diagonal mass (include mass
1177 – and moment of inertial) matrix and $V_i$ is a generalized velocity,
1178 – $V_i = V_i(v_i,\omega _i)$. The right side of Eq.
1179 – (\ref{LDGeneralizedForm}) consists of three generalized forces in
1180 – lab-fixed frame, systematic force $F_{s,i}$, dissipative force
1181 – $F_{f,i}$ and stochastic force $F_{r,i}$. While the evolution of the
1182 – system in Newtownian mechanics typically refers to lab-fixed frame,
1183 – it is also convenient to handle the rotation of rigid body in
1184 – body-fixed frame. Thus the friction and random forces are calculated
1185 – in body-fixed frame and converted back to lab-fixed frame by:
1186 – \[
1187 – \begin{array}{l}
1188 – F_{f,i}^l (t) = A^T F_{f,i}^b (t), \\
1189 – F_{r,i}^l (t) = A^T F_{r,i}^b (t) \\
1190 – \end{array}.
1191 – \]
1192 – Here, the body-fixed friction force $F_{r,i}^b$ is proportional to
1193 – the body-fixed velocity at center of resistance $v_{R,i}^b$ and
1194 – angular velocity $\omega _i$,
1195 – \begin{equation}
1196 – F_{r,i}^b (t) = \left( \begin{array}{l}
1197 – f_{r,i}^b (t) \\
1198 – \tau _{r,i}^b (t) \\
1199 – \end{array} \right) =  - \left( {\begin{array}{*{20}c}
1200 –   {\Xi _{R,t} } & {\Xi _{R,c}^T }  \\
1201 –   {\Xi _{R,c} } & {\Xi _{R,r} }  \\
1202 – \end{array}} \right)\left( \begin{array}{l}
1203 – v_{R,i}^b (t) \\
1204 – \omega _i (t) \\
1205 – \end{array} \right),
1206 – \end{equation}
1207 – while the random force $F_{r,i}^l$ is a Gaussian stochastic variable
1208 – with zero mean and variance
1209 – \begin{equation}
1210 – \left\langle {F_{r,i}^l (t)(F_{r,i}^l (t'))^T } \right\rangle  =
1211 – \left\langle {F_{r,i}^b (t)(F_{r,i}^b (t'))^T } \right\rangle  =
1212 – 2k_B T\Xi _R \delta (t - t').
1213 – \end{equation}
1214 – The equation of motion for $v_i$ can be written as
1215 – \begin{equation}
1216 – m\dot v_i (t) = f_{t,i} (t) = f_{s,i} (t) + f_{f,i}^l (t) +
1217 – f_{r,i}^l (t)
1218 – \end{equation}
1219 – Since the frictional force is applied at the center of resistance
1220 – which generally does not coincide with the center of mass, an extra
1221 – torque is exerted at the center of mass. Thus, the net body-fixed
1222 – frictional torque at the center of mass, $\tau _{n,i}^b (t)$, is
1223 – given by
1224 – \begin{equation}
1225 – \tau _{r,i}^b = \tau _{r,i}^b +r_{MR} \times f_{r,i}^b
1226 – \end{equation}
1227 – where $r_{MR}$ is the vector from the center of mass to the center
1228 – of the resistance. Instead of integrating angular velocity in
1229 – lab-fixed frame, we consider the equation of motion of angular
1230 – momentum in body-fixed frame
1231 – \begin{equation}
1232 – \dot \pi _i (t) = \tau _{t,i} (t) = \tau _{s,i} (t) + \tau _{f,i}^b
1233 – (t) + \tau _{r,i}^b(t)
1234 – \end{equation}
1235 –
1236 – Embedding the friction terms into force and torque, one can
1237 – integrate the langevin equations of motion for rigid body of
1238 – arbitrary shape in a velocity-Verlet style 2-part algorithm, where
1239 – $h= \delta t$:
1240 –
1241 – {\tt part one:}
1242 – \begin{align*}
1243 – v_i (t + h/2) &\leftarrow v_i (t) + \frac{{hf_{t,i}^l (t)}}{{2m_i }} \\
1244 – \pi _i (t + h/2) &\leftarrow \pi _i (t) + \frac{{h\tau _{t,i}^b (t)}}{2} \\
1245 – r_i (t + h) &\leftarrow r_i (t) + hv_i (t + h/2) \\
1246 – A_i (t + h) &\leftarrow rotate\left( {h\pi _i (t + h/2)I^{ - 1} } \right) \\
1247 – \end{align*}
1248 – In this context, the $\mathrm{rotate}$ function is the reversible
1249 – product of five consecutive body-fixed rotations,
1250 – \begin{equation}
1251 – \mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
1252 – \mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y
1253 – / 2) \cdot \mathsf{G}_x(a_x /2),
1254 – \end{equation}
1255 – where each rotational propagator, $\mathsf{G}_\alpha(\theta)$,
1256 – rotates both the rotation matrix ($\mathsf{A}$) and the body-fixed
1257 – angular momentum ($\pi$) by an angle $\theta$ around body-fixed axis
1258 – $\alpha$,
1259 – \begin{equation}
1260 – \mathsf{G}_\alpha( \theta ) = \left\{
1261 – \begin{array}{lcl}
1262 – \mathsf{A}(t) & \leftarrow & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^T, \\
1263 – {\bf j}(t) & \leftarrow & \mathsf{R}_\alpha(\theta) \cdot {\bf
1264 – j}(0).
1265 – \end{array}
1266 – \right.
1267 – \end{equation}
1268 – $\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis
1269 – rotation matrix.  For example, in the small-angle limit, the
1270 – rotation matrix around the body-fixed x-axis can be approximated as
1271 – \begin{equation}
1272 – \mathsf{R}_x(\theta) \approx \left(
1273 – \begin{array}{ccc}
1274 – 1 & 0 & 0 \\
1275 – 0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}  & -\frac{\theta}{1+
1276 – \theta^2 / 4} \\
1277 – 0 & \frac{\theta}{1+ \theta^2 / 4} & \frac{1-\theta^2 / 4}{1 +
1278 – \theta^2 / 4}
1279 – \end{array}
1280 – \right).
1281 – \end{equation}
1282 – All other rotations follow in a straightforward manner.
1283 –
1284 – After the first part of the propagation, the friction and random
1285 – forces are generated at the center of resistance in body-fixed frame
1286 – and converted back into lab-fixed frame
1287 – \[
1288 – f_{t,i}^l (t + h) =  - \left( {\frac{{\partial V}}{{\partial r_i }}}
1289 – \right)_{r_i (t + h)}  + A_i^T (t + h)[F_{f,i}^b (t + h) + F_{r,i}^b
1290 – (t + h)],
1291 – \]
1292 – while the system torque in lab-fixed frame is transformed into
1293 – body-fixed frame,
1294 – \[
1295 – \tau _{t,i}^b (t + h) = A\tau _{s,i}^l (t) + \tau _{n,i}^b (t) +
1296 – \tau _{r,i}^b (t).
1297 – \]
1298 – Once the forces and torques have been obtained at the new time step,
1299 – the velocities can be advanced to the same time value.
1300 –
1301 – {\tt part two:}
1302 – \begin{align*}
1303 – v_i (t) &\leftarrow v_i (t + h/2) + \frac{{hf_{t,i}^l (t + h)}}{{2m_i }} \\
1304 – \pi _i (t) &\leftarrow \pi _i (t + h/2) + \frac{{h\tau _{t,i}^b (t + h)}}{2} \\
1305 – \end{align*}
1306 –
1307 – \subsection{Results and discussion}

Diff Legend

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