| 20 |
|
\topmargin -21pt \headsep 10pt |
| 21 |
|
\textheight 9.0in \textwidth 6.5in |
| 22 |
|
\brokenpenalty=10000 |
| 23 |
< |
\renewcommand{\baselinestretch}{1.2} |
| 23 |
> |
%\renewcommand{\baselinestretch}{1.2} |
| 24 |
> |
\renewcommand{\baselinestretch}{2} |
| 25 |
|
\renewcommand\citemid{\ } % no comma in optional reference note |
| 26 |
+ |
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{2}} %doublespace captions |
| 27 |
+ |
\let\Caption\caption |
| 28 |
+ |
\renewcommand\caption[1]{% |
| 29 |
+ |
\Caption[#1]{}% |
| 30 |
+ |
} |
| 31 |
+ |
\setlength{\abovecaptionskip}{1.2in} |
| 32 |
+ |
\setlength{\belowcaptionskip}{1.2in} |
| 33 |
|
|
| 34 |
|
\begin{document} |
| 35 |
|
|
| 36 |
< |
\title{Is the Ewald Summation necessary? \\ |
| 37 |
< |
Pairwise alternatives to the accepted standard for \\ |
| 38 |
< |
long-range electrostatics} |
| 36 |
> |
\title{Is the Ewald summation still necessary? \\ |
| 37 |
> |
Pairwise alternatives to the accepted standard \\ |
| 38 |
> |
for long-range electrostatics} |
| 39 |
|
|
| 40 |
|
\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
| 41 |
|
gezelter@nd.edu} \\ |
| 46 |
|
\date{\today} |
| 47 |
|
|
| 48 |
|
\maketitle |
| 49 |
< |
\doublespacing |
| 49 |
> |
%\doublespacing |
| 50 |
|
|
| 43 |
– |
\nobibliography{} |
| 51 |
|
\begin{abstract} |
| 52 |
|
We investigate pairwise electrostatic interaction methods and show |
| 53 |
|
that there are viable and computationally efficient $(\mathscr{O}(N))$ |
| 54 |
|
alternatives to the Ewald summation for typical modern molecular |
| 55 |
|
simulations. These methods are extended from the damped and |
| 56 |
< |
cutoff-neutralized Coulombic sum originally proposed by Wolf |
| 57 |
< |
\textit{et al.} One of these, the damped shifted force method, shows |
| 56 |
> |
cutoff-neutralized Coulombic sum originally proposed by |
| 57 |
> |
[D. Wolf, P. Keblinski, S.~R. Phillpot, and J. Eggebrecht, {\it J. Chem. Phys.} {\bf 110}, 8255 (1999)] One of these, the damped shifted force method, shows |
| 58 |
|
a remarkable ability to reproduce the energetic and dynamic |
| 59 |
|
characteristics exhibited by simulations employing lattice summation |
| 60 |
|
techniques. Comparisons were performed with this and other pairwise |
| 61 |
< |
methods against the smooth particle mesh Ewald ({\sc spme}) summation to see |
| 62 |
< |
how well they reproduce the energetics and dynamics of a variety of |
| 63 |
< |
simulation types. |
| 61 |
> |
methods against the smooth particle mesh Ewald ({\sc spme}) summation |
| 62 |
> |
to see how well they reproduce the energetics and dynamics of a |
| 63 |
> |
variety of molecular simulations. |
| 64 |
|
\end{abstract} |
| 65 |
|
|
| 66 |
|
\newpage |
| 103 |
|
regarding possible artifacts caused by the inherent periodicity of the |
| 104 |
|
explicit Ewald summation.\cite{Tobias01} |
| 105 |
|
|
| 106 |
< |
In this paper, we focus on a new set of shifted methods devised by |
| 106 |
> |
In this paper, we focus on a new set of pairwise methods devised by |
| 107 |
|
Wolf {\it et al.},\cite{Wolf99} which we further extend. These |
| 108 |
|
methods along with a few other mixed methods (i.e. reaction field) are |
| 109 |
|
compared with the smooth particle mesh Ewald |
| 114 |
|
to the direct pairwise sum. They also lack the added periodicity of |
| 115 |
|
the Ewald sum, so they can be used for systems which are non-periodic |
| 116 |
|
or which have one- or two-dimensional periodicity. Below, these |
| 117 |
< |
methods are evaluated using a variety of model systems to establish |
| 118 |
< |
their usability in molecular simulations. |
| 117 |
> |
methods are evaluated using a variety of model systems to |
| 118 |
> |
establish their usability in molecular simulations. |
| 119 |
|
|
| 120 |
|
\subsection{The Ewald Sum} |
| 121 |
< |
The complete accumulation electrostatic interactions in a system with |
| 121 |
> |
The complete accumulation of the electrostatic interactions in a system with |
| 122 |
|
periodic boundary conditions (PBC) requires the consideration of the |
| 123 |
|
effect of all charges within a (cubic) simulation box as well as those |
| 124 |
|
in the periodic replicas, |
| 175 |
|
\begin{figure} |
| 176 |
|
\centering |
| 177 |
|
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 178 |
< |
\caption{The change in the application of the Ewald sum with |
| 179 |
< |
increasing computational power. Initially, only small systems could |
| 180 |
< |
be studied, and the Ewald sum replicated the simulation box to |
| 181 |
< |
convergence. Now, much larger systems of charges are investigated |
| 182 |
< |
with fixed-distance cutoffs.} |
| 178 |
> |
\caption{The change in the need for the Ewald sum with |
| 179 |
> |
increasing computational power. A:~Initially, only small systems |
| 180 |
> |
could be studied, and the Ewald sum replicated the simulation box to |
| 181 |
> |
convergence. B:~Now, radial cutoff methods should be able to reach |
| 182 |
> |
convergence for the larger systems of charges that are common today.} |
| 183 |
|
\label{fig:ewaldTime} |
| 184 |
|
\end{figure} |
| 185 |
|
|
| 208 |
|
conditions. However, in certain systems, such as vapor-liquid |
| 209 |
|
interfaces and membranes, the intrinsic three-dimensional periodicity |
| 210 |
|
can prove problematic. The Ewald sum has been reformulated to handle |
| 211 |
< |
2D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the |
| 212 |
< |
new methods are computationally expensive.\cite{Spohr97,Yeh99} |
| 213 |
< |
Inclusion of a correction term in the Ewald summation is a possible |
| 214 |
< |
direction for handling 2D systems while still enabling the use of the |
| 215 |
< |
modern optimizations.\cite{Yeh99} |
| 211 |
> |
2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} but these |
| 212 |
> |
methods are computationally expensive.\cite{Spohr97,Yeh99} More |
| 213 |
> |
recently, there have been several successful efforts toward reducing |
| 214 |
> |
the computational cost of 2-D lattice |
| 215 |
> |
summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} |
| 216 |
> |
bringing them more in line with the cost of the full 3-D summation. |
| 217 |
|
|
| 218 |
+ |
|
| 219 |
|
Several studies have recognized that the inherent periodicity in the |
| 220 |
|
Ewald sum can also have an effect on three-dimensional |
| 221 |
|
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 237 |
|
charge contained within the cutoff radius is crucial for potential |
| 238 |
|
stability. They devised a pairwise summation method that ensures |
| 239 |
|
charge neutrality and gives results similar to those obtained with the |
| 240 |
< |
Ewald summation. The resulting shifted Coulomb potential |
| 241 |
< |
(Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through |
| 242 |
< |
placement on the cutoff sphere and a distance-dependent damping |
| 243 |
< |
function (identical to that seen in the real-space portion of the |
| 235 |
< |
Ewald sum) to aid convergence |
| 240 |
> |
Ewald summation. The resulting shifted Coulomb potential includes |
| 241 |
> |
image-charges subtracted out through placement on the cutoff sphere |
| 242 |
> |
and a distance-dependent damping function (identical to that seen in |
| 243 |
> |
the real-space portion of the Ewald sum) to aid convergence |
| 244 |
|
\begin{equation} |
| 245 |
|
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\}. |
| 246 |
|
\label{eq:WolfPot} |
| 548 |
|
\label{fig:linearFit} |
| 549 |
|
\end{figure} |
| 550 |
|
|
| 551 |
< |
Each system type (detailed in section \ref{sec:RepSims}) was |
| 552 |
< |
represented using 500 independent configurations. Additionally, we |
| 553 |
< |
used seven different system types, so each of the alternative |
| 554 |
< |
(non-Ewald) electrostatic summation methods was evaluated using |
| 555 |
< |
873,250 configurational energy differences. |
| 551 |
> |
Each of the seven system types (detailed in section \ref{sec:RepSims}) |
| 552 |
> |
were represented using 500 independent configurations. Thus, each of |
| 553 |
> |
the alternative (non-Ewald) electrostatic summation methods was |
| 554 |
> |
evaluated using an accumulated 873,250 configurational energy |
| 555 |
> |
differences. |
| 556 |
|
|
| 557 |
|
Results and discussion for the individual analysis of each of the |
| 558 |
< |
system types appear in the supporting information, while the |
| 559 |
< |
cumulative results over all the investigated systems appears below in |
| 560 |
< |
section \ref{sec:EnergyResults}. |
| 558 |
> |
system types appear in the supporting information,\cite{EPAPSdeposit} |
| 559 |
> |
while the cumulative results over all the investigated systems appears |
| 560 |
> |
below in section \ref{sec:EnergyResults}. |
| 561 |
|
|
| 562 |
|
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
| 563 |
|
We evaluated the pairwise methods (outlined in section |
| 589 |
|
between two different electrostatic summation methods, there is no |
| 590 |
|
{\it a priori} reason for the profile to adhere to any specific |
| 591 |
|
shape. Thus, gaussian fits were used to measure the width of the |
| 592 |
< |
resulting distributions. |
| 593 |
< |
% |
| 594 |
< |
%\begin{figure} |
| 595 |
< |
%\centering |
| 596 |
< |
%\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 589 |
< |
%\caption{Sample fit of the angular distribution of the force vectors |
| 590 |
< |
%accumulated using all of the studied systems. Gaussian fits were used |
| 591 |
< |
%to obtain values for the variance in force and torque vectors.} |
| 592 |
< |
%\label{fig:gaussian} |
| 593 |
< |
%\end{figure} |
| 594 |
< |
% |
| 595 |
< |
%Figure \ref{fig:gaussian} shows an example distribution with applied |
| 596 |
< |
%non-linear fits. The solid line is a Gaussian profile, while the |
| 597 |
< |
%dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 598 |
< |
%Lorentzian. |
| 599 |
< |
%Since this distribution is a measure of angular error between two |
| 600 |
< |
%different electrostatic summation methods, there is no {\it a priori} |
| 601 |
< |
%reason for the profile to adhere to any specific shape. |
| 602 |
< |
%Gaussian fits was used to compare all the tested methods. |
| 603 |
< |
The variance ($\sigma^2$) was extracted from each of these fits and |
| 604 |
< |
was used to compare distribution widths. Values of $\sigma^2$ near |
| 605 |
< |
zero indicate vector directions indistinguishable from those |
| 606 |
< |
calculated when using the reference method ({\sc spme}). |
| 592 |
> |
resulting distributions. The variance ($\sigma^2$) was extracted from |
| 593 |
> |
each of these fits and was used to compare distribution widths. |
| 594 |
> |
Values of $\sigma^2$ near zero indicate vector directions |
| 595 |
> |
indistinguishable from those calculated when using the reference |
| 596 |
> |
method ({\sc spme}). |
| 597 |
|
|
| 598 |
|
\subsection{Short-time Dynamics} |
| 599 |
|
|
| 621 |
|
|
| 622 |
|
The effects of the same subset of alternative electrostatic methods on |
| 623 |
|
the {\it long-time} dynamics of charged systems were evaluated using |
| 624 |
< |
the same model system (NaCl crystals at 1000K). The power spectrum |
| 624 |
> |
the same model system (NaCl crystals at 1000~K). The power spectrum |
| 625 |
|
($I(\omega)$) was obtained via Fourier transform of the velocity |
| 626 |
|
autocorrelation function, \begin{equation} I(\omega) = |
| 627 |
|
\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 631 |
|
NaCl crystal is composed of two different atom types, the average of |
| 632 |
|
the two resulting power spectra was used for comparisons. Simulations |
| 633 |
|
were performed under the microcanonical ensemble, and velocity |
| 634 |
< |
information was saved every 5 fs over 100 ps trajectories. |
| 634 |
> |
information was saved every 5~fs over 100~ps trajectories. |
| 635 |
|
|
| 636 |
|
\subsection{Representative Simulations}\label{sec:RepSims} |
| 637 |
< |
A variety of representative simulations were analyzed to determine the |
| 638 |
< |
relative effectiveness of the pairwise summation techniques in |
| 639 |
< |
reproducing the energetics and dynamics exhibited by {\sc spme}. We wanted |
| 640 |
< |
to span the space of modern simulations (i.e. from liquids of neutral |
| 641 |
< |
molecules to ionic crystals), so the systems studied were: |
| 637 |
> |
A variety of representative molecular simulations were analyzed to |
| 638 |
> |
determine the relative effectiveness of the pairwise summation |
| 639 |
> |
techniques in reproducing the energetics and dynamics exhibited by |
| 640 |
> |
{\sc spme}. We wanted to span the space of typical molecular |
| 641 |
> |
simulations (i.e. from liquids of neutral molecules to ionic |
| 642 |
> |
crystals), so the systems studied were: |
| 643 |
|
\begin{enumerate} |
| 644 |
|
\item liquid water (SPC/E),\cite{Berendsen87} |
| 645 |
|
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), |
| 662 |
|
the crystal). The solid and liquid NaCl systems consisted of 500 |
| 663 |
|
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
| 664 |
|
these systems were selected and equilibrated in the same manner as the |
| 665 |
< |
water systems. The equilibrated temperatures were 1000~K for the NaCl |
| 666 |
< |
crystal and 7000~K for the liquid. The ionic solutions were made by |
| 667 |
< |
solvating 4 (or 40) ions in a periodic box containing 1000 SPC/E water |
| 668 |
< |
molecules. Ion and water positions were then randomly swapped, and |
| 669 |
< |
the resulting configurations were again equilibrated individually. |
| 670 |
< |
Finally, for the Argon / Water ``charge void'' systems, the identities |
| 671 |
< |
of all the SPC/E waters within 6 \AA\ of the center of the |
| 672 |
< |
equilibrated water configurations were converted to argon. |
| 673 |
< |
%(Fig. \ref{fig:argonSlice}). |
| 665 |
> |
water systems. In order to introduce measurable fluctuations in the |
| 666 |
> |
configuration energy differences, the crystalline simulations were |
| 667 |
> |
equilibrated at 1000~K, near the $T_\textrm{m}$ for NaCl. The liquid |
| 668 |
> |
NaCl configurations needed to represent a fully disordered array of |
| 669 |
> |
point charges, so the high temperature of 7000~K was selected for |
| 670 |
> |
equilibration. The ionic solutions were made by solvating 4 (or 40) |
| 671 |
> |
ions in a periodic box containing 1000 SPC/E water molecules. Ion and |
| 672 |
> |
water positions were then randomly swapped, and the resulting |
| 673 |
> |
configurations were again equilibrated individually. Finally, for the |
| 674 |
> |
Argon / Water ``charge void'' systems, the identities of all the SPC/E |
| 675 |
> |
waters within 6 \AA\ of the center of the equilibrated water |
| 676 |
> |
configurations were converted to argon. |
| 677 |
|
|
| 678 |
|
These procedures guaranteed us a set of representative configurations |
| 679 |
|
from chemically-relevant systems sampled from appropriate |
| 680 |
|
ensembles. Force field parameters for the ions and Argon were taken |
| 681 |
|
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 682 |
|
|
| 689 |
– |
%\begin{figure} |
| 690 |
– |
%\centering |
| 691 |
– |
%\includegraphics[width = \linewidth]{./slice.pdf} |
| 692 |
– |
%\caption{A slice from the center of a water box used in a charge void |
| 693 |
– |
%simulation. The darkened region represents the boundary sphere within |
| 694 |
– |
%which the water molecules were converted to argon atoms.} |
| 695 |
– |
%\label{fig:argonSlice} |
| 696 |
– |
%\end{figure} |
| 697 |
– |
|
| 683 |
|
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 684 |
|
We compared the following alternative summation methods with results |
| 685 |
|
from the reference method ({\sc spme}): |
| 701 |
|
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 702 |
|
manner across all systems and configurations. |
| 703 |
|
|
| 704 |
< |
The althernative methods were also evaluated with three different |
| 704 |
> |
The alternative methods were also evaluated with three different |
| 705 |
|
cutoff radii (9, 12, and 15 \AA). As noted previously, the |
| 706 |
|
convergence parameter ($\alpha$) plays a role in the balance of the |
| 707 |
|
real-space and reciprocal-space portions of the Ewald calculation. |
| 753 |
|
significant improvement using the group-switched cutoff because the |
| 754 |
|
salt and salt solution systems contain non-neutral groups. Interested |
| 755 |
|
readers can consult the accompanying supporting information for a |
| 756 |
< |
comparison where all groups are neutral. |
| 756 |
> |
comparison where all groups are neutral.\cite{EPAPSdeposit} |
| 757 |
|
|
| 758 |
|
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 759 |
|
the agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ |
| 901 |
|
particles in all seven systems, while torque vectors are only |
| 902 |
|
available for neutral molecular groups. Damping is more beneficial to |
| 903 |
|
charged bodies, and this observation is investigated further in the |
| 904 |
< |
accompanying supporting information. |
| 904 |
> |
accompanying supporting information.\cite{EPAPSdeposit} |
| 905 |
|
|
| 906 |
|
Although not discussed previously, group based cutoffs can be applied |
| 907 |
|
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 1081 |
|
\centering |
| 1082 |
|
\includegraphics[width = \linewidth]{./increasedDamping.pdf} |
| 1083 |
|
\caption{Effect of damping on the two lowest-frequency phonon modes in |
| 1084 |
< |
the NaCl crystal at 1000K. The undamped shifted force ({\sc sf}) |
| 1084 |
> |
the NaCl crystal at 1000~K. The undamped shifted force ({\sc sf}) |
| 1085 |
|
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1086 |
|
electrostatic damping to 0.25 \AA$^{-1}$ gives quantitative agreement |
| 1087 |
|
with the power spectrum obtained using the Ewald sum. Overdamping can |
