--- trunk/fennellDissertation/dissertation.tex	2006/07/03 14:51:28	2920
+++ trunk/fennellDissertation/dissertation.tex	2006/08/18 15:49:29	2973
@@ -1,12 +1,15 @@
-\documentclass[12pt]{ndthesis}
+\documentclass[11pt]{ndthesis}
 
 % some packages for things like equations and graphics
+\usepackage[tbtags]{amsmath}
 \usepackage{amsmath,bm}
 \usepackage{amssymb}
 \usepackage{mathrsfs}
 \usepackage{tabularx}
 \usepackage{graphicx}
 \usepackage{booktabs}
+\usepackage{cite}
+\usepackage{enumitem}
 
 \begin{document}
 
@@ -39,1093 +42,9 @@ STUDY OF WATER}     
 \mainmatter
 
 \chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND}
-
-\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION 
-TECHNIQUES}
-
-In molecular simulations, proper accumulation of the electrostatic
-interactions is essential and is one of the most
-computationally-demanding tasks.  The common molecular mechanics force
-fields represent atomic sites with full or partial charges protected
-by Lennard-Jones (short range) interactions.  This means that nearly
-every pair interaction involves a calculation of charge-charge forces.
-Coupled with relatively long-ranged $r^{-1}$ decay, the monopole
-interactions quickly become the most expensive part of molecular
-simulations.  Historically, the electrostatic pair interaction would
-not have decayed appreciably within the typical box lengths that could
-be feasibly simulated.  In the larger systems that are more typical of
-modern simulations, large cutoffs should be used to incorporate
-electrostatics correctly.
-
-There have been many efforts to address the proper and practical
-handling of electrostatic interactions, and these have resulted in a
-variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are
-typically classified as implicit methods (i.e., continuum dielectrics,
-static dipolar fields),\cite{Born20,Grossfield00} explicit methods
-(i.e., Ewald summations, interaction shifting or
-truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e.,
-reaction field type methods, fast multipole
-methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are
-often preferred because they physically incorporate solvent molecules
-in the system of interest, but these methods are sometimes difficult
-to utilize because of their high computational cost.\cite{Roux99} In
-addition to the computational cost, there have been some questions
-regarding possible artifacts caused by the inherent periodicity of the
-explicit Ewald summation.\cite{Tobias01}
-
-In this chapter, we focus on a new set of pairwise methods devised by
-Wolf {\it et al.},\cite{Wolf99} which we further extend.  These
-methods along with a few other mixed methods (i.e. reaction field) are
-compared with the smooth particle mesh Ewald
-sum,\cite{Onsager36,Essmann99} which is our reference method for
-handling long-range electrostatic interactions. The new methods for
-handling electrostatics have the potential to scale linearly with
-increasing system size since they involve only a simple modification
-to the direct pairwise sum.  They also lack the added periodicity of
-the Ewald sum, so they can be used for systems which are non-periodic
-or which have one- or two-dimensional periodicity.  Below, these
-methods are evaluated using a variety of model systems to
-establish their usability in molecular simulations.
-
-\section{The Ewald Sum}
-
-The complete accumulation of the electrostatic interactions in a system with
-periodic boundary conditions (PBC) requires the consideration of the
-effect of all charges within a (cubic) simulation box as well as those
-in the periodic replicas,
-\begin{equation}
-V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime 
-\left[ \sum_{i=1}^N\sum_{j=1}^N \phi
-\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right) 
-\right],
-\label{eq:PBCSum}
-\end{equation}
-where the sum over $\mathbf{n}$ is a sum over all periodic box
-replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the
-prime indicates $i = j$ are neglected for $\mathbf{n} =
-0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic
-particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is
-the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and
-$j$, and $\phi$ is the solution to Poisson's equation
-($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for
-charge-charge interactions). In the case of monopole electrostatics,
-eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for
-non-neutral systems.
-
-The electrostatic summation problem was originally studied by Ewald
-for the case of an infinite crystal.\cite{Ewald21}. The approach he
-took was to convert this conditionally convergent sum into two
-absolutely convergent summations: a short-ranged real-space summation
-and a long-ranged reciprocal-space summation,
-\begin{equation}
-\begin{split}
-V_\textrm{elec} = \frac{1}{2}& 
-\sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime 
-\frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)}
-{|\mathbf{r}_{ij}+\mathbf{n}|} \\ 
-&+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} 
-\exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) 
-\cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ 
-&- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 
-+ \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3}
-\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2,
-\end{split}
-\label{eq:EwaldSum}
-\end{equation}
-where $\alpha$ is the damping or convergence parameter with units of
-\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to
-$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric
-constant of the surrounding medium. The final two terms of
-eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term
-for interacting with a surrounding dielectric.\cite{Allen87} This
-dipolar term was neglected in early applications in molecular
-simulations,\cite{Brush66,Woodcock71} until it was introduced by de
-Leeuw {\it et al.} to address situations where the unit cell has a
-dipole moment which is magnified through replication of the periodic
-images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the
-system is said to be using conducting (or ``tin-foil'') boundary
-conditions, $\epsilon_{\rm S} = \infty$. Figure
-\ref{fig:ewaldTime} shows how the Ewald sum has been applied over
-time.  Initially, due to the small system sizes that could be
-simulated feasibly, the entire simulation box was replicated to
-convergence.  In more modern simulations, the systems have grown large
-enough that a real-space cutoff could potentially give convergent
-behavior.  Indeed, it has been observed that with the choice of a
-small $\alpha$, the reciprocal-space portion of the Ewald sum can be
-rapidly convergent and small relative to the real-space
-portion.\cite{Karasawa89,Kolafa92}
-
-\begin{figure}
-\includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf}
-\caption{The change in the need for the Ewald sum with
-increasing computational power.  A:~Initially, only small systems
-could be studied, and the Ewald sum replicated the simulation box to
-convergence.  B:~Now, radial cutoff methods should be able to reach
-convergence for the larger systems of charges that are common today.}
-\label{fig:ewaldTime}
-\end{figure}
-
-The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm.  The
-convergence parameter $(\alpha)$ plays an important role in balancing
-the computational cost between the direct and reciprocal-space
-portions of the summation.  The choice of this value allows one to
-select whether the real-space or reciprocal space portion of the
-summation is an $\mathscr{O}(N^2)$ calculation (with the other being
-$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of
-$\alpha$ and thoughtful algorithm development, this cost can be
-reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route
-taken to reduce the cost of the Ewald summation even further is to set
-$\alpha$ such that the real-space interactions decay rapidly, allowing
-for a short spherical cutoff. Then the reciprocal space summation is
-optimized.  These optimizations usually involve utilization of the
-fast Fourier transform (FFT),\cite{Hockney81} leading to the
-particle-particle particle-mesh (P3M) and particle mesh Ewald (PME)
-methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these
-methods, the cost of the reciprocal-space portion of the Ewald
-summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N
-\log N)$.
-
-These developments and optimizations have made the use of the Ewald
-summation routine in simulations with periodic boundary
-conditions. However, in certain systems, such as vapor-liquid
-interfaces and membranes, the intrinsic three-dimensional periodicity
-can prove problematic.  The Ewald sum has been reformulated to handle
-2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} but these
-methods are computationally expensive.\cite{Spohr97,Yeh99} More
-recently, there have been several successful efforts toward reducing
-the computational cost of 2-D lattice
-summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04}
-bringing them more in line with the cost of the full 3-D summation.
-
-Several studies have recognized that the inherent periodicity in the
-Ewald sum can also have an effect on three-dimensional
-systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
-Solvated proteins are essentially kept at high concentration due to
-the periodicity of the electrostatic summation method.  In these
-systems, the more compact folded states of a protein can be
-artificially stabilized by the periodic replicas introduced by the
-Ewald summation.\cite{Weber00} Thus, care must be taken when
-considering the use of the Ewald summation where the assumed
-periodicity would introduce spurious effects in the system dynamics.
-
-
-\section{The Wolf and Zahn Methods}
-
-In a recent paper by Wolf \textit{et al.}, a procedure was outlined
-for the accurate accumulation of electrostatic interactions in an
-efficient pairwise fashion.  This procedure lacks the inherent
-periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.}
-observed that the electrostatic interaction is effectively
-short-ranged in condensed phase systems and that neutralization of the
-charge contained within the cutoff radius is crucial for potential
-stability. They devised a pairwise summation method that ensures
-charge neutrality and gives results similar to those obtained with the
-Ewald summation.  The resulting shifted Coulomb potential includes
-image-charges subtracted out through placement on the cutoff sphere
-and a distance-dependent damping function (identical to that seen in
-the real-space portion of the Ewald sum) to aid convergence
-\begin{equation}
-V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}
-- \lim_{r_{ij}\rightarrow R_\textrm{c}}
-\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}.
-\label{eq:WolfPot}
-\end{equation}
-Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted
-potential.  However, neutralizing the charge contained within each
-cutoff sphere requires the placement of a self-image charge on the
-surface of the cutoff sphere.  This additional self-term in the total
-potential enabled Wolf {\it et al.}  to obtain excellent estimates of
-Madelung energies for many crystals.
-
-In order to use their charge-neutralized potential in molecular
-dynamics simulations, Wolf \textit{et al.} suggested taking the
-derivative of this potential prior to evaluation of the limit.  This
-procedure gives an expression for the forces,
-\begin{equation}
-\begin{split}
-F_{\textrm{Wolf}}(r_{ij}) = q_i q_j \Biggr\{& 
-\Biggr[\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r^2_{ij}} 
-+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{r_{ij}}
-\Biggr]\\ 
-&-\Biggr[
-\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2} 
-+ \frac{2\alpha}{\pi^{1/2}}
-\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}
-\Biggr]\Biggr\},
-\end{split}
-\label{eq:WolfForces}
-\end{equation}
-that incorporates both image charges and damping of the electrostatic
-interaction. 
-
-More recently, Zahn \textit{et al.} investigated these potential and
-force expressions for use in simulations involving water.\cite{Zahn02}
-In their work, they pointed out that the forces and derivative of
-the potential are not commensurate.  Attempts to use both
-eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead
-to poor energy conservation.  They correctly observed that taking the
-limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the
-derivatives gives forces for a different potential energy function
-than the one shown in eq. (\ref{eq:WolfPot}).
-
-Zahn \textit{et al.} introduced a modified form of this summation
-method as a way to use the technique in Molecular Dynamics
-simulations.  They proposed a new damped Coulomb potential,
-\begin{equation}
-\begin{split}
-V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\Biggr\{&
-\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} \\
-&- \left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}
-+ \frac{2\alpha}{\pi^{1/2}}
-\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}
-\right]\left(r_{ij}-R_\mathrm{c}\right)\Biggr\},
-\end{split}
-\label{eq:ZahnPot}
-\end{equation}
-and showed that this potential does fairly well at capturing the
-structural and dynamic properties of water compared the same
-properties obtained using the Ewald sum.
-
-\section{Simple Forms for Pairwise Electrostatics}
-
-The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et
-al.} are constructed using two different (and separable) computational
-tricks: 
-
-\begin{enumerate}
-\item shifting through the use of image charges, and 
-\item damping the electrostatic interaction.
-\end{enumerate}  
-Wolf \textit{et al.} treated the
-development of their summation method as a progressive application of
-these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded
-their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the
-post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using
-both techniques.  It is possible, however, to separate these
-tricks and study their effects independently. 
-
-Starting with the original observation that the effective range of the
-electrostatic interaction in condensed phases is considerably less
-than $r^{-1}$, either the cutoff sphere neutralization or the
-distance-dependent damping technique could be used as a foundation for
-a new pairwise summation method.  Wolf \textit{et al.} made the
-observation that charge neutralization within the cutoff sphere plays
-a significant role in energy convergence; therefore we will begin our
-analysis with the various shifted forms that maintain this charge
-neutralization.  We can evaluate the methods of Wolf
-\textit{et al.}  and Zahn \textit{et al.} by considering the standard
-shifted potential,
-\begin{equation}
-V_\textrm{SP}(r) = 	\begin{cases}
-v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r >
-R_\textrm{c}  
-\end{cases},
-\label{eq:shiftingPotForm}
-\end{equation}
-and shifted force,
-\begin{equation}
-V_\textrm{SF}(r) = \begin{cases}
-v(r) - v_\textrm{c}
-- \left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c})
-&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} 
-\end{cases},
-\label{eq:shiftingForm}
-\end{equation}
-functions where $v(r)$ is the unshifted form of the potential, and
-$v_c$ is $v(R_\textrm{c})$.  The Shifted Force ({\sc sf}) form ensures
-that both the potential and the forces goes to zero at the cutoff
-radius, while the Shifted Potential ({\sc sp}) form only ensures the
-potential is smooth at the cutoff radius
-($R_\textrm{c}$).\cite{Allen87}
 
