--- trunk/chrisDissertation/Introduction.tex	2006/09/07 20:42:27	3001
+++ 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,26 +93,164 @@ introduction to molecular dynamics and not Monte Carlo
 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,N).
+\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 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.
+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
+inventions.
 
 Moving a rigid body involves determination of both the force and
 torque applied by the surroundings, which directly affect the
@@ -74,7 +270,7 @@ position of, and torque on the component particles.
 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,
@@ -99,7 +295,7 @@ four quaternions ($q_w, q_x, q_y,$ and $q_z$),
 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) \\
@@ -110,25 +306,26 @@ small systems.\cite{Evans77} This integration methods 
 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
@@ -141,7 +338,7 @@ velocity-Verlet style two-part algorithm.\cite{Swope82
 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,
@@ -203,7 +400,7 @@ $\mathcal{O}\left(\delta t^4\right)$ for time steps of
 {\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
@@ -290,20 +487,8 @@ time. This demonstrates the computational advantage of
 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}.
-
-\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}
+time. This demonstrates the computational advantage of the {\sc dlm}
+integration scheme.
 
 \section{Accumulating Interactions}
 
@@ -312,13 +497,14 @@ interactions increases by $\mathcal{O}(N(N-1)/2)$ if y
 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
@@ -331,23 +517,24 @@ explicitly consider interactions between local particl
 \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
@@ -366,15 +553,16 @@ As particle move in and out of $R_\textrm{c}$, there w
 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}, \\
@@ -383,15 +571,16 @@ V_\textrm{LJ} = 4\epsilon\left[\left(\frac{\sigma}{r_{
 \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}) - 
@@ -401,11 +590,12 @@ of the derivative term distorts the potential even fur
 \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}, \\
@@ -423,89 +613,156 @@ 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}.
 
-\section{Measuring Properties}
+\subsection{\label{sec:periodicBoundaries}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.