ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/electrostaticMethodsPaper/electrostaticMethods.tex
(Generate patch)

Comparing trunk/electrostaticMethodsPaper/electrostaticMethods.tex (file contents):
Revision 2619 by chrisfen, Wed Mar 15 00:21:39 2006 UTC vs.
Revision 2639 by chrisfen, Mon Mar 20 02:00:26 2006 UTC

# Line 2 | Line 2
2   %\documentclass[aps,prb,preprint]{revtex4}
3   \documentclass[11pt]{article}
4   \usepackage{endfloat}
5 < \usepackage{amsmath}
5 > \usepackage{amsmath,bm}
6   \usepackage{amssymb}
7   \usepackage{epsf}
8   \usepackage{times}
# Line 58 | Line 58 | real-space portion of the lattice summation.
58  
59   %\narrowtext
60  
61 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61 > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62   %                              BODY OF TEXT
63 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64  
65   \section{Introduction}
66  
67   In molecular simulations, proper accumulation of the electrostatic
68   interactions is considered one of the most essential and
69 < computationally demanding tasks.
69 > computationally demanding tasks.  The common molecular mechanics force
70 > fields are founded on representation of the atomic sites centered on
71 > full or partial charges shielded by Lennard-Jones type interactions.
72 > This means that nearly every pair interaction involves an
73 > charge-charge calculation.  Coupled with $r^{-1}$ decay, the monopole
74 > interactions quickly become a burden for molecular systems of all
75 > sizes.  For example, in small systems, the electrostatic pair
76 > interaction may not have decayed appreciably within the box length
77 > leading to an effect excluded from the pair interactions within a unit
78 > box.  In large systems, excessively large cutoffs need to be used to
79 > accurately incorporate their effect, and since the computational cost
80 > increases proportionally with the cutoff sphere, it quickly becomes
81 > very time-consuming to perform these calculations.
82  
83 + There have been many efforts to address this issue of both proper and
84 + practical handling of electrostatic interactions, and these have
85 + resulted in the availability of a variety of
86 + techniques.\cite{Roux99,Sagui99,Tobias01} These are typically
87 + classified as implicit methods (i.e., continuum dielectrics, static
88 + dipolar fields),\cite{Born20,Grossfield00} explicit methods (i.e.,
89 + Ewald summations, interaction shifting or
90 + trucation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e.,
91 + reaction field type methods, fast multipole
92 + methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are
93 + often preferred because they incorporate dynamic solvent molecules in
94 + the system of interest, but these methods are sometimes difficult to
95 + utilize because of their high computational cost.\cite{Roux99} In
96 + addition to this cost, there has been some question of the inherent
97 + periodicity of the explicit Ewald summation artificially influencing
98 + systems dynamics.\cite{Tobias01}
99 +
100 + In this paper, we focus on the common mixed and explicit methods of
101 + reaction filed and smooth particle mesh
102 + Ewald\cite{Onsager36,Essmann99} and a new set of shifted methods
103 + devised by Wolf {\it et al.} which we further extend.\cite{Wolf99}
104 + These new methods for handling electrostatics are quite
105 + computationally efficient, since they involve only a simple
106 + modification to the direct pairwise sum, and they lack the added
107 + periodicity of the Ewald sum. Below, these methods are evaluated using
108 + a variety of model systems and comparison methodologies to establish
109 + their useability in molecular simulations.
110 +
111   \subsection{The Ewald Sum}
112 < blah blah blah Ewald Sum Important blah blah blah
112 > The complete accumulation electrostatic interactions in a system with
113 > periodic boundary conditions (PBC) requires the consideration of the
114 > effect of all charges within a simulation box, as well as those in the
115 > periodic replicas,
116 > \begin{equation}
117 > 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],
118 > \label{eq:PBCSum}
119 > \end{equation}
120 > where the sum over $\mathbf{n}$ is a sum over all periodic box
121 > replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the
122 > prime indicates $i = j$ are neglected for $\mathbf{n} =
123 > 0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic
124 > particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is
125 > the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and
126 > $j$, and $\phi$ is Poisson's equation ($\phi(\mathbf{r}_{ij}) = q_i
127 > q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge interactions). In the
128 > case of monopole electrostatics, eq. (\ref{eq:PBCSum}) is
129 > conditionally convergent and is discontiuous for non-neutral systems.
130  
131 + This electrostatic summation problem was originally studied by Ewald
132 + for the case of an infinite crystal.\cite{Ewald21}. The approach he
133 + took was to convert this conditionally convergent sum into two
134 + absolutely convergent summations: a short-ranged real-space summation
135 + and a long-ranged reciprocal-space summation,
136 + \begin{equation}
137 + \begin{split}
138 + 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,
139 + \end{split}
140 + \label{eq:EwaldSum}
141 + \end{equation}
142 + where $\alpha$ is a damping parameter, or separation constant, with
143 + units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and equal
144 + $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric
145 + constant of the encompassing medium. The final two terms of
146 + eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term
147 + for interacting with a surrounding dielectric.\cite{Allen87} This
148 + dipolar term was neglected in early applications in molecular
149 + simulations,\cite{Brush66,Woodcock71} until it was introduced by de
150 + Leeuw {\it et al.} to address situations where the unit cell has a
151 + dipole moment and this dipole moment gets magnified through
152 + replication of the periodic images.\cite{deLeeuw80,Smith81} If this
153 + term is taken to be zero, the system is using conducting boundary
154 + conditions, $\epsilon_{\rm S} = \infty$. Figure
155 + \ref{fig:ewaldTime} shows how the Ewald sum has been applied over
156 + time.  Initially, due to the small size of systems, the entire
157 + simulation box was replicated to convergence.  Currently, we balance a
158 + spherical real-space cutoff with the reciprocal sum and consider the
159 + surrounding dielectric.
160   \begin{figure}
161   \centering
162   \includegraphics[width = \linewidth]{./ewaldProgression.pdf}
# Line 84 | Line 170 | a surrounding dielectric term is included.}
170   \label{fig:ewaldTime}
171   \end{figure}
172  
173 + The Ewald summation in the straight-forward form is an
174 + $\mathscr{O}(N^2)$ algorithm.  The separation constant $(\alpha)$
175 + plays an important role in the computational cost balance between the
176 + direct and reciprocal-space portions of the summation.  The choice of
177 + the magnitude of this value allows one to select whether the
178 + real-space or reciprocal space portion of the summation is an
179 + $\mathscr{O}(N^2)$ calcualtion (with the other being
180 + $\mathscr{O}(N)$).\cite{Sagui99} With appropriate choice of $\alpha$
181 + and thoughtful algorithm development, this cost can be brought down to
182 + $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to
183 + reduce the cost of the Ewald summation further is to set $\alpha$ such
184 + that the real-space interactions decay rapidly, allowing for a short
185 + spherical cutoff, and then optimize the reciprocal space summation.
186 + These optimizations usually involve the utilization of the fast
187 + Fourier transform (FFT),\cite{Hockney81} leading to the
188 + particle-particle particle-mesh (P3M) and particle mesh Ewald (PME)
189 + methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these
190 + methods, the cost of the reciprocal-space portion of the Ewald
191 + summation is from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N \log N)$.
192 +
193 + These developments and optimizations have led the use of the Ewald
194 + summation to become routine in simulations with periodic boundary
195 + conditions. However, in certain systems the intrinsic three
196 + dimensional periodicity can prove to be problematic, such as two
197 + dimensional surfaces and membranes.  The Ewald sum has been
198 + reformulated to handle 2D
199 + systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the new
200 + methods have been found to be computationally
201 + expensive.\cite{Spohr97,Yeh99} Inclusion of a correction term in the
202 + full Ewald summation is a possible direction for enabling the handling
203 + of 2D systems and the inclusion of the optimizations described
204 + previously.\cite{Yeh99}
205 +
206 + Several studies have recognized that the inherent periodicity in the
207 + Ewald sum can also have an effect on systems that have the same
208 + dimensionality.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
209 + Good examples are solvated proteins kept at high relative
210 + concentration due to the periodicity of the electrostatics.  In these
211 + systems, the more compact folded states of a protein can be
212 + artificially stabilized by the periodic replicas introduced by the
213 + Ewald summation.\cite{Weber00} Thus, care ought to be taken when
214 + considering the use of the Ewald summation where the intrinsic
215 + perodicity may negatively affect the system dynamics.
216 +
217 +
218   \subsection{The Wolf and Zahn Methods}
219   In a recent paper by Wolf \textit{et al.}, a procedure was outlined
220 < for an accurate accumulation of electrostatic interactions in an
221 < efficient pairwise fashion.\cite{Wolf99} Wolf \textit{et al.} observed
222 < that the electrostatic interaction is effectively short-ranged in
223 < condensed phase systems and that neutralization of the charge
224 < contained within the cutoff radius is crucial for potential
225 < stability. They devised a pairwise summation method that ensures
226 < charge neutrality and gives results similar to those obtained with
227 < the Ewald summation.  The resulting shifted Coulomb potential
228 < (Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through
229 < placement on the cutoff sphere and a distance-dependent damping
230 < function (identical to that seen in the real-space portion of the
231 < Ewald sum) to aid convergence
220 > for the accurate accumulation of electrostatic interactions in an
221 > efficient pairwise fashion and lacks the inherent periodicity of the
222 > Ewald summation.\cite{Wolf99} Wolf \textit{et al.} observed that the
223 > electrostatic interaction is effectively short-ranged in condensed
224 > phase systems and that neutralization of the charge contained within
225 > the cutoff radius is crucial for potential stability. They devised a
226 > pairwise summation method that ensures charge neutrality and gives
227 > results similar to those obtained with the Ewald summation.  The
228 > resulting shifted Coulomb potential (Eq. \ref{eq:WolfPot}) includes
229 > image-charges subtracted out through placement on the cutoff sphere
230 > and a distance-dependent damping function (identical to that seen in
231 > the real-space portion of the Ewald sum) to aid convergence
232   \begin{equation}
233 < 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\}.
233 > 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\}.
234   \label{eq:WolfPot}
235   \end{equation}
236   Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted
237   potential.  However, neutralizing the charge contained within each
238   cutoff sphere requires the placement of a self-image charge on the
239   surface of the cutoff sphere.  This additional self-term in the total
240 < potential enables Wolf {\it et al.}  to obtain excellent estimates of
240 > potential enabled Wolf {\it et al.}  to obtain excellent estimates of
241   Madelung energies for many crystals.
242  
243   In order to use their charge-neutralized potential in molecular
# Line 114 | Line 245 | procedure gives an expression for the forces,
245   derivative of this potential prior to evaluation of the limit.  This
246   procedure gives an expression for the forces,
247   \begin{equation}
248 < 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\},
248 > 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\},
249   \label{eq:WolfForces}
250   \end{equation}
251   that incorporates both image charges and damping of the electrostatic
# Line 122 | Line 253 | force expressions for use in simulations involving wat
253  
254   More recently, Zahn \textit{et al.} investigated these potential and
255   force expressions for use in simulations involving water.\cite{Zahn02}
256 < In their work, they pointed out that the method that the forces and
257 < derivative of the potential are not commensurate.  Attempts to use
258 < both Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will
259 < lead to poor energy conservation.  They correctly observed that taking
260 < the limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating
261 < the derivatives is mathematically invalid.
256 > In their work, they pointed out that the forces and derivative of
257 > the potential are not commensurate.  Attempts to use both
258 > Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead
259 > to poor energy conservation.  They correctly observed that taking the
260 > limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the
261 > derivatives gives forces for a different potential energy function
262 > than the one shown in Eq. (\ref{eq:WolfPot}).
263  
264   Zahn \textit{et al.} proposed a modified form of this ``Wolf summation
265   method'' as a way to use this technique in Molecular Dynamics
# Line 135 | Line 267 | potential,
267   \ref{eq:WolfForces}, they proposed a new damped Coulomb
268   potential,
269   \begin{equation}
270 < 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\}.
270 > 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\}.
271   \label{eq:ZahnPot}
272   \end{equation}
273   They showed that this potential does fairly well at capturing the
# Line 146 | Line 278 | al.} are constructed using two different (and separabl
278  
279   The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et
280   al.} are constructed using two different (and separable) computational
281 < tricks: \begin{itemize}
281 > tricks: \begin{enumerate}
282   \item shifting through the use of image charges, and
283   \item damping the electrostatic interaction.
284 < \end{itemize}  Wolf \textit{et al.} treated the
284 > \end{enumerate}  Wolf \textit{et al.} treated the
285   development of their summation method as a progressive application of
286   these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded
287   their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the
# Line 169 | Line 301 | shifted potential,
301   \textit{et al.}  and Zahn \textit{et al.} by considering the standard
302   shifted potential,
303   \begin{equation}
304 < v^\textrm{SP}(r) =      \begin{cases}
304 > v_\textrm{SP}(r) =      \begin{cases}
305   v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r >
306   R_\textrm{c}  
307   \end{cases},
# Line 177 | Line 309 | and shifted force,
309   \end{equation}
310   and shifted force,
311   \begin{equation}
312 < v^\textrm{SF}(r) =      \begin{cases}
313 < v(r)-v_\textrm{c}-\left(\frac{\textrm{d}v(r)}{\textrm{d}r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c})
312 > v_\textrm{SF}(r) =      \begin{cases}
313 > v(r)-v_\textrm{c}-\left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c})
314   &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c}
315                                                  \end{cases},
316   \label{eq:shiftingForm}
# Line 190 | Line 322 | potential is smooth at the cutoff radius
322   potential is smooth at the cutoff radius
323   ($R_\textrm{c}$).\cite{Allen87}
324  
325 <
326 <
195 <
196 < 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}
325 > The forces associated with the shifted potential are simply the forces
326 > of the unshifted potential itself (when inside the cutoff sphere),
327   \begin{equation}
328 < V^\textrm{WSP}(r_{ij}) =        \begin{cases} q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c}
199 <                                                \end{cases}.
200 < \label{eq:WolfSP}
328 > f_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right),
329   \end{equation}
330 < The forces associated with this potential are obtained by taking the derivative, resulting in the following,
330 > and are zero outside.  Inside the cutoff sphere, the forces associated
331 > with the shifted force form can be written,
332   \begin{equation}
333 < F^\textrm{WSP}(r_{ij}) =        \begin{cases} q_iq_j\left(-\frac{1}{r_{ij}^2}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c}
334 <                                                \end{cases}.
206 < \label{eq:FWolfSP}
333 > f_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d
334 > v(r)}{dr} \right)_{r=R_\textrm{c}}.
335   \end{equation}
336 < 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}$.  
337 <
210 < If the derivative term in equation \ref{eq:shiftingForm} is evaluated, we obtain an hitherto undiscussed shifted force Coulomb potential,
336 >
337 > If the potential ($v(r)$) is taken to be the normal Coulomb potential,
338   \begin{equation}
339 < V^\textrm{SF}(r_{ij}) =         \begin{cases} 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} \\ 0 &\quad r_{ij}>R_\textrm{c}
340 <                                                \end{cases}.
341 < \label{eq:SFPot}
342 < \end{equation}
343 < Taking the derivative of this shifted force potential gives the following forces,
217 < \begin{equation}
218 < F^\textrm{SF}(r_{ij}) =         \begin{cases} q_iq_j\left(-\frac{1}{r_{ij}^2}+\frac{1}{R_\textrm{c}^2}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c}
219 <                                                \end{cases}.
220 < \label{eq:SFForces}
221 < \end{equation}
222 < 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.
223 <
224 < 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,
339 > v(r) = \frac{q_i q_j}{r},
340 > \label{eq:Coulomb}
341 > \end{equation}
342 > then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et
343 > al.}'s undamped prescription:
344   \begin{equation}
345 < v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}},
345 > v_\textrm{SP}(r) =
346 > q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad
347 > r\leqslant R_\textrm{c},
348 > \label{eq:SPPot}
349 > \end{equation}
350 > with associated forces,
351 > \begin{equation}
352 > f_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}.
353 > \label{eq:SPForces}
354 > \end{equation}
355 > These forces are identical to the forces of the standard Coulomb
356 > interaction, and cutting these off at $R_c$ was addressed by Wolf
357 > \textit{et al.} as undesirable.  They pointed out that the effect of
358 > the image charges is neglected in the forces when this form is
359 > used,\cite{Wolf99} thereby eliminating any benefit from the method in
360 > molecular dynamics.  Additionally, there is a discontinuity in the
361 > forces at the cutoff radius which results in energy drift during MD
362 > simulations.
363 >
364 > The shifted force ({\sc sf}) form using the normal Coulomb potential
365 > will give,
366 > \begin{equation}
367 > 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}.
368 > \label{eq:SFPot}
369 > \end{equation}
370 > with associated forces,
371 > \begin{equation}
372 > 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}.
373 > \label{eq:SFForces}
374 > \end{equation}
375 > This formulation has the benefits that there are no discontinuities at
376 > the cutoff distance, while the neutralizing image charges are present
377 > in both the energy and force expressions.  It would be simple to add
378 > the self-neutralizing term back when computing the total energy of the
379 > system, thereby maintaining the agreement with the Madelung energies.
380 > A side effect of this treatment is the alteration in the shape of the
381 > potential that comes from the derivative term.  Thus, a degree of
382 > clarity about agreement with the empirical potential is lost in order
383 > to gain functionality in dynamics simulations.
384 >
385 > Wolf \textit{et al.} originally discussed the energetics of the
386 > shifted Coulomb potential (Eq. \ref{eq:SPPot}), and they found that
387 > it was still insufficient for accurate determination of the energy
388 > with reasonable cutoff distances.  The calculated Madelung energies
389 > fluctuate around the expected value with increasing cutoff radius, but
390 > the oscillations converge toward the correct value.\cite{Wolf99} A
391 > damping function was incorporated to accelerate the convergence; and
392 > though alternative functional forms could be
393 > used,\cite{Jones56,Heyes81} the complimentary error function was
394 > chosen to mirror the effective screening used in the Ewald summation.
395 > Incorporating this error function damping into the simple Coulomb
396 > potential,
397 > \begin{equation}
398 > v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
399   \label{eq:dampCoulomb}
400   \end{equation}
401 < the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived \textit{via} equation \ref{eq:shiftingForm},
401 > the shifted potential (Eq. (\ref{eq:SPPot})) can be reacquired using
402 > eq. (\ref{eq:shiftingForm}),
403   \begin{equation}
404 < V^{\textrm{DSP}}(r_{ij}) = \begin{cases} 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} \\ 0 &\quad r_{ij}>R_\textrm{c}
232 < \end{cases}.
404 > 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},
405   \label{eq:DSPPot}
406   \end{equation}
407 < The derivative of this Shifted-Potential can be taken to obtain forces for use in MD,
407 > with associated forces,
408   \begin{equation}
409 < F^{\textrm{DSP}}(r_{ij}) = \begin{cases} 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} \\ 0 &\quad r_{ij}>R_\textrm{c}
238 < \end{cases}.
409 > 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}.
410   \label{eq:DSPForces}
411   \end{equation}
412 < 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,
412 > Again, this damped shifted potential suffers from a discontinuity and
413 > a lack of the image charges in the forces.  To remedy these concerns,
414 > one may derive a {\sc sf} variant by including  the derivative
415 > term in eq. (\ref{eq:shiftingForm}),
416   \begin{equation}
417 < V^\mathrm{DSF}(r_{ij}) = \begin{cases} q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}}\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\} &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c}
418 < \end{cases}.
417 > \begin{split}
418 > 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}.
419   \label{eq:DSFPot}
420 + \end{split}
421   \end{equation}
422 < The derivative of the above potential gives the following forces,
422 > The derivative of the above potential will lead to the following forces,
423   \begin{equation}
424 < F^\mathrm{DSF}(r_{ij}) = \begin{cases} q_iq_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}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2R_{\textrm{c}}^2)}}{R_{\textrm{c}}}\right]\right\} &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c}
425 < \end{cases}.
424 > \begin{split}
425 > 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}.
426   \label{eq:DSFForces}
427 + \end{split}
428   \end{equation}
429 + If the damping parameter $(\alpha)$ is chosen to be zero, the undamped
430 + case, eqs. (\ref{eq:SPPot}-\ref{eq:SFForces}) are correctly recovered
431 + from eqs. (\ref{eq:DSPPot}-\ref{eq:DSFForces}).
432  
433 < 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}
433 > This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot}
434 > derived by Zahn \textit{et al.}; however, there are two important
435 > differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from
436 > eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb})
437 > with $r$ replaced by $R_\textrm{c}$.  This term is {\it not} present
438 > in the Zahn potential, resulting in a potential discontinuity as
439 > particles cross $R_\textrm{c}$.  Second, the sign of the derivative
440 > portion is different.  The missing $v_\textrm{c}$ term would not
441 > affect molecular dynamics simulations (although the computed energy
442 > would be expected to have sudden jumps as particle distances crossed
443 > $R_c$).  The sign problem would be a potential source of errors,
444 > however.  In fact, it introduces a discontinuity in the forces at the
445 > cutoff, because the force function is shifted in the wrong direction
446 > and doesn't cross zero at $R_\textrm{c}$.
447  
448 + Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an
449 + electrostatic summation method that is continuous in both the
450 + potential and forces and which incorporates the damping function
451 + proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of this
452 + paper, we will evaluate exactly how good these methods ({\sc sp}, {\sc
453 + sf}, damping) are at reproducing the correct electrostatic summation
454 + performed by the Ewald sum.
455 +
456 + \subsection{Other alternatives}
457 + In addition to the methods described above, we will consider some
458 + other techniques that commonly get used in molecular simulations.  The
459 + simplest of these is group-based cutoffs.  Though of little use for
460 + non-neutral molecules, collecting atoms into neutral groups takes
461 + advantage of the observation that the electrostatic interactions decay
462 + faster than those for monopolar pairs.\cite{Steinbach94} When
463 + considering these molecules as groups, an orientational aspect is
464 + introduced to the interactions.  Consequently, as these molecular
465 + particles move through $R_\textrm{c}$, the energy will drift upward
466 + due to the anisotropy of the net molecular dipole
467 + interactions.\cite{Rahman71} To maintain good energy conservation,
468 + both the potential and derivative need to be smoothly switched to zero
469 + at $R_\textrm{c}$.\cite{Adams79} This is accomplished using a
470 + switching function,
471 + \begin{equation}
472 + S(r) = \begin{cases} 1 &\quad r\leqslant r_\textrm{sw} \\
473 + \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} \\
474 + 0 &\quad r>R_\textrm{c}
475 + \end{cases},
476 + \end{equation}
477 + where the above form is for a cubic function.  If a smooth second
478 + derivative is desired, a fifth (or higher) order polynomial can be
479 + used.\cite{Andrea83}
480 +
481 + Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$,
482 + and to incorporate their effect, a method like Reaction Field ({\sc
483 + rf}) can be used.  The original theory for {\sc rf} was originally
484 + developed by Onsager,\cite{Onsager36} and it was applied in
485 + simulations for the study of water by Barker and Watts.\cite{Barker73}
486 + In application, it is simply an extension of the group-based cutoff
487 + method where the net dipole within the cutoff sphere polarizes an
488 + external dielectric, which reacts back on the central dipole.  The
489 + same switching function considerations for group-based cutoffs need to
490 + made for {\sc rf}, with the additional pre-specification of a
491 + dielectric constant.
492 +
493   \section{Methods}
494  
495 < \subsection{What Qualities are Important?}\label{sec:Qualities}
496 < 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.  
495 > In classical molecular mechanics simulations, there are two primary
496 > techniques utilized to obtain information about the system of
497 > interest: Monte Carlo (MC) and Molecular Dynamics (MD).  Both of these
498 > techniques utilize pairwise summations of interactions between
499 > particle sites, but they use these summations in different ways.
500  
501 < 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.
501 > In MC, the potential energy difference between two subsequent
502 > configurations dictates the progression of MC sampling.  Going back to
503 > the origins of this method, the acceptance criterion for the canonical
504 > ensemble laid out by Metropolis \textit{et al.} states that a
505 > subsequent configuration is accepted if $\Delta E < 0$ or if $\xi <
506 > \exp(-\Delta E/kT)$, where $\xi$ is a random number between 0 and
507 > 1.\cite{Metropolis53} Maintaining the correct $\Delta E$ when using an
508 > alternate method for handling the long-range electrostatics will
509 > ensure proper sampling from the ensemble.
510  
511 < In MD, the derivative of the potential directs 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.
511 > In MD, the derivative of the potential governs how the system will
512 > progress in time.  Consequently, the force and torque vectors on each
513 > body in the system dictate how the system evolves.  If the magnitude
514 > and direction of these vectors are similar when using alternate
515 > electrostatic summation techniques, the dynamics in the short term
516 > will be indistinguishable.  Because error in MD calculations is
517 > cumulative, one should expect greater deviation at longer times,
518 > although methods which have large differences in the force and torque
519 > vectors will diverge from each other more rapidly.
520  
521   \subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
522 < 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}).
522 > The pairwise summation techniques (outlined in section
523 > \ref{sec:ESMethods}) were evaluated for use in MC simulations by
524 > studying the energy differences between conformations.  We took the
525 > SPME-computed energy difference between two conformations to be the
526 > correct behavior. An ideal performance by an alternative method would
527 > reproduce these energy differences exactly.  Since none of the methods
528 > provide exact energy differences, we used linear least squares
529 > regressions of the $\Delta E$ values between configurations using SPME
530 > against $\Delta E$ values using tested methods provides a quantitative
531 > comparison of this agreement.  Unitary results for both the
532 > correlation and correlation coefficient for these regressions indicate
533 > equivalent energetic results between the method under consideration
534 > and electrostatics handled using SPME.  Sample correlation plots for
535 > two alternate methods are shown in Fig. \ref{fig:linearFit}.
536  
537   \begin{figure}
538   \centering
# Line 272 | Line 541 | Evaluation of the pairwise summation techniques (outli
541   \label{fig:linearFit}
542   \end{figure}
543  
544 < Each system type (detailed in section \ref{sec:RepSims}) studied
545 < consisted of 500 independent configurations, each equilibrated from
546 < higher temperature trajectories. Thus, 124,750 $\Delta E$ data points
547 < are used in a regression of a single system type.  Results and
548 < discussion for the individual analysis of each of the system types
280 < appear in the supporting information, while the cumulative results
281 < over all the investigated systems appears below in section
282 < \ref{sec:EnergyResults}.
544 > Each system type (detailed in section \ref{sec:RepSims}) was
545 > represented using 500 independent configurations.  Additionally, we
546 > used seven different system types, so each of the alternate
547 > (non-Ewald) electrostatic summation methods was evaluated using
548 > 873,250 configurational energy differences.
549  
550 < \subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods}
551 < 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.
550 > Results and discussion for the individual analysis of each of the
551 > system types appear in the supporting information, while the
552 > cumulative results over all the investigated systems appears below in
553 > section \ref{sec:EnergyResults}.
554  
555 < 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
555 > \subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods}
556 > We evaluated the pairwise methods (outlined in section
557 > \ref{sec:ESMethods}) for use in MD simulations by
558 > comparing the force and torque vectors with those obtained using the
559 > reference Ewald summation (SPME).  Both the magnitude and the
560 > direction of these vectors on each of the bodies in the system were
561 > analyzed.  For the magnitude of these vectors, linear least squares
562 > regression analyses were performed as described previously for
563 > comparing $\Delta E$ values.  Instead of a single energy difference
564 > between two system configurations, we compared the magnitudes of the
565 > forces (and torques) on each molecule in each configuration.  For a
566 > system of 1000 water molecules and 40 ions, there are 1040 force
567 > vectors and 1000 torque vectors.  With 500 configurations, this
568 > results in 520,000 force and 500,000 torque vector comparisons.
569 > Additionally, data from seven different system types was aggregated
570 > before the comparison was made.
571 >
572 > The {\it directionality} of the force and torque vectors was
573 > investigated through measurement of the angle ($\theta$) formed
574 > between those computed from the particular method and those from SPME,
575   \begin{equation}
576 < \theta_F = \frac{\vec{F}_\textrm{SPME}}{|\vec{F}_\textrm{SPME}|}\cdot\frac{\vec{F}_\textrm{Method}}{|\vec{F}_\textrm{Method}|}.
576 > \theta_f = \cos^{-1} \left(\hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}\right),
577   \end{equation}
578 < 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.
578 > where $\hat{f}_\textrm{M}$ is the unit vector pointing along the
579 > force vector computed using method $M$.  
580  
581 + Each of these $\theta$ values was accumulated in a distribution
582 + function, weighted by the area on the unit sphere.  Non-linear
583 + Gaussian fits were used to measure the width of the resulting
584 + distributions.
585 +
586   \begin{figure}
587   \centering
588   \includegraphics[width = \linewidth]{./gaussFit.pdf}
# Line 297 | Line 590 | Each of these $\theta$ values was accumulated in a dis
590   \label{fig:gaussian}
591   \end{figure}
592  
593 < 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 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.
593 > Figure \ref{fig:gaussian} shows an example distribution with applied
594 > non-linear fits.  The solid line is a Gaussian profile, while the
595 > dotted line is a Voigt profile, a convolution of a Gaussian and a
596 > Lorentzian.  Since this distribution is a measure of angular error
597 > between two different electrostatic summation methods, there is no
598 > {\it a priori} reason for the profile to adhere to any specific shape.
599 > Gaussian fits was used to compare all the tested methods.  The
600 > variance ($\sigma^2$) was extracted from each of these fits and was
601 > used to compare distribution widths.  Values of $\sigma^2$ near zero
602 > indicate vector directions indistinguishable from those calculated
603 > when using the reference method (SPME).
604  
605 < \subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
606 < Evaluation of the long-time dynamics of charged systems was performed by considering the NaCl crystal system while using a subset of the best performing pairwise methods.  The NaCl crystal was chosen to avoid possible complications involving the propagation techniques of orientational motion in molecular systems.  To enhance the atomic motion, these crystals were equilibrated at 1000 K, near the experimental $T_m$ for NaCl.  Simulations were performed under the microcanonical ensemble, and velocity autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each of the trajectories,
605 > \subsection{Short-time Dynamics}
606 > Evaluation of the short-time dynamics of charged systems was performed
607 > by considering the 1000 K NaCl crystal system while using a subset of the
608 > best performing pairwise methods.  The NaCl crystal was chosen to
609 > avoid possible complications involving the propagation techniques of
610 > orientational motion in molecular systems.  All systems were started
611 > with the same initial positions and velocities.  Simulations were
612 > performed under the microcanonical ensemble, and velocity
613 > autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each
614 > of the trajectories,
615   \begin{equation}
616 < C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle.
616 > C_v(t) = \frac{\langle v_i(0)\cdot v_i(t)\rangle}{\langle v_i(0)\cdot v_i(0)\rangle}.
617   \label{eq:vCorr}
618   \end{equation}
619 < Velocity autocorrelation functions require detailed short time data and long trajectories for good statistics, thus velocity information was saved every 5 fs over 100 ps trajectories.  The power spectrum ($I(\omega)$) is obtained via Fourier transform of the autocorrelation function
619 > Velocity autocorrelation functions require detailed short time data,
620 > thus velocity information was saved every 2 fs over 10 ps
621 > trajectories. Because the NaCl crystal is composed of two different
622 > atom types, the average of the two resulting velocity autocorrelation
623 > functions was used for comparisons.
624 >
625 > \subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
626 > Evaluation of the long-time dynamics of charged systems was performed
627 > by considering the NaCl crystal system, again while using a subset of
628 > the best performing pairwise methods.  To enhance the atomic motion,
629 > these crystals were equilibrated at 1000 K, near the experimental
630 > $T_m$ for NaCl.  Simulations were performed under the microcanonical
631 > ensemble, and velocity information was saved every 5 fs over 100 ps
632 > trajectories.  The power spectrum ($I(\omega)$) was obtained via
633 > Fourier transform of the velocity autocorrelation function
634   \begin{equation}
635   I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt,
636   \label{eq:powerSpec}
637   \end{equation}
638 < where the frequency, $\omega=0,\ 1,\ ...,\ N-1$.
638 > where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the
639 > NaCl crystal is composed of two different atom types, the average of
640 > the two resulting power spectra was used for comparisons.
641  
642   \subsection{Representative Simulations}\label{sec:RepSims}
643 < A variety of common and representative simulations were analyzed to determine the relative effectiveness of the pairwise summation techniques in reproducing the energetics and dynamics exhibited by SPME.  The studied systems were as follows:
643 > A variety of common and representative simulations were analyzed to
644 > determine the relative effectiveness of the pairwise summation
645 > techniques in reproducing the energetics and dynamics exhibited by
646 > SPME.  The studied systems were as follows:
647   \begin{enumerate}
648   \item Liquid Water
649   \item Crystalline Water (Ice I$_\textrm{c}$)
# Line 323 | Line 653 | A variety of common and representative simulations wer
653   \item High Ionic Strength Solution of NaCl in Water
654   \item 6 \AA\  Radius Sphere of Argon in Water
655   \end{enumerate}
656 < 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 can comment on possible system dependence and/or universal applicability of the techniques.
656 > By utilizing the pairwise techniques (outlined in section
657 > \ref{sec:ESMethods}) in systems composed entirely of neutral groups,
658 > charged particles, and mixtures of the two, we can comment on possible
659 > system dependence and/or universal applicability of the techniques.
660  
661 < Generation of the system configurations was dependent on the system type.  For the solid and liquid water configurations, configuration snapshots were taken at regular intervals from higher temperature 1000 SPC/E water molecule trajectories and each equilibrated individually.  The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- ions and were selected and equilibrated in the same fashion as the water systems.  For the low and high ionic strength NaCl solutions, 4 and 40 ions were first solvated in a 1000 water molecule boxes respectively.  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}).
661 > Generation of the system configurations was dependent on the system
662 > type.  For the solid and liquid water configurations, configuration
663 > snapshots were taken at regular intervals from higher temperature 1000
664 > SPC/E water molecule trajectories and each equilibrated individually.
665 > The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl-
666 > ions and were selected and equilibrated in the same fashion as the
667 > water systems.  For the low and high ionic strength NaCl solutions, 4
668 > and 40 ions were first solvated in a 1000 water molecule boxes
669 > respectively.  Ion and water positions were then randomly swapped, and
670 > the resulting configurations were again equilibrated individually.
671 > Finally, for the Argon/Water "charge void" systems, the identities of
672 > all the SPC/E waters within 6 \AA\ of the center of the equilibrated
673 > water configurations were converted to argon
674 > (Fig. \ref{fig:argonSlice}).
675  
676   \begin{figure}
677   \centering
# Line 335 | Line 681 | Generation of the system configurations was dependent
681   \end{figure}
682  
683   \subsection{Electrostatic Summation Methods}\label{sec:ESMethods}
684 < Electrostatic summation method comparisons were performed using SPME, the Shifted-Potential and Shifted-Force methods - both with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak, moderate, and strong damping respectively), reaction field with an infinite dielectric constant, and an unmodified cutoff.  Group-based cutoffs with a fifth-order polynomial switching function were necessary for the reaction field simulations and were utilized in the SP, SF, and pure cutoff methods for comparison to the standard lack of group-based cutoffs with a hard truncation.  The SPME calculations 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}
684 > Electrostatic summation method comparisons were performed using SPME,
685 > the {\sc sp} and {\sc sf} methods - both with damping
686 > parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak,
687 > moderate, and strong damping respectively), reaction field with an
688 > infinite dielectric constant, and an unmodified cutoff.  Group-based
689 > cutoffs with a fifth-order polynomial switching function were
690 > necessary for the reaction field simulations and were utilized in the
691 > SP, SF, and pure cutoff methods for comparison to the standard lack of
692 > group-based cutoffs with a hard truncation.  The SPME calculations
693 > were performed using the TINKER implementation of SPME,\cite{Ponder87}
694 > while all other method calculations were performed using the OOPSE
695 > molecular mechanics package.\cite{Meineke05}
696  
697 < These methods were additionally evaluated with three different cutoff radii (9, 12, and 15 \AA) to investigate possible cutoff radius dependence.  It should be noted that the damping parameter chosen in SPME, or so called ``Ewald Coefficient", has a significant effect on the energies and forces calculated.  Typical molecular mechanics packages default 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 in the real-space portion of the summation.\cite{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.
697 > These methods were additionally evaluated with three different cutoff
698 > radii (9, 12, and 15 \AA) to investigate possible cutoff radius
699 > dependence.  It should be noted that the damping parameter chosen in
700 > SPME, or so called ``Ewald Coefficient", has a significant effect on
701 > the energies and forces calculated.  Typical molecular mechanics
702 > packages default this to a value dependent on the cutoff radius and a
703 > tolerance (typically less than $1 \times 10^{-4}$ kcal/mol).  Smaller
704 > tolerances are typically associated with increased accuracy, but this
705 > usually means more time spent calculating the reciprocal-space portion
706 > of the summation.\cite{Perram88,Essmann95} The default TINKER
707 > tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME
708 > calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and
709 > 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively.
710  
711   \section{Results and Discussion}
712  
713   \subsection{Configuration Energy Differences}\label{sec:EnergyResults}
714 < 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}.  
714 > In order to evaluate the performance of the pairwise electrostatic
715 > summation methods for Monte Carlo simulations, the energy differences
716 > between configurations were compared to the values obtained when using
717 > SPME.  The results for the subsequent regression analysis are shown in
718 > figure \ref{fig:delE}.
719  
720   \begin{figure}
721   \centering
# Line 351 | Line 724 | In order to evaluate the performance of the pairwise e
724   \label{fig:delE}
725   \end{figure}
726  
727 < In this figure, it is apparent that it is unreasonable to expect realistic results using an unmodified cutoff.  This is not all that surprising since this results in large energy fluctuations as atoms move in and out of the cutoff radius.  These fluctuations can be alleviated to some degree by using group based cutoffs with a switching function.\cite{Steinbach94}  The Group Switch Cutoff row doesn't show a significant improvement in this plot because the salt and salt solution systems contain non-neutral groups, see the accompanying supporting information for a comparison where all groups are neutral.  
727 > In this figure, it is apparent that it is unreasonable to expect
728 > realistic results using an unmodified cutoff.  This is not all that
729 > surprising since this results in large energy fluctuations as atoms
730 > move in and out of the cutoff radius.  These fluctuations can be
731 > alleviated to some degree by using group based cutoffs with a
732 > switching function.\cite{Steinbach94} The Group Switch Cutoff row
733 > doesn't show a significant improvement in this plot because the salt
734 > and salt solution systems contain non-neutral groups, see the
735 > accompanying supporting information for a comparison where all groups
736 > are neutral.
737  
738 < Correcting the resulting charged cutoff sphere is one of the purposes of the damped Coulomb summation proposed by Wolf \textit{et al.},\cite{Wolf99} and this correction indeed improves the results as seen in the Shifted-Potental rows.  While the undamped case of this method is a significant improvement over the pure cutoff, it still doesn't correlate that well with SPME.  Inclusion of potential damping improves the results, 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 Shifted-Force sets, increasing damping results in progressively poorer correlation.  Overall, the undamped case is the best performing set, as the correlation and quality of fits are consistently superior regardless of the cutoff distance.  This result is beneficial in that the undamped case is less computationally prohibitive do to the lack of complimentary error function calculation when performing the electrostatic pair interaction.  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.
738 > Correcting the resulting charged cutoff sphere is one of the purposes
739 > of the damped Coulomb summation proposed by Wolf \textit{et
740 > al.},\cite{Wolf99} and this correction indeed improves the results as
741 > seen in the Shifted-Potental rows.  While the undamped case of this
742 > method is a significant improvement over the pure cutoff, it still
743 > doesn't correlate that well with SPME.  Inclusion of potential damping
744 > improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows
745 > an excellent correlation and quality of fit with the SPME results,
746 > particularly with a cutoff radius greater than 12 \AA .  Use of a
747 > larger damping parameter is more helpful for the shortest cutoff
748 > shown, but it has a detrimental effect on simulations with larger
749 > cutoffs.  In the {\sc sf} sets, increasing damping results in
750 > progressively poorer correlation.  Overall, the undamped case is the
751 > best performing set, as the correlation and quality of fits are
752 > consistently superior regardless of the cutoff distance.  This result
753 > is beneficial in that the undamped case is less computationally
754 > prohibitive do to the lack of complimentary error function calculation
755 > when performing the electrostatic pair interaction.  The reaction
756 > field results illustrates some of that method's limitations, primarily
757 > that it was developed for use in homogenous systems; although it does
758 > provide results that are an improvement over those from an unmodified
759 > cutoff.
760  
761   \subsection{Magnitudes of the Force and Torque Vectors}
762  
763 < 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.
763 > Evaluation of pairwise methods for use in Molecular Dynamics
764 > simulations requires consideration of effects on the forces and
765 > torques.  Investigation of the force and torque vector magnitudes
766 > provides a measure of the strength of these values relative to SPME.
767 > Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the
768 > force and torque vector magnitude regression results for the
769 > accumulated analysis over all the system types.
770  
771   \begin{figure}
772   \centering
# Line 366 | Line 775 | Evaluation of pairwise methods for use in Molecular Dy
775   \label{fig:frcMag}
776   \end{figure}
777  
778 < 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 Shifted-Potential 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 Shifted-Force 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.
778 > Figure \ref{fig:frcMag}, for the most part, parallels the results seen
779 > in the previous $\Delta E$ section.  The unmodified cutoff results are
780 > poor, but using group based cutoffs and a switching function provides
781 > a improvement much more significant than what was seen with $\Delta
782 > E$.  Looking at the {\sc sp} sets, the slope and $R^2$
783 > improve with the use of damping to an optimal result of 0.2 \AA
784 > $^{-1}$ for the 12 and 15 \AA\ cutoffs.  Further increases in damping,
785 > while beneficial for simulations with a cutoff radius of 9 \AA\ , is
786 > detrimental to simulations with larger cutoff radii.  The undamped
787 > {\sc sf} method gives forces in line with those obtained using
788 > SPME, and use of a damping function results in minor improvement.  The
789 > reaction field results are surprisingly good, considering the poor
790 > quality of the fits for the $\Delta E$ results.  There is still a
791 > considerable degree of scatter in the data, but it correlates well in
792 > general.  To be fair, we again note that the reaction field
793 > calculations do not encompass NaCl crystal and melt systems, so these
794 > results are partly biased towards conditions in which the method
795 > performs more favorably.
796  
797   \begin{figure}
798   \centering
# Line 375 | Line 801 | Figure \ref{fig:frcMag}, for the most part, parallels
801   \label{fig:trqMag}
802   \end{figure}
803  
804 < 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 Shifted-Potential and helpful yet possibly unnecessary with the Shifted-Force 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.
804 > To evaluate the torque vector magnitudes, the data set from which
805 > values are drawn is limited to rigid molecules in the systems
806 > (i.e. water molecules).  In spite of this smaller sampling pool, the
807 > torque vector magnitude results in figure \ref{fig:trqMag} are still
808 > similar to those seen for the forces; however, they more clearly show
809 > the improved behavior that comes with increasing the cutoff radius.
810 > Moderate damping is beneficial to the {\sc sp} and helpful
811 > yet possibly unnecessary with the {\sc sf} method, and they also
812 > show that over-damping adversely effects all cutoff radii rather than
813 > showing an improvement for systems with short cutoffs.  The reaction
814 > field method performs well when calculating the torques, better than
815 > the Shifted Force method over this limited data set.
816  
817   \subsection{Directionality of the Force and Torque Vectors}
818  
819 < 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.  
819 > Having force and torque vectors with magnitudes that are well
820 > correlated to SPME is good, but if they are not pointing in the proper
821 > direction the results will be incorrect.  These vector directions were
822 > investigated through measurement of the angle formed between them and
823 > those from SPME.  The results (Fig. \ref{fig:frcTrqAng}) are compared
824 > through the variance ($\sigma^2$) of the Gaussian fits of the angle
825 > error distributions of the combined set over all system types.
826  
827   \begin{figure}
828   \centering
# Line 388 | Line 831 | Having force and torque vectors with magnitudes that a
831   \label{fig:frcTrqAng}
832   \end{figure}
833  
834 < 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 Shifted-Force, 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 Shifted-Potential and moderately for the Shifted-Force 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.  
834 > Both the force and torque $\sigma^2$ results from the analysis of the
835 > total accumulated system data are tabulated in figure
836 > \ref{fig:frcTrqAng}.  All of the sets, aside from the over-damped case
837 > show the improvement afforded by choosing a longer simulation cutoff.
838 > Increasing the cutoff from 9 to 12 \AA\ typically results in a halving
839 > of the distribution widths, with a similar improvement going from 12
840 > to 15 \AA .  The undamped {\sc sf}, Group Based Cutoff, and
841 > Reaction Field methods all do equivalently well at capturing the
842 > direction of both the force and torque vectors.  Using damping
843 > improves the angular behavior significantly for the {\sc sp}
844 > and moderately for the {\sc sf} methods.  Increasing the damping
845 > too far is destructive for both methods, particularly to the torque
846 > vectors.  Again it is important to recognize that the force vectors
847 > cover all particles in the systems, while torque vectors are only
848 > available for neutral molecular groups.  Damping appears to have a
849 > more beneficial effect on non-neutral bodies, and this observation is
850 > investigated further in the accompanying supporting information.
851  
852   \begin{table}[htbp]
853     \centering
# Line 423 | Line 882 | Both the force and torque $\sigma^2$ results from the
882     \label{tab:groupAngle}
883   \end{table}
884  
885 < Although not discussed previously, group based cutoffs can be applied to both the Shifted-Potential and Shifted-Force 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 Shifted-Potential shows much narrower angular distributions for both the force and torque vectors when using an $\alpha$ of 0.2 \AA$^{-1}$ or less, while Shifted-Force 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.  
885 > Although not discussed previously, group based cutoffs can be applied
886 > to both the {\sc sp} and {\sc sf} methods.  Use off a
887 > switching function corrects for the discontinuities that arise when
888 > atoms of a group exit the cutoff before the group's center of mass.
889 > Though there are no significant benefit or drawbacks observed in
890 > $\Delta E$ and vector magnitude results when doing this, there is a
891 > measurable improvement in the vector angle results.  Table
892 > \ref{tab:groupAngle} shows the angular variance values obtained using
893 > group based cutoffs and a switching function alongside the standard
894 > results seen in figure \ref{fig:frcTrqAng} for comparison purposes.
895 > The {\sc sp} shows much narrower angular distributions for
896 > both the force and torque vectors when using an $\alpha$ of 0.2
897 > \AA$^{-1}$ or less, while {\sc sf} shows improvements in the
898 > undamped and lightly damped cases.  Thus, by calculating the
899 > electrostatic interactions in terms of molecular pairs rather than
900 > atomic pairs, the direction of the force and torque vectors are
901 > determined more accurately.
902  
903 < One additional trend to recognize in table \ref{tab:groupAngle} is that the $\sigma^2$ values for both Shifted-Potential and Shifted-Force 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 Shifted-Force method.
903 > One additional trend to recognize in table \ref{tab:groupAngle} is
904 > that the $\sigma^2$ values for both {\sc sp} and
905 > {\sc sf} converge as $\alpha$ increases, something that is easier
906 > to see when using group based cutoffs.  Looking back on figures
907 > \ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this
908 > behavior clearly at large $\alpha$ and cutoff values.  The reason for
909 > this is that the complimentary error function inserted into the
910 > potential weakens the electrostatic interaction as $\alpha$ increases.
911 > Thus, at larger values of $\alpha$, both the summation method types
912 > progress toward non-interacting functions, so care is required in
913 > choosing large damping functions lest one generate an undesirable loss
914 > in the pair interaction.  Kast \textit{et al.}  developed a method for
915 > choosing appropriate $\alpha$ values for these types of electrostatic
916 > summation methods by fitting to $g(r)$ data, and their methods
917 > indicate optimal values of 0.34, 0.25, and 0.16 \AA$^{-1}$ for cutoff
918 > values of 9, 12, and 15 \AA\ respectively.\cite{Kast03} These appear
919 > to be reasonable choices to obtain proper MC behavior
920 > (Fig. \ref{fig:delE}); however, based on these findings, choices this
921 > high would introduce error in the molecular torques, particularly for
922 > the shorter cutoffs.  Based on the above findings, empirical damping
923 > up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably
924 > unnecessary when using the {\sc sf} method.
925  
926 + \subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}
927 +
928 + In the previous studies using a {\sc sf} variant of the damped
929 + Wolf coulomb potential, the structure and dynamics of water were
930 + investigated rather extensively.\cite{Zahn02,Kast03} Their results
931 + indicated that the damped {\sc sf} method results in properties
932 + very similar to those obtained when using the Ewald summation.
933 + Considering the statistical results shown above, the good performance
934 + of this method is not that surprising.  Rather than consider the same
935 + systems and simply recapitulate their results, we decided to look at
936 + the solid state dynamical behavior obtained using the best performing
937 + summation methods from the above results.
938 +
939 + \begin{figure}
940 + \centering
941 + \includegraphics[width = \linewidth]{./vCorrPlot.pdf}
942 + \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.}
943 + \label{fig:vCorrPlot}
944 + \end{figure}
945 +
946 + The short-time decays through the first collision are nearly identical
947 + in figure \ref{fig:vCorrPlot}, but the peaks and troughs of the
948 + functions show how the methods differ.  The undamped {\sc sf} method
949 + has deeper troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher
950 + peaks than any of the other methods.  As the damping function is
951 + increased, these peaks are smoothed out, and approach the SPME
952 + curve. The damping acts as a distance dependent Gaussian screening of
953 + the point charges for the pairwise summation methods; thus, the
954 + collisions are more elastic in the undamped {\sc sf} potental, and the
955 + stiffness of the potential is diminished as the electrostatic
956 + interactions are softened by the damping function.  With $\alpha$
957 + values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are
958 + nearly identical and track the SPME features quite well.  This is not
959 + too surprising in that the differences between the {\sc sf} and {\sc
960 + sp} potentials are mitigated with increased damping.  However, this
961 + appears to indicate that once damping is utilized, the form of the
962 + potential seems to play a lesser role in the crystal dynamics.
963 +
964   \subsection{Collective Motion: Power Spectra of NaCl Crystals}
965  
966 < In the previous studies using a Shifted-Force 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 Shifted-Force 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.
966 > The short time dynamics were extended to evaluate how the differences
967 > between the methods affect the collective long-time motion.  The same
968 > electrostatic summation methods were used as in the short time
969 > velocity autocorrelation function evaluation, but the trajectories
970 > were sampled over a much longer time. The power spectra of the
971 > resulting velocity autocorrelation functions were calculated and are
972 > displayed in figure \ref{fig:methodPS}.
973  
974   \begin{figure}
975   \centering
976   \includegraphics[width = \linewidth]{./spectraSquare.pdf}
977 < \caption{Power spectra obtained from the velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, Shifted-Force ($\alpha$ = 0, 0.1, \& 0.2), and Shifted-Potential ($\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.}
977 > \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.}
978   \label{fig:methodPS}
979   \end{figure}
980  
981 < Figure \ref{fig:methodPS} shows the power spectra for the NaCl crystals (from averaged Na and Cl ion velocity autocorrelation functions) using the stated electrostatic summation methods.  While high frequency peaks of all the spectra 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 Shifted-Force.  When using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the Shifted-Force and Shifted-Potential methods give near identical correlated motion behavior as the Ewald method (which has a damping value of 0.3119).  The damping acts as a distance dependent Gaussian screening of the point charges for the pairwise summation methods.  This weakening of the electrostatic interaction with distance 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 affects a simple real-space electrostatic potential,
981 > While high frequency peaks of the spectra in this figure overlap,
982 > showing the same general features, the low frequency region shows how
983 > the summation methods differ.  Considering the low-frequency inset
984 > (expanded in the upper frame of figure \ref{fig:dampInc}), at
985 > frequencies below 100 cm$^{-1}$, the correlated motions are
986 > blue-shifted when using undamped or weakly damped {\sc sf}.  When
987 > using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the {\sc sf}
988 > and {\sc sp} methods give near identical correlated motion behavior as
989 > the Ewald method (which has a damping value of 0.3119).  This
990 > weakening of the electrostatic interaction with increased damping
991 > explains why the long-ranged correlated motions are at lower
992 > frequencies for the moderately damped methods than for undamped or
993 > weakly damped methods.  To see this effect more clearly, we show how
994 > damping strength alone affects a simple real-space electrostatic
995 > potential,
996   \begin{equation}
997 < 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),
997 > V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r),
998   \end{equation}
999 < 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 Shifted-Force method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite respectable; however, it appears as though moderate damping is required for accurate reproduction of crystal dynamics.
999 > where $S(r)$ is a switching function that smoothly zeroes the
1000 > potential at the cutoff radius.  Figure \ref{fig:dampInc} shows how
1001 > the low frequency motions are dependent on the damping used in the
1002 > direct electrostatic sum.  As the damping increases, the peaks drop to
1003 > lower frequencies.  Incidentally, use of an $\alpha$ of 0.25
1004 > \AA$^{-1}$ on a simple electrostatic summation results in low
1005 > frequency correlated dynamics equivalent to a simulation using SPME.
1006 > When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks
1007 > shift to higher frequency in exponential fashion.  Though not shown,
1008 > the spectrum for the simple undamped electrostatic potential is
1009 > blue-shifted such that the lowest frequency peak resides near 325
1010 > cm$^{-1}$.  In light of these results, the undamped {\sc sf} method
1011 > producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite
1012 > respectable and shows that the shifted force procedure accounts for
1013 > most of the effect afforded through use of the Ewald summation.
1014 > However, it appears as though moderate damping is required for
1015 > accurate reproduction of crystal dynamics.
1016   \begin{figure}
1017   \centering
1018   \includegraphics[width = \linewidth]{./comboSquare.pdf}
1019 < \caption{Upper: Zoomed inset from figure \ref{fig:methodPS}.  As the damping value for the Shifted-Force potential increases, the low-frequency peaks red-shift.  Lower: Low-frequency correlated motions for NaCl crystals at 1000 K 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.}
1019 > \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.}
1020   \label{fig:dampInc}
1021   \end{figure}
1022  
1023   \section{Conclusions}
1024  
1025 < 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 Shifted-Force method, reformulated above as equation \ref{eq:SFPot}, 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 Shifted-Force and moderately damped Shifted-Potential methods produced results nearly identical to SPME.  Similarly for the dynamic features, the undamped or moderately damped Shifted-Force and moderately damped Shifted-Potential 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.
1025 > This investigation of pairwise electrostatic summation techniques
1026 > shows that there are viable and more computationally efficient
1027 > electrostatic summation techniques than the Ewald summation, chiefly
1028 > methods derived from the damped Coulombic sum originally proposed by
1029 > Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the
1030 > {\sc sf} method, reformulated above as eq. (\ref{eq:DSFPot}),
1031 > shows a remarkable ability to reproduce the energetic and dynamic
1032 > characteristics exhibited by simulations employing lattice summation
1033 > techniques.  The cumulative energy difference results showed the
1034 > undamped {\sc sf} and moderately damped {\sc sp} methods
1035 > produced results nearly identical to SPME.  Similarly for the dynamic
1036 > features, the undamped or moderately damped {\sc sf} and
1037 > moderately damped {\sc sp} methods produce force and torque
1038 > vector magnitude and directions very similar to the expected values.
1039 > These results translate into long-time dynamic behavior equivalent to
1040 > that produced in simulations using SPME.
1041  
1042 < 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.
1042 > Aside from the computational cost benefit, these techniques have
1043 > applicability in situations where the use of the Ewald sum can prove
1044 > problematic.  Primary among them is their use in interfacial systems,
1045 > where the unmodified lattice sum techniques artificially accentuate
1046 > the periodicity of the system in an undesirable manner.  There have
1047 > been alterations to the standard Ewald techniques, via corrections and
1048 > reformulations, to compensate for these systems; but the pairwise
1049 > techniques discussed here require no modifications, making them
1050 > natural tools to tackle these problems.  Additionally, this
1051 > transferability gives them benefits over other pairwise methods, like
1052 > reaction field, because estimations of physical properties (e.g. the
1053 > dielectric constant) are unnecessary.
1054  
1055 < 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 researcher have come to expect
1055 > We are not suggesting any flaw with the Ewald sum; in fact, it is the
1056 > standard by which these simple pairwise sums are judged.  However,
1057 > these results do suggest that in the typical simulations performed
1058 > today, the Ewald summation may no longer be required to obtain the
1059 > level of accuracy most researchers have come to expect
1060  
1061   \section{Acknowledgments}
1062   \newpage

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines