--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/01/28 22:42:46	2575
+++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/21 19:26:59	2652
@@ -1,13 +1,18 @@
 %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
-\documentclass[12pt]{article}
+%\documentclass[aps,prb,preprint]{revtex4}
+\documentclass[11pt]{article}
 \usepackage{endfloat}
-\usepackage{amsmath}
+\usepackage{amsmath,bm}
+\usepackage{amssymb}
 \usepackage{epsf}
 \usepackage{times}
-\usepackage{mathptm}
+\usepackage{mathptmx}
 \usepackage{setspace}
 \usepackage{tabularx}
 \usepackage{graphicx}
+\usepackage{booktabs}
+\usepackage{bibentry}
+\usepackage{mathrsfs}
 \usepackage[ref]{overcite}
 \pagestyle{plain}
 \pagenumbering{arabic}
@@ -20,9 +25,10 @@
 
 \begin{document}
 
-\title{On the necessity of the Ewald Summation in molecular simulations.}
+\title{Is the Ewald Summation necessary? Pairwise alternatives to the accepted standard for long-range electrostatics}
 
-\author{Christopher J. Fennell and J. Daniel Gezelter \\
+\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail:
+gezelter@nd.edu} \\
 Department of Chemistry and Biochemistry\\ 
 University of Notre Dame\\
 Notre Dame, Indiana 46556}
@@ -30,30 +36,1107 @@ Notre Dame, Indiana 46556}
 \date{\today}
 
 \maketitle
-%\doublespacing
+\doublespacing
 
+\nobibliography{}
 \begin{abstract}
+A new method for accumulating electrostatic interactions was derived
+from the previous efforts described in \bibentry{Wolf99} and
+\bibentry{Zahn02} as a possible replacement for lattice sum methods in
+molecular simulations.  Comparisons were performed with this and other
+pairwise electrostatic summation techniques against the smooth
+particle mesh Ewald (SPME) summation to see how well they reproduce
+the energetics and dynamics of a variety of simulation types.  The
+newly derived Shifted-Force technique shows a remarkable ability to
+reproduce the behavior exhibited in simulations using SPME with an
+$\mathscr{O}(N)$ computational cost, equivalent to merely the
+real-space portion of the lattice summation.
+
 \end{abstract}
 
+\newpage
+
 %\narrowtext 
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %	                       BODY OF TEXT
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 \section{Introduction}
+
+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 paper, we focus on a new set of shifted 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.
+
+\subsection{The Ewald Sum}
+The complete accumulation 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 sizes of the systems that could be
+feasibly simulated, the entire simulation box was replicated to
+convergence.  In more modern simulations, the simulation boxes have
+grown large enough that a real-space cutoff could potentially give
+convergent behavior.  Indeed, it has often been observed that the
+reciprocal-space portion of the Ewald sum can be small and rapidly
+convergent compared to the real-space portion with the choice of small
+$\alpha$.\cite{Karasawa89,Kolafa92}
+
+\begin{figure}
+\centering
+\includegraphics[width = \linewidth]{./ewaldProgression.pdf}
+\caption{How the application of the Ewald summation has changed with
+the increase in computer power.  Initially, only small numbers of
+particles could be studied, and the Ewald sum acted to replicate the
+unit cell charge distribution out to convergence.  Now, much larger
+systems of charges are investigated with fixed distance cutoffs.  The
+calculated structure factor is used to sum out to great distance, and
+a surrounding dielectric term is included.}
+\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
+2D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the
+new methods are computationally expensive.\cite{Spohr97,Yeh99}
+Inclusion of a correction term in the Ewald summation is a possible
+direction for handling 2D systems while still enabling the use of the
+modern optimizations.\cite{Yeh99}
+
+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.
+
+\subsection{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
+(Eq. \ref{eq:WolfPot}) 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}
+F_{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{\left[\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}}\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]\right\},
+\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}
+V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\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)\right\},
+\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.
+
+\subsection{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.+\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)\ \right] \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.-\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)\right] \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.
+
+\subsection{Other alternatives}
+In addition to the methods described above, we considered some other
+techniques that are commonly used in molecular simulations.  The
+simplest of these is group-based cutoffs.  Though of little use for
+charged molecules, collecting atoms into neutral groups takes
+advantage of the observation that the electrostatic interactions decay
+faster than those for monopolar pairs.\cite{Steinbach94} When
+considering these molecules as neutral groups, the relative
+orientations of the molecules control the strength of the interactions
+at the cutoff radius.  Consequently, as these molecular particles move
+through $R_\textrm{c}$, the energy will drift upward due to the
+anisotropy of the net molecular dipole interactions.\cite{Rahman71} To
+maintain good energy conservation, both the potential and derivative
+need to be smoothly switched to zero at $R_\textrm{c}$.\cite{Adams79}
+This is accomplished using a standard switching function.  If a smooth
+second derivative is desired, a fifth (or higher) order polynomial can
+be used.\cite{Andrea83}
+
+Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$,
+and to incorporate the effects of the surroundings, a method like
+Reaction Field ({\sc rf}) can be used.  The original theory for {\sc
+rf} was originally developed by Onsager,\cite{Onsager36} and it was
+applied in simulations for the study of water by Barker and
+Watts.\cite{Barker73} In modern simulation codes, {\sc rf} is simply
+an extension of the group-based cutoff method where the net dipole
+within the cutoff sphere polarizes an external dielectric, which
+reacts back on the central dipole.  The same switching function
+considerations for group-based cutoffs need to made for {\sc rf}, with
+the additional pre-specification of a dielectric constant.
+
 \section{Methods}
 
