--- trunk/chrisDissertation/Introduction.tex	2006/09/04 16:36:51	3000
+++ trunk/chrisDissertation/Introduction.tex	2006/09/07 20:42:27	3001
@@ -1,3 +1,511 @@
 \chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}
 
-\section{\label{sec:IntroIntegrate}Rigid Body Motion}
+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.
+ 
+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
+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.
+
+Depending on the system of interest, there are a variety of different
+computational techniques that can used to test and gather information
+from the developed models. In the study of classical systems, the two
+most commonly used techniques are Monte Carlo and molecular
+dynamics. Both of these methods operate by calculating interactions
+between particles of our model systems; however, the progression of
+the simulation under the different techniques is vastly
+different. Monte Carlo operates through random configuration changes
+that that follow rules adhering to a specific statistical mechanics
+ensemble, while molecular dynamics is chiefly concerned with solving
+the classic equations of motion to move between configurations within
+an ensemble. Thermodynamic properties can be calculated from both
+techniques; but because of the random nature of Monte Carlo, only
+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}
+
+\section{\label{sec:MolecularDynamics}Molecular Dynamics}
+
+\section{\label{sec:MovingParticles}Propagating Particle Motion}
+
+\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.
+
+Moving a rigid body involves determination of both the force and
+torque applied by the surroundings, which directly affect the
+translational and rotational motion in turn. In order to accumulate
+the total force on a rigid body, the external forces and torques must
+first be calculated for all the internal particles. The total force on
+the rigid body is simply the sum of these external forces.
+Accumulation of the total torque on the rigid body is more complex
+than the force because the torque is applied to the center of mass of
+the rigid body. The space-fixed torque on rigid body $i$ is
+\begin{equation}
+\boldsymbol{\tau}_i=
+	\sum_{a}\biggl[(\mathbf{r}_{ia}-\mathbf{r}_i)\times \mathbf{f}_{ia} 
+	+ \boldsymbol{\tau}_{ia}\biggr],
+\label{eq:torqueAccumulate}
+\end{equation}
+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.
+
+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,
+parameters describing the orientation must be maintained for each
+rigid body. At a minimum, the rotation matrix ($\mathsf{A}$) can be
+described by the three Euler angles ($\phi, \theta,$ and $\psi$),
+where the elements of $\mathsf{A}$ are composed of trigonometric
+operations involving $\phi, \theta,$ and $\psi$.\cite{Goldstein01}
+Direct propagation of the Euler angles has a known $1/\sin\theta$
+divergence in the equations of motion for $\phi$ and $\psi$, leading
+to numerical instabilities any time one of the directional atoms or
+rigid bodies has an orientation near $\theta=0$ or
+$\theta=\pi$.\cite{Allen87} One of the most practical ways to avoid
+this ``gimbal point'' is to switch to another angular set defining the
+orientation of the rigid body near this point.\cite{Barojas73} This
+procedure results in additional book-keeping and increased algorithm
+complexity. In the search for more elegant alternative methods, Evans
+proposed the use of quaternions to describe and propagate
+orientational motion.\cite{Evans77}
+
+The quaternion method for integration involves a four dimensional
+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$),
+\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) \\
+2(q_1q_2-q_0q_3) & q_0^2-q_1^2+q_2^2-q_3^2 & 2(q_2q_3+q_0q_1) \\
+2(q_1q_3+q_0q_2) & 2(q_2q_3-q_0q_1) & q_0^2-q_1^2-q_2^2+q_3^2 \\
+\end{array}\right).
+\end{equation}
+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
+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
+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}
+
+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$)
+ensemble.\cite{Frenkel02}
+
+The key aspect of the integration method proposed by Dullweber
+\emph{et al.} is that the entire $3 \times 3$ rotation matrix is
+propagated from one time step to the next. In the past, this would not
+have been as feasible, since the rotation matrix for a single body has
+nine elements compared with the more memory-efficient methods (using
+three Euler angles or four quaternions).  Computer memory has become
+much less costly in recent years, and this can be translated into
+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
+({\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,
+\begin{equation*}
+{\tt moveA} = \left\{\begin{array}{r@{\quad\leftarrow\quad}l}
+{\bf v}\left(t + \delta t / 2\right) & {\bf v}(t) 
+	+ \left( {\bf f}(t) / m \right)(\delta t/2), \\
+%
+{\bf r}(t + \delta t) & {\bf r}(t) 
+	+ {\bf v}\left(t + \delta t / 2 \right)\delta t, \\
+%
+{\bf j}\left(t + \delta t / 2 \right) & {\bf j}(t) 
+	+ \boldsymbol{\tau}^b(t)(\delta t/2), \\
+%
+\mathsf{A}(t + \delta t) & \mathrm{rotate}\left( {\bf j}
+	(t + \delta t / 2)\delta t \cdot 
+		\overleftrightarrow{\mathsf{I}}^{-1} \right),
+\end{array}\right.
+\end{equation*}
+where $\overleftrightarrow{\mathsf{I}}^{-1}$ is the inverse of the
+moment of inertia tensor. The $\mathrm{rotate}$ function is the
+product of rotations about the three body-fixed axes,
+\begin{equation}
+\mathrm{rotate}({\bf a}) = \mathsf{G}_x(a_x / 2) \cdot
+\mathsf{G}_y(a_y / 2) \cdot \mathsf{G}_z(a_z) \cdot \mathsf{G}_y(a_y /
+2) \cdot \mathsf{G}_x(a_x /2),
+\label{eq:dlmTrot}
+\end{equation}
+where each rotational propagator, $\mathsf{G}_\alpha(\theta)$, rotates
+both the rotation matrix ($\mathsf{A}$) and the body-fixed angular
+momentum (${\bf j}$) by an angle $\theta$ around body-fixed axis
+$\alpha$,
+\begin{equation}
+\mathsf{G}_\alpha( \theta ) = \left\{
+\begin{array}{l@{\quad\leftarrow\quad}l}
+\mathsf{A}(t) & \mathsf{A}(0) \cdot \mathsf{R}_\alpha(\theta)^\textrm{T},\\
+{\bf j}(t) & \mathsf{R}_\alpha(\theta) \cdot {\bf j}(0).
+\end{array}
+\right.
+\end{equation}
+$\mathsf{R}_\alpha$ is a quadratic approximation to the single-axis
+rotation matrix.  For example, in the small-angle limit, the rotation
+matrix around the body-fixed x-axis can be approximated as
+\begin{equation}
+\mathsf{R}_x(\theta) \approx \left(
+\begin{array}{ccc}
+1 & 0 & 0 \\
+0 & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4} & -\frac{\theta}{1+\theta^2 / 4} \\
+0 & \frac{\theta}{1+\theta^2 / 4} & \frac{1-\theta^2 / 4}{1 + \theta^2 / 4}
+\end{array}
+\right).
+\end{equation}
+The remaining rotations follow in a straightforward manner. As seen
+from the form of equation~(\ref{eq:dlmTrot}), the {\sc dlm} method
+uses a Trotter factorization of the orientational
+propagator.\cite{Trotter59} This has three effects:
+\begin{enumerate}
+\item the integrator is area-preserving in phase space (i.e. it is
+{\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$.
+\end{enumerate}
+
+After the initial half-step ({\tt moveA}), the forces and torques are
+evaluated for all of the particles. Once completed, the velocities can
+be advanced to complete the second-half of the two-part algorithm
+({\tt moveB}), resulting an a completed full step of both the
+positions and momenta,
+\begin{equation*}
+{\tt moveB} = \left\{\begin{array}{r@{\quad\leftarrow\quad}l}
+{\bf v}\left(t + \delta t \right) &
+	{\bf v}\left(t + \delta t / 2 \right) 
+	+ \left({\bf f}(t + \delta t) / m \right)(\delta t/2), \\
+%
+{\bf j}\left(t + \delta t \right) &
+	{\bf j}\left(t + \delta t / 2 \right) 
+	+ \boldsymbol{\tau}^b(t + \delta t)(\delta t/2).
+\end{array}\right.
+\end{equation*}
+
+The matrix rotations used in the {\sc dlm} method end up being more
+costly computationally than the simpler arithmetic quaternion
+propagation. With the same time step, a 1024-molecule water simulation
+incurs approximately a 10\% increase in computation time using the
+{\sc dlm} method in place of quaternions. This cost is more than
+justified when comparing the energy conservation achieved by the two
+methods. Figure \ref{fig:quatdlm} provides a comparative analysis of
+the {\sc dlm} method versus the traditional quaternion scheme.
+
+\begin{figure}
+\centering
+\includegraphics[width=3.5in]{./figures/dlmVsQuat.pdf}
+\caption[Energy conservation analysis of the {\sc dlm} and quaternion 
+integration methods]{Analysis of the energy conservation of the {\sc
+dlm} and quaternion integration methods.  $\delta \mathrm{E}_1$ is the
+linear drift in energy over time and $\delta \mathrm{E}_0$ is the
+standard deviation of energy fluctuations around this drift.  All
+simulations were of a 1024 SSD water system at 298 K starting from the
+same initial configuration. Note that the {\sc dlm} method provides
+more than an order-of-magnitude improvement in both the energy drift
+and the size of the energy fluctuations when compared with the
+quaternion method at any given time step.  At time steps larger than 4
+fs, the quaternion scheme resulted in rapidly rising energies which
+eventually lead to simulation failure.  Using the {\sc dlm} method,
+time steps up to 8 fs can be taken before this behavior is evident.}
+\label{fig:quatdlm}
+\end{figure}
+
+In figure \ref{fig:quatdlm}, $\delta \mbox{E}_1$ is a measure of the
+linear energy drift in units of $\mbox{kcal mol}^{-1}$ per particle
+over a nanosecond of simulation time, and $\delta \mbox{E}_0$ is the
+standard deviation of the energy fluctuations in units of $\mbox{kcal
+mol}^{-1}$ per particle. In the top plot, it is apparent that the
+energy drift is reduced by a significant amount (2 to 3 orders of
+magnitude improvement at all tested time steps) by choosing the {\sc
+dlm} method over the simple non-symplectic quaternion integration
+method.  In addition to this improvement in energy drift, the
+fluctuations in the total energy are also dampened by 1 to 2 orders of
+magnitude by utilizing the {\sc dlm} method.
+
+\begin{figure}
+\centering
+\includegraphics[width=\linewidth]{./figures/compCost.pdf}
+\caption[Energy drift as a function of required simulation run 
+time]{Energy drift as a function of required simulation run time.
+$\delta \mathrm{E}_1$ is the linear drift in energy over time.
+Simulations were performed on a single 2.5 GHz Pentium 4
+processor. Simulation time comparisons can be made by tracing
+horizontally from one curve to the other. For example, a simulation
+that takes 24 hours using the {\sc dlm} method will take roughly
+210 hours using the simple quaternion method if the same degree of
+energy conservation is desired.}
+\label{fig:cpuCost}
+\end{figure}
+Although the {\sc dlm} method is more computationally expensive than
+the traditional quaternion scheme for propagating of a time step,
+consideration of the computational cost for a long simulation with a
+particular level of energy conservation is in order.  A plot of energy
+drift versus computational cost was generated
+(Fig.~\ref{fig:cpuCost}). This figure provides an estimate of the CPU
+time required under the two integration schemes for 1 nanosecond of
+simulation time for the model 1024-molecule system.  By choosing a
+desired energy drift value it is possible to determine the CPU time
+required for both methods. If a $\delta \mbox{E}_1$ 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}
+
+\section{Accumulating Interactions}
+
+In the force calculation between {\tt moveA} and {\tt moveB} mentioned
+in section \ref{sec:IntroIntegrate}, we need to accumulate the
+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.
+
+\begin{figure}
+\centering
+\includegraphics[width=3.5in]{./figures/sphericalCut.pdf}
+\caption{When using a spherical cutoff, only particles within a chosen
+cutoff radius distance, $R_\textrm{c}$, of the central particle are
+included in the pairwise summation. This reduces a problem that scales
+by $\sim\mathcal{O}(N^2)$ to one that scales by $\sim\mathcal{O}(N)$.}
+\label{fig:sphereCut}
+\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
+\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
+regularly. Incorrect application of this technique leads to
+non-physical dynamics, such as the ``flying block of ice'' behavior
+for which improper neighbor list handling was identified a one of the
+possible causes.\cite{Harvey98,Sagui99}
+
+\subsection{Correcting Cutoff Discontinuities}
+\begin{figure}
+\centering
+\includegraphics[width=3.5in]{./figures/ljCutoffSquare.pdf}
+\caption{The common methods to smooth the potential discontinuity 
+introduced when using a cutoff include a shifted potential, a shifted
+force, and a switching function. The shifted potential and shifted
+force both lift the whole potential so that it zeroes at
+$R_\textrm{c}$, thereby reducing the strength of the interaction. The
+(cubic) switching function only alters the potential in the switching
+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
+\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}, \\
+	0 & r_{ij} \geqslant R_\textrm{c},
+\end{array}\right.
+\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].
+\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
+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
+\begin{equation}
+V_\textrm{SFLJ} = \left\{\begin{array}{l@{\quad\quad}l}
+	V_\textrm{LJ}({r_{ij}}) - V_\textrm{LJ}(R_\textrm{c}) - 
+		\left(\frac{d V(r_{ij})}{d r_{ij}}\right)_{r_{ij}=R_\textrm{c}}
+		(r_{ij} - R_\textrm{c}) & r_{ij} < R_\textrm{c}, \\
+	0 & r_{ij} \geqslant R_\textrm{c}.
+\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,
+\begin{equation}
+S(r) = \left\{\begin{array}{l@{\quad\quad}l}
+	1 & r_{ij} \leqslant R_\textrm{sw}, \\
+	\frac{(R_\textrm{c} + 2r_{ij} - 3R_\textrm{sw})
+		(R_\textrm{c} - r_{ij})^2}{(R_\textrm{c} - R_\textrm{sw})^3} 
+		& R_\textrm{sw} < r_{ij} \leqslant R_\textrm{c}, \\
+	0 & r_{ij} > R_\textrm{c},
+	\end{array}\right.
+\end{equation}
+is sufficient to smooth the potential (again, see figure
+\ref{fig:ljCutoff}) and the forces by only perturbing the potential in
+the switching region. If smooth second derivatives are required, a
+higher order polynomial switching function (i.e. fifth order
+polynomial) can be used.\cite{Andrea83,Leach01} It should be noted
+that the higher the order of the polynomial, the more abrupt the
+switching transition.
+
+\subsection{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
+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,
+\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.
+
+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
+chapter \ref{chap:electrostatics}.
+
+\section{Measuring Properties}
+
+\section{Application of the Techniques}
+
+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.
+
+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.