--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/20 15:43:13	2641
+++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/20 17:32:33	2643
@@ -65,54 +65,54 @@ In molecular simulations, proper accumulation of the e
 \section{Introduction}
 
 In molecular simulations, proper accumulation of the electrostatic
-interactions is considered one of the most essential and
-computationally demanding tasks.  The common molecular mechanics force
-fields are founded on representation of the atomic sites centered on
-full or partial charges shielded by Lennard-Jones type interactions.
-This means that nearly every pair interaction involves an
-charge-charge calculation.  Coupled with $r^{-1}$ decay, the monopole
-interactions quickly become a burden for molecular systems of all
-sizes.  For example, in small systems, the electrostatic pair
-interaction may not have decayed appreciably within the box length
-leading to an effect excluded from the pair interactions within a unit
-box.  In large systems, excessively large cutoffs need to be used to
-accurately incorporate their effect, and since the computational cost
-increases proportionally with the cutoff sphere, it quickly becomes
-very time-consuming to perform these calculations.
+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 this issue of both proper and
-practical handling of electrostatic interactions, and these have
-resulted in the availability of 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
+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 incorporate dynamic 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 this cost, there has been some question of the inherent
-periodicity of the explicit Ewald summation artificially influencing
-systems dynamics.\cite{Tobias01}
+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 the common mixed and explicit methods of
-reaction filed and smooth particle mesh
-Ewald\cite{Onsager36,Essmann99} and a new set of shifted methods
-devised by Wolf {\it et al.} which we further extend.\cite{Wolf99}
-These new methods for handling electrostatics are quite
-computationally efficient, since they involve only a simple
-modification to the direct pairwise sum, and they lack the added
-periodicity of the Ewald sum. Below, these methods are evaluated using
-a variety of model systems and comparison methodologies to establish
+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 simulation box, as well as those in the
-periodic replicas,
+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}
@@ -123,12 +123,13 @@ the cell length, $\bm{\Omega}_{i,j}$ are the Euler ang
 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 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 discontinuous for non-neutral systems.
+$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.
 
-This electrostatic summation problem was originally studied by Ewald
+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
@@ -140,23 +141,27 @@ where $\alpha$ is a damping parameter, or separation c
 \label{eq:EwaldSum}
 \end{equation}
 where $\alpha$ is a damping parameter, or separation constant, with
-units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and equal
-$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric
-constant of the encompassing medium. The final two terms of
+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 and this dipole moment gets magnified through
-replication of the periodic images.\cite{deLeeuw80,Smith81} If this
-term is taken to be zero, the system is using conducting boundary
+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 size of systems, the entire
-simulation box was replicated to convergence.  Currently, we balance a
-spherical real-space cutoff with the reciprocal sum and consider the
-surrounding dielectric.
+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 vanishingly
+small compared to the real-space portion.\cite{XXX}
+
 \begin{figure}
 \centering
 \includegraphics[width = \linewidth]{./ewaldProgression.pdf}
@@ -170,65 +175,63 @@ a surrounding dielectric term is included.}
 \label{fig:ewaldTime}
 \end{figure}
 
-The Ewald summation in the straight-forward form is an
-$\mathscr{O}(N^2)$ algorithm.  The separation constant $(\alpha)$
-plays an important role in the computational cost balance between the
-direct and reciprocal-space portions of the summation.  The choice of
-the magnitude 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 appropriate choice of $\alpha$
-and thoughtful algorithm development, this cost can be brought down to
-$\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to
-reduce the cost of the Ewald summation further is to set $\alpha$ such
-that the real-space interactions decay rapidly, allowing for a short
-spherical cutoff, and then optimize the reciprocal space summation.
-These optimizations usually involve the utilization of the fast
-Fourier transform (FFT),\cite{Hockney81} leading to the
+The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm.  The
+separation constant $(\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 from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N \log N)$.
+summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N
+\log N)$.
 
