--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/15 15:51:44	2623
+++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/15 17:09:09	2624
@@ -98,7 +98,7 @@ for an accurate accumulation of electrostatic interact
 
 \subsection{The Wolf and Zahn Methods}
 In a recent paper by Wolf \textit{et al.}, a procedure was outlined
-for an accurate accumulation of electrostatic interactions in an
+for the accurate accumulation of electrostatic interactions in an
 efficient pairwise fashion.\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
@@ -111,14 +111,14 @@ V^{\textrm{Wolf}}(r_{ij})= \frac{q_iq_j \textrm{erfc}(
 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_iq_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\}.
+V_{\textrm{Wolf}}(r_{ij})= \frac{q_iq_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 enables Wolf {\it et al.}  to obtain excellent estimates of
+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
@@ -126,7 +126,7 @@ F^{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{-\left[\frac
 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}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{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\},
+F_{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{-\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{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
@@ -134,12 +134,13 @@ In their work, they pointed out that the method that t
 
 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 method 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 is mathematically invalid.
+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.} proposed a modified form of this ``Wolf summation
 method'' as a way to use this technique in Molecular Dynamics
@@ -147,7 +148,7 @@ V^{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\frac{\mathrm
 \ref{eq:WolfForces}, 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\}.
+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}
 They showed that this potential does fairly well at capturing the
@@ -158,10 +159,10 @@ tricks: \begin{itemize}
 
 The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et
 al.} are constructed using two different (and separable) computational
-tricks: \begin{itemize}
+tricks: \begin{enumerate}
 \item shifting through the use of image charges, and 
 \item damping the electrostatic interaction.
-\end{itemize}  Wolf \textit{et al.} treated the
+\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
@@ -181,7 +182,7 @@ v^\textrm{SP}(r) = 	\begin{cases}
 \textit{et al.}  and Zahn \textit{et al.} by considering the standard
 shifted potential,
 \begin{equation}
-v^\textrm{SP}(r) = 	\begin{cases}
+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},
@@ -189,8 +190,8 @@ v^\textrm{SF}(r) = 	\begin{cases}
 \end{equation}
 and shifted force,
 \begin{equation}
-v^\textrm{SF}(r) = 	\begin{cases}
-v(r)-v_\textrm{c}-\left(\frac{\textrm{d}v(r)}{\textrm{d}r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c})
+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}
@@ -202,82 +203,107 @@ potential is smooth at the cutoff radius
 potential is smooth at the cutoff radius
 ($R_\textrm{c}$).\cite{Allen87}
 
-
-
-
-If the derivative term is taken to be zero, we are left with the shifted Coulomb potential devised by Wolf \textit{et al.},\cite{Wolf99}
+The forces associated with the shifted potential are simply the forces
+of the unshifted potential itself (when inside the cutoff sphere),
 \begin{equation}
-V^\textrm{SP}(r_{ij}) = q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}.				\label{eq:WolfSP}
+F_{\textrm{SP}} = \left( \frac{d v(r)}{dr} \right),
 \end{equation}
-The forces associated with this potential are obtained by taking the derivative, resulting in the following,
+and are zero outside.  Inside the cutoff sphere, the forces associated
+with the shifted force form can be written,
 \begin{equation}
-F^\textrm{SP}(r_{ij}) = q_iq_j\left(-\frac{1}{r_{ij}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}.
+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:WolfSP}
+\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:FWolfSP}
 \end{equation}
-These forces are identical to the forces of the standard electrostatic interaction, and this 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 they would expect there to be some pressure exerted due to their presence.\cite{Wolf99}  As a consequence the forces, though mathematically valid, may not be physically correct due to this missing component.  Additionally, there is a discontinuity in the forces.  This can be remedied with the use of a switching function to zero the potential and forces smoothly as particles near $R_\textrm{c}$.  
+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.
 
-If the derivative term in equation \ref{eq:shiftingForm} is evaluated, we obtain an hitherto undiscussed shifted force Coulomb potential,
+The shifted force ({\sc sf}) form using the normal Coulomb potential
+will give,
 \begin{equation}
-V^\textrm{SF}(r_{ij}) = q_iq_j\left[\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r_{ij}-R_\textrm{c})\right] \quad r_{ij}\leqslant R_\textrm{c}.
+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}
-Taking the derivative of this shifted force potential gives the
-following forces,
+with associated forces,
 \begin{equation}
