trunk/electrostaticMethodsPaper/electrostaticMethods.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
%\documentclass[aps,prb,preprint]{revtex4}
\documentclass[11pt]{article}
\usepackage{endfloat}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{epsf}
\usepackage{times}
\usepackage{mathptmx}
\usepackage{setspace}
\usepackage{tabularx}
\usepackage{graphicx}
\usepackage{booktabs}
\usepackage{bibentry}
\usepackage{mathrsfs}
\usepackage[ref]{overcite}
\pagestyle{plain}
\pagenumbering{arabic}
\oddsidemargin 0.0cm \evensidemargin 0.0cm
\topmargin -21pt \headsep 10pt
\textheight 9.0in \textwidth 6.5in
\brokenpenalty=10000
\renewcommand{\baselinestretch}{1.2}
\renewcommand\citemid{\ } % no comma in optional reference note

\begin{document}

\title{Is the Ewald Summation 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} \\
Department of Chemistry and Biochemistry\\ 
University of Notre Dame\\
Notre Dame, Indiana 46556}

\date{\today}

\maketitle
\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.

\end{abstract}

\newpage

%\narrowtext 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                              BODY OF TEXT
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

In molecular simulations, proper accumulation of the electrostatic
interactions is considered one of the most essential and
computationally demanding tasks.  The common molecular mechanics force
fields are founded on representation of the atomic sites centered on
full or partial charges shielded by Lennard-Jones type interactions.
This means that nearly every pair interaction involves an
charge-charge calculation.  Coupled with $r^{-1}$ decay, the monopole
interactions quickly become a burden for molecular systems of all
sizes.  For example, in small systems, the electrostatic pair
interaction may not have decayed appreciably within the box length
leading to an effect excluded from the pair interactions within a unit
box.  In large systems, excessively large cutoffs need to be used to
accurately incorporate their effect, and since the computational cost
increases proportionally with the cutoff sphere, it quickly becomes an
impractical task to perform these calculations.

\subsection{The Ewald Sum}
blah blah blah Ewald Sum Important blah blah blah

\begin{figure}
\centering
\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 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}
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_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}
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 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: \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{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}
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{SP}(r_{ij}) = q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}.                          \label{eq:WolfSP}
\end{equation}
The forces associated with this potential are obtained by taking the derivative, resulting in the following,
\begin{equation}
F^\textrm{SP}(r_{ij}) = q_iq_j\left(-\frac{1}{r_{ij}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}.
\label{eq:FWolfSP}
\end{equation}
These forces are identical to the forces of the standard electrostatic interaction, and this was addressed by Wolf \textit{et al.} as undesirable.  They pointed out that the effect of the image charges is neglected in the forces when they would expect there to be some pressure exerted due to their presence.\cite{Wolf99}  As a consequence the forces, though mathematically valid, may not be physically correct due to this missing component.  Additionally, there is a discontinuity in the forces.  This can be remedied with the use of a switching function to zero the potential and forces smoothly as particles near $R_\textrm{c}$.  

If the derivative term in equation \ref{eq:shiftingForm} is evaluated, we obtain an hitherto undiscussed shifted force Coulomb potential,
\begin{equation}
V^\textrm{SF}(r_{ij}) = q_iq_j\left[\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r_{ij}-R_\textrm{c})\right] \quad r_{ij}\leqslant R_\textrm{c}.
\label{eq:SFPot}
\end{equation}
Taking the derivative of this shifted force potential gives the
following forces,
\begin{equation}
F^\textrm{SF}(r_{ij} =  q_iq_j\left(-\frac{1}{r_{ij}^2}+\frac{1}{R_\textrm{c}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}.
\label{eq:SFForces}
\end{equation}
Using this formulation rather than the simple shifted potential
(Eq. \ref{eq:WolfSP}) means that there are no discontinuities in the
forces in addition to the potential.  This form also has the benefit
that the image charges have a force presence, addressing the concerns
about a missing physical component.  One side effect of this treatment
is a slight alteration in the shape of the potential that comes about
from the derivative term.  Thus, a degree of clarity about the
original formulation of the potential is lost in order to gain
functionality in dynamics simulations.

Wolf \textit{et al.} originally discussed the energetics of the
shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that
it was still insufficient for accurate determination of the energy.
The energy would fluctuate around the expected value with increasing
cutoff radius, but the oscillations appeared to be converging toward
the correct value.\cite{Wolf99} A damping function was incorporated to
accelerate convergence; and though alternative functional forms could
be used,\cite{Jones56,Heyes81} the complimentary error function was
chosen to draw parallels to the Ewald summation.  Incorporating
damping into the simple Coulomb potential,
\begin{equation}
v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}},
\label{eq:dampCoulomb}
\end{equation}
the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived
\textit{via} equation \ref{eq:shiftingForm},
\begin{equation}
V^{\textrm{DSP}}(r_{ij}) = q_iq_j\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}.
\label{eq:DSPPot}
\end{equation}
The derivative of this Shifted-Potential can be taken to obtain forces
for use in MD,
\begin{equation}
F^{\textrm{DSP}}(r_{ij}) = q_iq_j\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) \quad r_{ij}\leqslant R_\textrm{c}.
\label{eq:DSPForces}
\end{equation}
Again, this Shifted-Potential suffers from a discontinuity in the
forces, and a lack of an image-charge component in the forces.  To
remedy these concerns, a Shifted-Force variant is obtained by
inclusion of the derivative term in equation \ref{eq:shiftingForm} to
give,
\begin{equation}
\begin{split}
V^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ &\left.+\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] \quad r_{ij}\leqslant R_\textrm{c}.
\label{eq:DSFPot}
\end{split}
\end{equation}
The derivative of the above potential gives the following forces,
\begin{equation}
\begin{split}
F^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&-\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.+\left(\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right] \quad r_{ij}\leqslant R_\textrm{c}.
\label{eq:DSFForces}
\end{split}
\end{equation}

This new Shifted-Force potential is similar to equation
\ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are
two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term
from equation \ref{eq:shiftingForm} is equal to equation
\ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$.  This
term is not present in the Zahn potential, resulting in a
discontinuity as particles cross $R_\textrm{c}$.  Second, the sign of
the derivative portion is different.  The constant $v_\textrm{c}$ term
is not going to have a presence in the forces after performing the
derivative, but the negative sign does effect the derivative.  In
fact, it introduces a discontinuity in the forces at the cutoff,
because the force function is shifted in the wrong direction and
doesn't cross zero at $R_\textrm{c}$.  Thus, these alterations make
for an electrostatic summation method that is continuous in both the
potential and forces and incorporates the pairwise sum considerations
stressed by Wolf \textit{et al.}\cite{Wolf99}

\section{Methods}

\subsection{What Qualities are Important?}\label{sec:Qualities}
In classical molecular mechanics simulations, there are two primary
techniques utilized to obtain information about the system of
interest: Monte Carlo (MC) and Molecular Dynamics (MD).  Both of these
techniques utilize pairwise summations of interactions between
particle sites, but they use these summations in different ways.

In MC, the potential energy difference between two subsequent
configurations dictates the progression of MC sampling.  Going back to
the origins of this method, the Canonical ensemble acceptance criteria
laid out by Metropolis \textit{et al.} states that a subsequent
configuration is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta
E/kT)$, where $\xi$ is a random number between 0 and
1.\cite{Metropolis53} Maintaining a consistent $\Delta E$ when using
an alternate method for handling the long-range electrostatics ensures
proper sampling within the ensemble.

In MD, the derivative of the potential directs how the system will
progress in time.  Consequently, the force and torque vectors on each
body in the system dictate how it develops as a whole.  If the
magnitude and direction of these vectors are similar when using
alternate electrostatic summation techniques, the dynamics in the near
term will be indistinguishable.  Because error in MD calculations is
cumulative, one should expect greater deviation in the long term
trajectories with greater differences in these vectors between
configurations using different long-range electrostatics.

