--- trunk/chrisDissertation/Introduction.tex	2006/08/30 22:36:06	2987
+++ trunk/chrisDissertation/Introduction.tex	2006/09/26 16:02:25	3025
@@ -1,3 +1,768 @@
-\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}
+\chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND} 
 
-\section{\label{sec:IntroIntegrate}Rigid Body Motion}
+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. 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 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
+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
+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. 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}
+
+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
+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
+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 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,
+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_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) \\
+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 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 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 had also been observed by Kol {\it et al.} (see
+section~\ref{sec:waterSimMethods} for information on this
+investigation).\cite{Kol97}
+
+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
+\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^3\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 {\sc dlm}
+integration scheme.
+
+\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)$ 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
+\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 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}
+
+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
+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 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}, \\
+	0 & r_{ij} \geqslant R_\textrm{c},
+\end{array}\right.
+\end{equation}
+where
+\begin{equation}
+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 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 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}) - 
+		\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
+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}, \\
+	\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{\label{sec:LJCorrections}Long-Range Interaction Corrections}
+
+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 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 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 correction techniques, is presented in
+chapter \ref{chap:electrostatics}.
+
+\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 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 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.
+
+