-F^\textrm{SF}(r_{ij} =  q_iq_j\left(-\frac{1}{r_{ij}^2}+\frac{1}{R_\textrm{c}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}.
+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}
-Using this formulation rather than the simple shifted potential
-(Eq. \ref{eq:WolfSP}) means that there are no discontinuities in the
-forces in addition to the potential.  This form also has the benefit
-that the image charges have a force presence, addressing the concerns
-about a missing physical component.  One side effect of this treatment
-is a slight alteration in the shape of the potential that comes about
-from the derivative term.  Thus, a degree of clarity about the
-original formulation of the potential is lost in order to gain
-functionality in dynamics simulations.
+This formulation has the benefits that there are no discontinuities at
+the cutoff distance, 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:WolfSP}), and they found that
-it was still insufficient for accurate determination of the energy.
-The energy would fluctuate around the expected value with increasing
-cutoff radius, but the oscillations appeared to be converging toward
-the correct value.\cite{Wolf99} A damping function was incorporated to
-accelerate convergence; and though alternative functional forms could
-be used,\cite{Jones56,Heyes81} the complimentary error function was
-chosen to draw parallels to the Ewald summation.  Incorporating
-damping into the simple Coulomb potential,
+it was still insufficient for accurate determination of the energy
+with reasonable cutoff distances.  The calculated Madelung energies
+fluctuate around the expected value with increasing cutoff radius, but
+the oscillations converge toward the correct value.\cite{Wolf99} A
+damping function was incorporated to accelerate the convergence; and
+though alternative functional forms 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_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}},
+v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
 \label{eq:dampCoulomb}
 \end{equation}
-the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived
+the shifted potential (Eq. \ref{eq:WolfSP}) can be recovered
 \textit{via} equation \ref{eq:shiftingForm},
 \begin{equation}
-V^{\textrm{DSP}}(r_{ij}) = q_iq_j\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}.
+v_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}(\alpha r)}{r}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c}.
 \label{eq:DSPPot}
-\end{equation}
-The derivative of this Shifted-Potential can be taken to obtain forces
-for use in MD,
+\end{equation},
+with associated forces,
 \begin{equation}
-F^{\textrm{DSP}}(r_{ij}) = q_iq_j\left(-\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}-\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right) \quad r_{ij}\leqslant R_\textrm{c}.
+f_{\textrm{DSP}}(r) = q_iq_j\left(-\frac{\textrm{erfc}(\alpha r)}{r^2}-\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r^2)}}{r}\right) \quad r\leqslant R_\textrm{c}.
 \label{eq:DSPForces}
 \end{equation}
-Again, this Shifted-Potential suffers from a discontinuity in the
-forces, and a lack of an image-charge component in the forces.  To
-remedy these concerns, a Shifted-Force variant is obtained by
-inclusion of the derivative term in equation \ref{eq:shiftingForm} to
-give,
+Again, this damped shifted potential suffers from a discontinuity and
+a lack of the image charges in the forces.  To remedy these concerns,
+one may derive a Shifted-Force variant by including  the derivative
+term in equation \ref{eq:shiftingForm},
 \begin{equation}
 \begin{split}
-V^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\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_{ij}-R_\mathrm{c}\right)\ \right] \quad r_{ij}\leqslant R_\textrm{c}.
+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 gives the following forces,
 \begin{equation}
 \begin{split}
-F^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&-\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right) \\ &\left.+\left(\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{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_{ij}\leqslant R_\textrm{c}.
+f_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&-\left(\frac{\textrm{erfc}(\alpha r)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r^2)}}{r}\right) \\ &\left.+\left(\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{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}
@@ -285,23 +311,36 @@ from equation \ref{eq:shiftingForm} is equal to equati
 This new Shifted-Force 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 equation \ref{eq:shiftingForm} is equal to equation
-\ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$.  This
-term is not present in the Zahn potential, resulting in a
-discontinuity as particles cross $R_\textrm{c}$.  Second, the sign of
-the derivative portion is different.  The constant $v_\textrm{c}$ term
-is not going to have a presence in the forces after performing the
-derivative, but the negative sign does effect the derivative.  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}$.  Thus, these alterations make
-for an electrostatic summation method that is continuous in both the
-potential and forces and incorporates the pairwise sum considerations
-stressed by Wolf \textit{et al.}\cite{Wolf99}
+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 would be 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 that is continuous in both the
+potential and forces 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}
+
+Reaction Field blah
+
+Group-based cutoff blah
+
+
 \section{Methods}
 