-These developments and optimizations have led the use of the Ewald
-summation to become routine in simulations with periodic boundary
-conditions. However, in certain systems the intrinsic three
-dimensional periodicity can prove to be problematic, such as two
-dimensional surfaces and membranes.  The Ewald sum has been
-reformulated to handle 2D
-systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the new
-methods have been found to be computationally
-expensive.\cite{Spohr97,Yeh99} Inclusion of a correction term in the
-full Ewald summation is a possible direction for enabling the handling
-of 2D systems and the inclusion of the optimizations described
-previously.\cite{Yeh99}
+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 systems that have the same
-dimensionality.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
-Good examples are solvated proteins kept at high relative
-concentration due to the periodicity of the electrostatics.  In these
+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 ought to be taken when
-considering the use of the Ewald summation where the intrinsic
-periodicity may negatively affect the system dynamics.
+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 and 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
+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}
@@ -255,22 +258,20 @@ the potential are not commensurate.  Attempts to use b
 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
+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}).
+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
-simulations.  Taking the integral of the forces shown in equation
-\ref{eq:WolfForces}, they proposed a new damped Coulomb
-potential,
+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\}.
+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
+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.
 
@@ -301,7 +302,7 @@ shifted potential,
 \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},
@@ -309,7 +310,7 @@ and shifted force,
 \end{equation}
 and shifted force,
 \begin{equation}
-v_\textrm{SF}(r) = 	\begin{cases}
+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},
@@ -325,16 +326,16 @@ of the unshifted potential itself (when inside the cut
 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),
+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
+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,
+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}
@@ -342,14 +343,14 @@ al.}'s undamped prescription:
 then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et
 al.}'s undamped prescription:
 \begin{equation}
-v_\textrm{SP}(r) =
+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}.
+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
@@ -364,18 +365,18 @@ will give,
 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}.
+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}.
+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 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
+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
@@ -383,13 +384,13 @@ Wolf \textit{et al.} originally discussed the energeti
 to gain functionality in dynamics simulations.
 
 Wolf \textit{et al.} originally discussed the energetics of the
-shifted Coulomb potential (Eq. \ref{eq:SPPot}), and they found that
-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
+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 functional forms could be
+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
@@ -398,37 +399,37 @@ v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
 v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
 \label{eq:dampCoulomb}
 \end{equation}
-the shifted potential (Eq. (\ref{eq:SPPot})) can be reacquired using
-eq. (\ref{eq:shiftingForm}),
+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},
+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}.
+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 discontinuity and
-a lack of the image charges in the forces.  To remedy these concerns,
-one may derive a {\sc sf} variant by including  the derivative
-term in eq. (\ref{eq:shiftingForm}),
+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}.
+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}.
+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 chosen to be zero, the undamped
-case, eqs. (\ref{eq:SPPot}-\ref{eq:SFForces}) are correctly recovered
-from eqs. (\ref{eq:DSPPot}-\ref{eq:DSFForces}).
+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
@@ -440,55 +441,48 @@ would be expected to have sudden jumps as particle dis
 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}$.
+$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 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.
+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 will consider some
-other techniques that commonly get used in molecular simulations.  The
+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
-non-neutral molecules, collecting atoms into neutral groups takes
+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 groups, an orientational aspect is
-introduced to the interactions.  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
-switching function,
-\begin{equation}
-S(r) = \begin{cases} 1 &\quad r\leqslant r_\textrm{sw} \\
-\frac{(R_\textrm{c}+2r-3r_\textrm{sw})(R_\textrm{c}-r)^2}{(R_\textrm{c}-r_\textrm{sw})^3} &\quad r_\textrm{sw}<r\leqslant R_\textrm{c} \\
-0 &\quad r>R_\textrm{c}
-\end{cases},
-\end{equation}
-where the above form is for a cubic function.  If a smooth second
-derivative is desired, a fifth (or higher) order polynomial can be
-used.\cite{Andrea83} 
+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 their effect, 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 application, it 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.
+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}
 
