[ViewVC] Diff of: group/trunk/chrisDissertation/Introduction.tex

Comparing trunk/chrisDissertation/Introduction.tex (file contents):
Revision 3001 by chrisfen, Thu Sep 7 20:42:27 2006 UTC vs.
Revision 3023 by chrisfen, Tue Sep 26 03:07:59 2006 UTC

-<
+\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}
->
+\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}
-+
+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
-+
+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
-+
+outright in the later work; however, I feel that this organization is
-+
+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.
-+
-+
+In chapter \ref{chap:water}, we focus on simple water models,
-+
+specifically the single-point soft sticky dipole (SSD) family of water
-+
+models. These single-point models are more efficient than the common
-+
+multi-point partial charge models and, in many cases, better capture
-+
+the dynamic properties of water. We discuss improvements to these
-+
+models in regards to long-range electrostatic corrections and show
-+
+that these models can work well with the techniques discussed in
-+
+chapter \ref{chap:electrostatics}. By investigating and improving
-+
+simple water models, we are extending the limits of computational
-+
+efficiency for systems that depend heavily on water calculations.
-+
-+
+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
-+
+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
-+
+conditions. We also consider electrostatic corrections, again
-+
+including the techniques discussed in chapter
-+
+\ref{chap:electrostatics}, to see how they influence the free energy
-+
+results. This work, to some degree, addresses the appropriateness of
-+
+using these simplistic water models outside of the conditions for
-+
+which they were developed.
-+
-+
+Finally, in chapter \ref{chap:conclusion}, we summarize the work
-+
+presented in the previous chapters and connect ideas together with
-+
+some final comments. The supporting information follows in the
-+
+appendix, and it gives a more detailed look at systems discussed in
-+
+chapter \ref{chap:electrostatics}.
-+
-+
+\section{On Molecular Simulation}
-+
+In order to advance our understanding of natural chemical and physical
+processes, researchers develop explanations for events observed
+experimentally. These hypotheses, supported by a body of corroborating
+molecular dynamics can be used to investigate dynamical
+quantities. The research presented in the following chapters utilized
+molecular dynamics near exclusively, so we will present a general
-<
+introduction to molecular dynamics and not Monte Carlo. There are
-<
+several resources available for those desiring a more in-depth
-<
+presentation of either of these
-<
+techniques.\cite{Allen87,Frenkel02,Leach01}
->
+introduction to molecular dynamics. There are several resources
->
+available for those desiring a more in-depth presentation of either of
->
+these techniques.\cite{Allen87,Frenkel02,Leach01}
+\section{\label{sec:MolecularDynamics}Molecular Dynamics}
-<
+\section{\label{sec:MovingParticles}Propagating Particle Motion}
->
+As stated above, in molecular dynamics we focus on evolving
->
+configurations of molecules over time. In order to use this as a tool
->
+for understanding experiments and testing theories, we want the
->
+configuration to evolve in a manner that mimics real molecular
->
+systems. To do this, we start with clarifying what we know about a
->
+given configuration of particles at time $t_1$, basically that each
->
+particle in the configuration has a position ($\mathbf{q}$) and velocity
->
+($\dot{\mathbf{q}}$). We now want to know what the configuration will be at
->
+time $t_2$. To find out, we need the classical equations of
->
+motion, and one useful formulation of them is the Lagrangian form.
-+
+The Lagrangian ($L$) is a function of the positions and velocities that
-+
+takes the form,
-+
+\begin{equation}
-+
+L = K - V,
-+
+\label{eq:lagrangian}
-+
+\end{equation}
-+
+where $K$ is the kinetic energy and $V$ is the potential energy. We
-+
+can use Hamilton's principle, which states that the integral of the
-+
+Lagrangian over time has a stationary value for the correct path of
-+
+motion, to say that the variation of the integral of the Lagrangian
-+
+over time is zero,\cite{Tolman38}
-+
+\begin{equation}
-+
+\delta\int_{t_1}^{t_2}L(\mathbf{q},\dot{\mathbf{q}})dt = 0.
-+
+\end{equation}
-+
+The variation can be transferred to the variables that make up the
-+
+Lagrangian,
-+
+\begin{equation}
-+
+\int_{t_1}^{t_2}\sum_{i=1}^{3N}\left(
-+
+        \frac{\partial L}{\partial \mathbf{q}_i}\delta \mathbf{q}_i
-+
+        + \frac{\partial L}{\partial \dot{\mathbf{q}}_i}\delta
-+
+                \dot{\mathbf{q}}_i\right)dt = 0.
-+
+\end{equation}
-+
+Using the fact that $\dot{\mathbf{q}}$ is the derivative of
-+
+$\mathbf{q}$ with respect to time and integrating the second partial
-+
+derivative in the parenthesis by parts, this equation simplifies to
-+
+\begin{equation}
-+
+\int_{t_1}^{t_2}\sum_{i=1}^{3N}\left(
-+
+        \frac{d}{dt}\frac{\partial L}{\partial \dot{\mathbf{q}}_i}
-+
+        - \frac{\partial L}{\partial \mathbf{q}_i}\right)
-+
+                \delta {\mathbf{q}}_i dt = 0,
-+
+\end{equation}
-+
+and since each variable is independent, we can separate the
-+
+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).
-+
+\label{eq:formulation}
-+
+\end{equation}
-+
+To obtain the classical equations of motion for the particles, we can
-+
+substitute equation (\ref{eq:lagrangian}) into the above equation with
-+
+$m\dot{\mathbf{r}}^2/2$ for the kinetic energy, giving
-+
+\begin{equation}
-+
+\frac{d}{dt}(m\dot{\mathbf{r}})+\frac{dV}{d\mathbf{r}}=0,
-+
+\end{equation}
-+
+or more recognizably,
-+
+\begin{equation}
-+
+\mathbf{f} = m\mathbf{a},
-+
+\end{equation}
-+
+where $\mathbf{f} = -dV/d\mathbf{r}$ and $\mathbf{a} =
-+
+d^2\mathbf{r}/dt^2$. This Lagrangian formulation shown in equation
-+
+(\ref{eq:formulation}) is generalized, and it can be used to determine
-+
+equations of motion in forms outside of the typical Cartesian case
-+
+shown here.
-+
-+
+\subsection{\label{sec:Verlet}Verlet Integration}
-+
-+
+In order to perform molecular dynamics, we need an algorithm that
-+
+integrates the equations of motion described above. Ideal algorithms
-+
+are both simple in implementation and highly accurate. There have been
-+
+a large number of algorithms developed for this purpose; however, for
-+
+reasons discussed below, we are going to focus on the Verlet class of
-+
+integrators.\cite{Gear66,Beeman76,Berendsen86,Allen87,Verlet67,Swope82}
-+
-+
+In Verlet's original study of computer ``experiments'', he directly
-+
+integrated the Newtonian second order differential equation of motion,
-+
+\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:
-+
+\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,
-+
+\label{eq:verlet}
-+
+\end{equation}
-+
+where $\delta t$ is the time step of integration.\cite{Verlet67} It is
-+
+interesting to note that equation (\ref{eq:verlet}) does not include
-+
+velocities, and this makes sense since they are not present in the
-+
+differential equation. The velocities are necessary for the
-+
+calculation of the kinetic energy and can be calculated at each time
-+
+step with the equation:
-+
+\begin{equation}
-+
+\mathbf{v}_i(t) = \frac{\mathbf{r}_i(t+\delta t)-\mathbf{r}_i(t-\delta t)}
-+
+                       {2\delta t}.
-+
+\end{equation}
-+
-+
+Like the equation of motion it solves, the Verlet algorithm has the
-+
+beneficial property of being time-reversible, meaning that you can
-+
+integrate the configuration forward and then backward and end up at
-+
+the original configuration. Some other methods for integration, like
-+
+predictor-corrector methods, lack this property in that they require
-+
+higher order information that is discarded after integrating
-+
+steps. Another interesting property of this algorithm is that it is
-+
+symplectic, meaning that it preserves area in phase-space. Symplectic
-+
+algorithms keep the system evolving in the region of phase-space
-+
+dictated by the ensemble and enjoy a greater degree of energy
-+
+conservation.\cite{Frenkel02}
-+
-+
+While the error in the positions calculated using the Verlet algorithm
-+
+is small ($\mathcal{O}(\delta t^4)$), the error in the velocities is
-+
+substantially larger ($\mathcal{O}(\delta t^2)$).\cite{Allen87} Swope
-+
+{\it et al.}  developed a corrected version of this algorithm, the
-+
+`velocity Verlet' algorithm, which improves the error of the velocity
-+
+calculation and thus the energy conservation.\cite{Swope82} This
-+
+algorithm involves a full step of the positions,
-+
+\begin{equation}
-+
+\mathbf{r}(t+\delta t) = \mathbf{r}(t) + \delta t\mathbf{v}(t)
-+
+                                + \frac{1}{2}\delta t^2\mathbf{a}(t),
-+
+\end{equation}
-+
+and a half step of the velocities,
-+
+\begin{equation}
-+
+\mathbf{v}\left(t+\frac{1}{2}\delta t\right) = \mathbf{v}(t)
-+
+                                        + \frac{1}{2}\delta t\mathbf{a}(t).
-+
+\end{equation}
-+
+After forces are calculated at the new positions, the velocities can
-+
+be updated to a full step,
-+
+\begin{equation}
-+
+\mathbf{v}(t+\delta t) = \mathbf{v}\left(t+\frac{1}{2}\delta t\right)
-+
+                                + \frac{1}{2}\delta t\mathbf{a}(t+\delta t).
-+
+\end{equation}
-+
+By integrating in this manner, the error in the velocities reduces to
-+
+$\mathcal{O}(\delta t^3)$. It should be noted that the error in the
-+
+positions increases to $\mathcal{O}(\delta t^3)$, but the resulting
-+
+improvement in the energies coupled with the maintained simplicity,
-+
+time-reversibility, and symplectic nature make it an improvement over
-+
+the original form.
-+
+\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}, they are not included in
-<
+most simulation packages because of the algorithmic complexity
-<
+involved in propagating 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 inventions.
->
+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
->
+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
->
+inventions.
+Moving a rigid body involves determination of both the force and
+torque applied by the surroundings, which directly affect the
+representation of the orientation of a rigid
+body.\cite{Evans77,Evans77b,Allen87} Thus, the elements of
+$\mathsf{A}$ can be expressed as arithmetic operations involving the
-<
+four quaternions ($q_w, q_x, q_y,$ and $q_z$),
->
+four quaternions ($q_0, q_1, q_2,$ and $q_3$),
+\begin{equation}
+\mathsf{A} = \left( \begin{array}{l@{\quad}l@{\quad}l}
+q_0^2+q_1^2-q_2^2-q_3^2 & 2(q_1q_2+q_0q_3) & 2(q_1q_3-q_0q_2) \\
+Integration of the equations of motion involves a series of arithmetic
+operations involving the quaternions and angular momenta and leads to
+performance enhancements over Euler angles, particularly for very
-<
+small systems.\cite{Evans77} This integration methods works well for
->
+small systems.\cite{Evans77} This integration method works well for
+propagating orientational motion in the canonical ensemble ($NVT$);
+however, energy conservation concerns arise when using the simple
+quaternion technique under the microcanonical ($NVE$) ensemble.  An
-<
+earlier implementation of {\sc oopse} utilized quaternions for
->
+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
+{\it symplectic}),
+\item the integrator is time-{\it reversible}, and
+\item the error for a single time step is of order
-<
+$\mathcal{O}\left(\delta t^4\right)$ for time steps of length $\delta t$.
->
+$\mathcal{O}\left(\delta t^3\right)$ for time steps of length $\delta t$.
+\end{enumerate}
+After the initial half-step ({\tt moveA}), the forces and torques are
+.001~kcal~mol$^{-1}$ per particle is desired, a nanosecond of
+simulation time will require ~19 hours of CPU time with the {\sc dlm}
+integrator, while the quaternion scheme will require ~154 hours of CPU
-<
+time. This demonstrates the computational advantage of the integration
-<
+scheme utilized in {\sc oopse}.
->
+time. This demonstrates the computational advantage of the {\sc dlm}
->
+integration scheme.
-–
+\subsection{Periodic Boundary Conditions}
-–
-–
+\begin{figure}
-–
+\centering
-–
+\includegraphics[width=4.5in]{./figures/periodicImage.pdf}
-–
+\caption{With periodic boundary conditions imposed, when particles
-–
+move out of one side the simulation box, they wrap back in the
-–
+opposite side. In this manner, a finite system of particles behaves as
-–
+an infinite system.}
-–
+\label{fig:periodicImage}
-–
+\end{figure}
-–
+\section{Accumulating Interactions}
+In the force calculation between {\tt moveA} and {\tt moveB} mentioned
+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 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.
->
+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.
+\begin{figure}
+\centering
+\end{figure}
+With spherical cutoffs, rather than accumulating the full set of
+interactions between all particles in the simulation, we only
-<
+explicitly consider interactions between local particles out to a
-<
+specified cutoff radius distance, $R_\textrm{c}$, (see figure
->
+explicitly consider interactions between particles separated by less
->
+than a specified cutoff radius distance, $R_\textrm{c}$, (see figure
+\ref{fig:sphereCut}). This reduces the scaling of the interaction to
+$\mathcal{O}(N\cdot\textrm{c})$, where `c' is a value that depends on
+the size of $R_\textrm{c}$ (c $\approx R_\textrm{c}^3$). Determination
+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 of the particles within $R_\textrm{l}$ move a distance of
->
+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}
->
+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}
+region in order to smooth out the discontinuity.}
+\label{fig:ljCutoff}
+\end{figure}
-<
+As particle move in and out of $R_\textrm{c}$, there will be sudden
-<
+discontinuous jumps in the potential (and forces) due to their
-<
+appearance and disappearance. 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 function so that it
-<
+crosses $R_\textrm{c}$ in a continuous fashion, and the easiest
-<
+methods is shifting the potential. To shift the potential, we simply
-<
+subtract out the value we calculate at $R_\textrm{c}$ from the whole
-<
+potential. For the shifted form of the Lennard-Jones potential is
->
+As the distance between a pair of particles fluctuates around
->
+$R_\textrm{c}$, there will be sudden discontinuous jumps in the
->
+potential (and forces) due to their inclusion and exclusion from the
->
+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
->
+the potential. To shift the potential, we simply subtract out the
->
+value we calculate 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}
+        V_\textrm{LJ}(r_{ij}) - V_\textrm{LJ}(R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\
+\end{equation}
+where
+\begin{equation}
-<
+V_\textrm{LJ} = 4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12} -
-<
+                                \left(\frac{\sigma}{r_{ij}}\right)^6\right].
->
+V_\textrm{LJ}(r_{ij}) =
->
+\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12} -
->
+        \left(\frac{\sigma}{r_{ij}}\right)^6\right].
+\end{equation}
+In figure \ref{fig:ljCutoff}, the shifted form of the potential
+reaches zero at the cutoff radius at the expense of the correct
-<
+magnitude below $R_\textrm{c}$. This correction method also does
->
+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 rather than just the potential via
->
+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}) -
+that the higher the order of the polynomial, the more abrupt the
+switching transition.
-<
+\subsection{Long-Range Interaction Corrections}
->
+\subsection{\label{sec:LJCorrections}Long-Range Interaction Corrections}
-<
+While a good approximation, accumulating interaction only from the
-<
+nearby particles can lead to some issues, because the further away
-<
+surrounding particles 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 interaction that extends out to extremely long
->
+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 neglect by
->
+homogeneous Lennard-Jones systems, we can correct for this effect by
+assuming a uniform density 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,
+\end{equation}
+where $V_\textrm{c}$ is the truncated Lennard-Jones
-<
+potential.\cite{Allen87} Like the potential, other Lennard-Jones
-<
+properties can be corrected by integration over the relevant
-<
+function. Note that with heterogeneous systems, this correction begins
-<
+to break down because the density is no longer uniform.
->
+potential.\cite{Allen87} Like the potential, other properties can be
->
+corrected by integration over the relevant function. Note that with
->
+heterogeneous systems, this correction breaks down because the density
->
+is no longer uniform.
+Correcting long-range electrostatic interactions is a topic of great
+importance in the field of molecular simulations. There have been
+this problem, as well as useful corrective techniques, is presented in
+chapter \ref{chap:electrostatics}.
-<
+\section{Measuring Properties}
->
+\subsection{Periodic Boundary Conditions}
-<
+\section{Application of the Techniques}
->
+In typical molecular dynamics simulations there are no restrictions
->
+placed on the motion of particles outside of what the inter-particle
->
+interactions dictate. This means that if a particle collides with
->
+other particles, it is free to move away from the site of the
->
+collision. If we consider the entire system as a collection of
->
+particles, they are not confined by walls of the ``simulation box''
->
+and can freely move away from the other particles. With no boundary
->
+considerations, particles moving outside of the simulation box
->
+enter a vacuum. This is correct behavior for cluster simulations in a
->
+vacuum; however, if we want to simulate bulk or spatially infinite
->
+systems, we need to use periodic boundary conditions.
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=4.5in]{./figures/periodicImage.pdf}
-+
+\caption{With periodic boundary conditions imposed, when particles
-+
+move out of one side the simulation box, they wrap back in the
-+
+opposite side. In this manner, a finite system of particles behaves as
-+
+an infinite system.}
-+
+\label{fig:periodicImage}
-+
+\end{figure}
-+
+In periodic boundary conditions, as a particle moves outside one wall
-+
+of the simulation box, the coordinates are remapped such that the
-+
+particle enters the opposing side of the box. This process is easy to
-+
+visualize in two dimensions as shown in figure \ref{fig:periodicImage}
-+
+and can occur in three dimensions, though it is not as easy to
-+
+visualize. Remapping the actual coordinates of the particles can be
-+
+problematic in that the we are restricting the distance a particle can
-+
+move from it's point of origin to a diagonal of the simulation
-+
+box. Thus, even though we are not confining the system with hard
-+
+walls, we are confining the particle coordinates to a particular
-+
+region in space. To avoid this ``soft'' confinement, it is common
-+
+practice to allow the particle coordinates to expand in an
-+
+unrestricted fashion while calculating interactions using a wrapped
-+
+set of ``minimum image'' coordinates. These coordinates need not be
-+
+stored because they are easily calculated while determining particle
-+
+distances.
-+
-+
+\section{Calculating Properties}
-+
-+
+In order to use simulations to model experimental processes and
-+
+evaluate theories, we need to be able to extract properties from the
-+
+results. In experiments, we can measure thermodynamic properties such
-+
+as the pressure and free energy. In computer simulations, we can
-+
+calculate properties from the motion and configuration of particles in
-+
+the system and make connections between these properties and the
-+
+experimental thermodynamic properties through statistical mechanics.
-+
-+
+The work presented in the later chapters use the canonical ($NVT$),
-+
+isobaric-isothermal ($NPT$), and microcanonical ($NVE$) statistical
-+
+mechanical ensembles. The different ensembles lend themselves to more
-+
+effectively calculating specific properties. For instance, if we
-+
+concerned ourselves with the calculation of dynamic properties, which
-+
+are dependent upon the motion of the particles, it is better to choose
-+
+an ensemble that does not add artificial motions to keep properties
-+
+like the temperature or pressure constant. In this case, the $NVE$
-+
+ensemble would be the most appropriate choice. In chapter
-+
+\ref{chap:ice}, we discuss calculating free energies, which are
-+
+non-mechanical thermodynamic properties, and these calculations also
-+
+depend on the chosen ensemble.\cite{Allen87} The Helmholtz free energy
-+
+($A$) depends on the $NVT$ partition function ($Q_{NVT}$),
-+
+\begin{equation}
-+
+A = -k_\textrm{B}T\ln Q_{NVT},
-+
+\end{equation}
-+
+while the Gibbs free energy ($G$) depends on the $NPT$ partition
-+
+function ($Q_{NPT}$),
-+
+\begin{equation}
-+
+G = -k_\textrm{B}T\ln Q_{NPT}.
-+
+\end{equation}
-+
+It is also useful to note that the conserved quantities of the $NVT$
-+
+and $NPT$ ensembles are related to the Helmholtz and Gibbs free
-+
+energies respectively.\cite{Melchionna93}
-+
-+
+Integrating the equations of motion is a simple method to obtain a
-+
+sequence of configurations that sample the chosen ensemble. For each
-+
+of these configurations, we can calculate an instantaneous value for a
-+
+chosen property like the density in the $NPT$ ensemble, where the
-+
+volume is allowed to fluctuate. The density for the simulation is
-+
+calculated from an average over the densities for the individual
-+
+configurations. Since the configurations from the integration process
-+
+are related to one another by the time evolution of the interactions,
-+
+this average is technically a time average. In calculating
-+
+thermodynamic properties, we would actually prefer an ensemble average
-+
+that is representative of all accessible states of the system. We can
-+
+calculate thermodynamic properties from the time average by taking
-+
+advantage of the ergodic hypothesis, which states that for a
-+
+sufficiently chaotic system, and over a long enough period of time,
-+
+the time and ensemble averages are the same.
-+
-+
+In addition to the average, the fluctuations of a particular property
-+
+can be determined via the standard deviation. Not only are
-+
+fluctuations useful for determining the spread of values around the
-+
+average and the error in the calculation of the value, they are also
-+
+useful for measuring various thermodynamic properties in computer
-+
+simulations. In section \ref{sec:t5peThermo}, we use fluctuations in
-+
+properties like the enthalpy and volume to calculate other
-+
+thermodynamic properties, such as the constant pressure heat capacity
-+
+and the isothermal compressibility.
-+
-+
+\section{OOPSE}
-+
+In the following chapters, the above techniques are all utilized in
+the study of molecular systems. There are a number of excellent
+simulation packages available, both free and commercial, which
+incorporate many of these
+methods.\cite{Brooks83,MacKerell98,Pearlman95,Berendsen95,Lindahl01,Smith96,Ponder87}
+Because of our interest in rigid body dynamics, point multipoles, and
-<
+systems where the orientational degrees cannot be handled using the
-<
+common constraint procedures (like {\sc shake}), the majority of the
-<
+following work was performed using {\sc oopse}, the object-oriented
-<
+parallel simulation engine.\cite{Meineke05} The {\sc oopse} package
-<
+started out as a collection of separate programs written within our
-<
+group, and has developed into one of the few parallel molecular
-<
+dynamics packages capable of accurately integrating rigid bodies and
-<
+point multipoles.
->
+systems where the orientational degrees of freedom cannot be handled
->
+using the common constraint procedures (like {\sc shake}), the
->
+majority of the following work was performed using {\sc oopse}, the
->
+object-oriented parallel simulation engine.\cite{Meineke05} The {\sc
->
+oopse} package started out as a collection of separate programs
->
+written within our group, and has developed into one of the few
->
+parallel molecular dynamics packages capable of accurately integrating
->
+rigid bodies and point multipoles. This simulation package is
->
+open-source software, available to anyone interested in performing
->
+molecular dynamics simulations. More information about {\sc oopse} can
->
+be found in reference \cite{Meineke05} or at the {\tt
->
+http://oopse.org} website.
-–
+In chapter \ref{chap:electrostatics}, we investigate correction
-–
+techniques for proper handling of long-ranged electrostatic
-–
+interactions. In particular we develop and evaluate some new pairwise
-–
+methods which we have incorporated into {\sc oopse}. 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
-–
+models. These single-point models are more efficient than the common
-–
+multi-point partial charge models and, in many cases, better capture
-–
+the dynamic properties of water. We discuss improvements to these
-–
+models in regards to long-range electrostatic corrections and show
-–
+that these models can work well with the techniques discussed in
-–
+chapter \ref{chap:electrostatics}. By investigating and improving
-–
+simple water models, we are extending the limits of computational
-–
+efficiency for systems that heavily depend on water calculations.
-–
-–
+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
-–
+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
-–
+conditions. We also consider electrostatic corrections, again
-–
+including the techniques discussed in chapter
-–
+\ref{chap:electrostatics}, to see how they influence the free energy
-–
+results. This work, to some degree, addresses the appropriateness of
-–
+using these simplistic water models outside of the conditions for
-–
+which they were developed.
-–
-–
+In the final chapter we summarize the work presented previously. We
-–
+also close with some final comments before the supporting information
-–
+presented in the appendices.

Diff Legend

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

Comparing trunk/chrisDissertation/Introduction.tex (file contents): Revision 3001 by chrisfen, Thu Sep 7 20:42:27 2006 UTC vs. Revision 3023 by chrisfen, Tue Sep 26 03:07:59 2006 UTC

Diff Legend

Comparing trunk/chrisDissertation/Introduction.tex (file contents):
Revision 3001 by chrisfen, Thu Sep 7 20:42:27 2006 UTC vs.
Revision 3023 by chrisfen, Tue Sep 26 03:07:59 2006 UTC