-\subsection{What Qualities are Important?}\label{sec:Qualities}
 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
@@ -310,39 +349,39 @@ the origins of this method, the Canonical ensemble acc
 
 In MC, the potential energy difference between two subsequent
 configurations dictates the progression of MC sampling.  Going back to
-the origins of this method, the Canonical ensemble acceptance criteria
-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 a consistent $\Delta E$ when using
-an alternate method for handling the long-range electrostatics ensures
-proper sampling within the ensemble.
+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 directs how the system will
+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 it develops as a whole.  If the
-magnitude and direction of these vectors are similar when using
-alternate electrostatic summation techniques, the dynamics in the near
-term will be indistinguishable.  Because error in MD calculations is
-cumulative, one should expect greater deviation in the long term
-trajectories with greater differences in these vectors between
-configurations using different long-range electrostatics.
+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}
-Evaluation of the pairwise summation techniques (outlined in section
-\ref{sec:ESMethods}) for use in MC simulations was performed through
-study of the energy differences between conformations.  Considering
-the SPME results to be the correct or desired behavior, ideal
-performance of a tested method was taken to be agreement between the
-energy differences calculated.  Linear least squares regression of the
-$\Delta E$ values between configurations using SPME against $\Delta E$
-values using tested methods provides a quantitative comparison of this
-agreement.  Unitary results for both the correlation and correlation
-coefficient for these regressions indicate equivalent energetic
-results between the methods.  The correlation is the slope of the
-plotted data while the correlation coefficient ($R^2$) is a measure of
-the of the data scatter around the fitted line and tells about the
-quality of the fit (Fig. \ref{fig:linearFit}).
+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.  Since none of the methods
+provide exact energy differences, we used linear least squares
+regressions of the $\Delta E$ values between configurations using SPME
+against $\Delta E$ values using tested methods provides a quantitative
+comparison of this agreement.  Unitary results for both the
+correlation and correlation coefficient for these regressions indicate
+equivalent energetic results between the method under consideration
+and electrostatics handled using SPME.  Sample correlation plots for
+two alternate methods are shown in Fig. \ref{fig:linearFit}.
 
 \begin{figure}
 \centering
@@ -351,38 +390,47 @@ Each system type (detailed in section \ref{sec:RepSims
 \label{fig:linearFit}
 \end{figure}
 
-Each system type (detailed in section \ref{sec:RepSims}) studied
-consisted of 500 independent configurations, each equilibrated from
-higher temperature trajectories. Thus, 124,750 $\Delta E$ data points
-are used in a regression of a single system type.  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}.
+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 alternate
+(non-Ewald) electrostatic summation methods was evaluated using
+873,250 configurational energy differences.
 
-\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods}
-Evaluation of the pairwise methods (outlined in section
-\ref{sec:ESMethods}) for use in MD simulations was performed through
-comparison of the force and torque vectors obtained with those from
-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 analysis can be performed as
-described previously for comparing $\Delta E$ values. Instead of a
-single value between two system configurations, there is a value for
-each particle 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 samples for each system type.
+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}.
 
-The force and torque vector directions were investigated through
-measurement of the angle ($\theta$) formed between those from the
-particular method and those from SPME
+\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 = \frac{\vec{F}_\textrm{SPME}}{|\vec{F}_\textrm{SPME}|}\cdot\frac{\vec{F}_\textrm{Method}}{|\vec{F}_\textrm{Method}|}.
+\theta_f = \cos^{-1} \hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method},
 \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, weighted by the area on the unit sphere.  Non-linear fits
-were used to measure the shape of the resulting distributions.
+function, weighted by the area on the unit sphere.  Non-linear
+Gaussian fits were used to measure the width of the resulting
+distributions.
 
 \begin{figure}
 \centering
@@ -395,14 +443,15 @@ between two different electrostatic summation methods,
 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
-particular reason for the profile to adhere to a specific shape.
-Because of this and the Gaussian profile's more statistically
-meaningful properties, 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 SPME.
+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}
 
 \subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
 Evaluation of the long-time dynamics of charged systems was performed
@@ -752,7 +801,7 @@ V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\
 clearly, we show how damping strength 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_{ij}})}{r_{ij}}\right]S(r),
+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
