--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/03/22 21:20:40 2659 +++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/03/24 16:54:13 2669 @@ -25,7 +25,7 @@ \begin{document} -\title{Is the Ewald Summation necessary? \\ +\title{Is the Ewald summation still necessary? \\ Pairwise alternatives to the accepted standard for \\ long-range electrostatics} @@ -40,20 +40,19 @@ Notre Dame, Indiana 46556} \maketitle \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 simulation types. \end{abstract} \newpage @@ -96,7 +95,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 @@ -111,7 +110,7 @@ their usability in molecular simulations. 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 +167,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} @@ -228,11 +227,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} @@ -581,29 +579,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 +611,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, @@ -671,29 +651,23 @@ 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 @@ -716,7 +690,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. @@ -1096,7 +1070,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