\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
Evaluation of the pairwise summation techniques (outlined in section
\ref{sec:ESMethods}) for use in MC simulations was performed through
study of the energy differences between conformations.  Considering
the SPME results to be the correct or desired behavior, ideal
performance of a tested method was taken to be agreement between the
energy differences calculated.  Linear least squares regression of the
$\Delta E$ values between configurations using SPME against $\Delta E$
values using tested methods provides a quantitative comparison of this
agreement.  Unitary results for both the correlation and correlation
coefficient for these regressions indicate equivalent energetic
results between the methods.  The correlation is the slope of the
plotted data while the correlation coefficient ($R^2$) is a measure of
the of the data scatter around the fitted line and tells about the
quality of the fit (Fig. \ref{fig:linearFit}).

\begin{figure}
\centering
\includegraphics[width = \linewidth]{./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 shows a data set with a good correlation coefficient.}
\label{fig:linearFit}
\end{figure}

Each system type (detailed in section \ref{sec:RepSims}) studied
consisted of 500 independent configurations, each equilibrated from
higher temperature trajectories. Thus, 124,750 $\Delta E$ data points
are used in a regression of a single system type.  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}.

\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods}
Evaluation of the pairwise methods (outlined in section
\ref{sec:ESMethods}) for use in MD simulations was performed through
comparison of the force and torque vectors obtained with those from
SPME.  Both the magnitude and the direction of these vectors on each
of the bodies in the system were analyzed.  For the magnitude of these
vectors, linear least squares regression analysis can be performed as
described previously for comparing $\Delta E$ values. Instead of a
single value between two system configurations, there is a value for
each particle in each configuration.  For a system of 1000 water
molecules and 40 ions, there are 1040 force vectors and 1000 torque
vectors.  With 500 configurations, this results in 520,000 force and
500,000 torque vector comparisons samples for each system type.

The force and torque vector directions were investigated through
measurement of the angle ($\theta$) formed between those from the
particular method and those from SPME
\begin{equation}
\theta_F = \frac{\vec{F}_\textrm{SPME}}{|\vec{F}_\textrm{SPME}|}\cdot\frac{\vec{F}_\textrm{Method}}{|\vec{F}_\textrm{Method}|}.
\end{equation}
Each of these $\theta$ values was accumulated in a distribution
function, weighted by the area on the unit sphere.  Non-linear fits
were used to measure the shape of the resulting distributions.

\begin{figure}
\centering
\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}

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
particular reason for the profile to adhere to a specific shape.
Because of this and the Gaussian profile's more statistically
meaningful properties, 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 SPME.

\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}
Evaluation of the long-time dynamics of charged systems was performed
by considering the NaCl crystal system while using a subset of the
best performing pairwise methods.  The NaCl crystal was chosen to
avoid possible complications involving the propagation techniques of
orientational motion in molecular systems.  To enhance the atomic
motion, these crystals were equilibrated at 1000 K, near the
experimental $T_m$ for NaCl.  Simulations were performed under the
microcanonical ensemble, and velocity autocorrelation functions
(Eq. \ref{eq:vCorr}) were computed for each of the trajectories,
\begin{equation}
C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle.
\label{eq:vCorr}
\end{equation}
Velocity autocorrelation functions require detailed short time data
and long trajectories for good statistics, thus velocity information
was saved every 5 fs over 100 ps trajectories.  The power spectrum
($I(\omega)$) is obtained via Fourier transform of the autocorrelation
function
\begin{equation}
I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, 
\label{eq:powerSpec}
\end{equation}
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$.

