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 2599 by chrisfen, Tue Feb 28 14:09:55 2006 UTC vs.
Revision 2620 by chrisfen, Wed Mar 15 04:34:51 2006 UTC

# Line 1 | Line 1
1   %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
2 < \documentclass[12pt]{article}
2 > %\documentclass[aps,prb,preprint]{revtex4}
3 > \documentclass[11pt]{article}
4   \usepackage{endfloat}
5   \usepackage{amsmath}
6   \usepackage{amssymb}
6 – %\usepackage{ifsym}
7   \usepackage{epsf}
8   \usepackage{times}
9 < \usepackage{mathptm}
9 > \usepackage{mathptmx}
10   \usepackage{setspace}
11   \usepackage{tabularx}
12   \usepackage{graphicx}
13   \usepackage{booktabs}
14 < %\usepackage{berkeley}
14 > \usepackage{bibentry}
15 > \usepackage{mathrsfs}
16   \usepackage[ref]{overcite}
17   \pagestyle{plain}
18   \pagenumbering{arabic}
# Line 24 | Line 25
25  
26   \begin{document}
27  
28 < \title{Is the Ewald Summation necessary in typical molecular simulations: Alternatives to the accepted standard of cutoff policies}
28 > \title{Is the Ewald Summation necessary? Pairwise alternatives to the accepted standard for long-range electrostatics}
29  
30 < \author{Christopher J. Fennell and J. Daniel Gezelter \\
30 > \author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail:
31 > gezelter@nd.edu} \\
32   Department of Chemistry and Biochemistry\\
33   University of Notre Dame\\
34   Notre Dame, Indiana 46556}
# Line 34 | Line 36 | Notre Dame, Indiana 46556}
36   \date{\today}
37  
38   \maketitle
39 < %\doublespacing
39 > \doublespacing
40  
41 + \nobibliography{}
42   \begin{abstract}
43 + A new method for accumulating electrostatic interactions was derived
44 + from the previous efforts described in \bibentry{Wolf99} and
45 + \bibentry{Zahn02} as a possible replacement for lattice sum methods in
46 + molecular simulations.  Comparisons were performed with this and other
47 + pairwise electrostatic summation techniques against the smooth
48 + particle mesh Ewald (SPME) summation to see how well they reproduce
49 + the energetics and dynamics of a variety of simulation types.  The
50 + newly derived Shifted-Force technique shows a remarkable ability to
51 + reproduce the behavior exhibited in simulations using SPME with an
52 + $\mathscr{O}(N)$ computational cost, equivalent to merely the
53 + real-space portion of the lattice summation.
54 +
55   \end{abstract}
56  
57 + \newpage
58 +
59   %\narrowtext
60  
61 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61 > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
62   %                              BODY OF TEXT
63 < %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63 > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
64  
65   \section{Introduction}
66  
67 < In this paper, a variety of simulation situations were analyzed to determine the relative effectiveness of the adapted Wolf spherical truncation schemes at reproducing the results obtained using a smooth particle mesh Ewald (SPME) summation technique.  In addition to the Shifted-Potential and Shifted-Force adapted Wolf methods, both reaction field and uncorrected cutoff methods were included for comparison purposes.  The general usability of these methods in both Monte Carlo and Molecular Dynamics calculations was assessed through statistical analysis over the combined results from all of the following studied systems:
68 < \begin{enumerate}
69 < \item Liquid Water
70 < \item Crystalline Water (Ice I$_\textrm{c}$)
71 < \item NaCl Crystal
72 < \item NaCl Melt
73 < \item Low Ionic Strength Solution of NaCl in Water
74 < \item High Ionic Strength Solution of NaCl in Water
75 < \item 6 \AA\  Radius Sphere of Argon in Water
76 < \end{enumerate}
77 < Additional discussion on the results from the individual systems was also performed to identify limitations of the considered methods in specific systems.
67 > In molecular simulations, proper accumulation of the electrostatic
68 > interactions is considered one of the most essential and
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 an
81 > impractical task to perform these calculations.
82  
83 + \subsection{The Ewald Sum}
84 + blah blah blah Ewald Sum Important blah blah blah
85 +
86 + \begin{figure}
87 + \centering
88 + \includegraphics[width = \linewidth]{./ewaldProgression.pdf}
89 + \caption{How the application of the Ewald summation has changed with
90 + the increase in computer power.  Initially, only small numbers of
91 + particles could be studied, and the Ewald sum acted to replicate the
92 + unit cell charge distribution out to convergence.  Now, much larger
93 + systems of charges are investigated with fixed distance cutoffs.  The
94 + calculated structure factor is used to sum out to great distance, and
95 + a surrounding dielectric term is included.}
96 + \label{fig:ewaldTime}
97 + \end{figure}
98 +
99 + \subsection{The Wolf and Zahn Methods}
100 + In a recent paper by Wolf \textit{et al.}, a procedure was outlined
101 + for an accurate accumulation of electrostatic interactions in an
102 + efficient pairwise fashion.\cite{Wolf99} Wolf \textit{et al.} observed
103 + that the electrostatic interaction is effectively short-ranged in
104 + condensed phase systems and that neutralization of the charge
105 + contained within the cutoff radius is crucial for potential
106 + stability. They devised a pairwise summation method that ensures
107 + charge neutrality and gives results similar to those obtained with
108 + the Ewald summation.  The resulting shifted Coulomb potential
109 + (Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through
110 + placement on the cutoff sphere and a distance-dependent damping
111 + function (identical to that seen in the real-space portion of the
112 + Ewald sum) to aid convergence
113 + \begin{equation}
114 + 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\}.
115 + \label{eq:WolfPot}
116 + \end{equation}
117 + Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted
118 + potential.  However, neutralizing the charge contained within each
119 + cutoff sphere requires the placement of a self-image charge on the
120 + surface of the cutoff sphere.  This additional self-term in the total
121 + potential enables Wolf {\it et al.}  to obtain excellent estimates of
122 + Madelung energies for many crystals.
123 +
124 + In order to use their charge-neutralized potential in molecular
125 + dynamics simulations, Wolf \textit{et al.} suggested taking the
126 + derivative of this potential prior to evaluation of the limit.  This
127 + procedure gives an expression for the forces,
128 + \begin{equation}
129 + 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\},
130 + \label{eq:WolfForces}
131 + \end{equation}
132 + that incorporates both image charges and damping of the electrostatic
133 + interaction.
134 +
135 + More recently, Zahn \textit{et al.} investigated these potential and
136 + force expressions for use in simulations involving water.\cite{Zahn02}
137 + In their work, they pointed out that the method that the forces and
138 + derivative of the potential are not commensurate.  Attempts to use
139 + both Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will
140 + lead to poor energy conservation.  They correctly observed that taking
141 + the limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating
142 + the derivatives is mathematically invalid.
143 +
144 + Zahn \textit{et al.} proposed a modified form of this ``Wolf summation
145 + method'' as a way to use this technique in Molecular Dynamics
146 + simulations.  Taking the integral of the forces shown in equation
147 + \ref{eq:WolfForces}, they proposed a new damped Coulomb
148 + potential,
149 + \begin{equation}
150 + 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\}.
151 + \label{eq:ZahnPot}
152 + \end{equation}
153 + They showed that this potential does fairly well at capturing the
154 + structural and dynamic properties of water compared the same
155 + properties obtained using the Ewald sum.
156 +
157 + \subsection{Simple Forms for Pairwise Electrostatics}
158 +
159 + The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et
160 + al.} are constructed using two different (and separable) computational
161 + tricks: \begin{itemize}
162 + \item shifting through the use of image charges, and
163 + \item damping the electrostatic interaction.
164 + \end{itemize}  Wolf \textit{et al.} treated the
165 + development of their summation method as a progressive application of
166 + these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded
167 + their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the
168 + post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using
169 + both techniques.  It is possible, however, to separate these
170 + tricks and study their effects independently.
171 +
172 + Starting with the original observation that the effective range of the
173 + electrostatic interaction in condensed phases is considerably less
174 + than $r^{-1}$, either the cutoff sphere neutralization or the
175 + distance-dependent damping technique could be used as a foundation for
176 + a new pairwise summation method.  Wolf \textit{et al.} made the
177 + observation that charge neutralization within the cutoff sphere plays
178 + a significant role in energy convergence; therefore we will begin our
179 + analysis with the various shifted forms that maintain this charge
180 + neutralization.  We can evaluate the methods of Wolf
181 + \textit{et al.}  and Zahn \textit{et al.} by considering the standard
182 + shifted potential,
183 + \begin{equation}
184 + v^\textrm{SP}(r) =      \begin{cases}
185 + v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r >
186 + R_\textrm{c}  
187 + \end{cases},
188 + \label{eq:shiftingPotForm}
189 + \end{equation}
190 + and shifted force,
191 + \begin{equation}
192 + v^\textrm{SF}(r) =      \begin{cases}
193 + v(r)-v_\textrm{c}-\left(\frac{\textrm{d}v(r)}{\textrm{d}r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c})
194 + &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c}
195 +                                                \end{cases},
196 + \label{eq:shiftingForm}
197 + \end{equation}
198 + functions where $v(r)$ is the unshifted form of the potential, and
199 + $v_c$ is $v(R_\textrm{c})$.  The Shifted Force ({\sc sf}) form ensures
200 + that both the potential and the forces goes to zero at the cutoff
201 + radius, while the Shifted Potential ({\sc sp}) form only ensures the
202 + potential is smooth at the cutoff radius
203 + ($R_\textrm{c}$).\cite{Allen87}
204 +
205 +
206 +
207 +
208 + 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}
209 + \begin{equation}
210 + V^\textrm{SP}(r_{ij}) = q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}.                          \label{eq:WolfSP}
211 + \end{equation}
212 + The forces associated with this potential are obtained by taking the derivative, resulting in the following,
213 + \begin{equation}
214 + F^\textrm{SP}(r_{ij}) = q_iq_j\left(-\frac{1}{r_{ij}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}.
215 + \label{eq:FWolfSP}
216 + \end{equation}
217 + 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}$.  
218 +
219 + If the derivative term in equation \ref{eq:shiftingForm} is evaluated, we obtain an hitherto undiscussed shifted force Coulomb potential,
220 + \begin{equation}
221 + V^\textrm{SF}(r_{ij}) = q_iq_j\left[\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r_{ij}-R_\textrm{c})\right] \quad r_{ij}\leqslant R_\textrm{c}.
222 + \label{eq:SFPot}
223 + \end{equation}
224 + Taking the derivative of this shifted force potential gives the
225 + following forces,
226 + \begin{equation}
227 + F^\textrm{SF}(r_{ij} =  q_iq_j\left(-\frac{1}{r_{ij}^2}+\frac{1}{R_\textrm{c}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}.
228 + \label{eq:SFForces}
229 + \end{equation}
230 + Using this formulation rather than the simple shifted potential
231 + (Eq. \ref{eq:WolfSP}) means that there are no discontinuities in the
232 + forces in addition to the potential.  This form also has the benefit
233 + that the image charges have a force presence, addressing the concerns
234 + about a missing physical component.  One side effect of this treatment
235 + is a slight alteration in the shape of the potential that comes about
236 + from the derivative term.  Thus, a degree of clarity about the
237 + original formulation of the potential is lost in order to gain
238 + functionality in dynamics simulations.
239 +
240 + Wolf \textit{et al.} originally discussed the energetics of the
241 + shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that
242 + it was still insufficient for accurate determination of the energy.
243 + The energy would fluctuate around the expected value with increasing
244 + cutoff radius, but the oscillations appeared to be converging toward
245 + the correct value.\cite{Wolf99} A damping function was incorporated to
246 + accelerate convergence; and though alternative functional forms could
247 + be used,\cite{Jones56,Heyes81} the complimentary error function was
248 + chosen to draw parallels to the Ewald summation.  Incorporating
249 + damping into the simple Coulomb potential,
250 + \begin{equation}
251 + v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}},
252 + \label{eq:dampCoulomb}
253 + \end{equation}
254 + the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived
255 + \textit{via} equation \ref{eq:shiftingForm},
256 + \begin{equation}
257 + V^{\textrm{DSP}}(r_{ij}) = q_iq_j\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}.
258 + \label{eq:DSPPot}
259 + \end{equation}
260 + The derivative of this Shifted-Potential can be taken to obtain forces
261 + for use in MD,
262 + \begin{equation}
263 + F^{\textrm{DSP}}(r_{ij}) = q_iq_j\left(-\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}-\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right) \quad r_{ij}\leqslant R_\textrm{c}.
264 + \label{eq:DSPForces}
265 + \end{equation}
266 + Again, this Shifted-Potential suffers from a discontinuity in the
267 + forces, and a lack of an image-charge component in the forces.  To
268 + remedy these concerns, a Shifted-Force variant is obtained by
269 + inclusion of the derivative term in equation \ref{eq:shiftingForm} to
270 + give,
271 + \begin{equation}
272 + \begin{split}
273 + V^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ &\left.+\left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right)\left(r_{ij}-R_\mathrm{c}\right)\ \right] \quad r_{ij}\leqslant R_\textrm{c}.
274 + \label{eq:DSFPot}
275 + \end{split}
276 + \end{equation}
277 + The derivative of the above potential gives the following forces,
278 + \begin{equation}
279 + \begin{split}
280 + F^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&-\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right) \\ &\left.+\left(\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right] \quad r_{ij}\leqslant R_\textrm{c}.
281 + \label{eq:DSFForces}
282 + \end{split}
283 + \end{equation}
284 +
285 + This new Shifted-Force potential is similar to equation
286 + \ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are
287 + two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term
288 + from equation \ref{eq:shiftingForm} is equal to equation
289 + \ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$.  This
290 + term is not present in the Zahn potential, resulting in a
291 + discontinuity as particles cross $R_\textrm{c}$.  Second, the sign of
292 + the derivative portion is different.  The constant $v_\textrm{c}$ term
293 + is not going to have a presence in the forces after performing the
294 + derivative, but the negative sign does effect the derivative.  In
295 + fact, it introduces a discontinuity in the forces at the cutoff,
296 + because the force function is shifted in the wrong direction and
297 + doesn't cross zero at $R_\textrm{c}$.  Thus, these alterations make
298 + for an electrostatic summation method that is continuous in both the
299 + potential and forces and incorporates the pairwise sum considerations
300 + stressed by Wolf \textit{et al.}\cite{Wolf99}
301 +
302   \section{Methods}
303  
304 < In each of the simulated systems, 500 distinct configurations were generated, and the electrostatic summation methods were compared via sequential application on each of these fixed configurations.  The methods compared include SPME, the aforementioned Shifted Potential and Shifted Force methods - both with damping parameters ($\alpha$) of 0, 0.1, 0.2, and 0.3 \AA$^{-1}$, 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.  
304 > \subsection{What Qualities are Important?}\label{sec:Qualities}
305 > In classical molecular mechanics simulations, there are two primary
306 > techniques utilized to obtain information about the system of
307 > interest: Monte Carlo (MC) and Molecular Dynamics (MD).  Both of these
308 > techniques utilize pairwise summations of interactions between
309 > particle sites, but they use these summations in different ways.
310  
311 < 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{argonSlice}).
311 > In MC, the potential energy difference between two subsequent
312 > configurations dictates the progression of MC sampling.  Going back to
313 > the origins of this method, the Canonical ensemble acceptance criteria
314 > laid out by Metropolis \textit{et al.} states that a subsequent
315 > configuration is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta
316 > E/kT)$, where $\xi$ is a random number between 0 and
317 > 1.\cite{Metropolis53} Maintaining a consistent $\Delta E$ when using
318 > an alternate method for handling the long-range electrostatics ensures
319 > proper sampling within the ensemble.
320  
321 + In MD, the derivative of the potential directs how the system will
322 + progress in time.  Consequently, the force and torque vectors on each
323 + body in the system dictate how it develops as a whole.  If the
324 + magnitude and direction of these vectors are similar when using
325 + alternate electrostatic summation techniques, the dynamics in the near
326 + term will be indistinguishable.  Because error in MD calculations is
327 + cumulative, one should expect greater deviation in the long term
328 + trajectories with greater differences in these vectors between
329 + configurations using different long-range electrostatics.
330 +
331 + \subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
332 + Evaluation of the pairwise summation techniques (outlined in section
333 + \ref{sec:ESMethods}) for use in MC simulations was performed through
334 + study of the energy differences between conformations.  Considering
335 + the SPME results to be the correct or desired behavior, ideal
336 + performance of a tested method was taken to be agreement between the
337 + energy differences calculated.  Linear least squares regression of the
338 + $\Delta E$ values between configurations using SPME against $\Delta E$
339 + values using tested methods provides a quantitative comparison of this
340 + agreement.  Unitary results for both the correlation and correlation
341 + coefficient for these regressions indicate equivalent energetic
342 + results between the methods.  The correlation is the slope of the
343 + plotted data while the correlation coefficient ($R^2$) is a measure of
344 + the of the data scatter around the fitted line and tells about the
345 + quality of the fit (Fig. \ref{fig:linearFit}).
346 +
347   \begin{figure}
348   \centering
349 < \includegraphics[width=3.25in]{./slice.pdf}
350 < \caption{A slice from the center of a water box used in a charge void simulation.  The darkened region represents the boundary sphere within which the water molecules were converted to argon atoms.}
351 < \label{argonSlice}
349 > \includegraphics[width = \linewidth]{./dualLinear.pdf}
350 > \caption{Example least squares regressions of the configuration energy differences for SPC/E water systems. The upper plot shows a data set with a poor correlation coefficient ($R^2$), while the lower plot shows a data set with a good correlation coefficient.}
351 > \label{fig:linearFit}
352   \end{figure}
353  
354 < All of these comparisons were performed with three different cutoff radii (9, 12, and 15 \AA) to investigate the cutoff radius dependence of the various techniques.  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^{-5}$ kcal/mol).  We chose a tolerance of $1 \times 10^{-8}$, resulting in Ewald Coefficients of 0.4200, 0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively.
354 > Each system type (detailed in section \ref{sec:RepSims}) studied
355 > consisted of 500 independent configurations, each equilibrated from
356 > higher temperature trajectories. Thus, 124,750 $\Delta E$ data points
357 > are used in a regression of a single system type.  Results and
358 > discussion for the individual analysis of each of the system types
359 > appear in the supporting information, while the cumulative results
360 > over all the investigated systems appears below in section
361 > \ref{sec:EnergyResults}.
362  
363 < \section{Results and Discussion}
363 > \subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods}
364 > Evaluation of the pairwise methods (outlined in section
365 > \ref{sec:ESMethods}) for use in MD simulations was performed through
366 > comparison of the force and torque vectors obtained with those from
367 > SPME.  Both the magnitude and the direction of these vectors on each
368 > of the bodies in the system were analyzed.  For the magnitude of these
369 > vectors, linear least squares regression analysis can be performed as
370 > described previously for comparing $\Delta E$ values. Instead of a
371 > single value between two system configurations, there is a value for
372 > each particle in each configuration.  For a system of 1000 water
373 > molecules and 40 ions, there are 1040 force vectors and 1000 torque
374 > vectors.  With 500 configurations, this results in 520,000 force and
375 > 500,000 torque vector comparisons samples for each system type.
376  
377 < \subsection{$\Delta E$ Comparison}
378 < In order to evaluate the performance of the adapted Wolf Shifted Potential and Shifted Force electrostatic summation methods for Monte Carlo simulations, the energy differences between configurations need to be compared to the results using SPME.  Considering the SPME results to be the correct or desired behavior, ideal performance of a tested method is 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 gives an idea of the quality of the fit (Fig. \ref{linearFit}).
377 > The force and torque vector directions were investigated through
378 > measurement of the angle ($\theta$) formed between those from the
379 > particular method and those from SPME
380 > \begin{equation}
381 > \theta_F = \frac{\vec{F}_\textrm{SPME}}{|\vec{F}_\textrm{SPME}|}\cdot\frac{\vec{F}_\textrm{Method}}{|\vec{F}_\textrm{Method}|}.
382 > \end{equation}
383 > Each of these $\theta$ values was accumulated in a distribution
384 > function, weighted by the area on the unit sphere.  Non-linear fits
385 > were used to measure the shape of the resulting distributions.
386  
387   \begin{figure}
388   \centering
389 < \includegraphics[width=3.25in]{./linearFit.pdf}
390 < \caption{Example least squares regression of the $\Delta E$ between configurations for the SF method against SPME in the pure water system.  }
391 < \label{linearFit}
389 > \includegraphics[width = \linewidth]{./gaussFit.pdf}
390 > \caption{Sample fit of the angular distribution of the force vectors over all of the studied systems.  Gaussian fits were used to obtain values for the variance in force and torque vectors used in the following figure.}
391 > \label{fig:gaussian}
392   \end{figure}
393  
394 < With 500 independent configurations, 124,750 $\Delta E$ data points are used in a regression of a single system.  Results and discussion for the individual analysis of each of the system types appear in the supporting information.  To probe the applicability of each method in the general case, all the different system types were included in a single regression.  The results for this regression are shown in figure \ref{delE}.  
394 > Figure \ref{fig:gaussian} shows an example distribution with applied
395 > non-linear fits.  The solid line is a Gaussian profile, while the
396 > dotted line is a Voigt profile, a convolution of a Gaussian and a
397 > Lorentzian.  Since this distribution is a measure of angular error
398 > between two different electrostatic summation methods, there is
399 > particular reason for the profile to adhere to a specific shape.
400 > Because of this and the Gaussian profile's more statistically
401 > meaningful properties, Gaussian fits was used to compare all the
402 > tested methods.  The variance ($\sigma^2$) was extracted from each of
403 > these fits and was used to compare distribution widths.  Values of
404 > $\sigma^2$ near zero indicate vector directions indistinguishable from
405 > those calculated when using SPME.
406  
407 + \subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
408 + Evaluation of the long-time dynamics of charged systems was performed
409 + by considering the NaCl crystal system while using a subset of the
410 + best performing pairwise methods.  The NaCl crystal was chosen to
411 + avoid possible complications involving the propagation techniques of
412 + orientational motion in molecular systems.  To enhance the atomic
413 + motion, these crystals were equilibrated at 1000 K, near the
414 + experimental $T_m$ for NaCl.  Simulations were performed under the
415 + microcanonical ensemble, and velocity autocorrelation functions
416 + (Eq. \ref{eq:vCorr}) were computed for each of the trajectories,
417 + \begin{equation}
418 + C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle.
419 + \label{eq:vCorr}
420 + \end{equation}
421 + Velocity autocorrelation functions require detailed short time data
422 + and long trajectories for good statistics, thus velocity information
423 + was saved every 5 fs over 100 ps trajectories.  The power spectrum
424 + ($I(\omega)$) is obtained via Fourier transform of the autocorrelation
425 + function
426 + \begin{equation}
427 + I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt,
428 + \label{eq:powerSpec}
429 + \end{equation}
430 + where the frequency, $\omega=0,\ 1,\ ...,\ N-1$.
431 +
432 + \subsection{Representative Simulations}\label{sec:RepSims}
433 + A variety of common and representative simulations were analyzed to
434 + determine the relative effectiveness of the pairwise summation
435 + techniques in reproducing the energetics and dynamics exhibited by
436 + SPME.  The studied systems were as follows:
437 + \begin{enumerate}
438 + \item Liquid Water
439 + \item Crystalline Water (Ice I$_\textrm{c}$)
440 + \item NaCl Crystal
441 + \item NaCl Melt
442 + \item Low Ionic Strength Solution of NaCl in Water
443 + \item High Ionic Strength Solution of NaCl in Water
444 + \item 6 \AA\  Radius Sphere of Argon in Water
445 + \end{enumerate}
446 + By utilizing the pairwise techniques (outlined in section
447 + \ref{sec:ESMethods}) in systems composed entirely of neutral groups,
448 + charged particles, and mixtures of the two, we can comment on possible
449 + system dependence and/or universal applicability of the techniques.
450 +
451 + Generation of the system configurations was dependent on the system
452 + type.  For the solid and liquid water configurations, configuration
453 + snapshots were taken at regular intervals from higher temperature 1000
454 + SPC/E water molecule trajectories and each equilibrated individually.
455 + The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl-
456 + ions and were selected and equilibrated in the same fashion as the
457 + water systems.  For the low and high ionic strength NaCl solutions, 4
458 + and 40 ions were first solvated in a 1000 water molecule boxes
459 + respectively.  Ion and water positions were then randomly swapped, and
460 + the resulting configurations were again equilibrated individually.
461 + Finally, for the Argon/Water "charge void" systems, the identities of
462 + all the SPC/E waters within 6 \AA\ of the center of the equilibrated
463 + water configurations were converted to argon
464 + (Fig. \ref{fig:argonSlice}).
465 +
466   \begin{figure}
467   \centering
468 < \includegraphics[width=3.25in]{./delEplot.pdf}
469 < \caption{The results from the statistical analysis of the $\Delta$E results for all the system types at 9 \AA\ (${\bullet}$), 12 \AA\ ($\blacksquare$), and 15 \AA\ ($\blacktriangledown$) cutoff radii.  Results close to a value of 1 (dashed line) indicate $\Delta E$ values from that particular method (listed on the left) are nearly indistinguishable from those obtained from SPME.  Reaction Field results do not include NaCl crystal or melt configurations.}
470 < \label{delE}
468 > \includegraphics[width = \linewidth]{./slice.pdf}
469 > \caption{A slice from the center of a water box used in a charge void simulation.  The darkened region represents the boundary sphere within which the water molecules were converted to argon atoms.}
470 > \label{fig:argonSlice}
471   \end{figure}
472  
473 < In figure \ref{delE}, 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.  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.  Correcting the resulting charged cutoff sphere is one of the purposes of the shifted potential proposed by Wolf \textit{et al.}, 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.  This trend is repeated in the Shifted Force rows, where increasing damping results in progressively poorer correlation; however, damping looks to be unnecessary with this method.  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.
473 > \subsection{Electrostatic Summation Methods}\label{sec:ESMethods}
474 > Electrostatic summation method comparisons were performed using SPME,
475 > the Shifted-Potential and Shifted-Force methods - both with damping
476 > parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak,
477 > moderate, and strong damping respectively), reaction field with an
478 > infinite dielectric constant, and an unmodified cutoff.  Group-based
479 > cutoffs with a fifth-order polynomial switching function were
480 > necessary for the reaction field simulations and were utilized in the
481 > SP, SF, and pure cutoff methods for comparison to the standard lack of
482 > group-based cutoffs with a hard truncation.  The SPME calculations
483 > were performed using the TINKER implementation of SPME,\cite{Ponder87}
484 > while all other method calculations were performed using the OOPSE
485 > molecular mechanics package.\cite{Meineke05}
486  
487 < \subsection{Force Magnitude Comparison}
487 > These methods were additionally evaluated with three different cutoff
488 > radii (9, 12, and 15 \AA) to investigate possible cutoff radius
489 > dependence.  It should be noted that the damping parameter chosen in
490 > SPME, or so called ``Ewald Coefficient", has a significant effect on
491 > the energies and forces calculated.  Typical molecular mechanics
492 > packages default this to a value dependent on the cutoff radius and a
493 > tolerance (typically less than $1 \times 10^{-4}$ kcal/mol).  Smaller
494 > tolerances are typically associated with increased accuracy in the
495 > real-space portion of the summation.\cite{Essmann95} The default
496 > TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME
497 > calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and
498 > 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively.
499  
500 < While studying the energy differences provides insight into how comparable these methods are energetically, if we want to use these methods in Molecular Dynamics simulations, we also need to consider their effect on forces and torques.  Both the magnitude and the direction of the force and torque vectors of each of the bodies in the system can be compared to those observed while using SPME.  Analysis of the magnitude of these vectors can be performed in the manner described previously for comparing $\Delta E$ values, only 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 excess of 500,000 data samples for each system type.  Figures \ref{frcMag} and \ref{trqMag} respectively show the force and torque vector magnitude results for the accumulated analysis over all the system types.
500 > \section{Results and Discussion}
501  
502 + \subsection{Configuration Energy Differences}\label{sec:EnergyResults}
503 + In order to evaluate the performance of the pairwise electrostatic
504 + summation methods for Monte Carlo simulations, the energy differences
505 + between configurations were compared to the values obtained when using
506 + SPME.  The results for the subsequent regression analysis are shown in
507 + figure \ref{fig:delE}.
508 +
509   \begin{figure}
510   \centering
511 < \includegraphics[width=3.25in]{./frcMagplot.pdf}
512 < \caption{The results from the statistical analysis of the force vector magnitude results for all the system types at 9 \AA\ (${\bullet}$), 12 \AA\ ($\blacksquare$), and 15 \AA\ ($\blacktriangledown$) cutoff radii.  Results close to a value of 1 (dashed line) indicate force vector magnitude values from that particular method (listed on the left) are nearly indistinguishable from those obtained from SPME.}
513 < \label{frcMag}
511 > \includegraphics[width=5.5in]{./delEplot.pdf}
512 > \caption{Statistical analysis of the quality of configurational energy differences for a given electrostatic method compared with the reference Ewald sum.  Results with a value equal to 1 (dashed line) indicate $\Delta E$ values indistinguishable from those obtained using SPME.  Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).}
513 > \label{fig:delE}
514   \end{figure}
515  
516 < The results in figure \ref{frcMag} for the most part parallel those seen in the previous look at the $\Delta E$ results.  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 gives little to no gain.  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.
516 > In this figure, it is apparent that it is unreasonable to expect
517 > realistic results using an unmodified cutoff.  This is not all that
518 > surprising since this results in large energy fluctuations as atoms
519 > move in and out of the cutoff radius.  These fluctuations can be
520 > alleviated to some degree by using group based cutoffs with a
521 > switching function.\cite{Steinbach94} The Group Switch Cutoff row
522 > doesn't show a significant improvement in this plot because the salt
523 > and salt solution systems contain non-neutral groups, see the
524 > accompanying supporting information for a comparison where all groups
525 > are neutral.
526  
527 < \subsection{Torque Magnitude Comparison}
527 > Correcting the resulting charged cutoff sphere is one of the purposes
528 > of the damped Coulomb summation proposed by Wolf \textit{et
529 > al.},\cite{Wolf99} and this correction indeed improves the results as
530 > seen in the Shifted-Potental rows.  While the undamped case of this
531 > method is a significant improvement over the pure cutoff, it still
532 > doesn't correlate that well with SPME.  Inclusion of potential damping
533 > improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows
534 > an excellent correlation and quality of fit with the SPME results,
535 > particularly with a cutoff radius greater than 12 \AA .  Use of a
536 > larger damping parameter is more helpful for the shortest cutoff
537 > shown, but it has a detrimental effect on simulations with larger
538 > cutoffs.  In the Shifted-Force sets, increasing damping results in
539 > progressively poorer correlation.  Overall, the undamped case is the
540 > best performing set, as the correlation and quality of fits are
541 > consistently superior regardless of the cutoff distance.  This result
542 > is beneficial in that the undamped case is less computationally
543 > prohibitive do to the lack of complimentary error function calculation
544 > when performing the electrostatic pair interaction.  The reaction
545 > field results illustrates some of that method's limitations, primarily
546 > that it was developed for use in homogenous systems; although it does
547 > provide results that are an improvement over those from an unmodified
548 > cutoff.
549  
550 + \subsection{Magnitudes of the Force and Torque Vectors}
551 +
552 + Evaluation of pairwise methods for use in Molecular Dynamics
553 + simulations requires consideration of effects on the forces and
554 + torques.  Investigation of the force and torque vector magnitudes
555 + provides a measure of the strength of these values relative to SPME.
556 + Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the
557 + force and torque vector magnitude regression results for the
558 + accumulated analysis over all the system types.
559 +
560   \begin{figure}
561   \centering
562 < \includegraphics[width=3.25in]{./trqMagplot.pdf}
563 < \caption{The results from the statistical analysis of the torque vector magnitude results for all the system types at 9 \AA\ (${\bullet}$), 12 \AA\ ($\blacksquare$), and 15 \AA\ ($\blacktriangledown$) cutoff radii.  Results close to a value of 1 (dashed line) indicate torque vector magnitude values from that particular method (listed on the left) are nearly indistinguishable from those obtained from SPME.  Torques are only accumulated on the rigid water molecules, so these results exclude NaCl the systems.}
564 < \label{trqMag}
562 > \includegraphics[width=5.5in]{./frcMagplot.pdf}
563 > \caption{Statistical analysis of the quality of the force vector magnitudes for a given electrostatic method compared with the reference Ewald sum.  Results with a value equal to 1 (dashed line) indicate force magnitude values indistinguishable from those obtained using SPME.  Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).}
564 > \label{fig:frcMag}
565   \end{figure}
566  
567 < The torque vector magnitude results in figure \ref{trqMag} are similar to those seen for the forces, but more clearly show the improved behavior with increasing cutoff radius.  Moderate damping is beneficial to the Shifted Potential and 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.
567 > Figure \ref{fig:frcMag}, for the most part, parallels the results seen
568 > in the previous $\Delta E$ section.  The unmodified cutoff results are
569 > poor, but using group based cutoffs and a switching function provides
570 > a improvement much more significant than what was seen with $\Delta
571 > E$.  Looking at the Shifted-Potential sets, the slope and $R^2$
572 > improve with the use of damping to an optimal result of 0.2 \AA
573 > $^{-1}$ for the 12 and 15 \AA\ cutoffs.  Further increases in damping,
574 > while beneficial for simulations with a cutoff radius of 9 \AA\ , is
575 > detrimental to simulations with larger cutoff radii.  The undamped
576 > Shifted-Force method gives forces in line with those obtained using
577 > SPME, and use of a damping function results in minor improvement.  The
578 > reaction field results are surprisingly good, considering the poor
579 > quality of the fits for the $\Delta E$ results.  There is still a
580 > considerable degree of scatter in the data, but it correlates well in
581 > general.  To be fair, we again note that the reaction field
582 > calculations do not encompass NaCl crystal and melt systems, so these
583 > results are partly biased towards conditions in which the method
584 > performs more favorably.
585  
124 – \subsection{Force and Torque Direction Comparison}
125 –
126 – 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 dot product of these unit vectors provides a theta value that is accumulated in a distribution function, weighted by the area on the unit sphere.  Narrow distributions of theta values indicates similar to identical results between the tested method and SPME.  To measure the narrowness of the resulting distributions, non-linear Gaussian fits were performed.
127 –
586   \begin{figure}
587   \centering
588 < \includegraphics[width=3.25in]{./gaussFit.pdf}
589 < \caption{Example fitting of the angular distribution of the force vectors over all of the studied systems.  The solid and dotted lines show Gaussian and Voigt fits of the distribution data respectively.  Even though the Voigt profile make for a more accurate fit, the Gaussian was used due to more versatile statistical results.}
590 < \label{gaussian}
588 > \includegraphics[width=5.5in]{./trqMagplot.pdf}
589 > \caption{Statistical analysis of the quality of the torque vector magnitudes for a given electrostatic method compared with the reference Ewald sum.  Results with a value equal to 1 (dashed line) indicate torque magnitude values indistinguishable from those obtained using SPME.  Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).}
590 > \label{fig:trqMag}
591   \end{figure}
592  
593 < Figure \ref{gaussian} shows an example distribution and the non-linear fit applied.  The solid line is a Gaussian profile, while the dotted line is a Voigt profile, a convolution of a Gaussian and a Lorentzian profile.  Since this distribution is a measure of angular error between two different electrostatic summation methods, there is particular reason for it to adhere to a particular shape.  Because of this and the Gaussian profile's more statistically meaningful properties, Gaussian fitting was used to compare all the methods considered in this study.  The results (Fig. \ref{frcTrqAng}) are compared through the variance ($\sigma^2$) of these non-linear fits.  
593 > To evaluate the torque vector magnitudes, the data set from which
594 > values are drawn is limited to rigid molecules in the systems
595 > (i.e. water molecules).  In spite of this smaller sampling pool, the
596 > torque vector magnitude results in figure \ref{fig:trqMag} are still
597 > similar to those seen for the forces; however, they more clearly show
598 > the improved behavior that comes with increasing the cutoff radius.
599 > Moderate damping is beneficial to the Shifted-Potential and helpful
600 > yet possibly unnecessary with the Shifted-Force method, and they also
601 > show that over-damping adversely effects all cutoff radii rather than
602 > showing an improvement for systems with short cutoffs.  The reaction
603 > field method performs well when calculating the torques, better than
604 > the Shifted Force method over this limited data set.
605  
606 + \subsection{Directionality of the Force and Torque Vectors}
607 +
608 + Having force and torque vectors with magnitudes that are well
609 + correlated to SPME is good, but if they are not pointing in the proper
610 + direction the results will be incorrect.  These vector directions were
611 + investigated through measurement of the angle formed between them and
612 + those from SPME.  The results (Fig. \ref{fig:frcTrqAng}) are compared
613 + through the variance ($\sigma^2$) of the Gaussian fits of the angle
614 + error distributions of the combined set over all system types.
615 +
616   \begin{figure}
617   \centering
618 < \includegraphics[width=3.25in]{./frcTrqAngplot.pdf}
619 < \caption{The results from the statistical analysis of the force and torque vector angular distributions for all the system types at 9 \AA\ (${\bullet}$), 12 \AA\ ($\blacksquare$), and 15 \AA\ ($\blacktriangledown$) cutoff radii.  Plotted values are the variance ($\sigma^2$) of the Gaussian non-linear fits.  Results close to a value of 0 (dashed line) indicate force or torque vector directions from that particular method (listed on the left) are nearly indistinguishable from those obtained from SPME.  Torques are only accumulated on the rigid water molecules, so the torque vector angle results exclude NaCl the systems.}
620 < \label{frcTrqAng}
618 > \includegraphics[width=5.5in]{./frcTrqAngplot.pdf}
619 > \caption{Statistical analysis of the quality of the Gaussian fit of the force and torque vector angular distributions for a given electrostatic method compared with the reference Ewald sum.  Results with a variance ($\sigma^2$) equal to zero (dashed line) indicate force and torque directions indistinguishable from those obtained using SPME.  Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).}
620 > \label{fig:frcTrqAng}
621   \end{figure}
622  
623 < Both the force and torque $\sigma^2$ results from the analysis of the total accumulated system data are tabulated in figure \ref{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 $\sigma^2$, 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 non-neutral bodies, and this observation is investigated further in the accompanying supporting information.  
623 > Both the force and torque $\sigma^2$ results from the analysis of the
624 > total accumulated system data are tabulated in figure
625 > \ref{fig:frcTrqAng}.  All of the sets, aside from the over-damped case
626 > show the improvement afforded by choosing a longer simulation cutoff.
627 > Increasing the cutoff from 9 to 12 \AA\ typically results in a halving
628 > of the distribution widths, with a similar improvement going from 12
629 > to 15 \AA .  The undamped Shifted-Force, Group Based Cutoff, and
630 > Reaction Field methods all do equivalently well at capturing the
631 > direction of both the force and torque vectors.  Using damping
632 > improves the angular behavior significantly for the Shifted-Potential
633 > and moderately for the Shifted-Force methods.  Increasing the damping
634 > too far is destructive for both methods, particularly to the torque
635 > vectors.  Again it is important to recognize that the force vectors
636 > cover all particles in the systems, while torque vectors are only
637 > available for neutral molecular groups.  Damping appears to have a
638 > more beneficial effect on non-neutral bodies, and this observation is
639 > investigated further in the accompanying supporting information.
640  
641   \begin{table}[htbp]
642     \centering
# Line 173 | Line 668 | Both the force and torque $\sigma^2$ results from the
668  
669        \bottomrule
670     \end{tabular}
671 <   \label{groupAngle}
671 >   \label{tab:groupAngle}
672   \end{table}
673  
674 < Although not discussed previously, group based cutoffs can be applied to both the Shifted Potential and 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{groupAngle} shows the angular variance values obtained using group based cutoffs and a switching function alongside the standard results seen in figure \ref{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.  
674 > Although not discussed previously, group based cutoffs can be applied
675 > to both the Shifted-Potential and Shifted-Force methods.  Use off a
676 > switching function corrects for the discontinuities that arise when
677 > atoms of a group exit the cutoff before the group's center of mass.
678 > Though there are no significant benefit or drawbacks observed in
679 > $\Delta E$ and vector magnitude results when doing this, there is a
680 > measurable improvement in the vector angle results.  Table
681 > \ref{tab:groupAngle} shows the angular variance values obtained using
682 > group based cutoffs and a switching function alongside the standard
683 > results seen in figure \ref{fig:frcTrqAng} for comparison purposes.
684 > The Shifted-Potential shows much narrower angular distributions for
685 > both the force and torque vectors when using an $\alpha$ of 0.2
686 > \AA$^{-1}$ or less, while Shifted-Force shows improvements in the
687 > undamped and lightly damped cases.  Thus, by calculating the
688 > electrostatic interactions in terms of molecular pairs rather than
689 > atomic pairs, the direction of the force and torque vectors are
690 > determined more accurately.
691  
692 < One additional trend to recognize in table \ref{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{delE}, \ref{frcMag}, and \ref{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{delE}); however, based on these findings, choices this high would be introducing error in the molecular torques, particularly for the shorter cutoffs.  Based on the above findings, any empirical damping is arguably unnecessary with the choice of the Shifted Force method.
692 > One additional trend to recognize in table \ref{tab:groupAngle} is
693 > that the $\sigma^2$ values for both Shifted-Potential and
694 > Shifted-Force converge as $\alpha$ increases, something that is easier
695 > to see when using group based cutoffs.  Looking back on figures
696 > \ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this
697 > behavior clearly at large $\alpha$ and cutoff values.  The reason for
698 > this is that the complimentary error function inserted into the
699 > potential weakens the electrostatic interaction as $\alpha$ increases.
700 > Thus, at larger values of $\alpha$, both the summation method types
701 > progress toward non-interacting functions, so care is required in
702 > choosing large damping functions lest one generate an undesirable loss
703 > in the pair interaction.  Kast \textit{et al.}  developed a method for
704 > choosing appropriate $\alpha$ values for these types of electrostatic
705 > summation methods by fitting to $g(r)$ data, and their methods
706 > indicate optimal values of 0.34, 0.25, and 0.16 \AA$^{-1}$ for cutoff
707 > values of 9, 12, and 15 \AA\ respectively.\cite{Kast03} These appear
708 > to be reasonable choices to obtain proper MC behavior
709 > (Fig. \ref{fig:delE}); however, based on these findings, choices this
710 > high would introduce error in the molecular torques, particularly for
711 > the shorter cutoffs.  Based on the above findings, empirical damping
712 > up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably
713 > unnecessary when using the Shifted-Force method.
714  
715 + \subsection{Collective Motion: Power Spectra of NaCl Crystals}
716 +
717 + In the previous studies using a Shifted-Force variant of the damped
718 + Wolf coulomb potential, the structure and dynamics of water were
719 + investigated rather extensively.\cite{Zahn02,Kast03} Their results
720 + indicated that the damped Shifted-Force method results in properties
721 + very similar to those obtained when using the Ewald summation.
722 + Considering the statistical results shown above, the good performance
723 + of this method is not that surprising.  Rather than consider the same
724 + systems and simply recapitulate their results, we decided to look at
725 + the solid state dynamical behavior obtained using the best performing
726 + summation methods from the above results.
727 +
728 + \begin{figure}
729 + \centering
730 + \includegraphics[width = \linewidth]{./spectraSquare.pdf}
731 + \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.}
732 + \label{fig:methodPS}
733 + \end{figure}
734 +
735 + Figure \ref{fig:methodPS} shows the power spectra for the NaCl
736 + crystals (from averaged Na and Cl ion velocity autocorrelation
737 + functions) using the stated electrostatic summation methods.  While
738 + high frequency peaks of all the spectra overlap, showing the same
739 + general features, the low frequency region shows how the summation
740 + methods differ.  Considering the low-frequency inset (expanded in the
741 + upper frame of figure \ref{fig:dampInc}), at frequencies below 100
742 + cm$^{-1}$, the correlated motions are blue-shifted when using undamped
743 + or weakly damped Shifted-Force.  When using moderate damping ($\alpha
744 + = 0.2$ \AA$^{-1}$) both the Shifted-Force and Shifted-Potential
745 + methods give near identical correlated motion behavior as the Ewald
746 + method (which has a damping value of 0.3119).  The damping acts as a
747 + distance dependent Gaussian screening of the point charges for the
748 + pairwise summation methods.  This weakening of the electrostatic
749 + interaction with distance explains why the long-ranged correlated
750 + motions are at lower frequencies for the moderately damped methods
751 + than for undamped or weakly damped methods.  To see this effect more
752 + clearly, we show how damping strength affects a simple real-space
753 + electrostatic potential,
754 + \begin{equation}
755 + 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),
756 + \end{equation}
757 + where $S(r)$ is a switching function that smoothly zeroes the
758 + potential at the cutoff radius.  Figure \ref{fig:dampInc} shows how
759 + the low frequency motions are dependent on the damping used in the
760 + direct electrostatic sum.  As the damping increases, the peaks drop to
761 + lower frequencies.  Incidentally, use of an $\alpha$ of 0.25
762 + \AA$^{-1}$ on a simple electrostatic summation results in low
763 + frequency correlated dynamics equivalent to a simulation using SPME.
764 + When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks
765 + shift to higher frequency in exponential fashion.  Though not shown,
766 + the spectrum for the simple undamped electrostatic potential is
767 + blue-shifted such that the lowest frequency peak resides near 325
768 + cm$^{-1}$.  In light of these results, the undamped Shifted-Force
769 + method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is
770 + quite respectable; however, it appears as though moderate damping is
771 + required for accurate reproduction of crystal dynamics.
772 + \begin{figure}
773 + \centering
774 + \includegraphics[width = \linewidth]{./comboSquare.pdf}
775 + \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.}
776 + \label{fig:dampInc}
777 + \end{figure}
778 +
779   \section{Conclusions}
780  
781 < \section{Acknowledgments}
781 > This investigation of pairwise electrostatic summation techniques
782 > shows that there are viable and more computationally efficient
783 > electrostatic summation techniques than the Ewald summation, chiefly
784 > methods derived from the damped Coulombic sum originally proposed by
785 > Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the
786 > Shifted-Force method, reformulated above as equation \ref{eq:SFPot},
787 > shows a remarkable ability to reproduce the energetic and dynamic
788 > characteristics exhibited by simulations employing lattice summation
789 > techniques.  The cumulative energy difference results showed the
790 > undamped Shifted-Force and moderately damped Shifted-Potential methods
791 > produced results nearly identical to SPME.  Similarly for the dynamic
792 > features, the undamped or moderately damped Shifted-Force and
793 > moderately damped Shifted-Potential methods produce force and torque
794 > vector magnitude and directions very similar to the expected values.
795 > These results translate into long-time dynamic behavior equivalent to
796 > that produced in simulations using SPME.
797  
798 + Aside from the computational cost benefit, these techniques have
799 + applicability in situations where the use of the Ewald sum can prove
800 + problematic.  Primary among them is their use in interfacial systems,
801 + where the unmodified lattice sum techniques artificially accentuate
802 + the periodicity of the system in an undesirable manner.  There have
803 + been alterations to the standard Ewald techniques, via corrections and
804 + reformulations, to compensate for these systems; but the pairwise
805 + techniques discussed here require no modifications, making them
806 + natural tools to tackle these problems.  Additionally, this
807 + transferability gives them benefits over other pairwise methods, like
808 + reaction field, because estimations of physical properties (e.g. the
809 + dielectric constant) are unnecessary.
810 +
811 + We are not suggesting any flaw with the Ewald sum; in fact, it is the
812 + standard by which these simple pairwise sums are judged.  However,
813 + these results do suggest that in the typical simulations performed
814 + today, the Ewald summation may no longer be required to obtain the
815 + level of accuracy most researcher have come to expect
816 +
817 + \section{Acknowledgments}
818   \newpage
819  
820 < \bibliographystyle{achemso}
820 > \bibliographystyle{jcp2}
821   \bibliography{electrostaticMethods}
822  
823  

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines