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 2616 by chrisfen, Tue Mar 14 14:31:19 2006 UTC vs.
Revision 2617 by gezelter, Tue Mar 14 16:22:54 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}
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{bibentry}
15   \usepackage{mathrsfs}
15 – %\usepackage{berkeley}
16   \usepackage[ref]{overcite}
17   \pagestyle{plain}
18   \pagenumbering{arabic}
# Line 25 | Line 25
25  
26   \begin{document}
27  
28 < \title{Is the Ewald Summation necessary? : Pairwise alternatives to the accepted standard for long-range electrostatics}
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 35 | 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 from the previous efforts described in \bibentry{Wolf99} and \bibentry{Zahn02} as a possible replacement for lattice sum methods in molecular simulations.  Comparisons were performed with this and other pairwise electrostatic summation techniques against the smooth particle mesh Ewald (SPME) summation to see how well they reproduce the energetics and dynamics of a variety of simulation types.  The newly derived Shifted-Force technique shows a remarkable ability to reproduce the behavior exhibited in simulations using SPME with an $\mathscr{O}(N)$ computational cost, equivalent to merely the real-space portion of the lattice summation.  
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   \end{abstract}
55  
56 + \newpage
57 +
58   %\narrowtext
59  
60   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# Line 49 | Line 63 | In molecular simulations, proper accumulation of the e
63  
64   \section{Introduction}
65  
66 < In molecular simulations, proper accumulation of the electrostatic interactions is considered one of the most essential and computationally demanding tasks.  
66 > In molecular simulations, proper accumulation of the electrostatic
67 > interactions is considered one of the most essential and
68 > computationally demanding tasks.
69  
70   \subsection{The Ewald Sum}
71   blah blah blah Ewald Sum Important blah blah blah
72  
73   \begin{figure}
74   \centering
75 < \includegraphics[width = 3.25in]{./ewaldProgression.pdf}
76 < \caption{How the application of the Ewald summation has changed with the increase in computer power.  Initially, only small numbers of particles could be studied, and the Ewald sum acted to replicate the unit cell charge distribution out to convergence.  Now, much larger systems of charges are investigated with fixed distance cutoffs.  The calculated structure factor is used to sum out to great distance, and a surrounding dielectric term is included.}
75 > \includegraphics[width = \linewidth]{./ewaldProgression.pdf}
76 > \caption{How the application of the Ewald summation has changed with
77 > the increase in computer power.  Initially, only small numbers of
78 > particles could be studied, and the Ewald sum acted to replicate the
79 > unit cell charge distribution out to convergence.  Now, much larger
80 > systems of charges are investigated with fixed distance cutoffs.  The
81 > calculated structure factor is used to sum out to great distance, and
82 > a surrounding dielectric term is included.}
83   \label{fig:ewaldTime}
84   \end{figure}
85  
86   \subsection{The Wolf and Zahn Methods}
87 < In a recent paper by Wolf \textit{et al.}, a procedure was outlined for accumulation of electrostatic interactions in a simple pairwise fashion.\cite{Wolf99}  They took the observation that the effective electrostatic interaction is short-ranged in systems of charges and that charge neutralization within the cutoff spheres is crucial for potential stability. They devised a pairwise summation method that ensures charge neutrality and gives results similar to those obtained using the Ewald summation.  The resulting shifted Coulomb potential (Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through placement on the cutoff sphere and a distance-dependent damping function (identical to that seen in the real-space portion of the Ewald sum) to aid energetic convergence
87 > In a recent paper by Wolf \textit{et al.}, a procedure was outlined
88 > for an accurate accumulation of electrostatic interactions in an
89 > efficient pairwise fashion.\cite{Wolf99} Wolf \textit{et al.} observed
90 > that the electrostatic interaction is effectively short-ranged in
91 > condensed phase systems and that neutralization of the charge
92 > contained within the cutoff radius is crucial for potential
93 > stability. They devised a pairwise summation method that ensures
94 > charge neutrality and gives results similar to those obtained with
95 > the Ewald summation.  The resulting shifted Coulomb potential
96 > (Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through
97 > placement on the cutoff sphere and a distance-dependent damping
98 > function (identical to that seen in the real-space portion of the
99 > Ewald sum) to aid convergence
100   \begin{equation}
101   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\}.
102   \label{eq:WolfPot}
103   \end{equation}
104 < In order to use this potential in molecular dynamics simulations, Wolf \textit{et al.} suggested taking the derivative of this potential, followed by evaluation of the limit to give the following forces,
104 > Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted
105 > potential.  However, neutralizing the charge contained within each
106 > cutoff sphere requires the placement of a self-image charge on the
107 > surface of the cutoff sphere.  This additional self-term in the total
108 > potential enables Wolf {\it et al.}  to obtain excellent estimates of
109 > Madelung energies for many crystals.
110 >
111 > In order to use their charge-neutralized potential in molecular
112 > dynamics simulations, Wolf \textit{et al.} suggested taking the
113 > derivative of this potential prior to evaluation of the limit.  This
114 > procedure gives an expression for the forces,
115   \begin{equation}
116 < F^{\textrm{Wolf}}(r_{ij}) = q_iq_j\left\{\left[-\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right]+\left[\frac{\textrm{erfc}\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\}.
116 > 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\},
117   \label{eq:WolfForces}
118   \end{equation}
119 < More recently, Zahn \textit{et al.} investigated this electrostatic summation method for use in simulations involving water.\cite{Zahn02}  In their work, they point out that the method as proposed is problematic for use in Molecular Dynamics simulations, because the forces and derivative of the potential are not equivalent.  This comes about from the procedure of taking the limit shown in equation \ref{eq:WolfPot} after calculating the derivatives.\cite{Wolf99}  Zahn \textit{et al.} proposed a shifted force adaptation of this ``Wolf summation method" as a way to use this technique in Molecular Dynamics simulations.  Taking the integral of the forces shown in equation \ref{eq:WolfForces}, they obtained a new shifted damped Coulomb potential
119 > that incorporates both image charges and damping of the electrostatic
120 > interaction.
121 >
122 > More recently, Zahn \textit{et al.} investigated these potential and
123 > force expressions for use in simulations involving water.\cite{Zahn02}
124 > In their work, they pointed out that the method that the forces and
125 > derivative of the potential are not commensurate.  Attempts to use
126 > both Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will
127 > lead to poor energy conservation.  They correctly observed that taking
128 > the limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating
129 > the derivatives is mathematically invalid.
130 >
131 > Zahn \textit{et al.} proposed a modified form of this ``Wolf summation
132 > method'' as a way to use this technique in Molecular Dynamics
133 > simulations.  Taking the integral of the forces shown in equation
134 > \ref{eq:WolfForces}, they proposed a new damped Coulomb
135 > potential,
136   \begin{equation}
137   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\}.
138   \label{eq:ZahnPot}
139   \end{equation}
140 < They showed that this new potential does well in capturing the structural and dynamic properties present when using the Ewald sum with the models of water used in their simulations.
140 > They showed that this potential does fairly well at capturing the
141 > structural and dynamic properties of water compared the same
142 > properties obtained using the Ewald sum.
143  
144   \subsection{Simple Forms for Pairwise Electrostatics}
83 – The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et al.} are constructed using two different (and separable) computational tricks: shifting through use of image charges and damping of the electrostatic interaction.  Wolf \textit{et al.} treated the development of their summation method as a progressive application of these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded their shifted force adaptation \ref{eq:ZahnPot} on what they called "the formally incorrect prescription for the derivation of damped Coulomb pair forces".\cite{Zahn02}  Below, we consider the ideas encompassing these electrostatic summation method formulations and clarify their development.
145  
146 < Starting with the original observation that the effective range of the electrostatic interaction in condensed phases is considerably less than the $r^{-1}$ in vacuum, either the shifting or the distance-dependent damping technique could be used as a foundation for the summation method.  Wolf \textit{et al.} made the additional observation that charge neutralization within the cutoff sphere plays a significant role in energy convergence; thus, shifting through the use of image charges was taken as the initial step.  Using these image charges, the electrostatic summation is forced to converge at the cutoff radius.  We can incorporate the methods of Wolf \textit{et al.} and Zahn \textit{et al.} by considering the standard shifted force potential
146 > The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et
147 > al.} are constructed using two different (and separable) computational
148 > tricks: \begin{itemize}
149 > \item shifting through the use of image charges, and
150 > \item damping the electrostatic interaction.
151 > \end{itemize}  Wolf \textit{et al.} treated the
152 > development of their summation method as a progressive application of
153 > these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded
154 > their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the
155 > post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using
156 > both techniques.  It is possible, however, to separate these
157 > tricks and study their effects independently.
158 >
159 > Starting with the original observation that the effective range of the
160 > electrostatic interaction in condensed phases is considerably less
161 > than $r^{-1}$, either the cutoff sphere neutralization or the
162 > distance-dependent damping technique could be used as a foundation for
163 > a new pairwise summation method.  Wolf \textit{et al.} made the
164 > observation that charge neutralization within the cutoff sphere plays
165 > a significant role in energy convergence; therefore we will begin our
166 > analysis with the various shifted forms that maintain this charge
167 > neutralization.  We can evaluate the methods of Wolf
168 > \textit{et al.}  and Zahn \textit{et al.} by considering the standard
169 > shifted potential,
170   \begin{equation}
171 < V^\textrm{SF}(r_{ij}) =         \begin{cases} v(r_{ij})-v_\textrm{c}-\left[\frac{\textrm{d}v(r_{ij})}{\textrm{d}r_{ij}}\right]_{r_{ij}=R_\textrm{c}}(r_{ij}-R_\textrm{c}) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c}
171 > v^\textrm{SP}(r) =      \begin{cases}
172 > v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r >
173 > R_\textrm{c}  
174 > \end{cases},
175 > \label{eq:shiftingPotForm}
176 > \end{equation}
177 > and shifted force,
178 > \begin{equation}
179 > v^\textrm{SF}(r) =      \begin{cases}
180 > v(r)-v_\textrm{c}-\left(\frac{\textrm{d}v(r)}{\textrm{d}r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c})
181 > &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c}
182                                                  \end{cases},
183   \label{eq:shiftingForm}
184   \end{equation}
185 < where $v(r_{ij})$ is the unshifted form of the potential, and $v_c$ is $v(R_\textrm{c})$ and insures the potential goes to zero at the cutoff radius.\cite{Allen87}  If the derivative term is taken to be zero, we are left with the shifted Coulomb potential devised by Wolf \textit{et al.},\cite{Wolf99}
185 > functions where $v(r)$ is the unshifted form of the potential, and
186 > $v_c$ is $v(R_\textrm{c})$.  The Shifted Force ({\sc sf}) form ensures
187 > that both the potential and the forces goes to zero at the cutoff
188 > radius, while the Shifted Potential ({\sc sp}) form only ensures the
189 > potential is smooth at the cutoff radius
190 > ($R_\textrm{c}$).\cite{Allen87}
191 >
192 >
193 >
194 >
195 > 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}
196   \begin{equation}
197   V^\textrm{WSP}(r_{ij}) =        \begin{cases} q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c}
198                                                  \end{cases}.
# Line 162 | Line 266 | Evaluation of the pairwise summation techniques (outli
266  
267   \begin{figure}
268   \centering
269 < \includegraphics[width=3.25in]{./linearFit.pdf}
269 > \includegraphics[width = \linewidth]{./linearFit.pdf}
270   \caption{Example least squares regression of the $\Delta E$ between configurations for the SF method against SPME in the pure water system.  }
271   \label{fig:linearFit}
272   \end{figure}
# Line 180 | Line 284 | Each of these $\theta$ values was accumulated in a dis
284  
285   \begin{figure}
286   \centering
287 < \includegraphics[width=3.25in]{./gaussFit.pdf}
287 > \includegraphics[width = \linewidth]{./gaussFit.pdf}
288   \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.}
289   \label{fig:gaussian}
290   \end{figure}
# Line 217 | Line 321 | Generation of the system configurations was dependent
321  
322   \begin{figure}
323   \centering
324 < \includegraphics[width=3.25in]{./slice.pdf}
324 > \includegraphics[width = \linewidth]{./slice.pdf}
325   \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.}
326   \label{fig:argonSlice}
327   \end{figure}
# Line 234 | Line 338 | In order to evaluate the performance of the pairwise e
338  
339   \begin{figure}
340   \centering
341 < \includegraphics[width=3.25in]{./delEplot.pdf}
341 > \includegraphics[width=5.5in]{./delEplot.pdf}
342   \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).}
343   \label{fig:delE}
344   \end{figure}
# Line 249 | Line 353 | Evaluation of pairwise methods for use in Molecular Dy
353  
354   \begin{figure}
355   \centering
356 < \includegraphics[width=3.25in]{./frcMagplot.pdf}
356 > \includegraphics[width=5.5in]{./frcMagplot.pdf}
357   \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).}
358   \label{fig:frcMag}
359   \end{figure}
# Line 258 | Line 362 | Figure \ref{fig:frcMag}, for the most part, parallels
362  
363   \begin{figure}
364   \centering
365 < \includegraphics[width=3.25in]{./trqMagplot.pdf}
365 > \includegraphics[width=5.5in]{./trqMagplot.pdf}
366   \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).}
367   \label{fig:trqMag}
368   \end{figure}
# Line 271 | Line 375 | Having force and torque vectors with magnitudes that a
375  
376   \begin{figure}
377   \centering
378 < \includegraphics[width=3.25in]{./frcTrqAngplot.pdf}
378 > \includegraphics[width=5.5in]{./frcTrqAngplot.pdf}
379   \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).}
380   \label{fig:frcTrqAng}
381   \end{figure}
# Line 321 | Line 425 | In the previous studies using a Shifted-Force variant
425  
426   \begin{figure}
427   \centering
428 < \includegraphics[width = 3.25in]{./spectraSquare.pdf}
428 > \includegraphics[width = \linewidth]{./spectraSquare.pdf}
429   \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.}
430   \label{fig:methodPS}
431   \end{figure}
# Line 333 | Line 437 | where $S(r)$ is a switching function that smoothly zer
437   where $S(r)$ is a switching function that smoothly zeroes the potential at the cutoff radius.  Figure \ref{fig:dampInc} shows how the low frequency motions are dependent on the damping used in the direct electrostatic sum.  As the damping increases, the peaks drop to lower frequencies.  Incidentally, use of an $\alpha$ of 0.25 \AA$^{-1}$ on a simple electrostatic summation results in low frequency correlated dynamics equivalent to a simulation using SPME.  When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks shift to higher frequency in exponential fashion.  Though not shown, the spectrum for the simple undamped electrostatic potential is blue-shifted such that the lowest frequency peak resides near 325 cm$^{-1}$.  In light of these results, the undamped Shifted-Force method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite respectable; however, it appears as though moderate damping is required for accurate reproduction of crystal dynamics.
438   \begin{figure}
439   \centering
440 < \includegraphics[width = 3.25in]{./comboSquare.pdf}
440 > \includegraphics[width = \linewidth]{./comboSquare.pdf}
441   \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.}
442   \label{fig:dampInc}
443   \end{figure}
# Line 347 | Line 451 | We are not suggesting any flaw with the Ewald sum; in
451   We are not suggesting any flaw with the Ewald sum; in fact, it is the standard by which these simple pairwise sums are judged.  However, these results do suggest that in the typical simulations performed today, the Ewald summation may no longer be required to obtain the level of accuracy most researcher have come to expect
452  
453   \section{Acknowledgments}
350 –
454   \newpage
455  
456 < \bibliographystyle{achemso}
456 > \bibliographystyle{jcp2}
457   \bibliography{electrostaticMethods}
458  
459  

Diff Legend

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