--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/03/22 21:20:40 2659 +++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/05/02 13:11:41 2744 @@ -20,14 +20,22 @@ \topmargin -21pt \headsep 10pt \textheight 9.0in \textwidth 6.5in \brokenpenalty=10000 -\renewcommand{\baselinestretch}{1.2} +%\renewcommand{\baselinestretch}{1.2} +\renewcommand{\baselinestretch}{2} \renewcommand\citemid{\ } % no comma in optional reference note +\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{2}} %doublespace captions +\let\Caption\caption +\renewcommand\caption[1]{% + \Caption[#1]{}% +} +\setlength{\abovecaptionskip}{1.2in} +\setlength{\belowcaptionskip}{1.2in} \begin{document} -\title{Is the Ewald Summation necessary? \\ -Pairwise alternatives to the accepted standard for \\ -long-range electrostatics} +\title{Is the Ewald summation still necessary? \\ +Pairwise alternatives to the accepted standard \\ +for long-range electrostatics} \author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: gezelter@nd.edu} \\ @@ -38,22 +46,21 @@ Notre Dame, Indiana 46556} \date{\today} \maketitle -\doublespacing +%\doublespacing -\nobibliography{} \begin{abstract} We investigate pairwise electrostatic interaction methods and show that there are viable and computationally efficient $(\mathscr{O}(N))$ alternatives to the Ewald summation for typical modern molecular simulations. These methods are extended from the damped and -cutoff-neutralized Coulombic sum originally proposed by Wolf -\textit{et al.} One of these, the damped shifted force method, shows +cutoff-neutralized Coulombic sum originally proposed by +[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 a remarkable ability to reproduce the energetic and dynamic characteristics exhibited by simulations employing lattice summation techniques. Comparisons were performed with this and other pairwise -methods against the smooth particle mesh Ewald ({\sc spme}) summation to see -how well they reproduce the energetics and dynamics of a variety of -simulation types. +methods against the smooth particle mesh Ewald ({\sc spme}) summation +to see how well they reproduce the energetics and dynamics of a +variety of molecular simulations. \end{abstract} \newpage @@ -96,7 +103,7 @@ explicit Ewald summation.\cite{Tobias01} regarding possible artifacts caused by the inherent periodicity of the explicit Ewald summation.\cite{Tobias01} -In this paper, we focus on a new set of shifted methods devised by +In this paper, we focus on a new set of pairwise methods devised by Wolf {\it et al.},\cite{Wolf99} which we further extend. These methods along with a few other mixed methods (i.e. reaction field) are compared with the smooth particle mesh Ewald @@ -107,11 +114,11 @@ or which have one- or two-dimensional periodicity. Be to the direct pairwise sum. They also lack the added periodicity of the Ewald sum, so they can be used for systems which are non-periodic or which have one- or two-dimensional periodicity. Below, these -methods are evaluated using a variety of model systems to establish -their usability in molecular simulations. +methods are evaluated using a variety of model systems to +establish their usability in molecular simulations. \subsection{The Ewald Sum} -The complete accumulation electrostatic interactions in a system with +The complete accumulation of the electrostatic interactions in a system with periodic boundary conditions (PBC) requires the consideration of the effect of all charges within a (cubic) simulation box as well as those in the periodic replicas, @@ -168,11 +175,11 @@ portion.\cite{Karasawa89,Kolafa92} \begin{figure} \centering \includegraphics[width = \linewidth]{./ewaldProgression.pdf} -\caption{The change in the application of the Ewald sum with -increasing computational power. Initially, only small systems could -be studied, and the Ewald sum replicated the simulation box to -convergence. Now, much larger systems of charges are investigated -with fixed-distance cutoffs.} +\caption{The change in the need for the Ewald sum with +increasing computational power. A:~Initially, only small systems +could be studied, and the Ewald sum replicated the simulation box to +convergence. B:~Now, radial cutoff methods should be able to reach +convergence for the larger systems of charges that are common today.} \label{fig:ewaldTime} \end{figure} @@ -201,12 +208,14 @@ can prove problematic. The Ewald sum has been reformu conditions. However, in certain systems, such as vapor-liquid interfaces and membranes, the intrinsic three-dimensional periodicity can prove problematic. The Ewald sum has been reformulated to handle -2D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the -new methods are computationally expensive.\cite{Spohr97,Yeh99} -Inclusion of a correction term in the Ewald summation is a possible -direction for handling 2D systems while still enabling the use of the -modern optimizations.\cite{Yeh99} +2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} but these +methods are computationally expensive.\cite{Spohr97,Yeh99} More +recently, there have been several successful efforts toward reducing +the computational cost of 2-D lattice +summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} +bringing them more in line with the cost of the full 3-D summation. + Several studies have recognized that the inherent periodicity in the Ewald sum can also have an effect on three-dimensional systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} @@ -228,11 +237,10 @@ charge neutrality and gives results similar to those o charge contained within the cutoff radius is crucial for potential stability. They devised a pairwise summation method that ensures charge neutrality and gives results similar to those obtained with 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 convergence +Ewald summation. The resulting shifted Coulomb potential 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 convergence \begin{equation} 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\}. \label{eq:WolfPot} @@ -540,16 +548,16 @@ shows a data set with a good correlation coefficient.} \label{fig:linearFit} \end{figure} -Each system type (detailed in section \ref{sec:RepSims}) was -represented using 500 independent configurations. Additionally, we -used seven different system types, so each of the alternative -(non-Ewald) electrostatic summation methods was evaluated using -873,250 configurational energy differences. +Each of the seven system types (detailed in section \ref{sec:RepSims}) +were represented using 500 independent configurations. Thus, each of +the alternative (non-Ewald) electrostatic summation methods was +evaluated using an accumulated 873,250 configurational energy +differences. Results and discussion for the individual analysis of each of the -system types appear in the supporting information, while the -cumulative results over all the investigated systems appears below in -section \ref{sec:EnergyResults}. +system types appear in the supporting information,\cite{EPAPSdeposit} +while the cumulative results over all the investigated systems appears +below in section \ref{sec:EnergyResults}. \subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} We evaluated the pairwise methods (outlined in section @@ -581,29 +589,11 @@ shape. Thus, gaussian fits were used to measure the wi between two different electrostatic summation methods, there is no {\it a priori} reason for the profile to adhere to any specific shape. Thus, gaussian fits were used to measure the width of the -resulting distributions. -% -%\begin{figure} -%\centering -%\includegraphics[width = \linewidth]{./gaussFit.pdf} -%\caption{Sample fit of the angular distribution of the force vectors -%accumulated using all of the studied systems. Gaussian fits were used -%to obtain values for the variance in force and torque vectors.} -%\label{fig:gaussian} -%\end{figure} -% -%Figure \ref{fig:gaussian} shows an example distribution with applied -%non-linear fits. The solid line is a Gaussian profile, while the -%dotted line is a Voigt profile, a convolution of a Gaussian and a -%Lorentzian. -%Since this distribution is a measure of angular error between two -%different electrostatic summation methods, there is no {\it a priori} -%reason for the profile to adhere to any specific shape. -%Gaussian fits was used to compare all the tested methods. -The variance ($\sigma^2$) was extracted from each of these fits and -was used to compare distribution widths. Values of $\sigma^2$ near -zero indicate vector directions indistinguishable from those -calculated when using the reference method ({\sc spme}). +resulting distributions. The variance ($\sigma^2$) was extracted from +each of these fits and was used to compare distribution widths. +Values of $\sigma^2$ near zero indicate vector directions +indistinguishable from those calculated when using the reference +method ({\sc spme}). \subsection{Short-time Dynamics} @@ -631,7 +621,7 @@ the {\it long-time} dynamics of charged systems were e The effects of the same subset of alternative electrostatic methods on the {\it long-time} dynamics of charged systems were evaluated using -the same model system (NaCl crystals at 1000K). The power spectrum +the same model system (NaCl crystals at 1000~K). The power spectrum ($I(\omega)$) was obtained via Fourier transform of the velocity autocorrelation function, \begin{equation} I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, @@ -641,14 +631,15 @@ were performed under the microcanonical ensemble, and NaCl crystal is composed of two different atom types, the average of the two resulting power spectra was used for comparisons. Simulations were performed under the microcanonical ensemble, and velocity -information was saved every 5 fs over 100 ps trajectories. +information was saved every 5~fs over 100~ps trajectories. \subsection{Representative Simulations}\label{sec:RepSims} -A variety of representative simulations were analyzed to determine the -relative effectiveness of the pairwise summation techniques in -reproducing the energetics and dynamics exhibited by {\sc spme}. We wanted -to span the space of modern simulations (i.e. from liquids of neutral -molecules to ionic crystals), so the systems studied were: +A variety of representative molecular simulations were analyzed to +determine the relative effectiveness of the pairwise summation +techniques in reproducing the energetics and dynamics exhibited by +{\sc spme}. We wanted to span the space of typical molecular +simulations (i.e. from liquids of neutral molecules to ionic +crystals), so the systems studied were: \begin{enumerate} \item liquid water (SPC/E),\cite{Berendsen87} \item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), @@ -671,30 +662,24 @@ these systems were selected and equilibrated in the sa the crystal). The solid and liquid NaCl systems consisted of 500 $\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for these systems were selected and equilibrated in the same manner as the -water systems. The equilibrated temperatures were 1000~K for the NaCl -crystal and 7000~K for the liquid. The ionic solutions were made by -solvating 4 (or 40) ions in a periodic box containing 1000 SPC/E water -molecules. Ion and water positions were then randomly swapped, and -the resulting configurations were again equilibrated individually. -Finally, for the Argon / Water ``charge void'' systems, the identities -of all the SPC/E waters within 6 \AA\ of the center of the -equilibrated water configurations were converted to argon. -%(Fig. \ref{fig:argonSlice}). +water systems. In order to introduce measurable fluctuations in the +configuration energy differences, the crystalline simulations were +equilibrated at 1000~K, near the $T_\textrm{m}$ for NaCl. The liquid +NaCl configurations needed to represent a fully disordered array of +point charges, so the high temperature of 7000~K was selected for +equilibration. The ionic solutions were made by solvating 4 (or 40) +ions in a periodic box containing 1000 SPC/E water molecules. Ion and +water positions were then randomly swapped, and the resulting +configurations were again equilibrated individually. Finally, for the +Argon / Water ``charge void'' systems, the identities of all the SPC/E +waters within 6 \AA\ of the center of the equilibrated water +configurations were converted to argon. These procedures guaranteed us a set of representative configurations from chemically-relevant systems sampled from appropriate ensembles. Force field parameters for the ions and Argon were taken from the force field utilized by {\sc oopse}.\cite{Meineke05} -%\begin{figure} -%\centering -%\includegraphics[width = \linewidth]{./slice.pdf} -%\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.} -%\label{fig:argonSlice} -%\end{figure} - \subsection{Comparison of Summation Methods}\label{sec:ESMethods} We compared the following alternative summation methods with results from the reference method ({\sc spme}): @@ -716,7 +701,7 @@ manner across all systems and configurations. (i.e. Lennard-Jones interactions) were handled in exactly the same manner across all systems and configurations. -The althernative methods were also evaluated with three different +The alternative methods were also evaluated with three different cutoff radii (9, 12, and 15 \AA). As noted previously, the convergence parameter ($\alpha$) plays a role in the balance of the real-space and reciprocal-space portions of the Ewald calculation. @@ -768,7 +753,7 @@ readers can consult the accompanying supporting inform significant improvement using the group-switched cutoff because the salt and salt solution systems contain non-neutral groups. Interested readers can consult the accompanying supporting information for a -comparison where all groups are neutral. +comparison where all groups are neutral.\cite{EPAPSdeposit} For the {\sc sp} method, inclusion of electrostatic damping improves the agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ @@ -916,7 +901,7 @@ charged bodies, and this observation is investigated f particles in all seven systems, while torque vectors are only available for neutral molecular groups. Damping is more beneficial to charged bodies, and this observation is investigated further in the -accompanying supporting information. +accompanying supporting information.\cite{EPAPSdeposit} Although not discussed previously, group based cutoffs can be applied to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs @@ -1096,7 +1081,7 @@ cutoff distance. \centering \includegraphics[width = \linewidth]{./increasedDamping.pdf} \caption{Effect of damping on the two lowest-frequency phonon modes in -the NaCl crystal at 1000K. The undamped shifted force ({\sc sf}) +the NaCl crystal at 1000~K. The undamped shifted force ({\sc sf}) method is off by less than 10 cm$^{-1}$, and increasing the electrostatic damping to 0.25 \AA$^{-1}$ gives quantitative agreement with the power spectrum obtained using the Ewald sum. Overdamping can