--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/14 13:52:06	2613
+++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex	2006/03/14 16:22:54	2617
@@ -1,18 +1,18 @@
 %\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
-\documentclass[12pt]{article}
+%\documentclass[aps,prb,preprint]{revtex4}
+\documentclass[11pt]{article}
 \usepackage{endfloat}
 \usepackage{amsmath}
 \usepackage{amssymb}
 \usepackage{epsf}
 \usepackage{times}
-\usepackage{mathptm}
+\usepackage{mathptmx}
 \usepackage{setspace}
 \usepackage{tabularx}
 \usepackage{graphicx}
 \usepackage{booktabs}
 \usepackage{bibentry}
 \usepackage{mathrsfs}
-%\usepackage{berkeley}
 \usepackage[ref]{overcite}
 \pagestyle{plain}
 \pagenumbering{arabic}
@@ -25,9 +25,10 @@
 
 \begin{document}
 
-\title{Is the Ewald Summation necessary? : Pairwise alternatives to the accepted standard for long-range electrostatics}
+\title{Is the Ewald Summation necessary? Pairwise alternatives to the accepted standard for long-range electrostatics}
 
-\author{Christopher J. Fennell and J. Daniel Gezelter \\
+\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail:
+gezelter@nd.edu} \\
 Department of Chemistry and Biochemistry\\ 
 University of Notre Dame\\
 Notre Dame, Indiana 46556}
@@ -35,12 +36,25 @@ Notre Dame, Indiana 46556}
 \date{\today}
 
 \maketitle
-%\doublespacing
+\doublespacing
+
 \nobibliography{}
 \begin{abstract}
-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.  
+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.
 \end{abstract}
 
+\newpage
+
 %\narrowtext 
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -49,46 +63,136 @@ A new method for accumulating electrostatic interactio
 
 \section{Introduction}
 
-In molecular simulations, proper accumulation of the electrostatic interactions is considered one of the most essential and computationally demanding tasks.  
+In molecular simulations, proper accumulation of the electrostatic
+interactions is considered one of the most essential and
+computationally demanding tasks.
 
 \subsection{The Ewald Sum}
 blah blah blah Ewald Sum Important blah blah blah
 
 \begin{figure}
 \centering
-\includegraphics[width = 3.25in]{./ewaldProgression.pdf}
-\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.}
+\includegraphics[width = \linewidth]{./ewaldProgression.pdf}
+\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.}
 \label{fig:ewaldTime}
 \end{figure}
 
 \subsection{The Wolf and Zahn Methods}
-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
+In a recent paper by Wolf \textit{et al.}, a procedure was outlined
+for an accurate accumulation of electrostatic interactions in an
+efficient pairwise fashion.\cite{Wolf99} Wolf \textit{et al.} observed
+that the electrostatic interaction is effectively short-ranged in
+condensed phase systems and that neutralization of the 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
 \begin{equation}
 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\}.
 \label{eq:WolfPot}
 \end{equation}
-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, 
+Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted
+potential.  However, neutralizing the charge contained within each
+cutoff sphere requires the placement of a self-image charge on the
+surface of the cutoff sphere.  This additional self-term in the total
+potential enables Wolf {\it et al.}  to obtain excellent estimates of
+Madelung energies for many crystals.
+
+In order to use their charge-neutralized potential in molecular
+dynamics simulations, Wolf \textit{et al.} suggested taking the
+derivative of this potential prior to evaluation of the limit.  This
+procedure gives an expression for the forces,
 \begin{equation}
-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\}.
+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\},
 \label{eq:WolfForces}
 \end{equation}
-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
+that incorporates both image charges and damping of the electrostatic
+interaction. 
+
+More recently, Zahn \textit{et al.} investigated these potential and
+force expressions for use in simulations involving water.\cite{Zahn02}
+In their work, they pointed out that the method that the forces and
+derivative of the potential are not commensurate.  Attempts to use
+both Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will
+lead to poor energy conservation.  They correctly observed that taking
+the limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating
+the derivatives is mathematically invalid.
+
+Zahn \textit{et al.} proposed a modified form 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 proposed a new damped Coulomb
+potential,
 \begin{equation}
 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\}.
 \label{eq:ZahnPot}
 \end{equation}
-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.
+They showed that this potential does fairly well at capturing the
+structural and dynamic properties of water compared the same
+properties obtained using the Ewald sum.
 
 \subsection{Simple Forms for Pairwise Electrostatics}
-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.
 
-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
+The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et
+al.} are constructed using two different (and separable) computational
+tricks: \begin{itemize}
+\item shifting through the use of image charges, and 
+\item damping the electrostatic interaction.
+\end{itemize}  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 damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the
+post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using
+both techniques.  It is possible, however, to separate these
+tricks and study their effects independently. 
+
+Starting with the original observation that the effective range of the
+electrostatic interaction in condensed phases is considerably less
+than $r^{-1}$, either the cutoff sphere neutralization or the
+distance-dependent damping technique could be used as a foundation for
+a new pairwise summation method.  Wolf \textit{et al.} made the
+observation that charge neutralization within the cutoff sphere plays
+a significant role in energy convergence; therefore we will begin our
+analysis with the various shifted forms that maintain this charge
+neutralization.  We can evaluate the methods of Wolf
+\textit{et al.}  and Zahn \textit{et al.} by considering the standard
+shifted potential,
 \begin{equation}
-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} 
+v^\textrm{SP}(r) = 	\begin{cases}
+v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r >
+R_\textrm{c}  
+\end{cases},
+\label{eq:shiftingPotForm}
+\end{equation}
+and shifted force,
+\begin{equation}
+v^\textrm{SF}(r) = 	\begin{cases}
+v(r)-v_\textrm{c}-\left(\frac{\textrm{d}v(r)}{\textrm{d}r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c})
+&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} 
 						\end{cases},
 \label{eq:shiftingForm}
 \end{equation}
-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}
+functions where $v(r)$ is the unshifted form of the potential, and
+$v_c$ is $v(R_\textrm{c})$.  The Shifted Force ({\sc sf}) form ensures
+that both the potential and the forces goes to zero at the cutoff
+radius, while the Shifted Potential ({\sc sp}) form only ensures the
+potential is smooth at the cutoff radius
+($R_\textrm{c}$).\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}
 \begin{equation}
 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} 
 						\end{cases}.
@@ -162,7 +266,7 @@ Evaluation of the pairwise summation techniques (outli
 
 \begin{figure}
 \centering
-\includegraphics[width=3.25in]{./linearFit.pdf}
+\includegraphics[width = \linewidth]{./linearFit.pdf}
 \caption{Example least squares regression of the $\Delta E$ between configurations for the SF method against SPME in the pure water system.  }
 \label{fig:linearFit}
 \end{figure}
@@ -180,7 +284,7 @@ Each of these $\theta$ values was accumulated in a dis
 
 \begin{figure}
 \centering
-\includegraphics[width=3.25in]{./gaussFit.pdf}
+\includegraphics[width = \linewidth]{./gaussFit.pdf}
 \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.}
 \label{fig:gaussian}
 \end{figure}
@@ -217,7 +321,7 @@ Generation of the system configurations was dependent 
 
 \begin{figure}
 \centering
-\includegraphics[width=3.25in]{./slice.pdf}
+\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}
@@ -234,7 +338,7 @@ In order to evaluate the performance of the pairwise e
 
 \begin{figure}
 \centering
-\includegraphics[width=3.25in]{./delEplot.pdf}
+\includegraphics[width=5.5in]{./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) 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).}
 \label{fig:delE}
 \end{figure}
@@ -249,7 +353,7 @@ Evaluation of pairwise methods for use in Molecular Dy
 
 \begin{figure}
 \centering
-\includegraphics[width=3.25in]{./frcMagplot.pdf}
+\includegraphics[width=5.5in]{./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) 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).}
 \label{fig:frcMag}
 \end{figure}
@@ -258,7 +362,7 @@ Figure \ref{fig:frcMag}, for the most part, parallels 
 
 \begin{figure}
 \centering
-\includegraphics[width=3.25in]{./trqMagplot.pdf}
+\includegraphics[width=5.5in]{./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) 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).}
 \label{fig:trqMag}
 \end{figure}
@@ -271,7 +375,7 @@ Having force and torque vectors with magnitudes that a
 
 \begin{figure}
 \centering
-\includegraphics[width=3.25in]{./frcTrqAngplot.pdf}
+\includegraphics[width=5.5in]{./frcTrqAngplot.pdf}
 \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).}
 \label{fig:frcTrqAng}
 \end{figure}
@@ -321,8 +425,8 @@ In the previous studies using a Shifted-Force variant 
 
 \begin{figure}
 \centering
-\includegraphics[width = 3.25in]{./spectraSquare.pdf}
-\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 differentiate.}
+\includegraphics[width = \linewidth]{./spectraSquare.pdf}
+\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.}
 \label{fig:methodPS}
 \end{figure}
 
@@ -333,8 +437,8 @@ where $S(r)$ is a switching function that smoothly zer
 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.
 \begin{figure}
 \centering
-\includegraphics[width = 3.25in]{./comboSquare.pdf}
-\caption{Normal modes for an NaCl crystal 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.}
+\includegraphics[width = \linewidth]{./comboSquare.pdf}
+\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.}
 \label{fig:dampInc}
 \end{figure}
 
@@ -347,10 +451,9 @@ We are not suggesting any flaw with the Ewald sum; in 
 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
 
 \section{Acknowledgments}
-
 \newpage 
 
-\bibliographystyle{achemso}
+\bibliographystyle{jcp2}
 \bibliography{electrostaticMethods}
 
 