\subsection{Representative Simulations}\label{sec:RepSims}
A variety of common and representative simulations were analyzed to
determine the relative effectiveness of the pairwise summation
techniques in reproducing the energetics and dynamics exhibited by
SPME.  The studied systems were as follows:
\begin{enumerate}
\item Liquid Water
\item Crystalline Water (Ice I$_\textrm{c}$)
\item NaCl Crystal
\item NaCl Melt
\item Low Ionic Strength Solution of NaCl in Water 
\item High Ionic Strength Solution of NaCl in Water 
\item 6 \AA\  Radius Sphere of Argon in Water 
\end{enumerate}
By utilizing the pairwise techniques (outlined in section
\ref{sec:ESMethods}) in systems composed entirely of neutral groups,
charged particles, and mixtures of the two, we can comment on possible
system dependence and/or universal applicability of the techniques.

Generation of the system configurations was dependent on the system
type.  For the solid and liquid water configurations, configuration
snapshots were taken at regular intervals from higher temperature 1000
SPC/E water molecule trajectories and each equilibrated individually.
The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl-
ions and were selected and equilibrated in the same fashion as the
water systems.  For the low and high ionic strength NaCl solutions, 4
and 40 ions were first solvated in a 1000 water molecule boxes
respectively.  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}).

\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{Electrostatic Summation Methods}\label{sec:ESMethods}
Electrostatic summation method comparisons were performed using SPME,
the Shifted-Potential and Shifted-Force methods - both with damping
parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak,
moderate, and strong damping respectively), reaction field with an
infinite dielectric constant, and an unmodified cutoff.  Group-based
cutoffs with a fifth-order polynomial switching function were
necessary for the reaction field simulations and were utilized in the
SP, SF, and pure cutoff methods for comparison to the standard lack of
group-based cutoffs with a hard truncation.  The SPME calculations
were performed using the TINKER implementation of SPME,\cite{Ponder87}
while all other method calculations were performed using the OOPSE
molecular mechanics package.\cite{Meineke05}

These methods were additionally evaluated with three different cutoff
radii (9, 12, and 15 \AA) to investigate possible cutoff radius
dependence.  It should be noted that the damping parameter chosen in
SPME, or so called ``Ewald Coefficient", has a significant effect on
the energies and forces calculated.  Typical molecular mechanics
packages default this to a value dependent on the cutoff radius and a
tolerance (typically less than $1 \times 10^{-4}$ kcal/mol).  Smaller
tolerances are typically associated with increased accuracy in the
real-space portion of the summation.\cite{Essmann95} The default
TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME
calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and
0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively.

\section{Results and Discussion}

\subsection{Configuration Energy Differences}\label{sec:EnergyResults}
In order to evaluate the performance of the pairwise electrostatic
summation methods for Monte Carlo simulations, the energy differences
between configurations were compared to the values obtained when using
SPME.  The results for the subsequent regression analysis are shown in
figure \ref{fig:delE}.

\begin{figure}
\centering
\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}

In this figure, it is apparent that it is unreasonable to expect
realistic results using an unmodified cutoff.  This is not all that
surprising since this results in large energy fluctuations as atoms
move in and out of the cutoff radius.  These fluctuations can be
alleviated to some degree by using group based cutoffs with a
switching function.\cite{Steinbach94} The Group Switch Cutoff row
doesn't show a significant improvement in this plot because the salt
and salt solution systems contain non-neutral groups, see the
accompanying supporting information for a comparison where all groups
are neutral.

Correcting the resulting charged cutoff sphere is one of the purposes
of the damped Coulomb summation proposed by Wolf \textit{et
al.},\cite{Wolf99} and this correction indeed improves the results as
seen in the Shifted-Potental rows.  While the undamped case of this
method is a significant improvement over the pure cutoff, it still
doesn't correlate that well with SPME.  Inclusion of potential damping
improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows
an excellent correlation and quality of fit with the SPME results,
particularly with a cutoff radius greater than 12 \AA .  Use of a
larger damping parameter is more helpful for the shortest cutoff
shown, but it has a detrimental effect on simulations with larger
cutoffs.  In the Shifted-Force sets, increasing damping results in
progressively poorer correlation.  Overall, the undamped case is the
best performing set, as the correlation and quality of fits are
consistently superior regardless of the cutoff distance.  This result
is beneficial in that the undamped case is less computationally
prohibitive do to the lack of complimentary error function calculation
when performing the electrostatic pair interaction.  The reaction
field results illustrates some of that method's limitations, primarily
that it was developed for use in homogenous systems; although it does
provide results that are an improvement over those from an unmodified
cutoff.

\subsection{Magnitudes of the Force and Torque Vectors}

Evaluation of pairwise methods for use in Molecular Dynamics
simulations requires consideration of effects on the forces and
torques.  Investigation of the force and torque vector magnitudes
provides a measure of the strength of these values relative to SPME.
Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the
force and torque vector magnitude regression results for the
accumulated analysis over all the system types.

\begin{figure}
\centering
\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}

Figure \ref{fig:frcMag}, for the most part, parallels the results seen
in the previous $\Delta E$ section.  The unmodified cutoff results are
poor, but using group based cutoffs and a switching function provides
a improvement much more significant than what was seen with $\Delta
E$.  Looking at the Shifted-Potential sets, the slope and $R^2$
improve with the use of damping to an optimal result of 0.2 \AA
$^{-1}$ for the 12 and 15 \AA\ cutoffs.  Further increases in damping,
while beneficial for simulations with a cutoff radius of 9 \AA\ , is
detrimental to simulations with larger cutoff radii.  The undamped
Shifted-Force method gives forces in line with those obtained using
SPME, and use of a damping function results in minor improvement.  The
reaction field results are surprisingly good, considering the poor
quality of the fits for the $\Delta E$ results.  There is still a
considerable degree of scatter in the data, but it correlates well in
general.  To be fair, we again note that the reaction field
calculations do not encompass NaCl crystal and melt systems, so these
results are partly biased towards conditions in which the method
performs more favorably.

\begin{figure}
\centering
\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}

To evaluate the torque vector magnitudes, the data set from which
values are drawn is limited to rigid molecules in the systems
(i.e. water molecules).  In spite of this smaller sampling pool, the
torque vector magnitude results in figure \ref{fig:trqMag} are still
similar to those seen for the forces; however, they more clearly show
the improved behavior that comes with increasing the cutoff radius.
Moderate damping is beneficial to the Shifted-Potential and helpful
yet possibly unnecessary with the Shifted-Force method, and they also
show that over-damping adversely effects all cutoff radii rather than
showing an improvement for systems with short cutoffs.  The reaction
field method performs well when calculating the torques, better than
the Shifted Force method over this limited data set.

\subsection{Directionality of the Force and Torque Vectors}

Having force and torque vectors with magnitudes that are well
correlated to SPME is good, but if they are not pointing in the proper
direction the results will be incorrect.  These vector directions were
investigated through measurement of the angle formed between them and
those from SPME.  The results (Fig. \ref{fig:frcTrqAng}) are compared
through the variance ($\sigma^2$) of the Gaussian fits of the angle
error distributions of the combined set over all system types.

\begin{figure}
\centering
\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}

Both the force and torque $\sigma^2$ results from the analysis of the
total accumulated system data are tabulated in figure
\ref{fig:frcTrqAng}.  All of the sets, aside from the over-damped case
show the improvement afforded by choosing a longer simulation cutoff.
Increasing the cutoff from 9 to 12 \AA\ typically results in a halving
of the distribution widths, with a similar improvement going from 12
to 15 \AA .  The undamped Shifted-Force, Group Based Cutoff, and
Reaction Field methods all do equivalently well at capturing the
direction of both the force and torque vectors.  Using damping
improves the angular behavior significantly for the Shifted-Potential
and moderately for the Shifted-Force methods.  Increasing the damping
too far is destructive for both methods, particularly to the torque
vectors.  Again it is important to recognize that the force vectors
cover all particles in the systems, while torque vectors are only
available for neutral molecular groups.  Damping appears to have a
more beneficial effect on non-neutral bodies, and this observation is
investigated further in the accompanying supporting information.

