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

Comparing trunk/chrisDissertation/Introduction.tex (file contents):
Revision 3023 by chrisfen, Tue Sep 26 03:07:59 2006 UTC vs.
Revision 3025 by chrisfen, Tue Sep 26 16:02:25 2006 UTC

# Line 3 | Line 3 | years. Rather than presenting the topics in a chronolo
3   The following dissertation presents the primary aspects of the
4   research I have performed and been involved with over the last several
5   years. Rather than presenting the topics in a chronological fashion,
6 < they were arranged to form a series where the later topics apply and
6 > they are arranged to form a series where the later topics apply and
7   extend the findings of the former topics. This layout does lead to
8   occasional situations where knowledge gleaned from earlier chapters
9   (particularly chapter \ref{chap:electrostatics}) is not applied
# Line 11 | Line 11 | efforts.
11   more instructive and provides a more cohesive progression of research
12   efforts.
13  
14 < This chapter gives a general overview of molecular simulations, with
15 < particular emphasis on considerations that need to be made in order to
16 < apply the technique properly. This leads quite naturally into chapter
17 < \ref{chap:electrostatics}, where we investigate correction techniques
18 < for proper handling of long-ranged electrostatic interactions. In
19 < particular we develop and evaluate some new simple pairwise
20 < methods. These techniques make an appearance in the later chapters, as
21 < they are applied to specific systems of interest, showing how it they
22 < can improve the quality of various molecular simulations.
14 > This introductory chapter gives a general overview of molecular
15 > simulations, with particular emphasis on considerations that need to
16 > be made in order to apply the technique properly. This leads quite
17 > naturally into chapter \ref{chap:electrostatics}, where we investigate
18 > correction techniques for proper handling of long-ranged electrostatic
19 > interactions in molecular simulations. In particular we develop and
20 > evaluate some new simple pairwise methods. These techniques make an
21 > appearance in the later chapters, as they are applied to specific
22 > systems of interest, showing how it they can improve the quality of
23 > various molecular simulations.
24  
25   In chapter \ref{chap:water}, we focus on simple water models,
26   specifically the single-point soft sticky dipole (SSD) family of water
# Line 35 | Line 36 | water models discussed in the previous chapter. This f
36   In chapter \ref{chap:ice}, we study a unique polymorph of ice that we
37   discovered while performing water simulations with the fast simple
38   water models discussed in the previous chapter. This form of ice,
39 < which we called ``imaginary ice'' (Ice-$i$), has a low-density
40 < structure which is different from any known polymorph from either
41 < experiment or other simulations. In this study, we perform a free
39 > which we call ``imaginary ice'' (Ice-$i$), has a low-density structure
40 > which is different from any known polymorph characterized in either
41 > experiment or other simulations. In this work, we perform a free
42   energy analysis and see that this structure is in fact the
43   thermodynamically preferred form of ice for both the single-point and
44   commonly used multi-point water models under the chosen simulation
# Line 62 | Line 63 | limits of these theories, is essential in developing a
63   observations, can develop into accepted theories for how these
64   processes occur. This validation process, as well as testing the
65   limits of these theories, is essential in developing a firm foundation
66 < for our knowledge. Theories involving molecular scale systems cause
67 < particular difficulties in this process because their size and often
68 < fast motions make them difficult to observe experimentally. One useful
69 < tool for addressing these difficulties is computer simulation.
66 > for our knowledge. Developing and validating theories involving
67 > molecular scale systems is particularly difficult because their size
68 > and often fast motions make them difficult to observe
69 > experimentally. One useful tool for addressing these difficulties is
70 > computer simulation.
71  
72   In computer simulations, we can develop models from either the theory
73 < or experimental knowledge and then test them in a controlled
74 < environment. Work done in this manner allows us to further refine
73 > or our experimental knowledge and then test them in controlled
74 > environments. Work done in this manner allows us to further refine
75   theories, more accurately represent what is happening in experimental
76   observations, and even make predictions about what one will see in
77 < experiments. Thus, computer simulations of molecular systems act as a
78 < bridge between theory and experiment.
77 > future experiments. Thus, computer simulations of molecular systems
78 > act as a bridge between theory and experiment.
79  
80   Depending on the system of interest, there are a variety of different
81   computational techniques that can used to test and gather information
# Line 144 | Line 146 | contribution from each of the variables:
146   \begin{equation}
147   \frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i}
148          - \frac{\partial L}{\partial \mathbf{q}_i} = 0
149 <                \quad\quad(i=1,2,\dots,3N).
149 >                \quad\quad(i=1,2,\dots,N).
150   \label{eq:formulation}
151   \end{equation}
152   To obtain the classical equations of motion for the particles, we can
# Line 177 | Line 179 | m\frac{d^2\mathbf{r}_i}{dt^2} = \sum_{j\ne i}\mathbf{f
179   \begin{equation}
180   m\frac{d^2\mathbf{r}_i}{dt^2} = \sum_{j\ne i}\mathbf{f}(r_{ij}),
181   \end{equation}
182 < with the following simple algorithm:
182 > with the following algorithm:
183   \begin{equation}
184   \mathbf{r}_i(t+\delta t) = -\mathbf{r}_i(t-\delta t) + 2\mathbf{r}_i(t)
185          + \sum_{j\ne i}\mathbf{f}(r_{ij}(t))\delta t^2,
# Line 238 | Line 240 | Rigid bodies are non-spherical particles or collection
240   \subsection{\label{sec:IntroIntegrate}Rigid Body Motion}
241  
242   Rigid bodies are non-spherical particles or collections of particles
243 < (e.g. $\mbox{C}_{60}$) that have a constant internal potential and
244 < move collectively.\cite{Goldstein01} Discounting iterative constraint
245 < procedures like {\sc shake} and {\sc rattle} for approximating rigid
246 < bodies, they are not included in most simulation packages because of
247 < the algorithmic complexity involved in propagating orientational
243 > that have a constant internal potential and move
244 > collectively.\cite{Goldstein01} To move these particles, the
245 > translational and rotational motion can be integrated
246 > independently. Discounting iterative constraint procedures like {\sc
247 > shake} and {\sc rattle} for approximating rigid body dynamics, rigid
248 > bodies are not included in most simulation packages because of the
249 > algorithmic complexity involved in propagating the orientational
250   degrees of freedom.\cite{Ryckaert77,Andersen83,Krautler01} Integrators
251   which propagate orientational motion with an acceptable level of
252   energy conservation for molecular dynamics are relatively new
# Line 266 | Line 270 | $\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are th
270   where $\boldsymbol{\tau}_i$ and $\mathbf{r}_i$ are the torque on and
271   position of the center of mass respectively, while $\mathbf{f}_{ia}$,
272   $\mathbf{r}_{ia}$, and $\boldsymbol{\tau}_{ia}$ are the force on,
273 < position of, and torque on the component particles.
273 > position of, and torque on the component particles of the rigid body.
274  
275   The summation of the total torque is done in the body fixed axis. In
276   order to move between the space fixed and body fixed coordinate axes,
# Line 309 | Line 313 | showed that they resulted in a steady drift in the tot
313   earlier implementation of our simulation code utilized quaternions for
314   propagation of rotational motion; however, a detailed investigation
315   showed that they resulted in a steady drift in the total energy,
316 < something that has been observed by Kol {\it et al.} (also see
317 < section~\ref{sec:waterSimMethods}).\cite{Kol97}
316 > something that had also been observed by Kol {\it et al.} (see
317 > section~\ref{sec:waterSimMethods} for information on this
318 > investigation).\cite{Kol97}
319  
320 < Because of the outlined issues involving integration of the
321 < orientational motion using both Euler angles and quaternions, we
322 < decided to focus on a relatively new scheme that propagates the entire
323 < nine parameter rotation matrix. This techniques is a velocity-Verlet
324 < version of the symplectic splitting method proposed by Dullweber,
325 < Leimkuhler and McLachlan ({\sc dlm}).\cite{Dullweber97} When there are
326 < no directional atoms or rigid bodies present in the simulation, this
327 < integrator becomes the standard velocity-Verlet integrator which is
328 < known to effectively sample the microcanonical ($NVE$)
320 > Because of these issues involving integration of the orientational
321 > motion using both Euler angles and quaternions, we decided to focus on
322 > a relatively new scheme that propagates the entire nine parameter
323 > rotation matrix. This techniques is a velocity-Verlet version of the
324 > symplectic splitting method proposed by Dullweber, Leimkuhler and
325 > McLachlan ({\sc dlm}).\cite{Dullweber97} When there are no directional
326 > atoms or rigid bodies present in the simulation, this {\sc dlm}
327 > integrator reduces to the standard velocity-Verlet integrator, which
328 > is known to effectively sample the constant energy $NVE$
329   ensemble.\cite{Frenkel02}
330  
331   The key aspect of the integration method proposed by Dullweber
# Line 333 | Line 338 | The integration of the equations of motion is carried
338   substantial benefits in energy conservation.
339  
340   The integration of the equations of motion is carried out in a
341 < velocity-Verlet style two-part algorithm.\cite{Swope82} The first-part
341 > velocity-Verlet style two-part algorithm.\cite{Swope82} The first part
342   ({\tt moveA}) consists of a half-step ($t + \delta t/2$) of the linear
343   velocity (${\bf v}$) and angular momenta ({\bf j}) and a full-step ($t
344   + \delta t$) of the positions ({\bf r}) and rotation matrix,
# Line 492 | Line 497 | become a cumbersome task for large systems since the n
497   potential and forces (and torques if the particle is a rigid body or
498   multipole) on each particle from their surroundings. This can quickly
499   become a cumbersome task for large systems since the number of pair
500 < interactions increases by $\mathcal{O}(N(N-1)/2)$ if you accumulate
500 > interactions increases by $\mathcal{O}(N(N-1)/2)$ when accumulating
501   interactions between all particles in the system. (Note the
502   utilization of Newton's third law to reduce the interaction number
503 < from $\mathcal{O}(N^2)$.) The case of periodic boundary conditions
504 < further complicates matters by turning the finite system into an
505 < infinitely repeating one. Fortunately, we can reduce the scale of this
506 < problem by using spherical cutoff methods.
503 > from $\mathcal{O}(N^2)$.) Using periodic boundary conditions (section
504 > \ref{sec:periodicBoundaries}) further complicates matters by turning
505 > the finite system into an infinitely repeating one. Fortunately, we
506 > can reduce the scale of this problem by using spherical cutoff
507 > methods.
508  
509   \begin{figure}
510   \centering
# Line 519 | Line 525 | the this expense, we can use neighbor lists.\cite{Verl
525   of which particles are within the cutoff is also an issue, because
526   this process requires a full loop over all $N(N-1)/2$ pairs. To reduce
527   the this expense, we can use neighbor lists.\cite{Verlet67,Thompson83}
522 – With neighbor lists, we have a second list of particles built from a
523 – list radius $R_\textrm{l}$, which is larger than $R_\textrm{c}$. Once
524 – any particle within $R_\textrm{l}$ moves half the distance of
525 – $R_\textrm{l}-R_\textrm{c}$ (the ``skin'' thickness), we rebuild the
526 – list with the full $N(N-1)/2$ loop.\cite{Verlet67} With an appropriate
527 – skin thickness, these updates are only performed every $\sim$20 time
528 – steps, significantly reducing the time spent on pair-list bookkeeping
529 – operations.\cite{Allen87} If these neighbor lists are utilized, it is
530 – important that these list updates occur regularly. Incorrect
531 – application of this technique leads to non-physical dynamics, such as
532 – the ``flying block of ice'' behavior for which improper neighbor list
533 – handling was identified a one of the possible
534 – causes.\cite{Harvey98,Sagui99}
528  
529 + When using neighbor lists, we build a second list of particles built
530 + from a list radius $R_\textrm{l}$, which is larger than
531 + $R_\textrm{c}$. Once any particle within $R_\textrm{l}$ moves half the
532 + distance of $R_\textrm{l}-R_\textrm{c}$ (the ``skin'' thickness), we
533 + rebuild the list with the full $N(N-1)/2$ loop.\cite{Verlet67} With an
534 + appropriate skin thickness, these updates are only performed every
535 + $\sim$20 time steps, significantly reducing the time spent on
536 + pair-list bookkeeping operations.\cite{Allen87} If these neighbor
537 + lists are utilized, it is important that these list updates occur
538 + regularly. Incorrect application of this technique leads to
539 + non-physical dynamics, such as the ``flying block of ice'' behavior
540 + for which improper neighbor list handling was identified a one of the
541 + possible causes.\cite{Harvey98,Sagui99}
542 +
543   \subsection{Correcting Cutoff Discontinuities}
544   \begin{figure}
545   \centering
# Line 552 | Line 559 | smoothed out. There are several ways to modify the pot
559   interaction loop. In order to prevent heating and poor energy
560   conservation in the simulations, this discontinuity needs to be
561   smoothed out. There are several ways to modify the potential function
562 < to eliminate these discontinuties, and the easiest methods is shifting
562 > to eliminate these discontinuities, and the easiest method is shifting
563   the potential. To shift the potential, we simply subtract out the
564 < value we calculate at $R_\textrm{c}$ from the whole potential. The
564 > value of the function at $R_\textrm{c}$ from the whole potential. The
565   shifted form of the Lennard-Jones potential is
566   \begin{equation}
567   V_\textrm{SLJ} = \left\{\begin{array}{l@{\quad\quad}l}
# Line 572 | Line 579 | nothing to correct the discontinuity in the forces. We
579   reaches zero at the cutoff radius at the expense of the correct
580   magnitude inside $R_\textrm{c}$. This correction method also does
581   nothing to correct the discontinuity in the forces. We can account for
582 < this force discontinuity by shifting in the by applying the shifting
583 < in the forces as well as in the potential via
582 > this force discontinuity by applying the shifting in the forces as
583 > well as in the potential via
584   \begin{equation}
585   V_\textrm{SFLJ} = \left\{\begin{array}{l@{\quad\quad}l}
586          V_\textrm{LJ}({r_{ij}}) - V_\textrm{LJ}(R_\textrm{c}) -
# Line 583 | Line 590 | The forces are continuous with this potential; however
590   \end{array}\right.
591   \end{equation}
592   The forces are continuous with this potential; however, the inclusion
593 < of the derivative term distorts the potential even further than the
594 < simple shifting as shown in figure \ref{fig:ljCutoff}. The method for
595 < correcting the discontinuity which results in the smallest
596 < perturbation in both the potential and the forces is the use of a
597 < switching function. The cubic switching function,
593 > of the derivative term distorts the potential even further than
594 > shifting only the potential as shown in figure \ref{fig:ljCutoff}.
595 >
596 > The method for correcting these discontinuities which results in the
597 > smallest perturbation in both the potential and the forces is the use
598 > of a switching function. The cubic switching function,
599   \begin{equation}
600   S(r) = \left\{\begin{array}{l@{\quad\quad}l}
601          1 & r_{ij} \leqslant R_\textrm{sw}, \\
# Line 609 | Line 617 | particles can lead to some issues, because particles a
617  
618   While a good approximation, accumulating interactions only from nearby
619   particles can lead to some issues, because particles at distances
620 < greater than $R_\textrm{c}$ do still have a small effect. For
621 < instance, while the strength of the Lennard-Jones interaction is quite
622 < weak at $R_\textrm{c}$ in figure \ref{fig:ljCutoff}, we are discarding
623 < all of the attractive interactions that extend out to extremely long
624 < distances. Thus, the potential is a little too high and the pressure
625 < on the central particle from the surroundings is a little too low. For
626 < homogeneous Lennard-Jones systems, we can correct for this effect by
627 < assuming a uniform density and integrating the missing part,
620 > greater than $R_\textrm{c}$ still have a small effect. For instance,
621 > while the strength of the Lennard-Jones interaction is quite weak at
622 > $R_\textrm{c}$ in figure \ref{fig:ljCutoff}, we are discarding all of
623 > the attractive interactions that extend out to extremely long
624 > distances. Thus, the potential will be a little too high and the
625 > pressure on the central particle from the surroundings will be a
626 > little too low. For homogeneous Lennard-Jones systems, we can correct
627 > for this effect by assuming a uniform density ($\rho$) and integrating
628 > the missing part,
629   \begin{equation}
630   V_\textrm{full}(r_{ij}) \approx V_\textrm{c}
631                  + 2\pi N\rho\int^\infty_{R_\textrm{c}}r^2V_\textrm{LJ}(r)dr,
# Line 631 | Line 640 | summation and the reaction field technique. An in-dept
640   importance in the field of molecular simulations. There have been
641   several techniques developed to address this issue, like the Ewald
642   summation and the reaction field technique. An in-depth analysis of
643 < this problem, as well as useful corrective techniques, is presented in
643 > this problem, as well as useful correction techniques, is presented in
644   chapter \ref{chap:electrostatics}.
645  
646 < \subsection{Periodic Boundary Conditions}
646 > \subsection{\label{sec:periodicBoundaries}Periodic Boundary Conditions}
647  
648   In typical molecular dynamics simulations there are no restrictions
649   placed on the motion of particles outside of what the inter-particle

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines