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

Diff Legend

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