\begin{table}[htbp]
   \centering
   \caption{Variance ($\sigma^2$) of the force (top set) and torque (bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods.  Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function.  The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.}  
   \begin{tabular}{@{} ccrrrrrrrr @{}}
      \\
      \toprule
      & & \multicolumn{4}{c}{Shifted Potential} & \multicolumn{4}{c}{Shifted Force} \\
      \cmidrule(lr){3-6} 
      \cmidrule(l){7-10}
            $R_\textrm{c}$ & Groups & $\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$ & $\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$\\
      \midrule
     
9 \AA   & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\
        & \textbf{Y} & \textbf{2.486} & \textbf{2.160} & \textbf{0.667} & \textbf{0.608} & \textbf{1.768} & \textbf{1.766} & \textbf{0.676} & \textbf{0.609} \\
12 \AA  & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\
        & \textbf{Y} & \textbf{0.515} & \textbf{0.288} & \textbf{0.127} & \textbf{0.586} & \textbf{0.308} & \textbf{0.249} & \textbf{0.127} & \textbf{0.586} \\
15 \AA  & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\
        & \textbf{Y} & \textbf{0.228} & \textbf{0.099} & \textbf{0.121} & \textbf{0.598} & \textbf{0.144} & \textbf{0.090} & \textbf{0.121} & \textbf{0.598} \\      

      \midrule
      
9 \AA   & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\
        & \textbf{Y} & \textbf{2.115} & \textbf{1.914} & \textbf{1.878} & \textbf{5.142} & \textbf{2.076} & \textbf{2.039} & \textbf{1.972} & \textbf{5.146} \\
12 \AA  & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\
        & \textbf{Y} & \textbf{0.810} & \textbf{0.685} & \textbf{1.352} & \textbf{5.082} & \textbf{0.765} & \textbf{0.714} & \textbf{1.360} & \textbf{5.082} \\
15 \AA  & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\
        & \textbf{Y} & \textbf{0.282} & \textbf{0.294} & \textbf{1.272} & \textbf{4.999} & \textbf{0.324} & \textbf{0.318} & \textbf{1.272} & \textbf{4.999} \\

      \bottomrule
   \end{tabular}
   \label{tab:groupAngle}
\end{table}

Although not discussed previously, group based cutoffs can be applied
to both the Shifted-Potential and Shifted-Force methods.  Use off a
switching function corrects for the discontinuities that arise when
atoms of a group exit the cutoff before the group's center of mass.
Though there are no significant benefit or drawbacks observed in
$\Delta E$ and vector magnitude results when doing this, there is a
measurable improvement in the vector angle results.  Table
\ref{tab:groupAngle} shows the angular variance values obtained using
group based cutoffs and a switching function alongside the standard
results seen in figure \ref{fig:frcTrqAng} for comparison purposes.
The Shifted-Potential shows much narrower angular distributions for
both the force and torque vectors when using an $\alpha$ of 0.2
\AA$^{-1}$ or less, while Shifted-Force shows improvements in the
undamped and lightly damped cases.  Thus, by calculating the
electrostatic interactions in terms of molecular pairs rather than
atomic pairs, the direction of the force and torque vectors are
determined more accurately.

One additional trend to recognize in table \ref{tab:groupAngle} is
that the $\sigma^2$ values for both Shifted-Potential and
Shifted-Force converge as $\alpha$ increases, something that is easier
to see when using group based cutoffs.  Looking back on figures
\ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this
behavior clearly at large $\alpha$ and cutoff values.  The reason for
this is that the complimentary error function inserted into the
potential weakens the electrostatic interaction as $\alpha$ increases.
Thus, at larger values of $\alpha$, both the summation method types
progress toward non-interacting functions, so care is required in
choosing large damping functions lest one generate an undesirable loss
in the pair interaction.  Kast \textit{et al.}  developed a method for
choosing appropriate $\alpha$ values for these types of electrostatic
summation methods by fitting to $g(r)$ data, and their methods
indicate optimal values of 0.34, 0.25, and 0.16 \AA$^{-1}$ for cutoff
values of 9, 12, and 15 \AA\ respectively.\cite{Kast03} These appear
to be reasonable choices to obtain proper MC behavior
(Fig. \ref{fig:delE}); however, based on these findings, choices this
high would introduce error in the molecular torques, particularly for
the shorter cutoffs.  Based on the above findings, empirical damping
up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably
unnecessary when using the Shifted-Force method.

\subsection{Collective Motion: Power Spectra of NaCl Crystals}

In the previous studies using a Shifted-Force variant of the damped
Wolf coulomb potential, the structure and dynamics of water were
investigated rather extensively.\cite{Zahn02,Kast03} Their results
indicated that the damped Shifted-Force method results in properties
very similar to those obtained when using the Ewald summation.
Considering the statistical results shown above, the good performance
of this method is not that surprising.  Rather than consider the same
systems and simply recapitulate their results, we decided to look at
the solid state dynamical behavior obtained using the best performing
summation methods from the above results.

\begin{figure}
\centering
\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}

