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

Comparing trunk/tengDissertation/Introduction.tex (file contents):
Revision 2706 by tim, Thu Apr 13 04:47:47 2006 UTC vs.
Revision 2719 by tim, Tue Apr 18 22:38:19 2006 UTC

# Line 570 | Line 570 | The free rigid body is an example of Poisson system (a
570   \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian}
571   \end{equation}
572   The most obvious change being that matrix $J$ now depends on $x$.
573 The free rigid body is an example of Poisson system (actually a
574 Lie-Poisson system) with Hamiltonian function of angular kinetic
575 energy.
576 \begin{equation}
577 J(\pi ) = \left( {\begin{array}{*{20}c}
578   0 & {\pi _3 } & { - \pi _2 }  \\
579   { - \pi _3 } & 0 & {\pi _1 }  \\
580   {\pi _2 } & { - \pi _1 } & 0  \\
581 \end{array}} \right)
582 \end{equation}
583
584 \begin{equation}
585 H = \frac{1}{2}\left( {\frac{{\pi _1^2 }}{{I_1 }} + \frac{{\pi _2^2
586 }}{{I_2 }} + \frac{{\pi _3^2 }}{{I_3 }}} \right)
587 \end{equation}
573  
574   \subsection{\label{introSection:exactFlow}Exact Flow}
575  
# Line 837 | Line 822 | q(\Delta t)} \right]. %
822   %
823   q(\Delta t) &= q(0) + \frac{{\Delta t}}{2}\left[ {\dot q(0) + \dot
824   q(\Delta t)} \right]. %
825 < \label{introEquation:positionVerlet1}
825 > \label{introEquation:positionVerlet2}
826   \end{align}
827  
828   \subsubsection{\label{introSection:errorAnalysis}Error Analysis and Higher Order Methods}
# Line 903 | Line 888 | dynamical information.
888   biological systems, providing structural, thermodynamic and
889   dynamical information.
890  
891 + One of the principal tools for modeling proteins, nucleic acids and
892 + their complexes. Stability of proteins Folding of proteins.
893 + Molecular recognition by:proteins, DNA, RNA, lipids, hormones STP,
894 + etc. Enzyme reactions Rational design of biologically active
895 + molecules (drug design) Small and large-scale conformational
896 + changes. determination and construction of 3D structures (homology,
897 + Xray diffraction, NMR) Dynamic processes such as ion transport in
898 + biological systems.
899 +
900 + Macroscopic properties are related to microscopic behavior.
901 +
902 + Time dependent (and independent) microscopic behavior of a molecule
903 + can be calculated by molecular dynamics simulations.
904 +
905   \subsection{\label{introSec:mdInit}Initialization}
906  
907   \subsection{\label{introSec:forceEvaluation}Force Evaluation}
# Line 942 | Line 941 | rotation matrix $A$ and re-formulating Hamiltonian's e
941   The break through in geometric literature suggests that, in order to
942   develop a long-term integration scheme, one should preserve the
943   symplectic structure of the flow. Introducing conjugate momentum to
944 < rotation matrix $A$ and re-formulating Hamiltonian's equation, a
944 > rotation matrix $Q$ and re-formulating Hamiltonian's equation, a
945   symplectic integrator, RSHAKE, was proposed to evolve the
946   Hamiltonian system in a constraint manifold by iteratively
947 < satisfying the orthogonality constraint $A_t A = 1$. An alternative
947 > satisfying the orthogonality constraint $Q_T Q = 1$. An alternative
948   method using quaternion representation was developed by Omelyan.
949   However, both of these methods are iterative and inefficient. In
950   this section, we will present a symplectic Lie-Poisson integrator
951 < for rigid body developed by Dullweber and his coworkers\cite{}.
951 > for rigid body developed by Dullweber and his
952 > coworkers\cite{Dullweber1997} in depth.
953  
954 \subsection{\label{introSection:lieAlgebra}Lie Algebra}
955
954   \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
955 <
955 > The motion of the rigid body is Hamiltonian with the Hamiltonian
956 > function
957   \begin{equation}
958   H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) +
959   V(q,Q) + \frac{1}{2}tr[(QQ^T  - 1)\Lambda ].
# Line 975 | Line 974 | Q^t PJ^{ - 1}  + J^{ - 1} P^t Q = 0 . \\
974   Differentiating \ref{introEquation:orthogonalConstraint} and using
975   Equation \ref{introEquation:RBMotionMomentum}, one may obtain,
976   \begin{equation}
977 < Q^t PJ^{ - 1}  + J^{ - 1} P^t Q = 0 . \\
977 > Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
978   \label{introEquation:RBFirstOrderConstraint}
979   \end{equation}
980  
# Line 987 | Line 986 | the equations of motion,
986   \frac{{dq}}{{dt}} = \frac{p}{m} \label{introEquation:RBMotionPosition}\\
987   \frac{{dp}}{{dt}} =  - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
988   \frac{{dQ}}{{dt}} = PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
989 < \frac{{dP}}{{dt}} =  - \nabla _q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}\\
989 > \frac{{dP}}{{dt}} =  - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}\\
990   \end{array}
991   \]
992  
993 <
993 > In general, there are two ways to satisfy the holonomic constraints.
994 > We can use constraint force provided by lagrange multiplier on the
995 > normal manifold to keep the motion on constraint space. Or we can
996 > simply evolve the system in constraint manifold. The two method are
997 > proved to be equivalent. The holonomic constraint and equations of
998 > motions define a constraint manifold for rigid body
999   \[
1000 < M = \left\{ {(Q,P):Q^T Q = 1,Q^t PJ^{ - 1}  + J^{ - 1} P^t Q = 0}
1001 < \right\} .
1000 > M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0}
1001 > \right\}.
1002   \]
1003  
1004 < \subsection{\label{introSection:DLMMotionEquation}The Euler Equations of Rigid Body Motion}
1005 <
1006 < \subsection{\label{introSection:symplecticDiscretizationRB}Symplectic Discretization of Euler Equations}
1007 <
1004 <
1005 < \section{\label{introSection:langevinDynamics}Langevin Dynamics}
1006 <
1007 < \subsection{\label{introSection:LDIntroduction}Introduction and application of Langevin Dynamics}
1008 <
1009 < \subsection{\label{introSection:generalizedLangevinDynamics}Generalized Langevin Dynamics}
1010 <
1011 < \begin{equation}
1012 < H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
1013 < \label{introEquation:bathGLE}
1014 < \end{equation}
1015 < where $H_B$ is harmonic bath Hamiltonian,
1004 > Unfortunately, this constraint manifold is not the cotangent bundle
1005 > $T_{\star}SO(3)$. However, it turns out that under symplectic
1006 > transformation, the cotangent space and the phase space are
1007 > diffeomorphic. Introducing
1008   \[
1009 < H_B  =\sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
1018 < }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  w_\alpha ^2 } \right\}}
1009 > \tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
1010   \]
1011 < and $\Delta U$ is bilinear system-bath coupling,
1011 > the mechanical system subject to a holonomic constraint manifold $M$
1012 > can be re-formulated as a Hamiltonian system on the cotangent space
1013   \[
1014 < \Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
1014 > T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \tilde Q =
1015 > 1,\tilde Q^T \tilde PJ^{ - 1}  + J^{ - 1} P^T \tilde Q = 0} \right\}
1016   \]
1017 < Completing the square,
1017 >
1018 > For a body fixed vector $X_i$ with respect to the center of mass of
1019 > the rigid body, its corresponding lab fixed vector $X_0^{lab}$  is
1020 > given as
1021 > \begin{equation}
1022 > X_i^{lab} = Q X_i + q.
1023 > \end{equation}
1024 > Therefore, potential energy $V(q,Q)$ is defined by
1025   \[
1026 < H_B  + \Delta U = \sum\limits_{\alpha  = 1}^N {\left\{
1027 < {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
1028 < w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
1029 < w_\alpha ^2 }}x} \right)^2 } \right\}}  - \sum\limits_{\alpha  =
1030 < 1}^N {\frac{{g_\alpha ^2 }}{{2m_\alpha  w_\alpha ^2 }}} x^2
1026 > V(q,Q) = V(Q X_0 + q).
1027   \]
1028 < and putting it back into Eq.~\ref{introEquation:bathGLE},
1028 > Hence, the force and torque are given by
1029   \[
1030 < H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
1035 < {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
1036 < w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
1037 < w_\alpha ^2 }}x} \right)^2 } \right\}}
1030 > \nabla _q V(q,Q) = F(q,Q) = \sum\limits_i {F_i (q,Q)},
1031   \]
1032 < where
1032 > and
1033   \[
1034 < W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
1042 < }}{{2m_\alpha  w_\alpha ^2 }}} x^2
1034 > \nabla _Q V(q,Q) = F(q,Q)X_i^t
1035   \]
1036 < Since the first two terms of the new Hamiltonian depend only on the
1045 < system coordinates, we can get the equations of motion for
1046 < Generalized Langevin Dynamics by Hamilton's equations
1047 < \ref{introEquation:motionHamiltonianCoordinate,
1048 < introEquation:motionHamiltonianMomentum},
1049 < \begin{align}
1050 < \dot p &=  - \frac{{\partial H}}{{\partial x}}
1051 <       &= m\ddot x
1052 <       &= - \frac{{\partial W(x)}}{{\partial x}} - \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)}
1053 < \label{introEquation:Lp5}
1054 < \end{align}
1055 < , and
1056 < \begin{align}
1057 < \dot p_\alpha   &=  - \frac{{\partial H}}{{\partial x_\alpha  }}
1058 <                &= m\ddot x_\alpha
1059 <                &= \- m_\alpha  w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha}}{{m_\alpha  w_\alpha ^2 }}x} \right)
1060 < \end{align}
1036 > respectively.
1037  
1038 < \subsection{\label{introSection:laplaceTransform}The Laplace Transform}
1039 <
1038 > As a common choice to describe the rotation dynamics of the rigid
1039 > body, angular momentum on body frame $\Pi  = Q^t P$ is introduced to
1040 > rewrite the equations of motion,
1041 > \begin{equation}
1042 > \begin{array}{l}
1043 > \mathop \Pi \limits^ \bullet   = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda  \\
1044 > \mathop Q\limits^{{\rm{   }} \bullet }  = Q\Pi {\rm{ }}J^{ - 1}  \\
1045 > \end{array}
1046 > \label{introEqaution:RBMotionPI}
1047 > \end{equation}
1048 > , as well as holonomic constraints,
1049   \[
1050 < L(x) = \int_0^\infty  {x(t)e^{ - pt} dt}
1050 > \begin{array}{l}
1051 > \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
1052 > Q^T Q = 1 \\
1053 > \end{array}
1054   \]
1055  
1056 + For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a matrix $\hat v \in
1057 + so(3)^ \star$, the hat-map isomorphism,
1058 + \begin{equation}
1059 + v(v_1 ,v_2 ,v_3 ) \Leftrightarrow \hat v = \left(
1060 + {\begin{array}{*{20}c}
1061 +   0 & { - v_3 } & {v_2 }  \\
1062 +   {v_3 } & 0 & { - v_1 }  \\
1063 +   { - v_2 } & {v_1 } & 0  \\
1064 + \end{array}} \right),
1065 + \label{introEquation:hatmapIsomorphism}
1066 + \end{equation}
1067 + will let us associate the matrix products with traditional vector
1068 + operations
1069   \[
1070 < L(x + y) = L(x) + L(y)
1070 > \hat vu = v \times u
1071   \]
1072  
1073 < \[
1074 < L(ax) = aL(x)
1075 < \]
1073 > Using \ref{introEqaution:RBMotionPI}, one can construct a skew
1074 > matrix,
1075 > \begin{equation}
1076 > (\mathop \Pi \limits^ \bullet   - \mathop \Pi \limits^ \bullet  ^T
1077 > ){\rm{ }} = {\rm{ }}(\Pi  - \Pi ^T ){\rm{ }}(J^{ - 1} \Pi  + \Pi J^{
1078 > - 1} ) + \sum\limits_i {[Q^T F_i (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]} -
1079 > (\Lambda  - \Lambda ^T ) . \label{introEquation:skewMatrixPI}
1080 > \end{equation}
1081 > Since $\Lambda$ is symmetric, the last term of Equation
1082 > \ref{introEquation:skewMatrixPI} is zero, which implies the Lagrange
1083 > multiplier $\Lambda$ is absent from the equations of motion. This
1084 > unique property eliminate the requirement of iterations which can
1085 > not be avoided in other methods\cite{}.
1086  
1087 + Applying hat-map isomorphism, we obtain the equation of motion for
1088 + angular momentum on body frame
1089 + \begin{equation}
1090 + \dot \pi  = \pi  \times I^{ - 1} \pi  + \sum\limits_i {\left( {Q^T
1091 + F_i (r,Q)} \right) \times X_i }.
1092 + \label{introEquation:bodyAngularMotion}
1093 + \end{equation}
1094 + In the same manner, the equation of motion for rotation matrix is
1095 + given by
1096   \[
1097 < L(\dot x) = pL(x) - px(0)
1097 > \dot Q = Qskew(I^{ - 1} \pi )
1098   \]
1099  
1100 < \[
1101 < L(\ddot x) = p^2 L(x) - px(0) - \dot x(0)
1082 < \]
1100 > \subsection{\label{introSection:SymplecticFreeRB}Symplectic
1101 > Lie-Poisson Integrator for Free Rigid Body}
1102  
1103 < \[
1104 < L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p)
1105 < \]
1103 > If there is not external forces exerted on the rigid body, the only
1104 > contribution to the rotational is from the kinetic potential (the
1105 > first term of \ref{ introEquation:bodyAngularMotion}). The free
1106 > rigid body is an example of Lie-Poisson system with Hamiltonian
1107 > function
1108 > \begin{equation}
1109 > T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 )
1110 > \label{introEquation:rotationalKineticRB}
1111 > \end{equation}
1112 > where $T_i^r (\pi _i ) = \frac{{\pi _i ^2 }}{{2I_i }}$ and
1113 > Lie-Poisson structure matrix,
1114 > \begin{equation}
1115 > J(\pi ) = \left( {\begin{array}{*{20}c}
1116 >   0 & {\pi _3 } & { - \pi _2 }  \\
1117 >   { - \pi _3 } & 0 & {\pi _1 }  \\
1118 >   {\pi _2 } & { - \pi _1 } & 0  \\
1119 > \end{array}} \right)
1120 > \end{equation}
1121 > Thus, the dynamics of free rigid body is governed by
1122 > \begin{equation}
1123 > \frac{d}{{dt}}\pi  = J(\pi )\nabla _\pi  T^r (\pi )
1124 > \end{equation}
1125  
1126 < Some relatively important transformation,
1127 < \[
1128 < L(\cos at) = \frac{p}{{p^2  + a^2 }}
1126 > One may notice that each $T_i^r$ in Equation
1127 > \ref{introEquation:rotationalKineticRB} can be solved exactly. For
1128 > instance, the equations of motion due to $T_1^r$ are given by
1129 > \begin{equation}
1130 > \frac{d}{{dt}}\pi  = R_1 \pi ,\frac{d}{{dt}}Q = QR_1
1131 > \label{introEqaution:RBMotionSingleTerm}
1132 > \end{equation}
1133 > where
1134 > \[ R_1  = \left( {\begin{array}{*{20}c}
1135 >   0 & 0 & 0  \\
1136 >   0 & 0 & {\pi _1 }  \\
1137 >   0 & { - \pi _1 } & 0  \\
1138 > \end{array}} \right).
1139   \]
1140 + The solutions of Equation \ref{introEqaution:RBMotionSingleTerm} is
1141 + \[
1142 + \pi (\Delta t) = e^{\Delta tR_1 } \pi (0),Q(\Delta t) =
1143 + Q(0)e^{\Delta tR_1 }
1144 + \]
1145 + with
1146 + \[
1147 + e^{\Delta tR_1 }  = \left( {\begin{array}{*{20}c}
1148 +   0 & 0 & 0  \\
1149 +   0 & {\cos \theta _1 } & {\sin \theta _1 }  \\
1150 +   0 & { - \sin \theta _1 } & {\cos \theta _1 }  \\
1151 + \end{array}} \right),\theta _1  = \frac{{\pi _1 }}{{I_1 }}\Delta t.
1152 + \]
1153 + To reduce the cost of computing expensive functions in $e^{\Delta
1154 + tR_1 }$, we can use Cayley transformation,
1155 + \[
1156 + e^{\Delta tR_1 }  \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1
1157 + )
1158 + \]
1159  
1160 + The flow maps for $T_2^r$ and $T_2^r$ can be found in the same
1161 + manner.
1162 +
1163 + In order to construct a second-order symplectic method, we split the
1164 + angular kinetic Hamiltonian function can into five terms
1165   \[
1166 < L(\sin at) = \frac{a}{{p^2  + a^2 }}
1166 > T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2
1167 > ) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r
1168 > (\pi _1 )
1169 > \].
1170 > Concatenating flows corresponding to these five terms, we can obtain
1171 > an symplectic integrator,
1172 > \[
1173 > \varphi _{\Delta t,T^r }  = \varphi _{\Delta t/2,\pi _1 }  \circ
1174 > \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }
1175 > \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
1176 > _1 }.
1177   \]
1178  
1179 + The non-canonical Lie-Poisson bracket ${F, G}$ of two function
1180 + $F(\pi )$ and $G(\pi )$ is defined by
1181   \[
1182 < L(1) = \frac{1}{p}
1182 > \{ F,G\} (\pi ) = [\nabla _\pi  F(\pi )]^T J(\pi )\nabla _\pi  G(\pi
1183 > )
1184   \]
1185 + If the Poisson bracket of a function $F$ with an arbitrary smooth
1186 + function $G$ is zero, $F$ is a \emph{Casimir}, which is the
1187 + conserved quantity in Poisson system. We can easily verify that the
1188 + norm of the angular momentum, $\parallel \pi
1189 + \parallel$, is a \emph{Casimir}. Let$ F(\pi ) = S(\frac{{\parallel
1190 + \pi \parallel ^2 }}{2})$ for an arbitrary function $ S:R \to R$ ,
1191 + then by the chain rule
1192 + \[
1193 + \nabla _\pi  F(\pi ) = S'(\frac{{\parallel \pi \parallel ^2
1194 + }}{2})\pi
1195 + \]
1196 + Thus $ [\nabla _\pi  F(\pi )]^T J(\pi ) =  - S'(\frac{{\parallel \pi
1197 + \parallel ^2 }}{2})\pi  \times \pi  = 0 $. This explicit
1198 + Lie-Poisson integrator is found to be extremely efficient and stable
1199 + which can be explained by the fact the small angle approximation is
1200 + used and the norm of the angular momentum is conserved.
1201  
1202 < First, the bath coordinates,
1202 > \subsection{\label{introSection:RBHamiltonianSplitting} Hamiltonian
1203 > Splitting for Rigid Body}
1204 >
1205 > The Hamiltonian of rigid body can be separated in terms of kinetic
1206 > energy and potential energy,
1207   \[
1208 < p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega
1104 < _\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha
1105 < }}L(x)
1208 > H = T(p,\pi ) + V(q,Q)
1209   \]
1210 + The equations of motion corresponding to potential energy and
1211 + kinetic energy are listed in the below table,
1212 + \begin{center}
1213 + \begin{tabular}{|l|l|}
1214 +  \hline
1215 +  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
1216 +  Potential & Kinetic \\
1217 +  $\frac{{dq}}{{dt}} = \frac{p}{m}$ & $\frac{d}{{dt}}q = p$ \\
1218 +  $\frac{d}{{dt}}p =  - \frac{{\partial V}}{{\partial q}}$ & $ \frac{d}{{dt}}p = 0$ \\
1219 +  $\frac{d}{{dt}}Q = 0$ & $ \frac{d}{{dt}}Q = Qskew(I^{ - 1} j)$ \\
1220 +  $ \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$\\
1221 +  \hline
1222 + \end{tabular}
1223 + \end{center}
1224 + A second-order symplectic method is now obtained by the composition
1225 + of the flow maps,
1226   \[
1227 < L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) +
1228 < px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }}
1227 > \varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
1228 > _{\Delta t,T}  \circ \varphi _{\Delta t/2,V}.
1229 > \]
1230 > Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two
1231 > sub-flows which corresponding to force and torque respectively,
1232 > \[
1233 > \varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
1234 > _{\Delta t/2,\tau }.
1235 > \]
1236 > Since the associated operators of $\varphi _{\Delta t/2,F} $ and
1237 > $\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition
1238 > order inside $\varphi _{\Delta t/2,V}$ does not matter.
1239 >
1240 > Furthermore, kinetic potential can be separated to translational
1241 > kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$,
1242 > \begin{equation}
1243 > T(p,\pi ) =T^t (p) + T^r (\pi ).
1244 > \end{equation}
1245 > where $ T^t (p) = \frac{1}{2}p^T m^{ - 1} p $ and $T^r (\pi )$ is
1246 > defined by \ref{introEquation:rotationalKineticRB}. Therefore, the
1247 > corresponding flow maps are given by
1248 > \[
1249 > \varphi _{\Delta t,T}  = \varphi _{\Delta t,T^t }  \circ \varphi
1250 > _{\Delta t,T^r }.
1251 > \]
1252 > Finally, we obtain the overall symplectic flow maps for free moving
1253 > rigid body
1254 > \begin{equation}
1255 > \begin{array}{c}
1256 > \varphi _{\Delta t}  = \varphi _{\Delta t/2,F}  \circ \varphi _{\Delta t/2,\tau }  \\
1257 >  \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 }  \\
1258 >  \circ \varphi _{\Delta t/2,\tau }  \circ \varphi _{\Delta t/2,F}  .\\
1259 > \end{array}
1260 > \label{introEquation:overallRBFlowMaps}
1261 > \end{equation}
1262 >
1263 > \section{\label{introSection:langevinDynamics}Langevin Dynamics}
1264 > As an alternative to newtonian dynamics, Langevin dynamics, which
1265 > mimics a simple heat bath with stochastic and dissipative forces,
1266 > has been applied in a variety of studies. This section will review
1267 > the theory of Langevin dynamics simulation. A brief derivation of
1268 > generalized Langevin equation will be given first. Follow that, we
1269 > will discuss the physical meaning of the terms appearing in the
1270 > equation as well as the calculation of friction tensor from
1271 > hydrodynamics theory.
1272 >
1273 > \subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation}
1274 >
1275 > Harmonic bath model, in which an effective set of harmonic
1276 > oscillators are used to mimic the effect of a linearly responding
1277 > environment, has been widely used in quantum chemistry and
1278 > statistical mechanics. One of the successful applications of
1279 > Harmonic bath model is the derivation of Deriving Generalized
1280 > Langevin Dynamics. Lets consider a system, in which the degree of
1281 > freedom $x$ is assumed to couple to the bath linearly, giving a
1282 > Hamiltonian of the form
1283 > \begin{equation}
1284 > H = \frac{{p^2 }}{{2m}} + U(x) + H_B  + \Delta U(x,x_1 , \ldots x_N)
1285 > \label{introEquation:bathGLE}.
1286 > \end{equation}
1287 > Here $p$ is a momentum conjugate to $q$, $m$ is the mass associated
1288 > with this degree of freedom, $H_B$ is harmonic bath Hamiltonian,
1289 > \[
1290 > H_B  = \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{p_\alpha ^2
1291 > }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha  \omega _\alpha ^2 }
1292 > \right\}}
1293 > \]
1294 > where the index $\alpha$ runs over all the bath degrees of freedom,
1295 > $\omega _\alpha$ are the harmonic bath frequencies, $m_\alpha$ are
1296 > the harmonic bath masses, and $\Delta U$ is bilinear system-bath
1297 > coupling,
1298 > \[
1299 > \Delta U =  - \sum\limits_{\alpha  = 1}^N {g_\alpha  x_\alpha  x}
1300 > \]
1301 > where $g_\alpha$ are the coupling constants between the bath and the
1302 > coordinate $x$. Introducing
1303 > \[
1304 > W(x) = U(x) - \sum\limits_{\alpha  = 1}^N {\frac{{g_\alpha ^2
1305 > }}{{2m_\alpha  w_\alpha ^2 }}} x^2
1306 > \] and combining the last two terms in Equation
1307 > \ref{introEquation:bathGLE}, we may rewrite the Harmonic bath
1308 > Hamiltonian as
1309 > \[
1310 > H = \frac{{p^2 }}{{2m}} + W(x) + \sum\limits_{\alpha  = 1}^N
1311 > {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha  }} + \frac{1}{2}m_\alpha
1312 > w_\alpha ^2 \left( {x_\alpha   - \frac{{g_\alpha  }}{{m_\alpha
1313 > w_\alpha ^2 }}x} \right)^2 } \right\}}
1314 > \]
1315 > Since the first two terms of the new Hamiltonian depend only on the
1316 > system coordinates, we can get the equations of motion for
1317 > Generalized Langevin Dynamics by Hamilton's equations
1318 > \ref{introEquation:motionHamiltonianCoordinate,
1319 > introEquation:motionHamiltonianMomentum},
1320 > \begin{equation}
1321 > m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} -
1322 > \sum\limits_{\alpha  = 1}^N {g_\alpha  \left( {x_\alpha   -
1323 > \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right)},
1324 > \label{introEquation:coorMotionGLE}
1325 > \end{equation}
1326 > and
1327 > \begin{equation}
1328 > m\ddot x_\alpha   =  - m_\alpha  w_\alpha ^2 \left( {x_\alpha   -
1329 > \frac{{g_\alpha  }}{{m_\alpha  w_\alpha ^2 }}x} \right).
1330 > \label{introEquation:bathMotionGLE}
1331 > \end{equation}
1332 >
1333 > In order to derive an equation for $x$, the dynamics of the bath
1334 > variables $x_\alpha$ must be solved exactly first. As an integral
1335 > transform which is particularly useful in solving linear ordinary
1336 > differential equations, Laplace transform is the appropriate tool to
1337 > solve this problem. The basic idea is to transform the difficult
1338 > differential equations into simple algebra problems which can be
1339 > solved easily. Then applying inverse Laplace transform, also known
1340 > as the Bromwich integral, we can retrieve the solutions of the
1341 > original problems.
1342 >
1343 > Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace
1344 > transform of f(t) is a new function defined as
1345 > \[
1346 > L(f(t)) \equiv F(p) = \int_0^\infty  {f(t)e^{ - pt} dt}
1347 > \]
1348 > where  $p$ is real and  $L$ is called the Laplace Transform
1349 > Operator. Below are some important properties of Laplace transform
1350 > \begin{equation}
1351 > \begin{array}{c}
1352 > L(x + y) = L(x) + L(y) \\
1353 > L(ax) = aL(x) \\
1354 > L(\dot x) = pL(x) - px(0) \\
1355 > L(\ddot x) = p^2 L(x) - px(0) - \dot x(0) \\
1356 > L\left( {\int_0^t {g(t - \tau )h(\tau )d\tau } } \right) = G(p)H(p) \\
1357 > \end{array}
1358 > \end{equation}
1359 >
1360 > Applying Laplace transform to the bath coordinates, we obtain
1361 > \[
1362 > \begin{array}{c}
1363 > p^2 L(x_\alpha  ) - px_\alpha  (0) - \dot x_\alpha  (0) =  - \omega _\alpha ^2 L(x_\alpha  ) + \frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) \\
1364 > L(x_\alpha  ) = \frac{{\frac{{g_\alpha  }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha  (0)}}{{p^2  + \omega _\alpha ^2 }} \\
1365 > \end{array}
1366   \]
1367 < Then, the system coordinates,
1368 < \begin{align}
1369 < mL(\ddot x) &=  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
1370 < \sum\limits_{\alpha  = 1}^N {\left\{ {\frac{{\frac{{g_\alpha
1371 < }}{{\omega _\alpha  }}L(x) + px_\alpha  (0) + \dot x_\alpha
1372 < (0)}}{{p^2  + \omega _\alpha ^2 }} - \frac{{g_\alpha ^2 }}{{m_\alpha
1373 < }}\omega _\alpha ^2 L(x)} \right\}}
1118 < %
1119 < &= - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} -
1120 < \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x)
1121 < - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0)
1122 < - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}
1123 < \end{align}
1124 < Then, the inverse transform,
1367 > By the same way, the system coordinates become
1368 > \[
1369 > \begin{array}{c}
1370 > mL(\ddot x) =  - \frac{1}{p}\frac{{\partial W(x)}}{{\partial x}} \\
1371 >  - \sum\limits_{\alpha  = 1}^N {\left\{ { - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}\frac{p}{{p^2  + \omega _\alpha ^2 }}pL(x) - \frac{p}{{p^2  + \omega _\alpha ^2 }}g_\alpha  x_\alpha  (0) - \frac{1}{{p^2  + \omega _\alpha ^2 }}g_\alpha  \dot x_\alpha  (0)} \right\}}  \\
1372 > \end{array}
1373 > \]
1374  
1375 + With the help of some relatively important inverse Laplace
1376 + transformations:
1377 + \[
1378 + \begin{array}{c}
1379 + L(\cos at) = \frac{p}{{p^2  + a^2 }} \\
1380 + L(\sin at) = \frac{a}{{p^2  + a^2 }} \\
1381 + L(1) = \frac{1}{p} \\
1382 + \end{array}
1383 + \]
1384 + , we obtain
1385   \begin{align}
1386   m\ddot x &=  - \frac{{\partial W(x)}}{{\partial x}} -
1387   \sum\limits_{\alpha  = 1}^N {\left\{ {\left( { - \frac{{g_\alpha ^2
# Line 1142 | Line 1401 | t)\dot x(t - \tau )d} \tau }  + \sum\limits_{\alpha  =
1401   (\omega _\alpha  t)} \right\}}
1402   \end{align}
1403  
1404 + Introducing a \emph{dynamic friction kernel}
1405   \begin{equation}
1406 + \xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1407 + }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
1408 + \label{introEquation:dynamicFrictionKernelDefinition}
1409 + \end{equation}
1410 + and \emph{a random force}
1411 + \begin{equation}
1412 + R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
1413 + - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
1414 + \right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
1415 + (0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t),
1416 + \label{introEquation:randomForceDefinition}
1417 + \end{equation}
1418 + the equation of motion can be rewritten as
1419 + \begin{equation}
1420   m\ddot x =  - \frac{{\partial W}}{{\partial x}} - \int_0^t {\xi
1421   (t)\dot x(t - \tau )d\tau }  + R(t)
1422   \label{introEuqation:GeneralizedLangevinDynamics}
1423   \end{equation}
1424 < %where $ {\xi (t)}$ is friction kernel, $R(t)$ is random force and
1425 < %$W$ is the potential of mean force. $W(x) =  - kT\ln p(x)$
1424 > which is known as the \emph{generalized Langevin equation}.
1425 >
1426 > \subsubsection{\label{introSection:randomForceDynamicFrictionKernel}Random Force and Dynamic Friction Kernel}
1427 >
1428 > One may notice that $R(t)$ depends only on initial conditions, which
1429 > implies it is completely deterministic within the context of a
1430 > harmonic bath. However, it is easy to verify that $R(t)$ is totally
1431 > uncorrelated to $x$ and $\dot x$,
1432   \[
1433 < \xi (t) = \sum\limits_{\alpha  = 1}^N {\left( { - \frac{{g_\alpha ^2
1434 < }}{{m_\alpha  \omega _\alpha ^2 }}} \right)\cos (\omega _\alpha  t)}
1433 > \begin{array}{l}
1434 > \left\langle {x(t)R(t)} \right\rangle  = 0, \\
1435 > \left\langle {\dot x(t)R(t)} \right\rangle  = 0. \\
1436 > \end{array}
1437   \]
1438 < For an infinite harmonic bath, we can use the spectral density and
1439 < an integral over frequencies.
1438 > This property is what we expect from a truly random process. As long
1439 > as the model, which is gaussian distribution in general, chosen for
1440 > $R(t)$ is a truly random process, the stochastic nature of the GLE
1441 > still remains.
1442  
1443 + %dynamic friction kernel
1444 + The convolution integral
1445   \[
1446 < R(t) = \sum\limits_{\alpha  = 1}^N {\left( {g_\alpha  x_\alpha  (0)
1161 < - \frac{{g_\alpha ^2 }}{{m_\alpha  \omega _\alpha ^2 }}x(0)}
1162 < \right)\cos (\omega _\alpha  t)}  + \frac{{\dot x_\alpha
1163 < (0)}}{{\omega _\alpha  }}\sin (\omega _\alpha  t)
1446 > \int_0^t {\xi (t)\dot x(t - \tau )d\tau }
1447   \]
1448 < The random forces depend only on initial conditions.
1448 > depends on the entire history of the evolution of $x$, which implies
1449 > that the bath retains memory of previous motions. In other words,
1450 > the bath requires a finite time to respond to change in the motion
1451 > of the system. For a sluggish bath which responds slowly to changes
1452 > in the system coordinate, we may regard $\xi(t)$ as a constant
1453 > $\xi(t) = \Xi_0$. Hence, the convolution integral becomes
1454 > \[
1455 > \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = \xi _0 (x(t) - x(0))
1456 > \]
1457 > and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1458 > \[
1459 > m\ddot x =  - \frac{\partial }{{\partial x}}\left( {W(x) +
1460 > \frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t),
1461 > \]
1462 > which can be used to describe dynamic caging effect. The other
1463 > extreme is the bath that responds infinitely quickly to motions in
1464 > the system. Thus, $\xi (t)$ can be taken as a $delta$ function in
1465 > time:
1466 > \[
1467 > \xi (t) = 2\xi _0 \delta (t)
1468 > \]
1469 > Hence, the convolution integral becomes
1470 > \[
1471 > \int_0^t {\xi (t)\dot x(t - \tau )d\tau }  = 2\xi _0 \int_0^t
1472 > {\delta (t)\dot x(t - \tau )d\tau }  = \xi _0 \dot x(t),
1473 > \]
1474 > and Equation \ref{introEuqation:GeneralizedLangevinDynamics} becomes
1475 > \begin{equation}
1476 > m\ddot x =  - \frac{{\partial W(x)}}{{\partial x}} - \xi _0 \dot
1477 > x(t) + R(t) \label{introEquation:LangevinEquation}
1478 > \end{equation}
1479 > which is known as the Langevin equation. The static friction
1480 > coefficient $\xi _0$ can either be calculated from spectral density
1481 > or be determined by Stokes' law for regular shaped particles.A
1482 > briefly review on calculating friction tensor for arbitrary shaped
1483 > particles is given in section \ref{introSection:frictionTensor}.
1484  
1485   \subsubsection{\label{introSection:secondFluctuationDissipation}The Second Fluctuation Dissipation Theorem}
1486 < So we can define a new set of coordinates,
1486 >
1487 > Defining a new set of coordinates,
1488   \[
1489   q_\alpha  (t) = x_\alpha  (t) - \frac{1}{{m_\alpha  \omega _\alpha
1490   ^2 }}x(0)
1491 < \]
1492 < This makes
1491 > \],
1492 > we can rewrite $R(T)$ as
1493   \[
1494 < R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}
1494 > R(t) = \sum\limits_{\alpha  = 1}^N {g_\alpha  q_\alpha  (t)}.
1495   \]
1496   And since the $q$ coordinates are harmonic oscillators,
1497   \[
1498 < \begin{array}{l}
1498 > \begin{array}{c}
1499 > \left\langle {q_\alpha ^2 } \right\rangle  = \frac{{kT}}{{m_\alpha  \omega _\alpha ^2 }} \\
1500   \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  = \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t) \\
1501   \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle  = \delta _{\alpha \beta } \left\langle {q_\alpha  (t)q_\alpha  (0)} \right\rangle  \\
1502 + \left\langle {R(t)R(0)} \right\rangle  = \sum\limits_\alpha  {\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha  (t)q_\beta  (0)} \right\rangle } }  \\
1503 +  = \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)} \right\rangle \cos (\omega _\alpha  t)}  \\
1504 +  = kT\xi (t) \\
1505   \end{array}
1506   \]
1507 <
1185 < \begin{align}
1186 < \left\langle {R(t)R(0)} \right\rangle  &= \sum\limits_\alpha
1187 < {\sum\limits_\beta  {g_\alpha  g_\beta  \left\langle {q_\alpha
1188 < (t)q_\beta  (0)} \right\rangle } }
1189 < %
1190 < &= \sum\limits_\alpha  {g_\alpha ^2 \left\langle {q_\alpha ^2 (0)}
1191 < \right\rangle \cos (\omega _\alpha  t)}
1192 < %
1193 < &= kT\xi (t)
1194 < \end{align}
1195 <
1507 > Thus, we recover the \emph{second fluctuation dissipation theorem}
1508   \begin{equation}
1509   \xi (t) = \left\langle {R(t)R(0)} \right\rangle
1510 < \label{introEquation:secondFluctuationDissipation}
1510 > \label{introEquation:secondFluctuationDissipation}.
1511   \end{equation}
1512 + In effect, it acts as a constraint on the possible ways in which one
1513 + can model the random force and friction kernel.
1514  
1201 \section{\label{introSection:hydroynamics}Hydrodynamics}
1202
1515   \subsection{\label{introSection:frictionTensor} Friction Tensor}
1516 < \subsection{\label{introSection:analyticalApproach}Analytical
1517 < Approach}
1516 > Theoretically, the friction kernel can be determined using velocity
1517 > autocorrelation function. However, this approach become impractical
1518 > when the system become more and more complicate. Instead, various
1519 > approaches based on hydrodynamics have been developed to calculate
1520 > the friction coefficients. The friction effect is isotropic in
1521 > Equation, \zeta can be taken as a scalar. In general, friction
1522 > tensor \Xi is a $6\times 6$ matrix given by
1523 > \[
1524 > \Xi  = \left( {\begin{array}{*{20}c}
1525 >   {\Xi _{}^{tt} } & {\Xi _{}^{rt} }  \\
1526 >   {\Xi _{}^{tr} } & {\Xi _{}^{rr} }  \\
1527 > \end{array}} \right).
1528 > \]
1529 > Here, $ {\Xi^{tt} }$ and $ {\Xi^{rr} }$ are translational friction
1530 > tensor and rotational resistance (friction) tensor respectively,
1531 > while ${\Xi^{tr} }$ is translation-rotation coupling tensor and $
1532 > {\Xi^{rt} }$ is rotation-translation coupling tensor. When a
1533 > particle moves in a fluid, it may experience friction force or
1534 > torque along the opposite direction of the velocity or angular
1535 > velocity,
1536 > \[
1537 > \left( \begin{array}{l}
1538 > F_R  \\
1539 > \tau _R  \\
1540 > \end{array} \right) =  - \left( {\begin{array}{*{20}c}
1541 >   {\Xi ^{tt} } & {\Xi ^{rt} }  \\
1542 >   {\Xi ^{tr} } & {\Xi ^{rr} }  \\
1543 > \end{array}} \right)\left( \begin{array}{l}
1544 > v \\
1545 > w \\
1546 > \end{array} \right)
1547 > \]
1548 > where $F_r$ is the friction force and $\tau _R$ is the friction
1549 > toque.
1550  
1551 < \subsection{\label{introSection:approximationApproach}Approximation
1208 < Approach}
1551 > \subsubsection{\label{introSection:resistanceTensorRegular}The Resistance Tensor for Regular Shape}
1552  
1553 < \subsection{\label{introSection:centersRigidBody}Centers of Rigid
1554 < Body}
1553 > For a spherical particle, the translational and rotational friction
1554 > constant can be calculated from Stoke's law,
1555 > \[
1556 > \Xi ^{tt}  = \left( {\begin{array}{*{20}c}
1557 >   {6\pi \eta R} & 0 & 0  \\
1558 >   0 & {6\pi \eta R} & 0  \\
1559 >   0 & 0 & {6\pi \eta R}  \\
1560 > \end{array}} \right)
1561 > \]
1562 > and
1563 > \[
1564 > \Xi ^{rr}  = \left( {\begin{array}{*{20}c}
1565 >   {8\pi \eta R^3 } & 0 & 0  \\
1566 >   0 & {8\pi \eta R^3 } & 0  \\
1567 >   0 & 0 & {8\pi \eta R^3 }  \\
1568 > \end{array}} \right)
1569 > \]
1570 > where $\eta$ is the viscosity of the solvent and $R$ is the
1571 > hydrodynamics radius.
1572  
1573 < \section{\label{introSection:correlationFunctions}Correlation Functions}
1573 > Other non-spherical shape, such as cylinder and ellipsoid
1574 > \textit{etc}, are widely used as reference for developing new
1575 > hydrodynamics theory, because their properties can be calculated
1576 > exactly. In 1936, Perrin extended Stokes's law to general ellipsoid,
1577 > also called a triaxial ellipsoid, which is given in Cartesian
1578 > coordinates by
1579 > \[
1580 > \frac{{x^2 }}{{a^2 }} + \frac{{y^2 }}{{b^2 }} + \frac{{z^2 }}{{c^2
1581 > }} = 1
1582 > \]
1583 > where the semi-axes are of lengths $a$, $b$, and $c$. Unfortunately,
1584 > due to the complexity of the elliptic integral, only the ellipsoid
1585 > with the restriction of two axes having to be equal, \textit{i.e.}
1586 > prolate($ a \ge b = c$) and oblate ($ a < b = c $), can be solved
1587 > exactly. Introducing an elliptic integral parameter $S$ for prolate,
1588 > \[
1589 > S = \frac{2}{{\sqrt {a^2  - b^2 } }}\ln \frac{{a + \sqrt {a^2  - b^2
1590 > } }}{b},
1591 > \]
1592 > and oblate,
1593 > \[
1594 > S = \frac{2}{{\sqrt {b^2  - a^2 } }}arctg\frac{{\sqrt {b^2  - a^2 }
1595 > }}{a}
1596 > \],
1597 > one can write down the translational and rotational resistance
1598 > tensors
1599 > \[
1600 > \begin{array}{l}
1601 > \Xi _a^{tt}  = 16\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - b^2 )S - 2a}} \\
1602 > \Xi _b^{tt}  = \Xi _c^{tt}  = 32\pi \eta \frac{{a^2  - b^2 }}{{(2a^2  - 3b^2 )S + 2a}} \\
1603 > \end{array},
1604 > \]
1605 > and
1606 > \[
1607 > \begin{array}{l}
1608 > \Xi _a^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^2  - b^2 )b^2 }}{{2a - b^2 S}} \\
1609 > \Xi _b^{rr}  = \Xi _c^{rr}  = \frac{{32\pi }}{3}\eta \frac{{(a^4  - b^4 )}}{{(2a^2  - b^2 )S - 2a}} \\
1610 > \end{array}.
1611 > \]
1612 >
1613 > \subsubsection{\label{introSection:resistanceTensorRegularArbitrary}The Resistance Tensor for Arbitrary Shape}
1614 >
1615 > Unlike spherical and other regular shaped molecules, there is not
1616 > analytical solution for friction tensor of any arbitrary shaped
1617 > rigid molecules. The ellipsoid of revolution model and general
1618 > triaxial ellipsoid model have been used to approximate the
1619 > hydrodynamic properties of rigid bodies. However, since the mapping
1620 > from all possible ellipsoidal space, $r$-space, to all possible
1621 > combination of rotational diffusion coefficients, $D$-space is not
1622 > unique\cite{Wegener79} as well as the intrinsic coupling between
1623 > translational and rotational motion of rigid body\cite{}, general
1624 > ellipsoid is not always suitable for modeling arbitrarily shaped
1625 > rigid molecule. A number of studies have been devoted to determine
1626 > the friction tensor for irregularly shaped rigid bodies using more
1627 > advanced method\cite{} where the molecule of interest was modeled by
1628 > combinations of spheres(beads)\cite{} and the hydrodynamics
1629 > properties of the molecule can be calculated using the hydrodynamic
1630 > interaction tensor. Let us consider a rigid assembly of $N$ beads
1631 > immersed in a continuous medium. Due to hydrodynamics interaction,
1632 > the ``net'' velocity of $i$th bead, $v'_i$ is different than its
1633 > unperturbed velocity $v_i$,
1634 > \[
1635 > v'_i  = v_i  - \sum\limits_{j \ne i} {T_{ij} F_j }
1636 > \]
1637 > where $F_i$ is the frictional force, and $T_{ij}$ is the
1638 > hydrodynamic interaction tensor. The friction force of $i$th bead is
1639 > proportional to its ``net'' velocity
1640 > \begin{equation}
1641 > F_i  = \zeta _i v_i  - \zeta _i \sum\limits_{j \ne i} {T_{ij} F_j }.
1642 > \label{introEquation:tensorExpression}
1643 > \end{equation}
1644 > This equation is the basis for deriving the hydrodynamic tensor. In
1645 > 1930, Oseen and Burgers gave a simple solution to Equation
1646 > \ref{introEquation:tensorExpression}
1647 > \begin{equation}
1648 > T_{ij}  = \frac{1}{{8\pi \eta r_{ij} }}\left( {I + \frac{{R_{ij}
1649 > R_{ij}^T }}{{R_{ij}^2 }}} \right).
1650 > \label{introEquation:oseenTensor}
1651 > \end{equation}
1652 > Here $R_{ij}$ is the distance vector between bead $i$ and bead $j$.
1653 > A second order expression for element of different size was
1654 > introduced by Rotne and Prager\cite{} and improved by Garc\'{i}a de
1655 > la Torre and Bloomfield,
1656 > \begin{equation}
1657 > T_{ij}  = \frac{1}{{8\pi \eta R_{ij} }}\left[ {\left( {I +
1658 > \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right) + R\frac{{\sigma
1659 > _i^2  + \sigma _j^2 }}{{r_{ij}^2 }}\left( {\frac{I}{3} -
1660 > \frac{{R_{ij} R_{ij}^T }}{{R_{ij}^2 }}} \right)} \right].
1661 > \label{introEquation:RPTensorNonOverlapped}
1662 > \end{equation}
1663 > Both of the Equation \ref{introEquation:oseenTensor} and Equation
1664 > \ref{introEquation:RPTensorNonOverlapped} have an assumption $R_{ij}
1665 > \ge \sigma _i  + \sigma _j$. An alternative expression for
1666 > overlapping beads with the same radius, $\sigma$, is given by
1667 > \begin{equation}
1668 > T_{ij}  = \frac{1}{{6\pi \eta R_{ij} }}\left[ {\left( {1 -
1669 > \frac{2}{{32}}\frac{{R_{ij} }}{\sigma }} \right)I +
1670 > \frac{2}{{32}}\frac{{R_{ij} R_{ij}^T }}{{R_{ij} \sigma }}} \right]
1671 > \label{introEquation:RPTensorOverlapped}
1672 > \end{equation}
1673 >
1674 > To calculate the resistance tensor at an arbitrary origin $O$, we
1675 > construct a $3N \times 3N$ matrix consisting of $N \times N$
1676 > $B_{ij}$ blocks
1677 > \begin{equation}
1678 > B = \left( {\begin{array}{*{20}c}
1679 >   {B_{11} } &  \ldots  & {B_{1N} }  \\
1680 >    \vdots  &  \ddots  &  \vdots   \\
1681 >   {B_{N1} } &  \cdots  & {B_{NN} }  \\
1682 > \end{array}} \right),
1683 > \end{equation}
1684 > where $B_{ij}$ is given by
1685 > \[
1686 > B_{ij}  = \delta _{ij} \frac{I}{{6\pi \eta R}} + (1 - \delta _{ij}
1687 > )T_{ij}
1688 > \]
1689 > where $\delta _{ij}$ is Kronecker delta function. Inverting matrix
1690 > $B$, we obtain
1691 >
1692 > \[
1693 > C = B^{ - 1}  = \left( {\begin{array}{*{20}c}
1694 >   {C_{11} } &  \ldots  & {C_{1N} }  \\
1695 >    \vdots  &  \ddots  &  \vdots   \\
1696 >   {C_{N1} } &  \cdots  & {C_{NN} }  \\
1697 > \end{array}} \right)
1698 > \]
1699 > , which can be partitioned into $N \times N$ $3 \times 3$ block
1700 > $C_{ij}$. With the help of $C_{ij}$ and skew matrix $U_i$
1701 > \[
1702 > U_i  = \left( {\begin{array}{*{20}c}
1703 >   0 & { - z_i } & {y_i }  \\
1704 >   {z_i } & 0 & { - x_i }  \\
1705 >   { - y_i } & {x_i } & 0  \\
1706 > \end{array}} \right)
1707 > \]
1708 > where $x_i$, $y_i$, $z_i$ are the components of the vector joining
1709 > bead $i$ and origin $O$. Hence, the elements of resistance tensor at
1710 > arbitrary origin $O$ can be written as
1711 > \begin{equation}
1712 > \begin{array}{l}
1713 > \Xi _{}^{tt}  = \sum\limits_i {\sum\limits_j {C_{ij} } } , \\
1714 > \Xi _{}^{tr}  = \Xi _{}^{rt}  = \sum\limits_i {\sum\limits_j {U_i C_{ij} } } , \\
1715 > \Xi _{}^{rr}  =  - \sum\limits_i {\sum\limits_j {U_i C_{ij} } } U_j  \\
1716 > \end{array}
1717 > \label{introEquation:ResistanceTensorArbitraryOrigin}
1718 > \end{equation}
1719 >
1720 > The resistance tensor depends on the origin to which they refer. The
1721 > proper location for applying friction force is the center of
1722 > resistance (reaction), at which the trace of rotational resistance
1723 > tensor, $ \Xi ^{rr}$ reaches minimum. Mathematically, the center of
1724 > resistance is defined as an unique point of the rigid body at which
1725 > the translation-rotation coupling tensor are symmetric,
1726 > \begin{equation}
1727 > \Xi^{tr}  = \left( {\Xi^{tr} } \right)^T
1728 > \label{introEquation:definitionCR}
1729 > \end{equation}
1730 > Form Equation \ref{introEquation:ResistanceTensorArbitraryOrigin},
1731 > we can easily find out that the translational resistance tensor is
1732 > origin independent, while the rotational resistance tensor and
1733 > translation-rotation coupling resistance tensor depend on the
1734 > origin. Given resistance tensor at an arbitrary origin $O$, and a
1735 > vector ,$r_{OP}(x_{OP}, y_{OP}, z_{OP})$, from $O$ to $P$, we can
1736 > obtain the resistance tensor at $P$ by
1737 > \begin{equation}
1738 > \begin{array}{l}
1739 > \Xi _P^{tt}  = \Xi _O^{tt}  \\
1740 > \Xi _P^{tr}  = \Xi _P^{rt}  = \Xi _O^{tr}  - U_{OP} \Xi _O^{tt}  \\
1741 > \Xi _P^{rr}  = \Xi _O^{rr}  - U_{OP} \Xi _O^{tt} U_{OP}  + \Xi _O^{tr} U_{OP}  - U_{OP} \Xi _O^{tr} ^{^T }  \\
1742 > \end{array}
1743 > \label{introEquation:resistanceTensorTransformation}
1744 > \end{equation}
1745 > where
1746 > \[
1747 > U_{OP}  = \left( {\begin{array}{*{20}c}
1748 >   0 & { - z_{OP} } & {y_{OP} }  \\
1749 >   {z_i } & 0 & { - x_{OP} }  \\
1750 >   { - y_{OP} } & {x_{OP} } & 0  \\
1751 > \end{array}} \right)
1752 > \]
1753 > Using Equations \ref{introEquation:definitionCR} and
1754 > \ref{introEquation:resistanceTensorTransformation}, one can locate
1755 > the position of center of resistance,
1756 > \[
1757 > \left( \begin{array}{l}
1758 > x_{OR}  \\
1759 > y_{OR}  \\
1760 > z_{OR}  \\
1761 > \end{array} \right) = \left( {\begin{array}{*{20}c}
1762 >   {(\Xi _O^{rr} )_{yy}  + (\Xi _O^{rr} )_{zz} } & { - (\Xi _O^{rr} )_{xy} } & { - (\Xi _O^{rr} )_{xz} }  \\
1763 >   { - (\Xi _O^{rr} )_{xy} } & {(\Xi _O^{rr} )_{zz}  + (\Xi _O^{rr} )_{xx} } & { - (\Xi _O^{rr} )_{yz} }  \\
1764 >   { - (\Xi _O^{rr} )_{xz} } & { - (\Xi _O^{rr} )_{yz} } & {(\Xi _O^{rr} )_{xx}  + (\Xi _O^{rr} )_{yy} }  \\
1765 > \end{array}} \right)^{ - 1} \left( \begin{array}{l}
1766 > (\Xi _O^{tr} )_{yz}  - (\Xi _O^{tr} )_{zy}  \\
1767 > (\Xi _O^{tr} )_{zx}  - (\Xi _O^{tr} )_{xz}  \\
1768 > (\Xi _O^{tr} )_{xy}  - (\Xi _O^{tr} )_{yx}  \\
1769 > \end{array} \right).
1770 > \]
1771 > where $x_OR$, $y_OR$, $z_OR$ are the components of the vector
1772 > joining center of resistance $R$ and origin $O$.
1773 >
1774 > %\section{\label{introSection:correlationFunctions}Correlation Functions}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines