--- trunk/tengDissertation/Introduction.tex 2006/06/11 01:55:48 2850 +++ trunk/tengDissertation/Introduction.tex 2006/06/21 16:28:25 2872 @@ -496,7 +496,7 @@ Leimkuhler1999}. The velocity verlet method, which hap geometric integrators, which preserve various phase-flow invariants such as symplectic structure, volume and time reversal symmetry, are developed to address this issue\cite{Dullweber1997, McLachlan1998, -Leimkuhler1999}. The velocity verlet method, which happens to be a +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. @@ -591,18 +591,18 @@ Instead, we use a approximate map, $\psi_\tau$, which \end{equation} In most cases, it is not easy to find the exact flow $\varphi_\tau$. -Instead, we use a approximate map, $\psi_\tau$, which is usually +Instead, we use an approximate map, $\psi_\tau$, which is usually called integrator. The order of an integrator $\psi_\tau$ is $p$, if the Taylor series of $\psi_\tau$ agree to order $p$, \begin{equation} -\psi_tau(x) = x + \tau f(x) + O(\tau^{p+1}) +\psi_\tau(x) = x + \tau f(x) + O(\tau^{p+1}) \end{equation} \subsection{\label{introSection:geometricProperties}Geometric Properties} -The hidden geometric properties\cite{Budd1999, Marsden1998} of ODE -and its flow play important roles in numerical studies. Many of them -can be found in systems which occur naturally in applications. +The hidden geometric properties\cite{Budd1999, Marsden1998} of an +ODE and its flow play important roles in numerical studies. Many of +them can be found in systems which occur naturally in applications. Let $\varphi$ be the flow of Hamiltonian vector field, $\varphi$ is a \emph{symplectic} flow if it satisfies, @@ -617,15 +617,15 @@ is the property must be preserved by the integrator. \begin{equation} {\varphi '}^T J \varphi ' = J \circ \varphi \end{equation} -is the property must be preserved by the integrator. +is the property that must be preserved by the integrator. It is possible to construct a \emph{volume-preserving} flow for a -source free($ \nabla \cdot f = 0 $) ODE, if the flow satisfies $ +source free ODE ($ \nabla \cdot f = 0 $), if the flow satisfies $ \det d\varphi = 1$. One can show easily that a symplectic flow will be volume-preserving. -Changing the variables $y = h(x)$ in a ODE\ref{introEquation:ODE} -will result in a new system, +Changing the variables $y = h(x)$ in an ODE +(Eq.~\ref{introEquation:ODE}) will result in a new system, \[ \dot y = \tilde f(y) = ((dh \cdot f)h^{ - 1} )(y). \] @@ -675,7 +675,7 @@ constructed. The most famous example is the Verlet-lea A lot of well established and very effective numerical methods have been successful precisely because of their symplecticities even though this fact was not recognized when they were first -constructed. The most famous example is the Verlet-leapfrog methods +constructed. The most famous example is the Verlet-leapfrog method in molecular dynamics. In general, symplectic integrators can be constructed using one of four different methods. \begin{enumerate} @@ -754,14 +754,14 @@ to its symmetric property, \label{introEquation:timeReversible} \end{equation},appendixFig:architecture -\subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Example of Splitting Method}} +\subsubsection{\label{introSection:exampleSplittingMethod}\textbf{Examples of the Splitting Method}} The classical equation for a system consisting of interacting particles can be written in Hamiltonian form, \[ H = T + V \] where $T$ is the kinetic energy and $V$ is the potential energy. -Setting $H_1 = T, H_2 = V$ and applying Strang splitting, one +Setting $H_1 = T, H_2 = V$ and applying the Strang splitting, one obtains the following: \begin{align} q(\Delta t) &= q(0) + \dot{q}(0)\Delta t + @@ -788,7 +788,7 @@ q(\Delta t) &= q(0) + \Delta t\, \dot{q}\biggl (\frac{ \label{introEquation:Lp9b}\\% % \dot{q}(\Delta t) &= \dot{q}\biggl (\frac{\Delta t}{2}\biggr ) + - \frac{\Delta t}{2m}\, F[q(0)]. \label{introEquation:Lp9c} + \frac{\Delta t}{2m}\, F[q(t)]. \label{introEquation:Lp9c} \end{align} From the preceding splitting, one can see that the integration of the equations of motion would follow: @@ -797,13 +797,14 @@ the equations of motion would follow: \item Use the half step velocities to move positions one whole step, $\Delta t$. -\item Evaluate the forces at the new positions, $\mathbf{r}(\Delta t)$, and use the new forces to complete the velocity move. +\item Evaluate the forces at the new positions, $\mathbf{q}(\Delta t)$, and use the new forces to complete the velocity move. \item Repeat from step 1 with the new position, velocities, and forces assuming the roles of the initial values. \end{enumerate} -Simply switching the order of splitting and composing, a new -integrator, the \emph{position verlet} integrator, can be generated, +By simply switching the order of the propagators in the splitting +and composing a new integrator, the \emph{position verlet} +integrator, can be generated, \begin{align} \dot q(\Delta t) &= \dot q(0) + \Delta tF(q(0))\left[ {q(0) + \frac{{\Delta t}}{{2m}}\dot q(0)} \right], % @@ -816,10 +817,10 @@ Baker-Campbell-Hausdorff formula can be used to determ \subsubsection{\label{introSection:errorAnalysis}\textbf{Error Analysis and Higher Order Methods}} -Baker-Campbell-Hausdorff formula can be used to determine the local -error of splitting method in terms of commutator of the +The Baker-Campbell-Hausdorff formula can be used to determine the +local error of splitting method in terms of the commutator of the operators(\ref{introEquation:exponentialOperator}) associated with -the sub-flow. For operators $hX$ and $hY$ which are associate to +the sub-flow. 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) @@ -833,22 +834,22 @@ Applying Baker-Campbell-Hausdorff formula\cite{Varadar \[ [X,Y] = XY - YX . \] -Applying Baker-Campbell-Hausdorff formula\cite{Varadarajan1974} to -Sprang splitting, we can obtain +Applying the Baker-Campbell-Hausdorff formula\cite{Varadarajan1974} +to the Sprang splitting, we can obtain \begin{eqnarray*} \exp (h X/2)\exp (h Y)\exp (h X/2) & = & \exp (h X + h Y + h^2 [X,Y]/4 + h^2 [Y,X]/4 \\ & & \mbox{} + h^2 [X,X]/8 + h^2 [Y,Y]/8 \\ & & \mbox{} + h^3 [Y,[Y,X]]/12 - h^3[X,[X,Y]]/24 + \ldots ) \end{eqnarray*} -Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0\], the dominant local +Since \[ [X,Y] + [Y,X] = 0\] and \[ [X,X] = 0,\] the dominant local error of Spring splitting is proportional to $h^3$. The same -procedure can be applied to general splitting, of the form +procedure can be applied to a general splitting, of the form \begin{equation} \varphi _{b_m h}^2 \circ \varphi _{a_m h}^1 \circ \varphi _{b_{m - 1} h}^2 \circ \ldots \circ \varphi _{a_1 h}^1 . \end{equation} -Careful choice of coefficient $a_1 \ldots b_m$ will lead to higher -order method. Yoshida proposed an elegant way to compose higher +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 a symmetric second order base method $ \varphi _h^{(2)} $, a fourth-order symmetric method can be constructed by composing, @@ -861,9 +862,9 @@ _{\beta h}^{(2n)} \circ \varphi _{\alpha h}^{(2n)} integrator $ \varphi _h^{(2n + 2)}$ can be composed by \begin{equation} \varphi _h^{(2n + 2)} = \varphi _{\alpha h}^{(2n)} \circ \varphi -_{\beta h}^{(2n)} \circ \varphi _{\alpha h}^{(2n)} +_{\beta h}^{(2n)} \circ \varphi _{\alpha h}^{(2n)}, \end{equation} -, if the weights are chosen as +if the weights are chosen as \[ \alpha = - \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }},\beta = \frac{{2^{1/(2n + 1)} }}{{2 - 2^{1/(2n + 1)} }} . @@ -901,7 +902,7 @@ will discusses issues in production run. These three individual steps will be covered in the following sections. Sec.~\ref{introSec:initialSystemSettings} deals with the initialization of a simulation. Sec.~\ref{introSection:production} -will discusses issues in production run. +will discusse issues in production run. Sec.~\ref{introSection:Analysis} provides the theoretical tools for trajectory analysis. @@ -914,50 +915,52 @@ purification and crystallization. Even for the molecul databases, such as RCSB Protein Data Bank \textit{etc}. Although thousands of crystal structures of molecules are discovered every year, many more remain unknown due to the difficulties of -purification and crystallization. Even for the molecule with known -structure, some important information is missing. For example, the +purification and crystallization. Even for molecules with known +structure, some important information is missing. For example, a missing hydrogen atom which acts as donor in hydrogen bonding must be added. Moreover, in order to include electrostatic interaction, one may need to specify the partial charges for individual atoms. Under some circumstances, we may even need to prepare the system in -a special setup. For instance, when studying transport phenomenon in -membrane system, we may prepare the lipids in bilayer structure -instead of placing lipids randomly in solvent, since we are not -interested in self-aggregation and it takes a long time to happen. +a special configuration. For instance, when studying transport +phenomenon in membrane systems, we may prepare the lipids in a +bilayer structure instead of placing lipids randomly in solvent, +since we are not interested in the slow self-aggregation process. \subsubsection{\textbf{Minimization}} It is quite possible that some of molecules in the system from -preliminary preparation may be overlapped with each other. This -close proximity leads to high potential energy which consequently -jeopardizes any molecular dynamics simulations. To remove these -steric overlaps, one typically performs energy minimization to find -a more reasonable conformation. Several energy minimization methods -have been developed to exploit the energy surface and to locate the -local minimum. While converging slowly near the minimum, steepest -descent method is extremely robust when systems are far from -harmonic. Thus, it is often used to refine structure from -crystallographic data. Relied on the gradient or hessian, advanced -methods like conjugate gradient and Newton-Raphson converge rapidly -to a local minimum, while become unstable if the energy surface is -far from quadratic. Another factor must be taken into account, when +preliminary preparation may be overlapping with each other. This +close proximity leads to high initial potential energy which +consequently jeopardizes any molecular dynamics simulations. To +remove these steric overlaps, one typically performs energy +minimization to find a more reasonable conformation. Several energy +minimization methods have been developed to exploit the energy +surface and to locate the local minimum. While converging slowly +near the minimum, steepest descent method is extremely robust when +systems are strongly anharmonic. Thus, it is often used to refine +structure from crystallographic data. Relied on the gradient or +hessian, advanced methods like Newton-Raphson converge rapidly to a +local minimum, but become unstable if the energy surface is far from +quadratic. Another factor that must be taken into account, when choosing energy minimization method, is the size of the system. Steepest descent and conjugate gradient can deal with models of any -size. Because of the limit of computation power to calculate hessian -matrix and insufficient storage capacity to store them, most -Newton-Raphson methods can not be used with very large models. +size. Because of the limits on computer memory to store the hessian +matrix and the computing power needed to diagonalized these +matrices, most Newton-Raphson methods can not be used with very +large systems. \subsubsection{\textbf{Heating}} Typically, Heating is performed by assigning random velocities -according to a Gaussian distribution for a temperature. Beginning at -a lower temperature and gradually increasing the temperature by -assigning greater random velocities, we end up with setting the -temperature of the system to a final temperature at which the -simulation will be conducted. In heating phase, we should also keep -the system from drifting or rotating as a whole. Equivalently, the -net linear momentum and angular momentum of the system should be -shifted to zero. +according to a Maxwell-Boltzman distribution for a desired +temperature. Beginning at a lower temperature and gradually +increasing the temperature by assigning larger random velocities, we +end up with setting the temperature of the system to a final +temperature at which the simulation will be conducted. In heating +phase, we should also keep the system from drifting or rotating as a +whole. To do this, the net linear momentum and angular momentum of +the system is shifted to zero after each resampling from the Maxwell +-Boltzman distribution. \subsubsection{\textbf{Equilibration}} @@ -973,30 +976,32 @@ Production run is the most important step of the simul \subsection{\label{introSection:production}Production} -Production run is the most important step of the simulation, in +The production run is the most important step of the simulation, in which the equilibrated structure is used as a starting point and the motions of the molecules are collected for later analysis. In order to capture the macroscopic properties of the system, the molecular -dynamics simulation must be performed in correct and efficient way. +dynamics simulation must be performed by sampling correctly and +efficiently from the relevant thermodynamic ensemble. The most expensive part of a molecular dynamics simulation is the calculation of non-bonded forces, such as van der Waals force and Coulombic forces \textit{etc}. For a system of $N$ particles, the complexity of the algorithm for pair-wise interactions is $O(N^2 )$, which making large simulations prohibitive in the absence of any -computation saving techniques. +algorithmic tricks. -A natural approach to avoid system size issue is to represent the +A natural approach to avoid system size issues is to represent the bulk behavior by a finite number of the particles. However, this -approach will suffer from the surface effect. To offset this, -\textit{Periodic boundary condition} (see Fig.~\ref{introFig:pbc}) -is developed to simulate bulk properties with a relatively small -number of particles. In this method, the simulation box is -replicated throughout space to form an infinite lattice. During the -simulation, when a particle moves in the primary cell, its image in -other cells move in exactly the same direction with exactly the same -orientation. Thus, as a particle leaves the primary cell, one of its -images will enter through the opposite face. +approach will suffer from the surface effect at the edges of the +simulation. To offset this, \textit{Periodic boundary conditions} +(see Fig.~\ref{introFig:pbc}) is developed to simulate bulk +properties with a relatively small number of particles. In this +method, the simulation box is replicated throughout space to form an +infinite lattice. During the simulation, when a particle moves in +the primary cell, its image in other cells move in exactly the same +direction with exactly the same orientation. Thus, as a particle +leaves the primary cell, one of its images will enter through the +opposite face. \begin{figure} \centering \includegraphics[width=\linewidth]{pbc.eps} @@ -1008,32 +1013,32 @@ evaluation is to apply cutoff where particles farther %cutoff and minimum image convention Another important technique to improve the efficiency of force -evaluation is to apply cutoff where particles farther than a -predetermined distance, are not included in the calculation +evaluation is to apply spherical cutoff 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 -shift the potential to ensure the potential curve go smoothly to -zero at the cutoff radius. Cutoff strategy works pretty well for -Lennard-Jones interaction because of its short range nature. -However, simply truncating the electrostatic interaction with the -use of cutoff has been shown to lead to severe artifacts in -simulations. Ewald summation, in which the slowly conditionally -convergent 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 the -advantages of the fast Fourier transform (FFT) for calculating -discrete Fourier transforms, the particle mesh-based +shift 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 +nature. However, simply truncating the electrostatic interaction +with the use of cutoffs has been shown to lead to severe artifacts +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 the advantages +of the fast Fourier transform (FFT) 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 \emph{fast -multipole method}\cite{Greengard1987, Greengard1994}, which treats -Coulombic interaction exactly at short range, and approximate the -potential at long range through multipolar expansion. In spite of -their wide acceptances at the molecular simulation community, these -two methods are hard to be implemented 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 particle $i$ and particle $j$ at -distance $r_{rj}$ is given by: +$O(N^{3/2})$ to $O(N logN)$. An alternative approach is the +\emph{fast multipole method}\cite{Greengard1987, Greengard1994}, +which treats Coulombic interactions exactly at short range, and +approximate the potential at long range through multipolar +expansion. In spite of their wide acceptance at the molecular +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 +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 r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow @@ -1055,19 +1060,18 @@ Recently, advanced visualization technique are widely \subsection{\label{introSection:Analysis} Analysis} -Recently, advanced visualization technique are widely applied to +Recently, advanced visualization technique have become applied to monitor the motions of molecules. Although the dynamics of the system can be described qualitatively from animation, quantitative -trajectory analysis are more appreciable. According to the -principles of Statistical Mechanics, -Sec.~\ref{introSection:statisticalMechanics}, one can compute -thermodynamics properties, analyze fluctuations of structural -parameters, and investigate time-dependent processes of the molecule -from the trajectories. +trajectory analysis are more useful. According to the principles of +Statistical Mechanics, Sec.~\ref{introSection:statisticalMechanics}, +one can compute thermodynamic properties, analyze fluctuations of +structural parameters, and investigate time-dependent processes of +the molecule from the trajectories. -\subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamics Properties}} +\subsubsection{\label{introSection:thermodynamicsProperties}\textbf{Thermodynamic Properties}} -Thermodynamics properties, which can be expressed in terms of some +Thermodynamic properties, which can be expressed in terms of some function of the coordinates and momenta of all particles in the system, can be directly computed from molecular dynamics. The usual way to measure the pressure is based on virial theorem of Clausius @@ -1090,23 +1094,24 @@ distribution functions. Among these functions,\emph{pa \subsubsection{\label{introSection:structuralProperties}\textbf{Structural Properties}} Structural Properties of a simple fluid can be described by a set of -distribution functions. Among these functions,\emph{pair +distribution functions. Among these functions,the \emph{pair distribution function}, also known as \emph{radial distribution -function}, is of most fundamental importance to liquid-state theory. -Pair distribution function can be gathered by Fourier transforming -raw data from a series of neutron diffraction experiments and -integrating over the surface factor \cite{Powles1973}. The -experiment result can serve as a criterion to justify the -correctness of the theory. Moreover, various equilibrium -thermodynamic and structural properties can also be expressed in -terms of radial distribution function \cite{Allen1987}. +function}, is of most fundamental importance to liquid theory. +Experimentally, pair distribution function 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 radial distribution function \cite{Allen1987}. -A pair distribution functions $g(r)$ gives the probability that a +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 \[ g(r) = \frac{V}{{N^2 }}\left\langle {\sum\limits_i {\sum\limits_{j -\ne i} {\delta (r - r_{ij} )} } } \right\rangle. +\ne i} {\delta (r - r_{ij} )} } } \right\rangle = \fract{\rho +(r)}{\rho}. \] Note that the delta function can be replaced by a histogram in computer simulation. Figure @@ -1116,7 +1121,7 @@ locating particles at distances less than about 2.5{\A to find particles at certain radial values than at others. This is a result of the attractive interaction at such distances. Because of the strong repulsive forces at short distance, the probability of -locating particles at distances less than about 2.5{\AA} from each +locating particles at distances less than about 3.7{\AA} from each other is essentially zero. %\begin{figure} @@ -1131,25 +1136,25 @@ correlation function}, which correlates random variabl Properties}} Time-dependent properties are usually calculated using \emph{time -correlation function}, which correlates random variables $A$ and $B$ -at two different time +correlation functions}, which correlate random variables $A$ and $B$ +at two different times, \begin{equation} C_{AB} (t) = \left\langle {A(t)B(0)} \right\rangle. \label{introEquation:timeCorrelationFunction} \end{equation} If $A$ and $B$ refer to same variable, this kind of correlation -function is called \emph{auto correlation function}. One example of -auto correlation function is velocity auto-correlation function -which is directly related to transport properties of molecular -liquids: +function is called an \emph{autocorrelation function}. One example +of an auto correlation function is the velocity auto-correlation +function which is directly related to transport properties of +molecular liquids: \[ D = \frac{1}{3}\int\limits_0^\infty {\left\langle {v(t) \cdot v(0)} \right\rangle } dt \] -where $D$ is diffusion constant. Unlike velocity autocorrelation -function which is averaging over time origins and over all the -atoms, dipole autocorrelation are calculated for the entire system. -The dipole autocorrelation function is given by: +where $D$ is diffusion constant. Unlike the velocity autocorrelation +function, which is averaging over time origins and over all the +atoms, the dipole autocorrelation functions are calculated for the +entire system. The dipole autocorrelation function is given by: \[ c_{dipole} = \left\langle {u_{tot} (t) \cdot u_{tot} (t)} \right\rangle @@ -1175,46 +1180,46 @@ simulator is governed by the rigid body dynamics. In m areas, from engineering, physics, to chemistry. For example, missiles and vehicle are usually modeled by rigid bodies. The movement of the objects in 3D gaming engine or other physics -simulator is governed by the rigid body dynamics. In molecular -simulation, rigid body is used to simplify the model in -protein-protein docking study\cite{Gray2003}. +simulator is governed by rigid body dynamics. In molecular +simulations, rigid bodies are used to simplify protein-protein +docking studies\cite{Gray2003}. It is very important to develop stable and efficient methods to -integrate the equations of motion of orientational degrees of -freedom. Euler angles are the nature choice to describe the -rotational degrees of freedom. However, due to its singularity, the -numerical integration of corresponding equations of motion is very -inefficient and inaccurate. Although an alternative integrator using -different sets of Euler angles can overcome this -difficulty\cite{Barojas1973}, the computational penalty and the lost -of angular momentum conservation still remain. A singularity free -representation utilizing quaternions was developed by Evans in -1977\cite{Evans1977}. Unfortunately, this approach suffer from the -nonseparable Hamiltonian resulted from quaternion representation, -which prevents the symplectic algorithm to be utilized. Another -different approach is to apply holonomic constraints to the atoms -belonging to the rigid body. Each atom moves independently under the -normal forces deriving from potential energy and constraint forces -which are used to guarantee the rigidness. However, due to their -iterative nature, SHAKE and Rattle algorithm converge very slowly -when the number of constraint increases\cite{Ryckaert1977, -Andersen1983}. +integrate the equations of motion for orientational degrees of +freedom. Euler angles are the natural choice to describe the +rotational degrees of freedom. However, due to $\frac {1}{sin +\theta}$ singularities, the numerical integration of corresponding +equations of motion is very inefficient and inaccurate. Although an +alternative integrator using multiple sets of Euler 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}. Unfortunately, this +approach uses a nonseparable Hamiltonian resulting from the +quaternion representation, which prevents the symplectic algorithm +to be utilized. Another different approach is to apply holonomic +constraints to the atoms belonging to the rigid body. Each atom +moves independently under the normal forces 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}. -The break through in geometric literature suggests that, in order to +A break-through in geometric literature suggests that, in order to develop a long-term integration scheme, one should preserve the -symplectic structure of the flow. Introducing conjugate momentum to -rotation matrix $Q$ and re-formulating Hamiltonian's equation, a -symplectic integrator, RSHAKE\cite{Kol1997}, was proposed to evolve -the Hamiltonian system in a constraint manifold by iteratively -satisfying the orthogonality constraint $Q_T Q = 1$. An alternative -method using quaternion representation was developed by -Omelyan\cite{Omelyan1998}. However, both of these methods are -iterative and inefficient. In this section, we will present a +symplectic structure of the flow. By introducing a conjugate +momentum to the rotation matrix $Q$ and re-formulating Hamiltonian's +equation, a symplectic integrator, RSHAKE\cite{Kol1997}, was +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 +methods are iterative and inefficient. In this section, we descibe a symplectic Lie-Poisson integrator for rigid body developed by Dullweber and his coworkers\cite{Dullweber1997} in depth. -\subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Body} -The motion of the rigid body is Hamiltonian with the Hamiltonian +\subsection{\label{introSection:constrainedHamiltonianRB}Constrained Hamiltonian for Rigid Bodies} +The motion of a rigid body is Hamiltonian with the Hamiltonian function \begin{equation} H = \frac{1}{2}(p^T m^{ - 1} p) + \frac{1}{2}tr(PJ^{ - 1} P) + @@ -1228,13 +1233,14 @@ constrained Hamiltonian equation subjects to a holonom I_{ii}^{ - 1} = \frac{1}{2}\sum\limits_{i \ne j} {J_{jj}^{ - 1} } \] where $I_{ii}$ is the diagonal element of the inertia tensor. This -constrained Hamiltonian equation subjects to a holonomic constraint, +constrained Hamiltonian equation is subjected to a holonomic +constraint, \begin{equation} Q^T Q = 1, \label{introEquation:orthogonalConstraint} \end{equation} -which is used to ensure rotation matrix's orthogonality. -Differentiating \ref{introEquation:orthogonalConstraint} and using -Equation \ref{introEquation:RBMotionMomentum}, one may obtain, +which is used to ensure rotation matrix's unitarity. Differentiating +\ref{introEquation:orthogonalConstraint} and using Equation +\ref{introEquation:RBMotionMomentum}, one may obtain, \begin{equation} Q^T PJ^{ - 1} + J^{ - 1} P^T Q = 0 . \\ \label{introEquation:RBFirstOrderConstraint} @@ -1252,11 +1258,12 @@ We can use constraint force provided by lagrange multi \end{eqnarray} In general, there are two ways to satisfy the holonomic constraints. -We can use constraint force provided by lagrange multiplier on the -normal manifold to keep the motion on constraint space. Or we can -simply evolve the system in constraint manifold. These two methods -are proved to be equivalent. The holonomic constraint and equations -of motions define a constraint manifold for rigid body +We can use a constraint force provided by a Lagrange multiplier on +the normal manifold to keep the motion on constraint space. Or we +can simply evolve the system on the constraint manifold. These two +methods have been proved to be equivalent. The holonomic constraint +and equations of motions define a constraint manifold for rigid +bodies \[ M = \left\{ {(Q,P):Q^T Q = 1,Q^T PJ^{ - 1} + J^{ - 1} P^T Q = 0} \right\}. @@ -1265,7 +1272,7 @@ diffeomorphic. Introducing Unfortunately, this constraint manifold is not the cotangent bundle $T_{\star}SO(3)$. However, it turns out that under symplectic transformation, the cotangent space and the phase space are -diffeomorphic. Introducing +diffeomorphic. By introducing \[ \tilde Q = Q,\tilde P = \frac{1}{2}\left( {P - QP^T Q} \right), \] @@ -1297,8 +1304,8 @@ body, angular momentum on body frame $\Pi = Q^t P$ is respectively. As a common choice to describe the rotation dynamics of the rigid -body, angular momentum on body frame $\Pi = Q^t P$ is introduced to -rewrite the equations of motion, +body, the angular momentum on the body fixed frame $\Pi = Q^t P$ is +introduced to rewrite the equations of motion, \begin{equation} \begin{array}{l} \mathop \Pi \limits^ \bullet = J^{ - 1} \Pi ^T \Pi + Q^T \sum\limits_i {F_i (q,Q)X_i^T } - \Lambda \\ @@ -1341,11 +1348,11 @@ unique property eliminate the requirement of iteration Since $\Lambda$ is symmetric, the last term of Equation \ref{introEquation:skewMatrixPI} is zero, which implies the Lagrange multiplier $\Lambda$ is absent from the equations of motion. This -unique property eliminate the requirement of iterations which can +unique property eliminates the requirement of iterations which can not be avoided in other methods\cite{Kol1997, Omelyan1998}. -Applying hat-map isomorphism, we obtain the equation of motion for -angular momentum on body frame +Applying the hat-map isomorphism, we obtain the equation of motion +for angular momentum on body frame \begin{equation} \dot \pi = \pi \times I^{ - 1} \pi + \sum\limits_i {\left( {Q^T F_i (r,Q)} \right) \times X_i }. @@ -1360,10 +1367,11 @@ If there is not external forces exerted on the rigid b \subsection{\label{introSection:SymplecticFreeRB}Symplectic Lie-Poisson Integrator for Free Rigid Body} -If there is not external forces exerted on the rigid body, the only -contribution to the rotational is from the kinetic potential (the -first term of \ref{introEquation:bodyAngularMotion}). The free rigid -body is an example of Lie-Poisson system with Hamiltonian function +If there are no external forces exerted on the rigid body, the only +contribution to the rotational motion is from the kinetic energy +(the first term of \ref{introEquation:bodyAngularMotion}). The free +rigid body is an example of a Lie-Poisson system with Hamiltonian +function \begin{equation} T^r (\pi ) = T_1 ^r (\pi _1 ) + T_2^r (\pi _2 ) + T_3^r (\pi _3 ) \label{introEquation:rotationalKineticRB} @@ -1410,23 +1418,22 @@ tR_1 }$, we can use Cayley transformation, \end{array}} \right),\theta _1 = \frac{{\pi _1 }}{{I_1 }}\Delta t. \] To reduce the cost of computing expensive functions in $e^{\Delta -tR_1 }$, we can use Cayley transformation, +tR_1 }$, we can use Cayley transformation to obtain a single-aixs +propagator, \[ e^{\Delta tR_1 } \approx (1 - \Delta tR_1 )^{ - 1} (1 + \Delta tR_1 ) \] The flow maps for $T_2^r$ and $T_3^r$ can be found in the same -manner. - -In order to construct a second-order symplectic method, we split the -angular kinetic Hamiltonian function can into five terms +manner. In order to construct a second-order symplectic method, we +split the angular kinetic Hamiltonian function can into five terms \[ T^r (\pi ) = \frac{1}{2}T_1 ^r (\pi _1 ) + \frac{1}{2}T_2^r (\pi _2 ) + T_3^r (\pi _3 ) + \frac{1}{2}T_2^r (\pi _2 ) + \frac{1}{2}T_1 ^r -(\pi _1 ) -\]. -Concatenating flows corresponding to these five terms, we can obtain -an symplectic integrator, +(\pi _1 ). +\] +By concatenating the propagators corresponding to these five terms, +we can obtain an symplectic integrator, \[ \varphi _{\Delta t,T^r } = \varphi _{\Delta t/2,\pi _1 } \circ \varphi _{\Delta t/2,\pi _2 } \circ \varphi _{\Delta t,\pi _3 } @@ -1453,9 +1460,10 @@ Lie-Poisson integrator is found to be extremely effici \] Thus $ [\nabla _\pi F(\pi )]^T J(\pi ) = - S'(\frac{{\parallel \pi \parallel ^2 }}{2})\pi \times \pi = 0 $. This explicit -Lie-Poisson integrator is found to be extremely efficient and stable -which can be explained by the fact the small angle approximation is -used and the norm of the angular momentum is conserved. +Lie-Poisson integrator is found to be both extremely efficient and +stable. These properties can be explained by the fact the small +angle approximation is used and the norm of the angular momentum is +conserved. \subsection{\label{introSection:RBHamiltonianSplitting} Hamiltonian Splitting for Rigid Body} @@ -1482,36 +1490,36 @@ A second-order symplectic method is now obtained by th \end{tabular} \end{center} \end{table} -A second-order symplectic method is now obtained by the -composition of the flow maps, +A second-order symplectic method is now obtained by the composition +of the position and velocity propagators, \[ \varphi _{\Delta t} = \varphi _{\Delta t/2,V} \circ \varphi _{\Delta t,T} \circ \varphi _{\Delta t/2,V}. \] Moreover, $\varphi _{\Delta t/2,V}$ can be divided into two -sub-flows which corresponding to force and torque respectively, +sub-propagators which corresponding to force and torque +respectively, \[ \varphi _{\Delta t/2,V} = \varphi _{\Delta t/2,F} \circ \varphi _{\Delta t/2,\tau }. \] Since the associated operators of $\varphi _{\Delta t/2,F} $ and -$\circ \varphi _{\Delta t/2,\tau }$ are commuted, the composition -order inside $\varphi _{\Delta t/2,V}$ does not matter. - -Furthermore, kinetic potential can be separated to translational -kinetic term, $T^t (p)$, and rotational kinetic term, $T^r (\pi )$, +$\circ \varphi _{\Delta t/2,\tau }$ commute, the composition order +inside $\varphi _{\Delta t/2,V}$ does not matter. Furthermore, the +kinetic energy can be separated to translational kinetic term, $T^t +(p)$, and rotational kinetic term, $T^r (\pi )$, \begin{equation} T(p,\pi ) =T^t (p) + T^r (\pi ). \end{equation} where $ T^t (p) = \frac{1}{2}p^T m^{ - 1} p $ and $T^r (\pi )$ is defined by \ref{introEquation:rotationalKineticRB}. Therefore, the -corresponding flow maps are given by +corresponding propagators are given by \[ \varphi _{\Delta t,T} = \varphi _{\Delta t,T^t } \circ \varphi _{\Delta t,T^r }. \] -Finally, we obtain the overall symplectic flow maps for free moving -rigid body +Finally, we obtain the overall symplectic propagators for freely +moving rigid bodies \begin{equation} \begin{array}{c} \varphi _{\Delta t} = \varphi _{\Delta t/2,F} \circ \varphi _{\Delta t/2,\tau } \\ @@ -1525,28 +1533,28 @@ the theory of Langevin dynamics simulation. A brief de As an alternative to newtonian dynamics, Langevin dynamics, which mimics a simple heat bath with stochastic and dissipative forces, has been applied in a variety of studies. This section will review -the theory of Langevin dynamics simulation. A brief derivation of -generalized Langevin equation will be given first. Follow that, we -will discuss the physical meaning of the terms appearing in the -equation as well as the calculation of friction tensor from -hydrodynamics theory. +the theory of Langevin dynamics. A brief derivation of generalized +Langevin equation will be given first. Following that, we will +discuss the physical meaning of the terms appearing in the equation +as well as the calculation of friction tensor from hydrodynamics +theory. \subsection{\label{introSection:generalizedLangevinDynamics}Derivation of Generalized Langevin Equation} -Harmonic bath model, in which an effective set of harmonic +A harmonic bath model, in which an effective set of harmonic oscillators are used to mimic the effect of a linearly responding environment, has been widely used in quantum chemistry and statistical mechanics. One of the successful applications of -Harmonic bath model is the derivation of Deriving Generalized -Langevin Dynamics. Lets consider a system, in which the degree of +Harmonic bath model is the derivation of the Generalized Langevin +Dynamics (GLE). Lets consider a system, in which the degree of freedom $x$ is assumed to couple to the bath linearly, giving a Hamiltonian of the form \begin{equation} H = \frac{{p^2 }}{{2m}} + U(x) + H_B + \Delta U(x,x_1 , \ldots x_N) \label{introEquation:bathGLE}. \end{equation} -Here $p$ is a momentum conjugate to $q$, $m$ is the mass associated -with this degree of freedom, $H_B$ is harmonic bath Hamiltonian, +Here $p$ is a momentum conjugate to $x$, $m$ is the mass associated +with this degree of freedom, $H_B$ is a harmonic bath Hamiltonian, \[ H_B = \sum\limits_{\alpha = 1}^N {\left\{ {\frac{{p_\alpha ^2 }}{{2m_\alpha }} + \frac{1}{2}m_\alpha \omega _\alpha ^2 } @@ -1554,13 +1562,14 @@ the harmonic bath masses, and $\Delta U$ is bilinear s \] where the index $\alpha$ runs over all the bath degrees of freedom, $\omega _\alpha$ are the harmonic bath frequencies, $m_\alpha$ are -the harmonic bath masses, and $\Delta U$ is bilinear system-bath +the harmonic bath masses, and $\Delta U$ is a bilinear system-bath coupling, \[ \Delta U = - \sum\limits_{\alpha = 1}^N {g_\alpha x_\alpha x} \] -where $g_\alpha$ are the coupling constants between the bath and the -coordinate $x$. Introducing +where $g_\alpha$ are the coupling constants between the bath +coordinates ($x_ \apha$) and the system coordinate ($x$). +Introducing \[ W(x) = U(x) - \sum\limits_{\alpha = 1}^N {\frac{{g_\alpha ^2 }}{{2m_\alpha w_\alpha ^2 }}} x^2 @@ -1575,9 +1584,7 @@ Generalized Langevin Dynamics by Hamilton's equations \] Since the first two terms of the new Hamiltonian depend only on the system coordinates, we can get the equations of motion for -Generalized Langevin Dynamics by Hamilton's equations -\ref{introEquation:motionHamiltonianCoordinate, -introEquation:motionHamiltonianMomentum}, +Generalized Langevin Dynamics by Hamilton's equations, \begin{equation} m\ddot x = - \frac{{\partial W(x)}}{{\partial x}} - \sum\limits_{\alpha = 1}^N {g_\alpha \left( {x_\alpha - @@ -1594,11 +1601,11 @@ differential equations, Laplace transform is the appro In order to derive an equation for $x$, the dynamics of the bath variables $x_\alpha$ must be solved exactly first. As an integral transform which is particularly useful in solving linear ordinary -differential equations, Laplace transform is the appropriate tool to -solve this problem. The basic idea is to transform the difficult +differential equations,the Laplace transform is the appropriate tool +to solve this problem. The basic idea is to transform the difficult differential equations into simple algebra problems which can be -solved easily. Then applying inverse Laplace transform, also known -as the Bromwich integral, we can retrieve the solutions of the +solved easily. Then, by applying the inverse Laplace transform, also +known as the Bromwich integral, we can retrieve the solutions of the original problems. Let $f(t)$ be a function defined on $ [0,\infty ) $. The Laplace @@ -1618,7 +1625,7 @@ Applying Laplace transform to the bath coordinates, we \end{eqnarray*} -Applying Laplace transform to the bath coordinates, we obtain +Applying the Laplace transform to the bath coordinates, we obtain \begin{eqnarray*} 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) \\ L(x_\alpha ) & = & \frac{{\frac{{g_\alpha }}{{\omega _\alpha }}L(x) + px_\alpha (0) + \dot x_\alpha (0)}}{{p^2 + \omega _\alpha ^2 }} \\ @@ -1695,9 +1702,8 @@ as the model, which is gaussian distribution in genera \end{array} \] This property is what we expect from a truly random process. As long -as the model, which is gaussian distribution in general, chosen for -$R(t)$ is a truly random process, the stochastic nature of the GLE -still remains. +as the model chosen for $R(t)$ was a gaussian distribution in +general, the stochastic nature of the GLE still remains. %dynamic friction kernel The convolution integral @@ -1718,10 +1724,10 @@ which can be used to describe dynamic caging effect. T m\ddot x = - \frac{\partial }{{\partial x}}\left( {W(x) + \frac{1}{2}\xi _0 (x - x_0 )^2 } \right) + R(t), \] -which can be used to describe dynamic caging effect. The other -extreme is the bath that responds infinitely quickly to motions in -the system. Thus, $\xi (t)$ can be taken as a $delta$ function in -time: +which can be used to describe the effect of dynamic caging in +viscous solvents. The other extreme is the bath that responds +infinitely quickly to motions in the system. Thus, $\xi (t)$ can be +taken as a $delta$ function in time: \[ \xi (t) = 2\xi _0 \delta (t) \]