Figure \ref{fig:methodPS} shows the power spectra for the NaCl
crystals (from averaged Na and Cl ion velocity autocorrelation
functions) using the stated electrostatic summation methods.  While
high frequency peaks of all the spectra overlap, showing the same
general features, the low frequency region shows how the summation
methods differ.  Considering the low-frequency inset (expanded in the
upper frame of figure \ref{fig:dampInc}), at frequencies below 100
cm$^{-1}$, the correlated motions are blue-shifted when using undamped
or weakly damped Shifted-Force.  When using moderate damping ($\alpha
= 0.2$ \AA$^{-1}$) both the Shifted-Force and Shifted-Potential
methods give near identical correlated motion behavior as the Ewald
method (which has a damping value of 0.3119).  The damping acts as a
distance dependent Gaussian screening of the point charges for the
pairwise summation methods.  This weakening of the electrostatic
interaction with distance explains why the long-ranged correlated
motions are at lower frequencies for the moderately damped methods
than for undamped or weakly damped methods.  To see this effect more
clearly, we show how damping strength affects a simple real-space
electrostatic potential,
\begin{equation}
V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r_{ij}})}{r_{ij}}\right]S(r),
\end{equation}
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 = \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}

\section{Conclusions}

This investigation of pairwise electrostatic summation techniques
shows that there are viable and more computationally efficient
electrostatic summation techniques than the Ewald summation, chiefly
methods derived from the damped Coulombic sum originally proposed by
Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the
Shifted-Force method, reformulated above as equation \ref{eq:SFPot},
shows a remarkable ability to reproduce the energetic and dynamic
characteristics exhibited by simulations employing lattice summation
techniques.  The cumulative energy difference results showed the
undamped Shifted-Force and moderately damped Shifted-Potential methods
produced results nearly identical to SPME.  Similarly for the dynamic
features, the undamped or moderately damped Shifted-Force and
moderately damped Shifted-Potential methods produce force and torque
vector magnitude and directions very similar to the expected values.
These results translate into long-time dynamic behavior equivalent to
that produced in simulations using SPME.

Aside from the computational cost benefit, these techniques have
applicability in situations where the use of the Ewald sum can prove
problematic.  Primary among them is their use in interfacial systems,
where the unmodified lattice sum techniques artificially accentuate
the periodicity of the system in an undesirable manner.  There have
been alterations to the standard Ewald techniques, via corrections and
reformulations, to compensate for these systems; but the pairwise
techniques discussed here require no modifications, making them
natural tools to tackle these problems.  Additionally, this
transferability gives them benefits over other pairwise methods, like
reaction field, because estimations of physical properties (e.g. the
dielectric constant) are unnecessary.

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{jcp2}
\bibliography{electrostaticMethods}


\end{document}
Revision:	2620
Committed:	Wed Mar 15 04:34:51 2006 UTC (19 years, 1 month ago) by chrisfen
Content type:	application/x-tex
File size:	47092 byte(s)
Log Message:	changing equations and adding text