+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
+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 SPME.
+Sample correlation plots for two alternate methods are shown in
+Fig. \ref{fig:linearFit}.
+
+\begin{figure}
+\centering
+\includegraphics[width = \linewidth]{./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 system type (detailed in section \ref{sec:RepSims}) was
+represented using 500 independent configurations.  Additionally, we
+used seven different system types, so each of the alternative
+(non-Ewald) electrostatic summation methods was evaluated using
+873,250 configurational energy differences.
+
+Results and discussion for the individual analysis of each of the
+system types appear in the supporting information, while the
+cumulative results over all the investigated systems appears below in
+section \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 (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 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.
+%
+%\begin{figure}
+%\centering
+%\includegraphics[width = \linewidth]{./gaussFit.pdf}
+%\caption{Sample fit of the angular distribution of the force vectors
+%accumulated using all of the studied systems.  Gaussian fits were used
+%to obtain values for the variance in force and torque vectors.}
+%\label{fig:gaussian}
+%\end{figure}
+%
+%Figure \ref{fig:gaussian} shows an example distribution with applied
+%non-linear fits.  The solid line is a Gaussian profile, while the
+%dotted line is a Voigt profile, a convolution of a Gaussian and a
+%Lorentzian.  
+%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.
+%Gaussian fits was used to compare all the tested methods.  
+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 (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_i(0)\cdot v_i(t)\rangle}{\langle v_i(0)\cdot v_i(0)\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 simulations were analyzed to determine the
+relative effectiveness of the pairwise summation techniques in
+reproducing the energetics and dynamics exhibited by SPME.  We wanted
+to span the space of modern 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 (300~K for the liquid, 200~K 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.  The equilibrated temperatures were 1000~K for the NaCl
+crystal and 7000~K for the liquid. 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.
+%(Fig. \ref{fig:argonSlice}).
+
+These procedures guaranteed us a set of representative configurations
+from chemically-relevant systems sampled from an appropriate
+ensemble. Force field parameters for the ions and Argon were taken
+from the force field utilized by {\sc oopse}.\cite{Meineke05}
+
+%\begin{figure}
+%\centering
+%\includegraphics[width = \linewidth]{./slice.pdf}
+%\caption{A slice from the center of a water box used in a charge void
+%simulation.  The darkened region represents the boundary sphere within
+%which the water molecules were converted to argon atoms.} 
+%\label{fig:argonSlice}
+%\end{figure}
+
+\subsection{Comparison of Summation Methods}\label{sec:ESMethods}
+We compared the following alternative summation methods with results
+from the reference method (SPME):
+\begin{itemize}
+\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{itemize}
+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 SP, SF, and pure
+cutoff.  The SPME electrostatics were performed using the TINKER
+implementation of SPME,\cite{Ponder87} while all other method
+calculations were performed using the 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 althernative 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
+increased accuracy at the expense of increased time spent calculating
+the reciprocal-space portion of the
+summation.\cite{Perram88,Essmann95} The default TINKER tolerance of $1
+\times 10^{-8}$ kcal/mol was used in all 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{Results and Discussion}
+
+\subsection{Configuration Energy Differences}\label{sec:EnergyResults}
+In order to evaluate the performance of the pairwise electrostatic
+summation methods for Monte Carlo simulations, the energy differences
+between configurations were compared to the values obtained when using
+SPME.  The results for the subsequent regression analysis are shown in
+figure \ref{fig:delE}.
+
+\begin{figure}
+\centering
+\includegraphics[width=5.5in]{./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
+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.  Interested
+readers can consult the accompanying supporting information for a
+comparison where all groups are neutral.
+
+For the {\sc sp} method, inclusion of potential 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 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
+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.
+
+\subsection{Magnitudes of the Force and Torque Vectors}
+
+Evaluation of pairwise methods for use in Molecular Dynamics
+simulations requires consideration of effects on the forces and
+torques.  Investigation of the force and torque vector magnitudes
+provides a measure of the strength of these values relative to SPME.
+Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the
+force and torque vector magnitude regression results for the
+accumulated analysis over all the system types.
+
+\begin{figure}
+\centering
+\includegraphics[width=5.5in]{./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 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}
+
+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
+a improvement much more significant than what was seen with $\Delta
+E$.  Looking at the {\sc sp} sets, the slope and $R^2$
+improve with the use of damping to an optimal result of 0.2 \AA
+$^{-1}$ for the 12 and 15 \AA\ cutoffs.  Further increases in damping,
+while beneficial for simulations with a cutoff radius of 9 \AA\ , is
+detrimental to simulations with larger cutoff radii.  The undamped
+{\sc sf} method gives forces in line with those obtained using
+SPME, and use of a damping function results in minor improvement.  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 it correlates well in
+general.  To be fair, we again note that the reaction field
+calculations do not encompass NaCl crystal and melt systems, so these
+results are partly biased towards conditions in which the method
+performs more favorably.
+
+\begin{figure}
+\centering
+\includegraphics[width=5.5in]{./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 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}
+
+To evaluate the torque vector magnitudes, the data set from which
+values are drawn is limited to rigid molecules in the systems
+(i.e. water molecules).  In spite of this smaller sampling pool, the
+torque vector magnitude results in figure \ref{fig:trqMag} are still
+similar to those seen for the forces; however, they more clearly show
+the improved behavior that comes with increasing the cutoff radius.
+Moderate damping is beneficial to the {\sc sp} and helpful
+yet possibly unnecessary with the {\sc sf} method, and they also
+show that over-damping adversely effects all cutoff radii rather than
+showing an improvement for systems with short cutoffs.  The reaction
+field method performs well when calculating the torques, better than
+the Shifted Force method over this limited data set.
+
+\subsection{Directionality of the Force and Torque Vectors}
+
+Having force and torque vectors with magnitudes that are well
+correlated to SPME is good, but if they are not pointing in the proper
+direction the results will be incorrect.  These vector directions were
+investigated through measurement of the angle formed between them and
+those from 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=5.5in]{./frcTrqAngplot.pdf}
+\caption{Statistical analysis of the quality of the Gaussian fit of 
+the force and torque vector angular distributions for a given
+electrostatic method compared with 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 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}.  All of the sets, aside from the over-damped case
+show the improvement afforded by choosing a longer simulation cutoff.
+Increasing the cutoff from 9 to 12 \AA\ typically results in a halving
+of the distribution widths, with a similar improvement 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 damping
+improves the angular behavior significantly for the {\sc sp}
+and moderately for the {\sc sf} methods.  Increasing the damping
+too far is destructive for both methods, particularly to the torque
+vectors.  Again it is important to recognize that the force vectors
+cover all particles in the systems, while torque vectors are only
+available for neutral molecular groups.  Damping appears to have a
+more beneficial effect on non-neutral bodies, and this observation is
+investigated further in the accompanying supporting information.
+
+\begin{table}[htbp]
+   \centering
+   \caption{Variance ($\sigma^2$) of the force (top set) and torque 
+(bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods.  Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function.  The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.}	
+   \begin{tabular}{@{} ccrrrrrrrr @{}}
+      \\
+      \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}
+   \label{tab:groupAngle}
+\end{table}
+
+Although not discussed previously, group based cutoffs can be applied
+to both the {\sc sp} and {\sc sf} methods.  Use off a
+switching function corrects for the discontinuities that arise when
+atoms of a group exit the cutoff before the group's center of mass.
+Though there are no significant benefit or drawbacks observed in
+$\Delta E$ and vector magnitude results when doing this, there is a
+measurable improvement in the vector angle results.  Table
+\ref{tab:groupAngle} shows the angular variance values obtained using
+group based cutoffs and a switching function alongside the standard
+results seen in figure \ref{fig:frcTrqAng} for comparison purposes.
+The {\sc sp} shows much narrower angular distributions for
+both the force and torque vectors when using an $\alpha$ of 0.2
+\AA$^{-1}$ or less, while {\sc sf} shows improvements in the
+undamped and lightly damped cases.  Thus, by calculating the
+electrostatic interactions in terms of molecular pairs rather than
+atomic pairs, the direction of the force and torque vectors are
+determined more accurately.
+
+One additional trend to recognize 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 easier
+to see when using group based cutoffs.  Looking back on figures
+\ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this
+behavior clearly at large $\alpha$ and cutoff values.  The reason for
+this is that the complimentary error function inserted into the
+potential weakens the electrostatic interaction as $\alpha$ increases.
+Thus, at larger values of $\alpha$, both the summation method types
+progress toward non-interacting functions, so care is required in
+choosing large damping functions lest one generate an undesirable loss
+in the pair interaction.  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 the above findings, empirical damping
+up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably
+unnecessary when using the {\sc sf} method.
+
+\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}
+
+In the previous studies using a {\sc sf} variant of the damped
+Wolf coulomb potential, the structure and dynamics of water were
+investigated rather extensively.\cite{Zahn02,Kast03} Their results
+indicated that the damped {\sc sf} method results in properties
+very similar to those obtained when using the Ewald summation.
+Considering the statistical results shown above, the good performance
+of this method is not that surprising.  Rather than consider the same
+systems and simply recapitulate their results, we decided to look at
+the solid state dynamical behavior obtained using the best performing
+summation methods from the above results.
+
+\begin{figure}
+\centering
+\includegraphics[width = \linewidth]{./vCorrPlot.pdf}
+\caption{Velocity auto-correlation functions of NaCl crystals at 
+1000 K while using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and
+{\sc sp} ($\alpha$ = 0.2). The inset is a magnification of the first
+trough. The times to first collision are nearly identical, but the
+differences can be seen in the peaks and troughs, where the undamped
+to weakly damped methods are stiffer than the moderately damped and
+SPME methods.}
+\label{fig:vCorrPlot}
+\end{figure}
+
+The short-time decays 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 function is
+increased, these peaks are smoothed out, and approach the SPME
+curve. The damping acts as a distance dependent Gaussian screening of
+the point charges for the pairwise summation methods; thus, the
+collisions are more elastic in the undamped {\sc sf} potential, and the
+stiffness of the potential is diminished as the electrostatic
+interactions are softened by the damping function.  With $\alpha$
+values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are
+nearly identical and track the SPME features quite well.  This is not
+too surprising in that the differences between the {\sc sf} and {\sc
+sp} potentials are mitigated with increased damping.  However, this
+appears to indicate that once damping is utilized, the form of the
+potential seems to play a lesser role in the crystal dynamics.
+
+\subsection{Collective Motion: Power Spectra of NaCl Crystals}
+
+The short time dynamics were extended to evaluate how the differences
+between the methods affect the collective long-time motion.  The same
+electrostatic summation methods were used as in the short time
+velocity autocorrelation function evaluation, but the trajectories
+were sampled over a much longer time. The power spectra of the
+resulting velocity autocorrelation functions were calculated and are
+displayed in figure \ref{fig:methodPS}.
+
+\begin{figure}
+\centering
+\includegraphics[width = \linewidth]{./spectraSquare.pdf}
+\caption{Power spectra obtained from the velocity auto-correlation 
+functions of NaCl crystals at 1000 K while using SPME, {\sc sf}
+($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2).
+Apodization of the correlation functions via a cubic switching
+function between 40 and 50 ps was used to clear up the spectral noise
+resulting from data truncation, and had no noticeable effect on peak
+location or magnitude.  The inset shows the frequency region below 100
+cm$^{-1}$ to highlight where the spectra begin to differ.}
+\label{fig:methodPS}
+\end{figure}
+
+While high frequency peaks of the spectra in this figure overlap,
+showing the same general features, 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 near identical correlated motion behavior as
+the Ewald method (which has a damping value of 0.3119).  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 see this effect more clearly, we show how
+damping strength alone affects a simple real-space electrostatic
+potential,
+\begin{equation}
+V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r),
+\end{equation}
+where $S(r)$ is a switching function that smoothly zeroes the
+potential at the cutoff radius.  Figure \ref{fig:dampInc} shows how
+the low frequency motions are dependent on the damping used in the
+direct electrostatic sum.  As the damping increases, the peaks drop to
+lower frequencies.  Incidentally, use of an $\alpha$ of 0.25
+\AA$^{-1}$ on a simple electrostatic summation results in low
+frequency correlated dynamics equivalent to a simulation using SPME.
+When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks
+shift to higher frequency in exponential fashion.  Though not shown,
+the spectrum for the simple undamped electrostatic potential is
+blue-shifted such that the lowest frequency peak resides near 325
+cm$^{-1}$.  In light of these results, the undamped {\sc sf} method
+producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite
+respectable and shows that the shifted force procedure accounts for
+most of the effect afforded through use of the Ewald summation.
+However, it appears as though moderate damping is required for
+accurate reproduction of crystal dynamics.
+\begin{figure}
+\centering
+\includegraphics[width = \linewidth]{./comboSquare.pdf}
+\caption{Regions of spectra showing the low-frequency correlated 
+motions for NaCl crystals at 1000 K using various electrostatic
+summation methods.  The upper plot is a zoomed inset from figure
+\ref{fig:methodPS}.  As the damping value for the {\sc sf} potential
+increases, the low-frequency peaks red-shift.  The lower plot is of
+spectra when using SPME and a simple damped Coulombic sum with damping
+coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$.  As
+$\alpha$ increases, the peaks are red-shifted toward and eventually
+beyond the values given by SPME.  The larger $\alpha$ values weaken
+the real-space electrostatics, explaining this shift towards less
+strongly correlated motions in the crystal.}
+\label{fig:dampInc}
+\end{figure}
+
 \section{Conclusions}
 
-\section{Acknowledgments}
+This investigation of pairwise electrostatic summation techniques
+shows that there are viable and more computationally efficient
+electrostatic summation techniques than the Ewald summation, chiefly
+methods derived from the damped Coulombic sum originally proposed by
+Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the
+{\sc sf} method, reformulated above as eq. (\ref{eq:DSFPot}),
+shows a remarkable ability to reproduce the energetic and dynamic
+characteristics exhibited by simulations employing lattice summation
+techniques.  The cumulative energy difference results showed the
+undamped {\sc sf} and moderately damped {\sc sp} methods
+produced results nearly identical to SPME.  Similarly for the dynamic
+features, the undamped or moderately damped {\sc sf} and
+moderately damped {\sc sp} methods produce force and torque
+vector magnitude and directions very similar to the expected values.
+These results translate into long-time dynamic behavior equivalent to
+that produced in simulations using SPME.
 
-\newpage
+Aside from the computational cost benefit, these techniques have
+applicability in situations where the use of the Ewald sum can prove
+problematic.  Primary among them is their use in interfacial systems,
+where the unmodified lattice sum techniques artificially accentuate
+the periodicity of the system in an undesirable manner.  There have
+been alterations to the standard Ewald techniques, via corrections and
+reformulations, to compensate for these systems; but the pairwise
+techniques discussed here require no modifications, making them
+natural tools to tackle these problems.  Additionally, this
+transferability gives them benefits over other pairwise methods, like
+reaction field, because estimations of physical properties (e.g. the
+dielectric constant) are unnecessary.
 
-\bibliographystyle{achemso}
+We are not suggesting any flaw with the Ewald sum; in fact, it is the
+standard by which these simple pairwise sums are judged.  However,
+these results do suggest that in the typical simulations performed
+today, the Ewald summation may no longer be required to obtain the
+level of accuracy most researchers have come to expect
+
+\section{Acknowledgments}
+\newpage 
+
+\bibliographystyle{jcp2}
 \bibliography{electrostaticMethods}