--- trunk/tengDissertation/Introduction.tex 2006/07/17 18:13:26 2940 +++ trunk/tengDissertation/Introduction.tex 2006/07/17 20:01:05 2941 @@ -67,8 +67,8 @@ the quality of numerical integration schemes for rigid \begin{equation}E = T + V. \label{introEquation:energyConservation} \end{equation} All of these conserved quantities are important factors to determine -the quality of numerical integration schemes for rigid bodies -\cite{Dullweber1997}. +the quality of numerical integration schemes for rigid +bodies.\cite{Dullweber1997} \subsection{\label{introSection:lagrangian}Lagrangian Mechanics} @@ -178,7 +178,7 @@ known as the canonical equations of motions \cite{Gold where Eq.~\ref{introEquation:motionHamiltonianCoordinate} and Eq.~\ref{introEquation:motionHamiltonianMomentum} are Hamilton's equation of motion. Due to their symmetrical formula, they are also -known as the canonical equations of motions \cite{Goldstein2001}. +known as the canonical equations of motions.\cite{Goldstein2001} An important difference between Lagrangian approach and the Hamiltonian approach is that the Lagrangian is considered to be a @@ -188,7 +188,7 @@ only works with 1st-order differential equations\cite{ Hamiltonian Mechanics is more appropriate for application to statistical mechanics and quantum mechanics, since it treats the coordinate and its time derivative as independent variables and it -only works with 1st-order differential equations\cite{Marion1990}. +only works with 1st-order differential equations.\cite{Marion1990} In Newtonian Mechanics, a system described by conservative forces conserves the total energy (Eq.~\ref{introEquation:energyConservation}). It follows that @@ -416,7 +416,7 @@ statistical ensemble are identical \cite{Frenkel1996, many-body system in Statistical Mechanics. Fortunately, the Ergodic Hypothesis makes a connection between time average and the ensemble average. It states that the time average and average over the -statistical ensemble are identical \cite{Frenkel1996, Leach2001}: +statistical ensemble are identical:\cite{Frenkel1996, Leach2001} \begin{equation} \langle A(q , p) \rangle_t = \mathop {\lim }\limits_{t \to \infty } \frac{1}{t}\int\limits_0^t {A(q(t),p(t))dt = \int\limits_\Gamma @@ -434,21 +434,21 @@ choice\cite{Frenkel1996}. utilized. Or if the system lends itself to a time averaging approach, the Molecular Dynamics techniques in Sec.~\ref{introSection:molecularDynamics} will be the best -choice\cite{Frenkel1996}. +choice.\cite{Frenkel1996} \section{\label{introSection:geometricIntegratos}Geometric Integrators} A variety of numerical integrators have been proposed to simulate the motions of atoms in MD simulation. They usually begin with -initial conditions and move the objects in the direction governed -by the differential equations. However, most of them ignore the -hidden physical laws contained within the equations. Since 1990, -geometric integrators, which preserve various phase-flow invariants -such as symplectic structure, volume and time reversal symmetry, -were developed to address this issue\cite{Dullweber1997, -McLachlan1998, Leimkuhler1999}. The velocity Verlet method, which -happens to be a simple example of symplectic integrator, continues -to gain popularity in the molecular dynamics community. This fact -can be partly explained by its geometric nature. +initial conditions and move the objects in the direction governed by +the differential equations. However, most of them ignore the hidden +physical laws contained within the equations. Since 1990, geometric +integrators, which preserve various phase-flow invariants such as +symplectic structure, volume and time reversal symmetry, were +developed to address this issue.\cite{Dullweber1997, McLachlan1998, +Leimkuhler1999} The velocity Verlet method, which happens to be a +simple example of symplectic integrator, continues to gain +popularity in the molecular dynamics community. This fact can be +partly explained by its geometric nature. \subsection{\label{introSection:symplecticManifold}Symplectic Manifolds} A \emph{manifold} is an abstract mathematical space. It looks @@ -457,19 +457,19 @@ apply calculus\cite{Hirsch1997}. A \emph{symplectic ma surface of Earth. It seems to be flat locally, but it is round if viewed as a whole. A \emph{differentiable manifold} (also known as \emph{smooth manifold}) is a manifold on which it is possible to -apply calculus\cite{Hirsch1997}. A \emph{symplectic manifold} is +apply calculus.\cite{Hirsch1997} A \emph{symplectic manifold} is defined as a pair $(M, \omega)$ which consists of a \emph{differentiable manifold} $M$ and a close, non-degenerate, bilinear symplectic form, $\omega$. A symplectic form on a vector space $V$ is a function $\omega(x, y)$ which satisfies $\omega(\lambda_1x_1+\lambda_2x_2, y) = \lambda_1\omega(x_1, y)+ \lambda_2\omega(x_2, y)$, $\omega(x, y) = - \omega(y, x)$ and -$\omega(x, x) = 0$\cite{McDuff1998}. The cross product operation in +$\omega(x, x) = 0$.\cite{McDuff1998} The cross product operation in vector field is an example of symplectic form. One of the motivations to study \emph{symplectic manifolds} in Hamiltonian Mechanics is that a symplectic manifold can represent all possible configurations of the system and the phase space of the system can -be described by it's cotangent bundle\cite{Jost2002}. Every +be described by it's cotangent bundle.\cite{Jost2002} Every symplectic manifold is even dimensional. For instance, in Hamilton equations, coordinate and momentum always appear in pairs. @@ -496,11 +496,12 @@ Hamiltonian dynamics is Poisson Dynamics\cite{Olver198 \label{introEquation:compactHamiltonian} \end{equation}In this case, $f$ is called a \emph{Hamiltonian vector field}. Another generalization of -Hamiltonian dynamics is Poisson Dynamics\cite{Olver1986}, +Hamiltonian dynamics is Poisson Dynamics,\cite{Olver1986} \begin{equation} \dot x = J(x)\nabla _x H \label{introEquation:poissonHamiltonian} \end{equation} -The most obvious change being that matrix $J$ now depends on $x$. +where the most obvious change being that matrix $J$ now depends on +$x$. \subsection{\label{introSection:exactFlow}Exact Propagator} @@ -620,7 +621,7 @@ accurately\cite{Kane2000}. Since they are geometricall Generating functions\cite{Channell1990} tend to lead to methods which are cumbersome and difficult to use. In dissipative systems, variational methods can capture the decay of energy -accurately\cite{Kane2000}. Since they are geometrically unstable +accurately.\cite{Kane2000} Since they are geometrically unstable against non-Hamiltonian perturbations, ordinary implicit Runge-Kutta methods are not suitable for Hamiltonian system. Recently, various high-order explicit Runge-Kutta methods \cite{Owren1992,Chen2003} @@ -629,7 +630,7 @@ system\cite{Tuckerman1992, McLachlan1998}. methods, they have not attracted much attention from the Molecular Dynamics community. Instead, splitting methods have been widely accepted since they exploit natural decompositions of the -system\cite{Tuckerman1992, McLachlan1998}. +system.\cite{McLachlan1998, Tuckerman1992} \subsubsection{\label{introSection:splittingMethod}\textbf{Splitting Methods}} @@ -673,7 +674,7 @@ a second-order decomposition, The Lie-Trotter splitting(Eq.~\ref{introEquation:firstOrderSplitting}) introduces local errors proportional to $h^2$, while the Strang splitting gives -a second-order decomposition, +a second-order decomposition,\cite{Strang1968} \begin{equation} \varphi _h = \varphi _{1,h/2} \circ \varphi _{2,h} \circ \varphi _{1,h/2} , \label{introEquation:secondOrderSplitting} @@ -748,11 +749,13 @@ The Baker-Campbell-Hausdorff formula can be used to de \subsubsection{\label{introSection:errorAnalysis}\textbf{Error Analysis and Higher Order Methods}} -The Baker-Campbell-Hausdorff formula can be used to determine the -local error of a splitting method in terms of the commutator of the -operators(Eq.~\ref{introEquation:exponentialOperator}) associated with -the sub-propagator. For operators $hX$ and $hY$ which are associated -with $\varphi_1(t)$ and $\varphi_2(t)$ respectively , we have +The Baker-Campbell-Hausdorff formula\cite{Gilmore1974} can be used +to determine the local error of a splitting method in terms of the +commutator of the +operators(Eq.~\ref{introEquation:exponentialOperator}) associated +with the sub-propagator. For operators $hX$ and $hY$ which are +associated with $\varphi_1(t)$ and $\varphi_2(t)$ respectively , we +have \begin{equation} \exp (hX + hY) = \exp (hZ) \end{equation} @@ -782,7 +785,7 @@ order methods based on symmetric splitting\cite{Yoshid \end{equation} A careful choice of coefficient $a_1 \ldots b_m$ will lead to higher order methods. Yoshida proposed an elegant way to compose higher -order methods based on symmetric splitting\cite{Yoshida1990}. Given +order methods based on symmetric splitting.\cite{Yoshida1990} Given a symmetric second order base method $ \varphi _h^{(2)} $, a fourth-order symmetric method can be constructed by composing, \[ @@ -943,9 +946,9 @@ than a predetermined distance are not included in the %cutoff and minimum image convention Another important technique to improve the efficiency of force evaluation is to apply spherical cutoffs where particles farther -than a predetermined distance are not included in the calculation -\cite{Frenkel1996}. The use of a cutoff radius will cause a -discontinuity in the potential energy curve. Fortunately, one can +than a predetermined distance are not included in the +calculation.\cite{Frenkel1996} The use of a cutoff radius will cause +a discontinuity in the potential energy curve. Fortunately, one can shift a simple radial potential to ensure the potential curve go smoothly to zero at the cutoff radius. The cutoff strategy works well for Lennard-Jones interaction because of its short range @@ -954,9 +957,9 @@ periodicity artifacts in liquid simulations. Taking ad in simulations. The Ewald summation, in which the slowly decaying Coulomb potential is transformed into direct and reciprocal sums with rapid and absolute convergence, has proved to minimize the -periodicity artifacts in liquid simulations. Taking advantage -of fast Fourier transform (FFT) techniques for calculating discrete Fourier -transforms, the particle mesh-based +periodicity artifacts in liquid simulations. Taking advantage of +fast Fourier transform (FFT) techniques for calculating discrete +Fourier transforms, the particle mesh-based methods\cite{Hockney1981,Shimada1993, Luty1994} are accelerated from $O(N^{3/2})$ to $O(N logN)$. An alternative approach is the \emph{fast multipole method}\cite{Greengard1987, Greengard1994}, @@ -966,7 +969,7 @@ his coworkers\cite{Wolf1999}. The shifted Coulomb pote simulation community, these two methods are difficult to implement correctly and efficiently. Instead, we use a damped and charge-neutralized Coulomb potential method developed by Wolf and -his coworkers\cite{Wolf1999}. The shifted Coulomb potential for +his coworkers.\cite{Wolf1999} The shifted Coulomb potential for particle $i$ and particle $j$ at distance $r_{rj}$ is given by: \begin{equation} V(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha @@ -1029,14 +1032,14 @@ experiments and integrating over the surface factor function}, is of most fundamental importance to liquid theory. Experimentally, pair distribution functions can be gathered by Fourier transforming raw data from a series of neutron diffraction -experiments and integrating over the surface factor -\cite{Powles1973}. The experimental results can serve as a criterion -to justify the correctness of a liquid model. Moreover, various -equilibrium thermodynamic and structural properties can also be -expressed in terms of the radial distribution function -\cite{Allen1987}. The pair distribution functions $g(r)$ gives the -probability that a particle $i$ will be located at a distance $r$ -from a another particle $j$ in the system +experiments and integrating over the surface +factor.\cite{Powles1973} The experimental results can serve as a +criterion to justify the correctness of a liquid model. Moreover, +various equilibrium thermodynamic and structural properties can also +be expressed in terms of the radial distribution +function.\cite{Allen1987} The pair distribution functions $g(r)$ +gives the probability that a particle $i$ will be located at a +distance $r$ from a another particle $j$ in the system \begin{equation} g(r) = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j \ne i} {\delta (r - r_{ij} )} } } \right\rangle = \frac{\rho @@ -1097,7 +1100,7 @@ docking studies\cite{Gray2003}. movement of the objects in 3D gaming engines or other physics simulators is governed by rigid body dynamics. In molecular simulations, rigid bodies are used to simplify protein-protein -docking studies\cite{Gray2003}. +docking studies.\cite{Gray2003} It is very important to develop stable and efficient methods to integrate the equations of motion for orientational degrees of @@ -1109,7 +1112,7 @@ quaternions was developed by Evans in 1977\cite{Evans1 angles can overcome this difficulty\cite{Barojas1973}, the computational penalty and the loss of angular momentum conservation still remain. A singularity-free representation utilizing -quaternions was developed by Evans in 1977\cite{Evans1977}. +quaternions was developed by Evans in 1977.\cite{Evans1977} Unfortunately, this approach used a nonseparable Hamiltonian resulting from the quaternion representation, which prevented the symplectic algorithm from being utilized. Another different approach @@ -1118,7 +1121,7 @@ number of constraints increases\cite{Ryckaert1977, And deriving from potential energy and constraint forces which are used to guarantee the rigidness. However, due to their iterative nature, the SHAKE and Rattle algorithms also converge very slowly when the -number of constraints increases\cite{Ryckaert1977, Andersen1983}. +number of constraints increases.\cite{Ryckaert1977, Andersen1983} A break-through in geometric literature suggests that, in order to develop a long-term integration scheme, one should preserve the @@ -1128,13 +1131,13 @@ developed by Omelyan\cite{Omelyan1998}. However, both proposed to evolve the Hamiltonian system in a constraint manifold by iteratively satisfying the orthogonality constraint $Q^T Q = 1$. An alternative method using the quaternion representation was -developed by Omelyan\cite{Omelyan1998}. However, both of these +developed by Omelyan.\cite{Omelyan1998} However, both of these methods are iterative and inefficient. In this section, we descibe a symplectic Lie-Poisson integrator for rigid bodies developed by Dullweber and his coworkers\cite{Dullweber1997} in depth. \subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Bodies} -The Hamiltonian of a rigid body is given by +The Hamiltonian of a rigid body is given by \begin{equation} H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) + V(q,Q) + \frac{1}{2}tr[(QQ^T - 1)\Lambda ]. @@ -1248,8 +1251,8 @@ iterations which can not be avoided in other methods\c Eq.~\ref{introEquation:skewMatrixPI} is zero, which implies the Lagrange multiplier $\Lambda$ is absent from the equations of motion. This unique property eliminates the requirement of -iterations which can not be avoided in other methods\cite{Kol1997, -Omelyan1998}. Applying the hat-map isomorphism, we obtain the +iterations which can not be avoided in other methods.\cite{Kol1997, +Omelyan1998} Applying the hat-map isomorphism, we obtain the equation of motion for angular momentum in the body frame \begin{equation} \dot \pi = \pi \times I^{ - 1} \pi + \sum\limits_i {\left( {Q^T @@ -1355,7 +1358,7 @@ norm of the angular momentum, $\parallel \pi function $G$ is zero, $F$ is a \emph{Casimir}, which is the conserved quantity in Poisson system. We can easily verify that the norm of the angular momentum, $\parallel \pi -\parallel$, is a \emph{Casimir}\cite{McLachlan1993}. Let $F(\pi ) = S(\frac{{\parallel +\parallel$, is a \emph{Casimir}.\cite{McLachlan1993} Let $F(\pi ) = S(\frac{{\parallel \pi \parallel ^2 }}{2})$ for an arbitrary function $ S:R \to R$ , then by the chain rule \[ @@ -1376,7 +1379,7 @@ listed in Table~\ref{introTable:rbEquations}. The Hamiltonian of rigid body can be separated in terms of kinetic energy and potential energy, $H = T(p,\pi ) + V(q,Q)$. The equations of motion corresponding to potential energy and kinetic energy are -listed in Table~\ref{introTable:rbEquations}. +listed in Table~\ref{introTable:rbEquations}. \begin{table} \caption{EQUATIONS OF MOTION DUE TO POTENTIAL AND KINETIC ENERGIES} \label{introTable:rbEquations}