-The forces associated with the shifted potential are simply the forces
-of the unshifted potential itself (when inside the cutoff sphere),
-\begin{equation}
-F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right),
-\end{equation}
-and are zero outside.  Inside the cutoff sphere, the forces associated
-with the shifted force form can be written,
-\begin{equation}
-F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d
-v(r)}{dr} \right)_{r=R_\textrm{c}}.
-\end{equation}
- 
-If the potential, $v(r)$, is taken to be the normal Coulomb potential,
-\begin{equation}
-v(r) = \frac{q_i q_j}{r},
-\label{eq:Coulomb}
-\end{equation}
-then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et
-al.}'s undamped prescription:
-\begin{equation}
-V_\textrm{SP}(r) =
-q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad
-r\leqslant R_\textrm{c},
-\label{eq:SPPot}
-\end{equation}
-with associated forces,
-\begin{equation}
-F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) 
-\quad r\leqslant R_\textrm{c}.
-\label{eq:SPForces}
-\end{equation}
-These forces are identical to the forces of the standard Coulomb
-interaction, and cutting these off at $R_c$ was addressed by Wolf
-\textit{et al.} as undesirable.  They pointed out that the effect of
-the image charges is neglected in the forces when this form is
-used,\cite{Wolf99} thereby eliminating any benefit from the method in
-molecular dynamics.  Additionally, there is a discontinuity in the
-forces at the cutoff radius which results in energy drift during MD
-simulations.
-
-The shifted force ({\sc sf}) form using the normal Coulomb potential
-will give,
-\begin{equation}
-V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}} 
-+ \left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] 
-\quad r\leqslant R_\textrm{c}.
-\label{eq:SFPot}
-\end{equation}
-with associated forces,
-\begin{equation}
-F_\textrm{SF}(r) =  q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) 
-\quad r\leqslant R_\textrm{c}.
-\label{eq:SFForces}
-\end{equation}
-This formulation has the benefits that there are no discontinuities at
-the cutoff radius, while the neutralizing image charges are present in
-both the energy and force expressions.  It would be simple to add the
-self-neutralizing term back when computing the total energy of the
-system, thereby maintaining the agreement with the Madelung energies.
-A side effect of this treatment is the alteration in the shape of the
-potential that comes from the derivative term.  Thus, a degree of
-clarity about agreement with the empirical potential is lost in order
-to gain functionality in dynamics simulations.
-
-Wolf \textit{et al.} originally discussed the energetics of the
-shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was
-insufficient for accurate determination of the energy with reasonable
-cutoff distances.  The calculated Madelung energies fluctuated around
-the expected value as the cutoff radius was increased, but the
-oscillations converged toward the correct value.\cite{Wolf99} A
-damping function was incorporated to accelerate the convergence; and
-though alternative forms for the damping function could be
-used,\cite{Jones56,Heyes81} the complimentary error function was
-chosen to mirror the effective screening used in the Ewald summation.
-Incorporating this error function damping into the simple Coulomb
-potential,
-\begin{equation}
-v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
-\label{eq:dampCoulomb}
-\end{equation}
-the shifted potential (eq. (\ref{eq:SPPot})) becomes
-\begin{equation}
-V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}
-- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) 
-\quad r\leqslant R_\textrm{c},
-\label{eq:DSPPot}
-\end{equation}
-with associated forces,
-\begin{equation}
-F_{\textrm{DSP}}(r) = q_iq_j
-\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2} 
-+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) 
-\quad r\leqslant R_\textrm{c}.
-\label{eq:DSPForces}
-\end{equation}
-Again, this damped shifted potential suffers from a
-force-discontinuity at the cutoff radius, and the image charges play
-no role in the forces.  To remedy these concerns, one may derive a
-{\sc sf} variant by including the derivative term in
-eq. (\ref{eq:shiftingForm}),
-\begin{equation}
-\begin{split}
-V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&
-\frac{\mathrm{erfc}\left(\alpha r\right)}{r} 
-- \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ 
-&+ \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} 
-+ \frac{2\alpha}{\pi^{1/2}}
-\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}
-\right)\left(r-R_\mathrm{c}\right)\ \Biggr{]} 
-\quad r\leqslant R_\textrm{c}.
-\label{eq:DSFPot}
-\end{split}
-\end{equation}
-The derivative of the above potential will lead to the following forces,
-\begin{equation}
-\begin{split}
-F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&
-\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}
-+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ 
-&- \left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}
-{R_{\textrm{c}}^2}
-+ \frac{2\alpha}{\pi^{1/2}}
-\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}
-\right)\Biggr{]} 
-\quad r\leqslant R_\textrm{c}.
-\label{eq:DSFForces}
-\end{split}
-\end{equation}
-If the damping parameter $(\alpha)$ is set to zero, the undamped case,
-eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly
-recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}).
-
-This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot}
-derived by Zahn \textit{et al.}; however, there are two important
-differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from
-eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb})
-with $r$ replaced by $R_\textrm{c}$.  This term is {\it not} present
-in the Zahn potential, resulting in a potential discontinuity as
-particles cross $R_\textrm{c}$.  Second, the sign of the derivative
-portion is different.  The missing $v_\textrm{c}$ term would not
-affect molecular dynamics simulations (although the computed energy
-would be expected to have sudden jumps as particle distances crossed
-$R_c$).  The sign problem is a potential source of errors, however.
-In fact, it introduces a discontinuity in the forces at the cutoff,
-because the force function is shifted in the wrong direction and
-doesn't cross zero at $R_\textrm{c}$.
-
-Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an
-electrostatic summation method in which the potential and forces are
-continuous at the cutoff radius and which incorporates the damping
-function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of
-this paper, we will evaluate exactly how good these methods ({\sc sp},
-{\sc sf}, damping) are at reproducing the correct electrostatic
-summation performed by the Ewald sum.
-
-
-\section{Evaluating Pairwise Summation Techniques}
-
-In classical molecular mechanics simulations, there are two primary
-techniques utilized to obtain information about the system of
-interest: Monte Carlo (MC) and Molecular Dynamics (MD).  Both of these
-techniques utilize pairwise summations of interactions between
-particle sites, but they use these summations in different ways.
-
-In MC, the potential energy difference between configurations dictates
-the progression of MC sampling.  Going back to the origins of this
-method, the acceptance criterion for the canonical ensemble laid out
-by Metropolis \textit{et al.} states that a subsequent configuration
-is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where
-$\xi$ is a random number between 0 and 1.\cite{Metropolis53}
-Maintaining the correct $\Delta E$ when using an alternate method for
-handling the long-range electrostatics will ensure proper sampling
-from the ensemble.
-
-In MD, the derivative of the potential governs how the system will
-progress in time.  Consequently, the force and torque vectors on each
-body in the system dictate how the system evolves.  If the magnitude
-and direction of these vectors are similar when using alternate
-electrostatic summation techniques, the dynamics in the short term
-will be indistinguishable.  Because error in MD calculations is
-cumulative, one should expect greater deviation at longer times,
-although methods which have large differences in the force and torque
-vectors will diverge from each other more rapidly.
-
-\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
-
-The pairwise summation techniques (outlined in section
-\ref{sec:ESMethods}) were evaluated for use in MC simulations by 
-studying the energy differences between conformations.  We took the
-{\sc spme}-computed energy difference between two conformations to be the
-correct behavior. An ideal performance by an alternative method would
-reproduce these energy differences exactly (even if the absolute
-energies calculated by the methods are different).  Since none of the
-methods provide exact energy differences, we used linear least squares
-regressions of energy gap data to evaluate how closely the methods
-mimicked the Ewald energy gaps.  Unitary results for both the
-correlation (slope) and correlation coefficient for these regressions
-indicate perfect agreement between the alternative method and {\sc spme}.
-Sample correlation plots for two alternate methods are shown in
-Fig. \ref{fig:linearFit}.
-
-\begin{figure}
-\centering
-\includegraphics[width = \linewidth]{./figures/dualLinear.pdf}
-\caption{Example least squares regressions of the configuration energy
-differences for SPC/E water systems. The upper plot shows a data set
-with a poor correlation coefficient ($R^2$), while the lower plot
-shows a data set with a good correlation coefficient.} 
-\label{fig:linearFit} 
-\end{figure}
-
-Each of the seven system types (detailed in section \ref{sec:RepSims})
-were represented using 500 independent configurations.  Thus, each of
-the alternative (non-Ewald) electrostatic summation methods was
-evaluated using an accumulated 873,250 configurational energy
-differences.
-
-Results and discussion for the individual analysis of each of the
-system types appear in sections \ref{sec:SystemResults}, while the
-cumulative results over all the investigated systems appear below in
-sections \ref{sec:EnergyResults}.
-
-\subsection{Molecular Dynamics and the Force and Torque 
-Vectors}\label{sec:MDMethods} We evaluated the pairwise methods
-(outlined in section \ref{sec:ESMethods}) for use in MD simulations by
-comparing the force and torque vectors with those obtained using the
-reference Ewald summation ({\sc spme}).  Both the magnitude and the
-direction of these vectors on each of the bodies in the system were
-analyzed.  For the magnitude of these vectors, linear least squares
-regression analyses were performed as described previously for
-comparing $\Delta E$ values.  Instead of a single energy difference
-between two system configurations, we compared the magnitudes of the
-forces (and torques) on each molecule in each configuration.  For a
-system of 1000 water molecules and 40 ions, there are 1040 force
-vectors and 1000 torque vectors.  With 500 configurations, this
-results in 520,000 force and 500,000 torque vector comparisons.
-Additionally, data from seven different system types was aggregated
-before the comparison was made.
-
-The {\it directionality} of the force and torque vectors was
-investigated through measurement of the angle ($\theta$) formed
-between those computed from the particular method and those from {\sc spme},
-\begin{equation}
-\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} 
-\cdot \hat{F}_\textrm{M}\right),
-\end{equation}
-where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force
-vector computed using method M.  Each of these $\theta$ values was
-accumulated in a distribution function and weighted by the area on the
-unit sphere.  Since this distribution is a measure of angular error
-between two different electrostatic summation methods, there is no
-{\it a priori} reason for the profile to adhere to any specific
-shape. Thus, gaussian fits were used to measure the width of the
-resulting distributions. The variance ($\sigma^2$) was extracted from
-each of these fits and was used to compare distribution widths.
-Values of $\sigma^2$ near zero indicate vector directions
-indistinguishable from those calculated when using the reference
-method ({\sc spme}).
-
-\subsection{Short-time Dynamics}
-
-The effects of the alternative electrostatic summation methods on the
-short-time dynamics of charged systems were evaluated by considering a
-NaCl crystal at a temperature of 1000 K.  A subset of the best
-performing pairwise methods was used in this comparison.  The NaCl
-crystal was chosen to avoid possible complications from the treatment
-of orientational motion in molecular systems.  All systems were
-started with the same initial positions and velocities.  Simulations
-were performed under the microcanonical ensemble, and velocity
-autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each
-of the trajectories,
-\begin{equation}
-C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}.
-\label{eq:vCorr}
-\end{equation}
-Velocity autocorrelation functions require detailed short time data,
-thus velocity information was saved every 2 fs over 10 ps
-trajectories. Because the NaCl crystal is composed of two different
-atom types, the average of the two resulting velocity autocorrelation
-functions was used for comparisons.
-
-\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
-
-The effects of the same subset of alternative electrostatic methods on
-the {\it long-time} dynamics of charged systems were evaluated using
-the same model system (NaCl crystals at 1000K).  The power spectrum
-($I(\omega)$) was obtained via Fourier transform of the velocity
-autocorrelation function, \begin{equation} I(\omega) =
-\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt,
-\label{eq:powerSpec}
-\end{equation}
-where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the
-NaCl crystal is composed of two different atom types, the average of
-the two resulting power spectra was used for comparisons. Simulations
-were performed under the microcanonical ensemble, and velocity
-information was saved every 5~fs over 100~ps trajectories. 
-
-\subsection{Representative Simulations}\label{sec:RepSims}
-A variety of representative molecular simulations were analyzed to
-determine the relative effectiveness of the pairwise summation
-techniques in reproducing the energetics and dynamics exhibited by
-{\sc spme}.  We wanted to span the space of typical molecular
-simulations (i.e. from liquids of neutral molecules to ionic
-crystals), so the systems studied were:
-
-\begin{enumerate}
-\item liquid water (SPC/E),\cite{Berendsen87}
-\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E),
-\item NaCl crystals,
-\item NaCl melts,
-\item a low ionic strength solution of NaCl in water (0.11 M),
-\item a high ionic strength solution of NaCl in water (1.1 M), and
-\item a 6\AA\  radius sphere of Argon in water.
-\end{enumerate}
-
-By utilizing the pairwise techniques (outlined in section
-\ref{sec:ESMethods}) in systems composed entirely of neutral groups,
-charged particles, and mixtures of the two, we hope to discern under
-which conditions it will be possible to use one of the alternative
-summation methodologies instead of the Ewald sum.
-
-For the solid and liquid water configurations, configurations were
-taken at regular intervals from high temperature trajectories of 1000
-SPC/E water molecules.  Each configuration was equilibrated
-independently at a lower temperature (300K for the liquid, 200K for
-the crystal).  The solid and liquid NaCl systems consisted of 500
-$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions.  Configurations for
-these systems were selected and equilibrated in the same manner as the
-water systems. In order to introduce measurable fluctuations in the
-configuration energy differences, the crystalline simulations were
-equilibrated at 1000K, near the $T_\textrm{m}$ for NaCl. The liquid
-NaCl configurations needed to represent a fully disordered array of
-point charges, so the high temperature of 7000K was selected for
-equilibration. The ionic solutions were made by solvating 4 (or 40)
-ions in a periodic box containing 1000 SPC/E water molecules.  Ion and
-water positions were then randomly swapped, and the resulting
-configurations were again equilibrated individually.  Finally, for the
-Argon / Water ``charge void'' systems, the identities of all the SPC/E
-waters within 6\AA\ of the center of the equilibrated water
-configurations were converted to argon.
-
-These procedures guaranteed us a set of representative configurations
-from chemically-relevant systems sampled from appropriate
-ensembles. Force field parameters for the ions and Argon were taken
-from the force field utilized by {\sc oopse}.\cite{Meineke05}
-
-\subsection{Comparison of Summation Methods}\label{sec:ESMethods}
-We compared the following alternative summation methods with results
-from the reference method ({\sc spme}):
-
-\begin{enumerate}
-\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
-and 0.3\AA$^{-1}$,
-\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
-and 0.3\AA$^{-1}$,
-\item reaction field with an infinite dielectric constant, and 
-\item an unmodified cutoff.
-\end{enumerate}
-
-Group-based cutoffs with a fifth-order polynomial switching function
-were utilized for the reaction field simulations.  Additionally, we
-investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure
-cutoff.  The {\sc spme} electrostatics were performed using the {\sc tinker}
-implementation of {\sc spme},\cite{Ponder87} while all other calculations
-were performed using the {\sc oopse} molecular mechanics
-package.\cite{Meineke05} All other portions of the energy calculation
-(i.e. Lennard-Jones interactions) were handled in exactly the same
-manner across all systems and configurations.
-
-The alternative methods were also evaluated with three different
-cutoff radii (9, 12, and 15\AA).  As noted previously, the
-convergence parameter ($\alpha$) plays a role in the balance of the
-real-space and reciprocal-space portions of the Ewald calculation.
-Typical molecular mechanics packages set this to a value dependent on
-the cutoff radius and a tolerance (typically less than $1 \times
-10^{-4}$ kcal/mol).  Smaller tolerances are typically associated with
-increasing accuracy at the expense of computational time spent on the
-reciprocal-space portion of the summation.\cite{Perram88,Essmann95}
-The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used
-in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200,
-0.3119, and 0.2476\AA$^{-1}$ for cutoff radii of 9, 12, and 15\AA\
-respectively.
-
-\section{Combined Configuration Energy Difference Results}\label{sec:EnergyResults}
-In order to evaluate the performance of the pairwise electrostatic
-summation methods for Monte Carlo (MC) simulations, the energy
-differences between configurations were compared to the values
-obtained when using {\sc spme}.  The results for the combined
-regression analysis of all of the systems are shown in figure
-\ref{fig:delE}.
-
-\begin{figure}
-\centering
-\includegraphics[width=4.75in]{./figures/delEplot.pdf}
-\caption{Statistical analysis of the quality of configurational energy
-differences for a given electrostatic method compared with the
-reference Ewald sum.  Results with a value equal to 1 (dashed line)
-indicate $\Delta E$ values indistinguishable from those obtained using
-{\sc spme}.  Different values of the cutoff radius are indicated with
-different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ =
-inverted triangles).} 
-\label{fig:delE}
-\end{figure}
-
-The most striking feature of this plot is how well the Shifted Force
-({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy
-differences.  For the undamped {\sc sf} method, and the
-moderately-damped {\sc sp} methods, the results are nearly
-indistinguishable from the Ewald results.  The other common methods do
-significantly less well.
-
-The unmodified cutoff method is essentially unusable.  This is not
-surprising since hard cutoffs give large energy fluctuations as atoms
-or molecules move in and out of the cutoff
-radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to
-some degree by using group based cutoffs with a switching
-function.\cite{Adams79,Steinbach94,Leach01} However, we do not see
-significant improvement using the group-switched cutoff because the
-salt and salt solution systems contain non-neutral groups.  Section
-\ref{sec:SystemResults} includes results for systems comprised entirely 
-of neutral groups.
-
-For the {\sc sp} method, inclusion of electrostatic damping improves
-the agreement with Ewald, and using an $\alpha$ of 0.2\AA $^{-1}$
-shows an excellent correlation and quality of fit with the {\sc spme}
-results, particularly with a cutoff radius greater than 12
-\AA .  Use of a larger damping parameter is more helpful for the
-shortest cutoff shown, but it has a detrimental effect on simulations
-with larger cutoffs.  
-
-In the {\sc sf} sets, increasing damping results in progressively {\it
-worse} correlation with Ewald.  Overall, the undamped case is the best
-performing set, as the correlation and quality of fits are
-consistently superior regardless of the cutoff distance.  The undamped
-case is also less computationally demanding (because no evaluation of
-the complementary error function is required).
-
-The reaction field results illustrates some of that method's
-limitations, primarily that it was developed for use in homogenous
-systems; although it does provide results that are an improvement over
-those from an unmodified cutoff.
-
-\section{Magnitude of the Force and Torque Vector Results}
-
-Evaluation of pairwise methods for use in Molecular Dynamics
-simulations requires consideration of effects on the forces and
-torques.  Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the
-regression results for the force and torque vector magnitudes,
-respectively.  The data in these figures was generated from an
-accumulation of the statistics from all of the system types.
-
-\begin{figure}
-\centering
-\includegraphics[width=4.75in]{./figures/frcMagplot.pdf}
-\caption{Statistical analysis of the quality of the force vector 
-magnitudes for a given electrostatic method compared with the
-reference Ewald sum.  Results with a value equal to 1 (dashed line)
-indicate force magnitude values indistinguishable from those obtained
-using {\sc spme}.  Different values of the cutoff radius are indicated with
-different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ =
-inverted triangles).}
-\label{fig:frcMag}
-\end{figure}
-
-Again, it is striking how well the Shifted Potential and Shifted Force
-methods are doing at reproducing the {\sc spme} forces.  The undamped and
-weakly-damped {\sc sf} method gives the best agreement with Ewald.
-This is perhaps expected because this method explicitly incorporates a
-smooth transition in the forces at the cutoff radius as well as the
-neutralizing image charges.
-
-Figure \ref{fig:frcMag}, for the most part, parallels the results seen
-in the previous $\Delta E$ section.  The unmodified cutoff results are
-poor, but using group based cutoffs and a switching function provides
-an improvement much more significant than what was seen with $\Delta
-E$.
-
-With moderate damping and a large enough cutoff radius, the {\sc sp}
-method is generating usable forces.  Further increases in damping,
-while beneficial for simulations with a cutoff radius of 9\AA\ , is
-detrimental to simulations with larger cutoff radii.
-
-The reaction field results are surprisingly good, considering the poor
-quality of the fits for the $\Delta E$ results.  There is still a
-considerable degree of scatter in the data, but the forces correlate
-well with the Ewald forces in general.  We note that the reaction
-field calculations do not include the pure NaCl systems, so these
-results are partly biased towards conditions in which the method
-performs more favorably.
-
-\begin{figure}
-\centering
-\includegraphics[width=4.75in]{./figures/trqMagplot.pdf}
-\caption{Statistical analysis of the quality of the torque vector 
-magnitudes for a given electrostatic method compared with the
-reference Ewald sum.  Results with a value equal to 1 (dashed line)
-indicate torque magnitude values indistinguishable from those obtained
-using {\sc spme}.  Different values of the cutoff radius are indicated with
-different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ =
-inverted triangles).}
-\label{fig:trqMag}
-\end{figure}
-
-Molecular torques were only available from the systems which contained
-rigid molecules (i.e. the systems containing water).  The data in
-fig. \ref{fig:trqMag} is taken from this smaller sampling pool.
-
-Torques appear to be much more sensitive to charges at a longer
-distance.   The striking feature in comparing the new electrostatic
-methods with {\sc spme} is how much the agreement improves with increasing
-cutoff radius.  Again, the weakly damped and undamped {\sc sf} method
-appears to be reproducing the {\sc spme} torques most accurately.  
-
-Water molecules are dipolar, and the reaction field method reproduces
-the effect of the surrounding polarized medium on each of the
-molecular bodies. Therefore it is not surprising that reaction field
-performs best of all of the methods on molecular torques.
-
-\section{Directionality of the Force and Torque Vector Results}
+\input{electrostaticsChapter}
 
-It is clearly important that a new electrostatic method can reproduce
-the magnitudes of the force and torque vectors obtained via the Ewald
-sum. However, the {\it directionality} of these vectors will also be
-vital in calculating dynamical quantities accurately.  Force and
-torque directionalities were investigated by measuring the angles
-formed between these vectors and the same vectors calculated using
-{\sc spme}.  The results (Fig. \ref{fig:frcTrqAng}) are compared through the
-variance ($\sigma^2$) of the Gaussian fits of the angle error
-distributions of the combined set over all system types.
-
-\begin{figure}
-\centering
-\includegraphics[width=4.75in]{./figures/frcTrqAngplot.pdf}
-\caption{Statistical analysis of the width of the angular distribution
-that the force and torque vectors from a given electrostatic method
-make with their counterparts obtained using the reference Ewald sum.
-Results with a variance ($\sigma^2$) equal to zero (dashed line)
-indicate force and torque directions indistinguishable from those
-obtained using {\sc spme}.  Different values of the cutoff radius are
-indicated with different symbols (9\AA\ = circles, 12\AA\ = squares,
-and 15\AA\ = inverted triangles).}
-\label{fig:frcTrqAng}
-\end{figure}
-
-Both the force and torque $\sigma^2$ results from the analysis of the
-total accumulated system data are tabulated in figure
-\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc
-sp}) method would be essentially unusable for molecular dynamics
-unless the damping function is added.  The Shifted Force ({\sc sf})
-method, however, is generating force and torque vectors which are
-within a few degrees of the Ewald results even with weak (or no)
-damping.
-
-All of the sets (aside from the over-damped case) show the improvement
-afforded by choosing a larger cutoff radius.  Increasing the cutoff
-from 9 to 12\AA\ typically results in a halving of the width of the
-distribution, with a similar improvement when going from 12 to 15
-\AA . 
-
-The undamped {\sc sf}, group-based cutoff, and reaction field methods
-all do equivalently well at capturing the direction of both the force
-and torque vectors.  Using the electrostatic damping improves the
-angular behavior significantly for the {\sc sp} and moderately for the
-{\sc sf} methods.  Overdamping is detrimental to both methods.  Again
-it is important to recognize that the force vectors cover all
-particles in all seven systems, while torque vectors are only
-available for neutral molecular groups.  Damping is more beneficial to
-charged bodies, and this observation is investigated further in
-section \ref{SystemResults}.
-
-Although not discussed previously, group based cutoffs can be applied
-to both the {\sc sp} and {\sc sf} methods.  The group-based cutoffs
-will reintroduce small discontinuities at the cutoff radius, but the
-effects of these can be minimized by utilizing a switching function.
-Though there are no significant benefits or drawbacks observed in
-$\Delta E$ and the force and torque magnitudes when doing this, there
-is a measurable improvement in the directionality of the forces and
-torques. Table \ref{tab:groupAngle} shows the angular variances
-obtained both without (N) and with (Y) group based cutoffs and a
-switching function.  Note that the $\alpha$ values have units of
-\AA$^{-1}$ and the variance values have units of degrees$^2$.  The
-{\sc sp} (with an $\alpha$ of 0.2\AA$^{-1}$ or smaller) shows much
-narrower angular distributions when using group-based cutoffs. The
-{\sc sf} method likewise shows improvement in the undamped and lightly
-damped cases.
-
-\begin{table}
-\caption{STATISTICAL ANALYSIS OF THE ANGULAR DISTRIBUTIONS ($\sigma^2$) 
-THAT THE FORCE ({\it upper}) AND TORQUE ({\it lower}) VECTORS FROM A
-GIVEN ELECTROSTATIC METHOD MAKE WITH THEIR COUNTERPARTS OBTAINED USING
-THE REFERENCE EWALD SUMMATION}
-
-\footnotesize
-\begin{center}
-\begin{tabular}{@{} ccrrrrrrrr @{}} \\ 
-\toprule 
-\toprule 
-
-& &\multicolumn{4}{c}{Shifted Potential} & \multicolumn{4}{c}{Shifted
-Force} \\
-\cmidrule(lr){3-6} \cmidrule(l){7-10} $R_\textrm{c}$ & Groups & 
-$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$ &
-$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$\\
-
-\midrule
-     
-9\AA   & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\
-       & \textbf{Y} & \textbf{2.486} & \textbf{2.160} & \textbf{0.667} & \textbf{0.608} & \textbf{1.768} & \textbf{1.766} & \textbf{0.676} & \textbf{0.609} \\
-12\AA  & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\
-       & \textbf{Y} & \textbf{0.515} & \textbf{0.288} & \textbf{0.127} & \textbf{0.586} & \textbf{0.308} & \textbf{0.249} & \textbf{0.127} & \textbf{0.586} \\
-15\AA  & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\
-       & \textbf{Y} & \textbf{0.228} & \textbf{0.099} & \textbf{0.121} & \textbf{0.598} & \textbf{0.144} & \textbf{0.090} & \textbf{0.121} & \textbf{0.598} \\
-
-\midrule
-      
-9\AA   & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\
-       & \textbf{Y} & \textbf{2.115} & \textbf{1.914} & \textbf{1.878} & \textbf{5.142} & \textbf{2.076} & \textbf{2.039} & \textbf{1.972} & \textbf{5.146} \\
-12\AA  & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\
-       & \textbf{Y} & \textbf{0.810} & \textbf{0.685} & \textbf{1.352} & \textbf{5.082} & \textbf{0.765} & \textbf{0.714} & \textbf{1.360} & \textbf{5.082} \\
-15\AA  & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\
-       & \textbf{Y} & \textbf{0.282} & \textbf{0.294} & \textbf{1.272} & \textbf{4.999} & \textbf{0.324} & \textbf{0.318} & \textbf{1.272} & \textbf{4.999} \\
-
-\bottomrule 
-\end{tabular} 
-\end{center}
-\label{tab:groupAngle}
-\end{table}
-
-One additional trend in table \ref{tab:groupAngle} is that the
-$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$
-increases, something that is more obvious with group-based cutoffs.
-The complimentary error function inserted into the potential weakens
-the electrostatic interaction as the value of $\alpha$ is increased.
-However, at larger values of $\alpha$, it is possible to overdamp the
-electrostatic interaction and to remove it completely.  Kast
-\textit{et al.}  developed a method for choosing appropriate $\alpha$
-values for these types of electrostatic summation methods by fitting
-to $g(r)$ data, and their methods indicate optimal values of 0.34,
-0.25, and 0.16\AA$^{-1}$ for cutoff values of 9, 12, and 15\AA\
-respectively.\cite{Kast03} These appear to be reasonable choices to
-obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on
-these findings, choices this high would introduce error in the
-molecular torques, particularly for the shorter cutoffs.  Based on our
-observations, empirical damping up to 0.2\AA$^{-1}$ is beneficial,
-but damping may be unnecessary when using the {\sc sf} method.
-
-\section{Individual System Analysis Results}
-
-The combined results of the previous sections show how the pairwise
-methods compare to the Ewald summation in the general sense over all
-of the system types.  It is also useful to consider each of the
-studied systems in an individual fashion, so that we can identify
-conditions that are particularly difficult for a selected pairwise
-method to address. This allows us to further establish the limitations
-of these pairwise techniques.  Below, the energy difference, force
-vector, and torque vector analyses are presented on an individual
-system basis.
-
-\subsection{SPC/E Water Results}
-
-\subsection{SPC/E Ice I$_\textrm{c}$ Results}
-
-\subsection{NaCl Melt Results}
-
-\subsection{NaCl Crystal Results}
-
-\subsection{0.1M NaCl Solution Results}
-
-\subsection{1M NaCl Solution Results}
-
-\subsection{6\AA\ Argon Sphere in SPC/E Water Results}
-
-\section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}
-
-Zahn {\it et al.} investigated the structure and dynamics of water
-using eqs. (\ref{eq:ZahnPot}) and
-(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated
-that a method similar (but not identical with) the damped {\sc sf}
-method resulted in properties very similar to those obtained when
-using the Ewald summation.  The properties they studied (pair
-distribution functions, diffusion constants, and velocity and
-orientational correlation functions) may not be particularly sensitive
-to the long-range and collective behavior that governs the
-low-frequency behavior in crystalline systems.  Additionally, the
-ionic crystals are the worst case scenario for the pairwise methods
-because they lack the reciprocal space contribution contained in the
-Ewald summation.  
-
-We are using two separate measures to probe the effects of these
-alternative electrostatic methods on the dynamics in crystalline
-materials.  For short- and intermediate-time dynamics, we are
-computing the velocity autocorrelation function, and for long-time
-and large length-scale collective motions, we are looking at the
-low-frequency portion of the power spectrum.
-
-\begin{figure}
-\centering
-\includegraphics[width = \linewidth]{./figures/vCorrPlot.pdf}
-\caption{Velocity autocorrelation functions of NaCl crystals at 
-1000K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc
-sp} ($\alpha$ = 0.2). The inset is a magnification of the area around
-the first minimum.  The times to first collision are nearly identical,
-but differences can be seen in the peaks and troughs, where the
-undamped and weakly damped methods are stiffer than the moderately
-damped and {\sc spme} methods.}
-\label{fig:vCorrPlot}
-\end{figure}
-
-The short-time decay of the velocity autocorrelation function through
-the first collision are nearly identical in figure
-\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show
-how the methods differ.  The undamped {\sc sf} method has deeper
-troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than
-any of the other methods.  As the damping parameter ($\alpha$) is
-increased, these peaks are smoothed out, and the {\sc sf} method
-approaches the {\sc spme} results.  With $\alpha$ values of 0.2\AA$^{-1}$,
-the {\sc sf} and {\sc sp} functions are nearly identical and track the
-{\sc spme} features quite well.  This is not surprising because the {\sc sf}
-and {\sc sp} potentials become nearly identical with increased
-damping.  However, this appears to indicate that once damping is
-utilized, the details of the form of the potential (and forces)
-constructed out of the damped electrostatic interaction are less
-important.
-
-\section{Collective Motion: Power Spectra of NaCl Crystals}
-
-To evaluate how the differences between the methods affect the
-collective long-time motion, we computed power spectra from long-time
-traces of the velocity autocorrelation function. The power spectra for
-the best-performing alternative methods are shown in
-fig. \ref{fig:methodPS}.  Apodization of the correlation functions via
-a cubic switching function between 40 and 50 ps was used to reduce the
-ringing resulting from data truncation.  This procedure had no
-noticeable effect on peak location or magnitude.
-
-\begin{figure}
-\centering
-\includegraphics[width = \linewidth]{./figures/spectraSquare.pdf}
-\caption{Power spectra obtained from the velocity auto-correlation 
-functions of NaCl crystals at 1000K while using {\sc spme}, {\sc sf}
-($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2).  The inset
-shows the frequency region below 100 cm$^{-1}$ to highlight where the
-spectra differ.}
-\label{fig:methodPS}
-\end{figure}
-
-While the high frequency regions of the power spectra for the
-alternative methods are quantitatively identical with Ewald spectrum,
-the low frequency region shows how the summation methods differ.
-Considering the low-frequency inset (expanded in the upper frame of
-figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the
-correlated motions are blue-shifted when using undamped or weakly
-damped {\sc sf}.  When using moderate damping ($\alpha = 0.2$
-\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical
-correlated motion to the Ewald method (which has a convergence
-parameter of 0.3119\AA$^{-1}$).  This weakening of the electrostatic
-interaction with increased damping explains why the long-ranged
-correlated motions are at lower frequencies for the moderately damped
-methods than for undamped or weakly damped methods. 
-
-To isolate the role of the damping constant, we have computed the
-spectra for a single method ({\sc sf}) with a range of damping
-constants and compared this with the {\sc spme} spectrum.
-Fig. \ref{fig:dampInc} shows more clearly that increasing the
-electrostatic damping red-shifts the lowest frequency phonon modes.
-However, even without any electrostatic damping, the {\sc sf} method
-has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode.
-Without the {\sc sf} modifications, an undamped (pure cutoff) method
-would predict the lowest frequency peak near 325 cm$^{-1}$.  {\it
-Most} of the collective behavior in the crystal is accurately captured
-using the {\sc sf} method.  Quantitative agreement with Ewald can be
-obtained using moderate damping in addition to the shifting at the
-cutoff distance.
-
-\begin{figure}
-\centering
-\includegraphics[width = \linewidth]{./figures/increasedDamping.pdf}
-\caption{Effect of damping on the two lowest-frequency phonon modes in
-the NaCl crystal at 1000K.  The undamped shifted force ({\sc sf})
-method is off by less than 10 cm$^{-1}$, and increasing the
-electrostatic damping to 0.25\AA$^{-1}$ gives quantitative agreement
-with the power spectrum obtained using the Ewald sum.  Overdamping can
-result in underestimates of frequencies of the long-wavelength
-motions.}
-\label{fig:dampInc}
-\end{figure}
-
-
 \chapter{\label{chap:water}SIMPLE MODELS FOR WATER}
 
 \chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS}
@@ -1138,7 +57,7 @@ SIMULATIONS}
 \backmatter
  
 \bibliographystyle{ndthesis}
-\bibliography{dissertation}           
+\bibliography{dissertation}  
 
 \end{document}