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 2702 by tim, Wed Apr 12 03:37:46 2006 UTC vs.
Revision 2707 by tim, Thu Apr 13 22:17:13 2006 UTC

# Line 315 | Line 315 | partition function like,
315   isolated and conserve energy, Microcanonical ensemble(NVE) has a
316   partition function like,
317   \begin{equation}
318 < \Omega (N,V,E) = e^{\beta TS}
319 < \label{introEqaution:NVEPartition}.
318 > \Omega (N,V,E) = e^{\beta TS} \label{introEquation:NVEPartition}.
319   \end{equation}
320   A canonical ensemble(NVT)is an ensemble of systems, each of which
321   can share its energy with a large heat reservoir. The distribution
# Line 635 | Line 634 | a \emph{symplectic} flow if it satisfies,
634   Let $\varphi$ be the flow of Hamiltonian vector field, $\varphi$ is
635   a \emph{symplectic} flow if it satisfies,
636   \begin{equation}
637 < '\varphi^T J '\varphi = J.
637 > {\varphi '}^T J \varphi ' = J.
638   \end{equation}
639   According to Liouville's theorem, the symplectic volume is invariant
640   under a Hamiltonian flow, which is the basis for classical
# Line 643 | Line 642 | symplectomorphism. As to the Poisson system,
642   field on a symplectic manifold can be shown to be a
643   symplectomorphism. As to the Poisson system,
644   \begin{equation}
645 < '\varphi ^T J '\varphi  = J \circ \varphi
645 > {\varphi '}^T J \varphi ' = J \circ \varphi
646   \end{equation}
647   is the property must be preserved by the integrator.
648  
# Line 661 | Line 660 | When designing any numerical methods, one should alway
660   In other words, the flow of this vector field is reversible if and
661   only if $ h \circ \varphi ^{ - 1}  = \varphi  \circ h $.
662  
663 < When designing any numerical methods, one should always try to
663 > A \emph{first integral}, or conserved quantity of a general
664 > differential function is a function $ G:R^{2d}  \to R^d $ which is
665 > constant for all solutions of the ODE $\frac{{dx}}{{dt}} = f(x)$ ,
666 > \[
667 > \frac{{dG(x(t))}}{{dt}} = 0.
668 > \]
669 > Using chain rule, one may obtain,
670 > \[
671 > \sum\limits_i {\frac{{dG}}{{dx_i }}} f_i (x) = f \bullet \nabla G,
672 > \]
673 > which is the condition for conserving \emph{first integral}. For a
674 > canonical Hamiltonian system, the time evolution of an arbitrary
675 > smooth function $G$ is given by,
676 > \begin{equation}
677 > \begin{array}{c}
678 > \frac{{dG(x(t))}}{{dt}} = [\nabla _x G(x(t))]^T \dot x(t) \\
679 >  = [\nabla _x G(x(t))]^T J\nabla _x H(x(t)). \\
680 > \end{array}
681 > \label{introEquation:firstIntegral1}
682 > \end{equation}
683 > Using poisson bracket notion, Equation
684 > \ref{introEquation:firstIntegral1} can be rewritten as
685 > \[
686 > \frac{d}{{dt}}G(x(t)) = \left\{ {G,H} \right\}(x(t)).
687 > \]
688 > Therefore, the sufficient condition for $G$ to be the \emph{first
689 > integral} of a Hamiltonian system is
690 > \[
691 > \left\{ {G,H} \right\} = 0.
692 > \]
693 > As well known, the Hamiltonian (or energy) H of a Hamiltonian system
694 > is a \emph{first integral}, which is due to the fact $\{ H,H\}  =
695 > 0$.
696 >
697 >
698 > When designing any numerical methods, one should always try to
699   preserve the structural properties of the original ODE and its flow.
700  
701   \subsection{\label{introSection:constructionSymplectic}Construction of Symplectic Methods}
# Line 685 | Line 719 | instability \cite{}. However, due to computational pen
719   ordinary implicit Runge-Kutta methods are not suitable for
720   Hamiltonian system. Recently, various high-order explicit
721   Runge--Kutta methods have been developed to overcome this
722 < instability \cite{}. However, due to computational penalty involved
723 < in implementing the Runge-Kutta methods, they do not attract too
724 < much attention from Molecular Dynamics community. Instead, splitting
725 < have been widely accepted since they exploit natural decompositions
726 < of the system\cite{Tuckerman92}.
722 > instability. However, due to computational penalty involved in
723 > implementing the Runge-Kutta methods, they do not attract too much
724 > attention from Molecular Dynamics community. Instead, splitting have
725 > been widely accepted since they exploit natural decompositions of
726 > the system\cite{Tuckerman92}.
727  
728   \subsubsection{\label{introSection:splittingMethod}Splitting Method}
729  
# Line 736 | Line 770 | _{1,h/2} ,
770   splitting gives a second-order decomposition,
771   \begin{equation}
772   \varphi _h  = \varphi _{1,h/2}  \circ \varphi _{2,h}  \circ \varphi
773 < _{1,h/2} ,
740 < \label{introEqaution:secondOrderSplitting}
773 > _{1,h/2} , \label{introEquation:secondOrderSplitting}
774   \end{equation}
775   which has a local error proportional to $h^3$. Sprang splitting's
776   popularity in molecular simulation community attribute to its
777   symmetric property,
778   \begin{equation}
779   \varphi _h^{ - 1} = \varphi _{ - h}.
780 < \lable{introEquation:timeReversible}
780 > \label{introEquation:timeReversible}
781   \end{equation}
782  
783   \subsubsection{\label{introSection:exampleSplittingMethod}Example of Splitting Method}
# Line 802 | Line 835 | q(\Delta t) = q(0) + \frac{{\Delta t}}{2}\left[ {\dot
835   \frac{{\Delta t}}{{2m}}\dot q(0)} \right], %
836   \label{introEquation:positionVerlet1} \\%
837   %
838 < q(\Delta t) = q(0) + \frac{{\Delta t}}{2}\left[ {\dot q(0) + \dot
838 > q(\Delta t) &= q(0) + \frac{{\Delta t}}{2}\left[ {\dot q(0) + \dot
839   q(\Delta t)} \right]. %
840   \label{introEquation:positionVerlet1}
841   \end{align}
# Line 828 | Line 861 | can obtain
861   \]
862   Applying Baker-Campbell-Hausdorff formula to Sprang splitting, we
863   can obtain
864 < \begin{eqnarray}
864 > \begin{eqnarray*}
865   \exp (h X/2)\exp (h Y)\exp (h X/2) & = & \exp (h X + h Y + h^2
866 < [X,Y]/4 + h^2 [Y,X]/4 \\ & & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 +
867 < h^3 [Y,[Y,X]]/12 - h^3 [X,[X,Y]]/24 +  \ldots )
868 < \end{eqnarray}
866 > [X,Y]/4 + h^2 [Y,X]/4 \\ & & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\
867 > & & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3 [X,[X,Y]]/24 & & \mbox{} +
868 > \ldots )
869 > \end{eqnarray*}
870   Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0\], the dominant local
871   error of Spring splitting is proportional to $h^3$. The same
872   procedure can be applied to general splitting,  of the form
# Line 871 | Line 905 | dynamical information.
905  
906   \subsection{\label{introSec:mdInit}Initialization}
907  
908 + \subsection{\label{introSec:forceEvaluation}Force Evaluation}
909 +
910   \subsection{\label{introSection:mdIntegration} Integration of the Equations of Motion}
911  
912   \section{\label{introSection:rigidBody}Dynamics of Rigid Bodies}
913  
914 < A rigid body is a body in which the distance between any two given
915 < points of a rigid body remains constant regardless of external
916 < forces exerted on it. A rigid body therefore conserves its shape
917 < during its motion.
914 > Rigid bodies are frequently involved in the modeling of different
915 > areas, from engineering, physics, to chemistry. For example,
916 > missiles and vehicle are usually modeled by rigid bodies.  The
917 > movement of the objects in 3D gaming engine or other physics
918 > simulator is governed by the rigid body dynamics. In molecular
919 > simulation, rigid body is used to simplify the model in
920 > protein-protein docking study{\cite{Gray03}}.
921  
922 < Applications of dynamics of rigid bodies.
922 > It is very important to develop stable and efficient methods to
923 > integrate the equations of motion of orientational degrees of
924 > freedom. Euler angles are the nature choice to describe the
925 > rotational degrees of freedom. However, due to its singularity, the
926 > numerical integration of corresponding equations of motion is very
927 > inefficient and inaccurate. Although an alternative integrator using
928 > different sets of Euler angles can overcome this difficulty\cite{},
929 > the computational penalty and the lost of angular momentum
930 > conservation still remain. A singularity free representation
931 > utilizing quaternions was developed by Evans in 1977. Unfortunately,
932 > this approach suffer from the nonseparable Hamiltonian resulted from
933 > quaternion representation, which prevents the symplectic algorithm
934 > to be utilized. Another different approach is to apply holonomic
935 > constraints to the atoms belonging to the rigid body. Each atom
936 > moves independently under the normal forces deriving from potential
937 > energy and constraint forces which are used to guarantee the
938 > rigidness. However, due to their iterative nature, SHAKE and Rattle
939 > algorithm converge very slowly when the number of constraint
940 > increases.
941  
942 + The break through in geometric literature suggests that, in order to
943 + develop a long-term integration scheme, one should preserve the
944 + symplectic structure of the flow. Introducing conjugate momentum to
945 + rotation matrix $A$ and re-formulating Hamiltonian's equation, a
946 + symplectic integrator, RSHAKE, was proposed to evolve the
947 + Hamiltonian system in a constraint manifold by iteratively
948 + satisfying the orthogonality constraint $A_t A = 1$. An alternative
949 + method using quaternion representation was developed by Omelyan.
950 + However, both of these methods are iterative and inefficient. In
951 + this section, we will present a symplectic Lie-Poisson integrator
952 + for rigid body developed by Dullweber and his
953 + coworkers\cite{Dullweber1997}.
954 +
955   \subsection{\label{introSection:lieAlgebra}Lie Algebra}
956  
957 < \subsection{\label{introSection:DLMMotionEquation}The Euler Equations of Rigid Body Motion}
957 > \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body}
958  
959 < \subsection{\label{introSection:otherRBMotionEquation}Other Formulations for Rigid Body Motion}
959 > \begin{equation}
960 > H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) +
961 > V(q,Q) + \frac{1}{2}tr[(QQ^T  - 1)\Lambda ].
962 > \label{introEquation:RBHamiltonian}
963 > \end{equation}
964 > Here, $q$ and $Q$  are the position and rotation matrix for the
965 > rigid-body, $p$ and $P$  are conjugate momenta to $q$  and $Q$ , and
966 > $J$, a diagonal matrix, is defined by
967 > \[
968 > I_{ii}^{ - 1}  = \frac{1}{2}\sum\limits_{i \ne j} {J_{jj}^{ - 1} }
969 > \]
970 > where $I_{ii}$ is the diagonal element of the inertia tensor. This
971 > constrained Hamiltonian equation subjects to a holonomic constraint,
972 > \begin{equation}
973 > Q^T Q = 1$, \label{introEquation:orthogonalConstraint}
974 > \end{equation}
975 > which is used to ensure rotation matrix's orthogonality.
976 > Differentiating \ref{introEquation:orthogonalConstraint} and using
977 > Equation \ref{introEquation:RBMotionMomentum}, one may obtain,
978 > \begin{equation}
979 > Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0 . \\
980 > \label{introEquation:RBFirstOrderConstraint}
981 > \end{equation}
982  
983 < \section{\label{introSection:correlationFunctions}Correlation Functions}
983 > Using Equation (\ref{introEquation:motionHamiltonianCoordinate},
984 > \ref{introEquation:motionHamiltonianMomentum}), one can write down
985 > the equations of motion,
986 > \[
987 > \begin{array}{c}
988 > \frac{{dq}}{{dt}} = \frac{p}{m} \label{introEquation:RBMotionPosition}\\
989 > \frac{{dp}}{{dt}} =  - \nabla _q V(q,Q) \label{introEquation:RBMotionMomentum}\\
990 > \frac{{dQ}}{{dt}} = PJ^{ - 1}  \label{introEquation:RBMotionRotation}\\
991 > \frac{{dP}}{{dt}} =  - \nabla _Q V(q,Q) - 2Q\Lambda . \label{introEquation:RBMotionP}\\
992 > \end{array}
993 > \]
994  
995 + In general, there are two ways to satisfy the holonomic constraints.
996 + We can use constraint force provided by lagrange multiplier on the
997 + normal manifold to keep the motion on constraint space. Or we can
998 + simply evolve the system in constraint manifold. The two method are
999 + proved to be equivalent. The holonomic constraint and equations of
1000 + motions define a constraint manifold for rigid body
1001 + \[
1002 + M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1}  + J^{ - 1} P^T Q = 0}
1003 + \right\}.
1004 + \]
1005 +
1006 + Unfortunately, this constraint manifold is not the cotangent bundle
1007 + $T_{\star}SO(3)$. However, it turns out that under symplectic
1008 + transformation, the cotangent space and the phase space are
1009 + diffeomorphic. Introducing
1010 + \[
1011 + \tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right),
1012 + \]
1013 + the mechanical system subject to a holonomic constraint manifold $M$
1014 + can be re-formulated as a Hamiltonian system on the cotangent space
1015 + \[
1016 + T^* SO(3) = \left\{ {(\tilde Q,\tilde P):\tilde Q^T \tilde Q =
1017 + 1,\tilde Q^T \tilde PJ^{ - 1}  + J^{ - 1} P^T \tilde Q = 0} \right\}
1018 + \]
1019 +
1020 + For a body fixed vector $X_i$ with respect to the center of mass of
1021 + the rigid body, its corresponding lab fixed vector $X_0^{lab}$  is
1022 + given as
1023 + \begin{equation}
1024 + X_i^{lab} = Q X_i + q.
1025 + \end{equation}
1026 + Therefore, potential energy $V(q,Q)$ is defined by
1027 + \[
1028 + V(q,Q) = V(Q X_0 + q).
1029 + \]
1030 + Hence,
1031 + \[
1032 + \nabla _q V(q,Q) = F(q,Q) = \sum\limits_i {F_i (q,Q)}
1033 + \]
1034 +
1035 + \[
1036 + \nabla _Q V(q,Q) = F(q,Q)X_i^t
1037 + \]
1038 +
1039 + As a common choice to describe the rotation dynamics of the rigid
1040 + body, angular momentum on body frame $\Pi  = Q^t P$ is introduced to
1041 + rewrite the equations of motion,
1042 + \begin{equation}
1043 + \begin{array}{l}
1044 + \mathop \Pi \limits^ \bullet   = J^{ - 1} \Pi ^T \Pi  + Q^T \sum\limits_i {F_i (q,Q)X_i^T }  - \Lambda  \\
1045 + \mathop Q\limits^{{\rm{   }} \bullet }  = Q\Pi {\rm{ }}J^{ - 1}  \\
1046 + \end{array}
1047 + \label{introEqaution:RBMotionPI}
1048 + \end{equation}
1049 + , as well as holonomic constraints,
1050 + \[
1051 + \begin{array}{l}
1052 + \Pi J^{ - 1}  + J^{ - 1} \Pi ^t  = 0 \\
1053 + Q^T Q = 1 \\
1054 + \end{array}
1055 + \]
1056 +
1057 + For a vector $v(v_1 ,v_2 ,v_3 ) \in R^3$ and a matrix $\hat v \in
1058 + so(3)^ \star$, the hat-map isomorphism,
1059 + \begin{equation}
1060 + v(v_1 ,v_2 ,v_3 ) \Leftrightarrow \hat v = \left(
1061 + {\begin{array}{*{20}c}
1062 +   0 & { - v_3 } & {v_2 }  \\
1063 +   {v_3 } & 0 & { - v_1 }  \\
1064 +   { - v_2 } & {v_1 } & 0  \\
1065 + \end{array}} \right),
1066 + \label{introEquation:hatmapIsomorphism}
1067 + \end{equation}
1068 + will let us associate the matrix products with traditional vector
1069 + operations
1070 + \[
1071 + \hat vu = v \times u
1072 + \]
1073 +
1074 + Using \ref{introEqaution:RBMotionPI}, one can construct a skew
1075 + matrix,
1076 + \begin{equation}
1077 + (\mathop \Pi \limits^ \bullet   - \mathop \Pi \limits^ \bullet  ^T
1078 + ){\rm{ }} = {\rm{ }}(\Pi  - \Pi ^T ){\rm{ }}(J^{ - 1} \Pi  + \Pi J^{
1079 + - 1} ) + \sum\limits_i {[Q^T F_i (r,Q)X_i^T  - X_i F_i (r,Q)^T Q]} -
1080 + (\Lambda  - \Lambda ^T ) . \label{introEquation:skewMatrixPI}
1081 + \end{equation}
1082 + Since $\Lambda$ is symmetric, the last term of Equation
1083 + \ref{introEquation:skewMatrixPI}, which implies the Lagrange
1084 + multiplier $\Lambda$ is ignored in the integration.
1085 +
1086 + Hence, applying hat-map isomorphism, we obtain the equation of
1087 + motion for angular momentum on body frame
1088 + \[
1089 + \dot \pi  = \pi  \times I^{ - 1} \pi  + Q^T \sum\limits_i {F_i (r,Q)
1090 + \times X_i }
1091 + \]
1092 + In the same manner, the equation of motion for rotation matrix is
1093 + given by
1094 + \[
1095 + \dot Q = Qskew(M^{ - 1} \pi )
1096 + \]
1097 +
1098 + The free rigid body equation is an example of a non-canonical
1099 + Hamiltonian system.
1100 +
1101 + \subsection{\label{introSection:symplecticDiscretizationRB}Symplectic Integration of Euler Equations}
1102 +
1103 + \[
1104 + \varphi _{\Delta t}  = \varphi _{\Delta t/2,V}  \circ \varphi
1105 + _{\Delta t,T}  \circ \varphi _{\Delta t/2,V}
1106 + \]
1107 +
1108 + \[
1109 + \varphi _{\Delta t,T}  = \varphi _{\Delta t,R}  \circ \varphi
1110 + _{\Delta t,\pi }
1111 + \]
1112 +
1113 + \[
1114 + \varphi _{\Delta t,\pi }  = \varphi _{\Delta t/2,\pi _1 }  \circ
1115 + \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t,\pi _3 }
1116 + \circ \varphi _{\Delta t/2,\pi _2 }  \circ \varphi _{\Delta t/2,\pi
1117 + _1 }
1118 + \]
1119 +
1120 + \[
1121 + \varphi _{\Delta t/2,V}  = \varphi _{\Delta t/2,F}  \circ \varphi
1122 + _{\Delta t/2,\tau }
1123 + \]
1124 +
1125 +
1126   \section{\label{introSection:langevinDynamics}Langevin Dynamics}
1127  
1128   \subsection{\label{introSection:LDIntroduction}Introduction and application of Langevin Dynamics}
# Line 1097 | Line 1330 | Body}
1330  
1331   \subsection{\label{introSection:centersRigidBody}Centers of Rigid
1332   Body}
1333 +
1334 + \section{\label{introSection:correlationFunctions}Correlation Functions}

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines