--- trunk/chrisDissertation/Introduction.tex	2006/09/21 23:21:37	3016
+++ trunk/chrisDissertation/Introduction.tex	2006/09/26 16:02:25	3025
@@ -1,23 +1,81 @@
-\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 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
+outright in the later work; however, I feel that this organization is
+more instructive and provides a more cohesive progression of research
+efforts. 
+
+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
+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 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
+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
 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
@@ -35,10 +93,9 @@ molecular dynamics near exclusively, so we will presen
 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}
 
@@ -48,12 +105,12 @@ given configuration of particles at time $t_1$, basica
 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 ($q$) and velocity
-($\dot{q}$). We now want to know what the configuration will be at
+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 velocites that
+The Lagrangian ($L$) is a function of the positions and velocities that
 takes the form,
 \begin{equation}
 L = K - V,
@@ -65,25 +122,31 @@ over time is zero,\cite{Tolman38}
 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(q,\dot{q})dt = 0.
+\delta\int_{t_1}^{t_2}L(\mathbf{q},\dot{\mathbf{q}})dt = 0.
 \end{equation}
-This can be written as a summation of integrals to give
+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 q_i}\delta q_i
-	+ \frac{\partial L}{\partial \dot{q}_i}\delta \dot{q}_i\right)dt = 0.
+\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 fact that $\dot q$ is the derivative with respect to time of $q$
-and integrating the second partial derivative in the parenthesis by
-parts, this equation simplifies to
+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{q}_i}
-	- \frac{\partial L}{\partial q_i}\right)\delta {q}_i dt = 0,
+	\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 the above equation is only valid for
+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{q}_i}
-	- \frac{\partial L}{\partial q_i} = 0\quad\quad(i=1,2,\dots,3N).
+\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,N).
 \label{eq:formulation}
 \end{equation}
 To obtain the classical equations of motion for the particles, we can
@@ -116,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,
@@ -126,8 +189,8 @@ differential equation. The velocities are necessary fo
 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 they can be calculated at each
-time step with the following equation:
+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}.
@@ -141,14 +204,14 @@ symplectic, meaning that it preserves area in phase-sp
 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 system stays in the region of phase-space dictated by the
-ensemble and enjoy greater long-time energy
-conservations.\cite{Frenkel02}
+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
-quite a bit larger ($\mathcal{O}(\delta t^2)$).\cite{Allen87} Swope
-{\it et al.}  developed a corrected for of this algorithm, the
+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,
@@ -161,8 +224,8 @@ and a half step of the velocities,
 \mathbf{v}\left(t+\frac{1}{2}\delta t\right) = \mathbf{v}(t) 
 					+ \frac{1}{2}\delta t\mathbf{a}(t).
 \end{equation}
-After calculating new forces, the velocities can be updated to a full
-step,
+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).
@@ -177,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
@@ -205,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,
@@ -230,7 +295,7 @@ $\mathsf{A}$ can be expressed as arithmetic operations
 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) \\
@@ -241,25 +306,26 @@ performance enhancements over Euler angles, particular
 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
-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
@@ -272,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,
@@ -334,7 +400,7 @@ propagator.\cite{Trotter59} This has three effects:
 {\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
@@ -421,8 +487,8 @@ integrator, while the quaternion scheme will require ~
 0.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.
 
 \section{Accumulating Interactions}
 
@@ -431,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 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 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)$.) 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
@@ -450,23 +517,24 @@ interactions between all particles in the simulation, 
 \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
 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 of the particles within $R_\textrm{l}$ move a 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
+
+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
@@ -485,15 +553,16 @@ region in order to smooth out the discontinuity.}
 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 discontinuities, and the easiest method is shifting
+the potential. To shift the potential, we simply subtract out 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}
 	V_\textrm{LJ}(r_{ij}) - V_\textrm{LJ}(R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\
@@ -502,15 +571,16 @@ where
 \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}) =
+	4\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
+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}) - 
@@ -520,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}, \\
@@ -542,39 +613,40 @@ switching transition.
 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 a good approximation, accumulating interactions only from nearby
+particles can lead to some issues, because particles at distances
+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 interaction that extends 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
-assuming a uniform density and integrating the missing part,
+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,
 \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
 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 interparticle
+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
@@ -607,14 +679,14 @@ unrestricted fashion while calculating interactions us
 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 minimum image coordinates
-need not be stored because they are easily calculated on the fly when
-determining particle distances.
+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 process and evaluate
-theories, we need to be able to extract properties from the
+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
@@ -623,17 +695,17 @@ isobaric-isothermal ($NPT$), and microcanonical ($NVE$
 
 The work presented in the later chapters use the canonical ($NVT$),
 isobaric-isothermal ($NPT$), and microcanonical ($NVE$) statistical
-mechanics ensembles. The different ensembles lend themselves to more
+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 system 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}$),
+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}
@@ -652,25 +724,27 @@ calculated from an average over the densities for the 
 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. Being that 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 over a long period of time, the time average and the ensemble
-average are the same.
+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 that particular
-property can be determined via the standard deviation. Fluctuations
-are useful for measuring various thermodynamic properties in computer
+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
-propeties like the enthalpy and volume to calculate various
+properties like the enthalpy and volume to calculate other
 thermodynamic properties, such as the constant pressure heat capacity
 and the isothermal compressibility.
 
-\section{Application of the Techniques}
+\section{OOPSE}
 
 In the following chapters, the above techniques are all utilized in
 the study of molecular systems. There are a number of excellent
@@ -678,50 +752,17 @@ Because of our interest in rigid body dynamics, point 
 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.