--- trunk/chrisDissertation/Introduction.tex 2006/09/26 03:07:59 3023 +++ trunk/chrisDissertation/Introduction.tex 2006/09/26 16:02:25 3025 @@ -3,7 +3,7 @@ years. Rather than presenting the topics in a chronolo The following dissertation presents the primary aspects of the research I have performed and been involved with over the last several years. Rather than presenting the topics in a chronological fashion, -they were arranged to form a series where the later topics apply and +they are arranged to form a series where the later topics apply and extend the findings of the former topics. This layout does lead to occasional situations where knowledge gleaned from earlier chapters (particularly chapter \ref{chap:electrostatics}) is not applied @@ -11,15 +11,16 @@ efforts. more instructive and provides a more cohesive progression of research efforts. -This chapter gives a general overview of molecular simulations, with -particular emphasis on considerations that need to be made in order to -apply the technique properly. This leads quite naturally into chapter -\ref{chap:electrostatics}, where we investigate correction techniques -for proper handling of long-ranged electrostatic interactions. In -particular we develop and evaluate some new simple pairwise -methods. These techniques make an appearance in the later chapters, as -they are applied to specific systems of interest, showing how it they -can improve the quality of various molecular simulations. +This introductory chapter gives a general overview of molecular +simulations, with particular emphasis on considerations that need to +be made in order to apply the technique properly. This leads quite +naturally into chapter \ref{chap:electrostatics}, where we investigate +correction techniques for proper handling of long-ranged electrostatic +interactions in molecular simulations. In particular we develop and +evaluate some new simple pairwise methods. These techniques make an +appearance in the later chapters, as they are applied to specific +systems of interest, showing how it they can improve the quality of +various molecular simulations. In chapter \ref{chap:water}, we focus on simple water models, specifically the single-point soft sticky dipole (SSD) family of water @@ -35,9 +36,9 @@ water models discussed in the previous chapter. This f In chapter \ref{chap:ice}, we study a unique polymorph of ice that we discovered while performing water simulations with the fast simple water models discussed in the previous chapter. This form of ice, -which we called ``imaginary ice'' (Ice-$i$), has a low-density -structure which is different from any known polymorph from either -experiment or other simulations. In this study, we perform a free +which we call ``imaginary ice'' (Ice-$i$), has a low-density structure +which is different from any known polymorph characterized in either +experiment or other simulations. In this work, we perform a free energy analysis and see that this structure is in fact the thermodynamically preferred form of ice for both the single-point and commonly used multi-point water models under the chosen simulation @@ -62,18 +63,19 @@ limits of these theories, is essential in developing a observations, can develop into accepted theories for how these processes occur. This validation process, as well as testing the limits of these theories, is essential in developing a firm foundation -for our knowledge. Theories involving molecular scale systems cause -particular difficulties in this process because their size and often -fast motions make them difficult to observe experimentally. One useful -tool for addressing these difficulties is computer simulation. +for our knowledge. Developing and validating theories involving +molecular scale systems is particularly difficult because their size +and often fast motions make them difficult to observe +experimentally. One useful tool for addressing these difficulties is +computer simulation. In computer simulations, we can develop models from either the theory -or experimental knowledge and then test them in a controlled -environment. Work done in this manner allows us to further refine +or our experimental knowledge and then test them in controlled +environments. Work done in this manner allows us to further refine theories, more accurately represent what is happening in experimental observations, and even make predictions about what one will see in -experiments. Thus, computer simulations of molecular systems act as a -bridge between theory and experiment. +future experiments. Thus, computer simulations of molecular systems +act as a bridge between theory and experiment. Depending on the system of interest, there are a variety of different computational techniques that can used to test and gather information @@ -144,7 +146,7 @@ contribution from each of the variables: \begin{equation} \frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i} - \frac{\partial L}{\partial \mathbf{q}_i} = 0 - \quad\quad(i=1,2,\dots,3N). + \quad\quad(i=1,2,\dots,N). \label{eq:formulation} \end{equation} To obtain the classical equations of motion for the particles, we can @@ -177,7 +179,7 @@ m\frac{d^2\mathbf{r}_i}{dt^2} = \sum_{j\ne i}\mathbf{f \begin{equation} m\frac{d^2\mathbf{r}_i}{dt^2} = \sum_{j\ne i}\mathbf{f}(r_{ij}), \end{equation} -with the following simple algorithm: +with the following algorithm: \begin{equation} \mathbf{r}_i(t+\delta t) = -\mathbf{r}_i(t-\delta t) + 2\mathbf{r}_i(t) + \sum_{j\ne i}\mathbf{f}(r_{ij}(t))\delta t^2, @@ -238,11 +240,13 @@ Rigid bodies are non-spherical particles or collection \subsection{\label{sec:IntroIntegrate}Rigid Body Motion} Rigid bodies are non-spherical particles or collections of particles -(e.g. $\mbox{C}_{60}$) that have a constant internal potential and -move collectively.\cite{Goldstein01} Discounting iterative constraint -procedures like {\sc shake} and {\sc rattle} for approximating rigid -bodies, they are not included in most simulation packages because of -the algorithmic complexity involved in propagating orientational +that have a constant internal potential and move +collectively.\cite{Goldstein01} To move these particles, the +translational and rotational motion can be integrated +independently. Discounting iterative constraint procedures like {\sc +shake} and {\sc rattle} for approximating rigid body dynamics, rigid +bodies are not included in most simulation packages because of the +algorithmic complexity involved in propagating the orientational degrees of freedom.\cite{Ryckaert77,Andersen83,Krautler01} Integrators which propagate orientational motion with an acceptable level of energy conservation for molecular dynamics are relatively new @@ -266,7 +270,7 @@ $\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are th where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and position of the center of mass respectively, while $\mathbf{f}_{ia}$, $\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on, -position of, and torque on the component particles. +position of, and torque on the component particles of the rigid body. The summation of the total torque is done in the body fixed axis. In order to move between the space fixed and body fixed coordinate axes, @@ -309,18 +313,19 @@ showed that they resulted in a steady drift in the tot earlier implementation of our simulation code utilized quaternions for propagation of rotational motion; however, a detailed investigation showed that they resulted in a steady drift in the total energy, -something that has been observed by Kol {\it et al.} (also see -section~\ref{sec:waterSimMethods}).\cite{Kol97} +something that had also been observed by Kol {\it et al.} (see +section~\ref{sec:waterSimMethods} for information on this +investigation).\cite{Kol97} -Because of the outlined issues involving integration of the -orientational motion using both Euler angles and quaternions, we -decided to focus on a relatively new scheme that propagates the entire -nine parameter rotation matrix. This techniques is a velocity-Verlet -version of the symplectic splitting method proposed by Dullweber, -Leimkuhler and McLachlan ({\sc dlm}).\cite{Dullweber97} When there are -no directional atoms or rigid bodies present in the simulation, this -integrator becomes the standard velocity-Verlet integrator which is -known to effectively sample the microcanonical ($NVE$) +Because of these issues involving integration of the orientational +motion using both Euler angles and quaternions, we decided to focus on +a relatively new scheme that propagates the entire nine parameter +rotation matrix. This techniques is a velocity-Verlet version of the +symplectic splitting method proposed by Dullweber, Leimkuhler and +McLachlan ({\sc dlm}).\cite{Dullweber97} When there are no directional +atoms or rigid bodies present in the simulation, this {\sc dlm} +integrator reduces to the standard velocity-Verlet integrator, which +is known to effectively sample the constant energy $NVE$ ensemble.\cite{Frenkel02} The key aspect of the integration method proposed by Dullweber @@ -333,7 +338,7 @@ The integration of the equations of motion is carried substantial benefits in energy conservation. The integration of the equations of motion is carried out in a -velocity-Verlet style two-part algorithm.\cite{Swope82} The first-part +velocity-Verlet style two-part algorithm.\cite{Swope82} The first part ({\tt moveA}) consists of a half-step ($t + \delta t/2$) of the linear velocity (${\bf v}$) and angular momenta ({\bf j}) and a full-step ($t + \delta t$) of the positions ({\bf r}) and rotation matrix, @@ -492,13 +497,14 @@ become a cumbersome task for large systems since the n potential and forces (and torques if the particle is a rigid body or multipole) on each particle from their surroundings. This can quickly become a cumbersome task for large systems since the number of pair -interactions increases by $\mathcal{O}(N(N-1)/2)$ if you accumulate +interactions increases by $\mathcal{O}(N(N-1)/2)$ when accumulating interactions between all particles in the system. (Note the utilization of Newton's third law to reduce the interaction number -from $\mathcal{O}(N^2)$.) The case of periodic boundary conditions -further complicates matters by turning the finite system into an -infinitely repeating one. Fortunately, we can reduce the scale of this -problem by using spherical cutoff methods. +from $\mathcal{O}(N^2)$.) Using periodic boundary conditions (section +\ref{sec:periodicBoundaries}) further complicates matters by turning +the finite system into an infinitely repeating one. Fortunately, we +can reduce the scale of this problem by using spherical cutoff +methods. \begin{figure} \centering @@ -519,20 +525,21 @@ the this expense, we can use neighbor lists.\cite{Verl of which particles are within the cutoff is also an issue, because this process requires a full loop over all $N(N-1)/2$ pairs. To reduce the this expense, we can use neighbor lists.\cite{Verlet67,Thompson83} -With neighbor lists, we have a second list of particles built from a -list radius $R_\textrm{l}$, which is larger than $R_\textrm{c}$. Once -any particle within $R_\textrm{l}$ moves half the distance of -$R_\textrm{l}-R_\textrm{c}$ (the ``skin'' thickness), we rebuild the -list with the full $N(N-1)/2$ loop.\cite{Verlet67} With an appropriate -skin thickness, these updates are only performed every $\sim$20 time -steps, significantly reducing the time spent on pair-list bookkeeping -operations.\cite{Allen87} If these neighbor lists are utilized, it is -important that these list updates occur regularly. Incorrect -application of this technique leads to non-physical dynamics, such as -the ``flying block of ice'' behavior for which improper neighbor list -handling was identified a one of the possible -causes.\cite{Harvey98,Sagui99} +When using neighbor lists, we build a second list of particles built +from a list radius $R_\textrm{l}$, which is larger than +$R_\textrm{c}$. Once any particle within $R_\textrm{l}$ moves half the +distance of $R_\textrm{l}-R_\textrm{c}$ (the ``skin'' thickness), we +rebuild the list with the full $N(N-1)/2$ loop.\cite{Verlet67} With an +appropriate skin thickness, these updates are only performed every +$\sim$20 time steps, significantly reducing the time spent on +pair-list bookkeeping operations.\cite{Allen87} If these neighbor +lists are utilized, it is important that these list updates occur +regularly. Incorrect application of this technique leads to +non-physical dynamics, such as the ``flying block of ice'' behavior +for which improper neighbor list handling was identified a one of the +possible causes.\cite{Harvey98,Sagui99} + \subsection{Correcting Cutoff Discontinuities} \begin{figure} \centering @@ -552,9 +559,9 @@ smoothed out. There are several ways to modify the pot interaction loop. In order to prevent heating and poor energy conservation in the simulations, this discontinuity needs to be smoothed out. There are several ways to modify the potential function -to eliminate these discontinuties, and the easiest methods is shifting +to eliminate these discontinuities, and the easiest method is shifting the potential. To shift the potential, we simply subtract out the -value we calculate at $R_\textrm{c}$ from the whole potential. The +value of the function at $R_\textrm{c}$ from the whole potential. The shifted form of the Lennard-Jones potential is \begin{equation} V_\textrm{SLJ} = \left\{\begin{array}{l@{\quad\quad}l} @@ -572,8 +579,8 @@ nothing to correct the discontinuity in the forces. We reaches zero at the cutoff radius at the expense of the correct magnitude inside $R_\textrm{c}$. This correction method also does nothing to correct the discontinuity in the forces. We can account for -this force discontinuity by shifting in the by applying the shifting -in the forces as well as in the potential via +this force discontinuity by applying the shifting in the forces as +well as in the potential via \begin{equation} V_\textrm{SFLJ} = \left\{\begin{array}{l@{\quad\quad}l} V_\textrm{LJ}({r_{ij}}) - V_\textrm{LJ}(R_\textrm{c}) - @@ -583,11 +590,12 @@ The forces are continuous with this potential; however \end{array}\right. \end{equation} The forces are continuous with this potential; however, the inclusion -of the derivative term distorts the potential even further than the -simple shifting as shown in figure \ref{fig:ljCutoff}. The method for -correcting the discontinuity which results in the smallest -perturbation in both the potential and the forces is the use of a -switching function. The cubic switching function, +of the derivative term distorts the potential even further than +shifting only the potential as shown in figure \ref{fig:ljCutoff}. + +The method for correcting these discontinuities which results in the +smallest perturbation in both the potential and the forces is the use +of a switching function. The cubic switching function, \begin{equation} S(r) = \left\{\begin{array}{l@{\quad\quad}l} 1 & r_{ij} \leqslant R_\textrm{sw}, \\ @@ -609,14 +617,15 @@ particles can lead to some issues, because particles a While a good approximation, accumulating interactions only from nearby particles can lead to some issues, because particles at distances -greater than $R_\textrm{c}$ do still have a small effect. For -instance, while the strength of the Lennard-Jones interaction is quite -weak at $R_\textrm{c}$ in figure \ref{fig:ljCutoff}, we are discarding -all of the attractive interactions that extend out to extremely long -distances. Thus, the potential is a little too high and the pressure -on the central particle from the surroundings is a little too low. For -homogeneous Lennard-Jones systems, we can correct for this effect by -assuming a uniform density and integrating the missing part, +greater than $R_\textrm{c}$ still have a small effect. For instance, +while the strength of the Lennard-Jones interaction is quite weak at +$R_\textrm{c}$ in figure \ref{fig:ljCutoff}, we are discarding all of +the attractive interactions that extend out to extremely long +distances. Thus, the potential will be a little too high and the +pressure on the central particle from the surroundings will be a +little too low. For homogeneous Lennard-Jones systems, we can correct +for this effect by assuming a uniform density ($\rho$) and integrating +the missing part, \begin{equation} V_\textrm{full}(r_{ij}) \approx V_\textrm{c} + 2\pi N\rho\int^\infty_{R_\textrm{c}}r^2V_\textrm{LJ}(r)dr, @@ -631,10 +640,10 @@ summation and the reaction field technique. An in-dept importance in the field of molecular simulations. There have been several techniques developed to address this issue, like the Ewald summation and the reaction field technique. An in-depth analysis of -this problem, as well as useful corrective techniques, is presented in +this problem, as well as useful correction techniques, is presented in chapter \ref{chap:electrostatics}. -\subsection{Periodic Boundary Conditions} +\subsection{\label{sec:periodicBoundaries}Periodic Boundary Conditions} In typical molecular dynamics simulations there are no restrictions placed on the motion of particles outside of what the inter-particle