--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/22 20:49:53	2656
+++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/04/26 14:20:17	2740
@@ -1,7 +1,7 @@
 %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
 %\documentclass[aps,prb,preprint]{revtex4}
-\documentclass[11pt]{article}
-\usepackage{endfloat}
+\documentclass[10pt]{article}
+%\usepackage{endfloat}
 \usepackage{amsmath,bm}
 \usepackage{amssymb}
 \usepackage{epsf}
@@ -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}
 
@@ -38,22 +38,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 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
@@ -107,11 +106,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 +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}
 
@@ -202,10 +201,11 @@ can prove problematic.  The Ewald sum has been reformu
 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}
+new methods are computationally expensive.\cite{Spohr97,Yeh99} More
+recently, there have been several successful efforts toward reducing
+the computational cost of 2D lattice summations, often enabling the
+use of the mentioned
+optimizations.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04}
 
 Several studies have recognized that the inherent periodicity in the
 Ewald sum can also have an effect on three-dimensional
@@ -228,11 +228,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}
@@ -532,7 +531,7 @@ Fig. \ref{fig:linearFit}.
 
 \begin{figure}
 \centering
-\includegraphics[width = \linewidth]{./dualLinear.pdf}
+\includegraphics[width = 3.25in]{./dualLinear.pdf}
 \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
@@ -540,11 +539,11 @@ 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
@@ -581,29 +580,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 +612,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 +622,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 +653,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 +692,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.
@@ -741,7 +717,7 @@ figure \ref{fig:delE}.
 
 \begin{figure}
 \centering
-\includegraphics[width=5.5in]{./delEplot.pdf}
+\includegraphics[width=3.25in]{./delEplot.pdf}
 \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)
@@ -801,7 +777,7 @@ accumulation of the statistics from all of the system 
 
 \begin{figure}
 \centering
-\includegraphics[width=5.5in]{./frcMagplot.pdf}
+\includegraphics[width=3.25in]{./frcMagplot.pdf}
 \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)
@@ -840,7 +816,7 @@ performs more favorably.
 
 \begin{figure}
 \centering
-\includegraphics[width=5.5in]{./trqMagplot.pdf}
+\includegraphics[width=3.25in]{./trqMagplot.pdf}
 \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)
@@ -880,7 +856,7 @@ distributions of the combined set over all system type
 
 \begin{figure}
 \centering
-\includegraphics[width=5.5in]{./frcTrqAngplot.pdf}
+\includegraphics[width=3.25in]{./frcTrqAngplot.pdf}
 \caption{Statistical analysis of the width of the angular distribution
 that the force and torque vectors from a given electrostatic method
 make with their counterparts obtained using the reference Ewald sum.
@@ -1015,7 +991,7 @@ low-frequency portion of the power spectrum.
 
 \begin{figure}
 \centering
-\includegraphics[width = \linewidth]{./vCorrPlot.pdf}
+\includegraphics[width = 3.25in]{./vCorrPlot.pdf}
 \caption{Velocity autocorrelation functions of NaCl crystals at 
 1000 K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc
 sp} ($\alpha$ = 0.2). The inset is a magnification of the area around
@@ -1055,7 +1031,7 @@ noticeable effect on peak location or magnitude.
 
 \begin{figure}
 \centering
-\includegraphics[width = \linewidth]{./spectraSquare.pdf}
+\includegraphics[width = 3.25in]{./spectraSquare.pdf}
 \caption{Power spectra obtained from the velocity auto-correlation 
 functions of NaCl crystals at 1000 K while using {\sc spme}, {\sc sf}
 ($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2).  The inset
@@ -1094,9 +1070,9 @@ cutoff distance.
 
 \begin{figure}
 \centering
-\includegraphics[width = \linewidth]{./comboSquare.pdf}
+\includegraphics[width = 3.25in]{./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
