trunk/electrostaticMethodsPaper/electrostaticMethods.tex

%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
%\documentclass[aps,prb,preprint]{revtex4}
\documentclass[11pt]{article}
\usepackage{endfloat}
\usepackage{amsmath,bm}
\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 essential and is one of the most
computationally-demanding tasks.  The common molecular mechanics force
fields represent atomic sites with full or partial charges protected
by Lennard-Jones (short range) interactions.  This means that nearly
every pair interaction involves a calculation of charge-charge forces.
Coupled with relatively long-ranged $r^{-1}$ decay, the monopole
interactions quickly become the most expensive part of molecular
simulations.  Historically, the electrostatic pair interaction would
not have decayed appreciably within the typical box lengths that could
be feasibly simulated.  In the larger systems that are more typical of
modern simulations, large cutoffs should be used to incorporate
electrostatics correctly.

There have been many efforts to address the proper and practical
handling of electrostatic interactions, and these have resulted in a
variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are
typically classified as implicit methods (i.e., continuum dielectrics,
static dipolar fields),\cite{Born20,Grossfield00} explicit methods
(i.e., Ewald summations, interaction shifting or
truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e.,
reaction field type methods, fast multipole
methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are
often preferred because they physically incorporate solvent molecules
in the system of interest, but these methods are sometimes difficult
to utilize because of their high computational cost.\cite{Roux99} In
addition to the computational cost, there have been some questions
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
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
sum,\cite{Onsager36,Essmann99} which is our reference method for
handling long-range electrostatic interactions. The new methods for
handling electrostatics have the potential to scale linearly with
increasing system size since they involve only a simple modification
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.

\subsection{The Ewald Sum}
The complete accumulation 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,
\begin{equation}
V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime \left[ \sum_{i=1}^N\sum_{j=1}^N \phi\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right) \right],
\label{eq:PBCSum}
\end{equation}
where the sum over $\mathbf{n}$ is a sum over all periodic box
replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the
prime indicates $i = j$ are neglected for $\mathbf{n} =
0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic
particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is
the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and
$j$, and $\phi$ is the solution to Poisson's equation
($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for
charge-charge interactions). In the case of monopole electrostatics,
eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for
non-neutral systems.

The electrostatic summation problem was originally studied by Ewald
for the case of an infinite crystal.\cite{Ewald21}. The approach he
took was to convert this conditionally convergent sum into two
absolutely convergent summations: a short-ranged real-space summation
and a long-ranged reciprocal-space summation,
\begin{equation}
\begin{split}
V_\textrm{elec} = \frac{1}{2}& \sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime \frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)}{|\mathbf{r}_{ij}+\mathbf{n}|} \\ &+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} \exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) \cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ &- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 + \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3}\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2,
\end{split}
\label{eq:EwaldSum}
\end{equation}
where $\alpha$ is the damping or convergence parameter with units of
\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric
constant of the surrounding medium. The final two terms of
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term
for interacting with a surrounding dielectric.\cite{Allen87} This
dipolar term was neglected in early applications in molecular
simulations,\cite{Brush66,Woodcock71} until it was introduced by de
Leeuw {\it et al.} to address situations where the unit cell has a
dipole moment which is magnified through replication of the periodic
images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the
system is said to be using conducting (or ``tin-foil'') boundary
conditions, $\epsilon_{\rm S} = \infty$. Figure
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over
time.  Initially, due to the small sizes of the systems that could be
feasibly simulated, the entire simulation box was replicated to
convergence.  In more modern simulations, the simulation boxes have
grown large enough that a real-space cutoff could potentially give
convergent behavior.  Indeed, it has often been observed that the
reciprocal-space portion of the Ewald sum can be small and rapidly
convergent compared to the real-space portion with the choice of small
$\alpha$.\cite{Karasawa89,Kolafa92}

\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}

The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm.  The
convergence parameter $(\alpha)$ plays an important role in balancing
the computational cost between the direct and reciprocal-space
portions of the summation.  The choice of this value allows one to
select whether the real-space or reciprocal space portion of the
summation is an $\mathscr{O}(N^2)$ calculation (with the other being
$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of
$\alpha$ and thoughtful algorithm development, this cost can be
reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route
taken to reduce the cost of the Ewald summation even further is to set
$\alpha$ such that the real-space interactions decay rapidly, allowing
for a short spherical cutoff. Then the reciprocal space summation is
optimized.  These optimizations usually involve utilization of the
fast Fourier transform (FFT),\cite{Hockney81} leading to the
particle-particle particle-mesh (P3M) and particle mesh Ewald (PME)
methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these
methods, the cost of the reciprocal-space portion of the Ewald
summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N
\log N)$.

These developments and optimizations have made the use of the Ewald
summation routine in simulations with periodic boundary
conditions. However, in certain systems, such as vapor-liquid
interfaces and membranes, the intrinsic three-dimensional periodicity
can prove problematic.  The Ewald sum has been reformulated to handle
2D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the
new methods are computationally expensive.\cite{Spohr97,Yeh99}
Inclusion of a correction term in the Ewald summation is a possible
direction for handling 2D systems while still enabling the use of the
modern optimizations.\cite{Yeh99}

Several studies have recognized that the inherent periodicity in the
Ewald sum can also have an effect on three-dimensional
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
Solvated proteins are essentially kept at high concentration due to
the periodicity of the electrostatic summation method.  In these
systems, the more compact folded states of a protein can be
artificially stabilized by the periodic replicas introduced by the
Ewald summation.\cite{Weber00} Thus, care must be taken when
considering the use of the Ewald summation where the assumed
periodicity would introduce spurious effects in the system dynamics.

\subsection{The Wolf and Zahn Methods}
In a recent paper by Wolf \textit{et al.}, a procedure was outlined
for the accurate accumulation of electrostatic interactions in an
efficient pairwise fashion.  This procedure lacks the inherent
periodicity of the Ewald summation.\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_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}
\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 enabled 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}\left(\alpha r_{ij}\right)}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{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 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 gives forces for a different potential energy function
than the one shown in eq. (\ref{eq:WolfPot}).

Zahn \textit{et al.} introduced a modified form of this summation
method as a way to use the technique in Molecular Dynamics
simulations.  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}
and 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{enumerate}
\item shifting through the use of image charges, and 
\item damping the electrostatic interaction.
\end{enumerate}  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{d v(r)}{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}

The forces associated with the shifted potential are simply the forces
of the unshifted potential itself (when inside the cutoff sphere),
\begin{equation}
F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right),
\end{equation}
and are zero outside.  Inside the cutoff sphere, the forces associated
with the shifted force form can be written,
\begin{equation}
F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d
v(r)}{dr} \right)_{r=R_\textrm{c}}.
\end{equation}
 
If the potential, $v(r)$, is taken to be the normal Coulomb potential,
\begin{equation}
v(r) = \frac{q_i q_j}{r},
\label{eq:Coulomb}
\end{equation}
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et
al.}'s undamped prescription:
\begin{equation}
V_\textrm{SP}(r) =
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad
r\leqslant R_\textrm{c},
\label{eq:SPPot}
\end{equation}
with associated forces,
\begin{equation}
F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}.
\label{eq:SPForces}
\end{equation}
These forces are identical to the forces of the standard Coulomb
interaction, and cutting these off at $R_c$ 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 this form is
used,\cite{Wolf99} thereby eliminating any benefit from the method in
molecular dynamics.  Additionally, there is a discontinuity in the
forces at the cutoff radius which results in energy drift during MD
simulations.

The shifted force ({\sc sf}) form using the normal Coulomb potential
will give,
\begin{equation}
V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] \quad r\leqslant R_\textrm{c}.
\label{eq:SFPot}
\end{equation}
with associated forces,
\begin{equation}
F_\textrm{SF}(r) =  q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) \quad r\leqslant R_\textrm{c}.
\label{eq:SFForces}
\end{equation}
This formulation has the benefits that there are no discontinuities at
the cutoff radius, while the neutralizing image charges are present in
both the energy and force expressions.  It would be simple to add the
self-neutralizing term back when computing the total energy of the
system, thereby maintaining the agreement with the Madelung energies.
A side effect of this treatment is the alteration in the shape of the
potential that comes from the derivative term.  Thus, a degree of
clarity about agreement with the empirical 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:SPPot}) and found that it was
insufficient for accurate determination of the energy with reasonable
cutoff distances.  The calculated Madelung energies fluctuated around
the expected value as the cutoff radius was increased, but the
oscillations converged toward the correct value.\cite{Wolf99} A
damping function was incorporated to accelerate the convergence; and
though alternative forms for the damping function could be
used,\cite{Jones56,Heyes81} the complimentary error function was
chosen to mirror the effective screening used in the Ewald summation.
Incorporating this error function damping into the simple Coulomb
potential,
\begin{equation}
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
\label{eq:dampCoulomb}
\end{equation}
the shifted potential (eq. (\ref{eq:SPPot})) becomes
\begin{equation}
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}-\frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c},
\label{eq:DSPPot}
\end{equation}
with associated forces,
\begin{equation}
F_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \quad r\leqslant R_\textrm{c}.
\label{eq:DSPForces}
\end{equation}
Again, this damped shifted potential suffers from a
force-discontinuity at the cutoff radius, and the image charges play
no role in the forces.  To remedy these concerns, one may derive a
{\sc sf} variant by including the derivative term in
eq. (\ref{eq:shiftingForm}),
\begin{equation}
\begin{split}
V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r\right)}{r}-\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-R_\mathrm{c}\right)\ \right] \quad r\leqslant R_\textrm{c}.
\label{eq:DSFPot}
\end{split}
\end{equation}
The derivative of the above potential will lead to the following forces,
\begin{equation}
\begin{split}
F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ &\left.-\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] \quad r\leqslant R_\textrm{c}.
\label{eq:DSFForces}
\end{split}
\end{equation}
If the damping parameter $(\alpha)$ is set to zero, the undamped case,
eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly
recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}).

This new {\sc sf} 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
eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb})
with $r$ replaced by $R_\textrm{c}$.  This term is {\it not} present
in the Zahn potential, resulting in a potential discontinuity as
particles cross $R_\textrm{c}$.  Second, the sign of the derivative
portion is different.  The missing $v_\textrm{c}$ term would not
affect molecular dynamics simulations (although the computed energy
would be expected to have sudden jumps as particle distances crossed
$R_c$).  The sign problem is a potential source of errors, however.
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}$.

Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an
electrostatic summation method in which the potential and forces are
continuous at the cutoff radius and which incorporates the damping
function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of
this paper, we will evaluate exactly how good these methods ({\sc sp},
{\sc sf}, damping) are at reproducing the correct electrostatic
summation performed by the Ewald sum.

\subsection{Other alternatives}
In addition to the methods described above, we considered some other
techniques that are commonly used in molecular simulations.  The
simplest of these is group-based cutoffs.  Though of little use for
charged molecules, collecting atoms into neutral groups takes
advantage of the observation that the electrostatic interactions decay
faster than those for monopolar pairs.\cite{Steinbach94} When
considering these molecules as neutral groups, the relative
orientations of the molecules control the strength of the interactions
at the cutoff radius.  Consequently, as these molecular particles move
through $R_\textrm{c}$, the energy will drift upward due to the
anisotropy of the net molecular dipole interactions.\cite{Rahman71} To
maintain good energy conservation, both the potential and derivative
need to be smoothly switched to zero at $R_\textrm{c}$.\cite{Adams79}
This is accomplished using a standard switching function.  If a smooth
second derivative is desired, a fifth (or higher) order polynomial can
be used.\cite{Andrea83}

Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$,
and to incorporate the effects of the surroundings, a method like
Reaction Field ({\sc rf}) can be used.  The original theory for {\sc
rf} was originally developed by Onsager,\cite{Onsager36} and it was
applied in simulations for the study of water by Barker and
Watts.\cite{Barker73} In modern simulation codes, {\sc rf} is simply
an extension of the group-based cutoff method where the net dipole
within the cutoff sphere polarizes an external dielectric, which
reacts back on the central dipole.  The same switching function
considerations for group-based cutoffs need to made for {\sc rf}, with
the additional pre-specification of a dielectric constant.

\section{Methods}

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 configurations dictates
the progression of MC sampling.  Going back to the origins of this
method, the acceptance criterion for the canonical ensemble 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 the correct $\Delta E$ when using an alternate method for
handling the long-range electrostatics will ensure proper sampling
from the ensemble.

In MD, the derivative of the potential governs how the system will
progress in time.  Consequently, the force and torque vectors on each
body in the system dictate how the system evolves.  If the magnitude
and direction of these vectors are similar when using alternate
electrostatic summation techniques, the dynamics in the short term
will be indistinguishable.  Because error in MD calculations is
cumulative, one should expect greater deviation at longer times,
although methods which have large differences in the force and torque
vectors will diverge from each other more rapidly.

\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}

The pairwise summation techniques (outlined in section
\ref{sec:ESMethods}) were evaluated for use in MC simulations by 
studying the energy differences between conformations.  We took the
SPME-computed energy difference between two conformations to be the
correct behavior. An ideal performance by an alternative method would
reproduce these energy differences exactly (even if the absolute
energies calculated by the methods are different).  Since none of the
methods provide exact energy differences, we used linear least squares
regressions of energy gap data to evaluate how closely the methods
mimicked the Ewald energy gaps.  Unitary results for both the
correlation (slope) and correlation coefficient for these regressions
indicate perfect agreement between the alternative method and SPME.
Sample correlation plots for two alternate methods are shown in
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}) 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.

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}
We evaluated the pairwise methods (outlined in section
\ref{sec:ESMethods}) for use in MD simulations by 
comparing the force and torque vectors with those obtained using the
reference Ewald summation (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 analyses were performed as described previously for
comparing $\Delta E$ values.  Instead of a single energy difference
between two system configurations, we compared the magnitudes of the
forces (and torques) on each molecule 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.
Additionally, data from seven different system types was aggregated
before the comparison was made.

The {\it directionality} of the force and torque vectors was
investigated through measurement of the angle ($\theta$) formed
between those computed from the particular method and those from SPME,
\begin{equation}
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right),
\end{equation}
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force
vector computed using method M.

Each of these $\theta$ values was accumulated in a distribution
function and weighted by the area on the unit sphere.  Non-linear
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 (SPME).

\subsection{Short-time Dynamics}

The effects of the alternative electrostatic summation methods on the
short-time dynamics of charged systems were evaluated by considering a
NaCl crystal at a temperature of 1000 K.  A subset of the best
performing pairwise methods was used in this comparison.  The NaCl
crystal was chosen to avoid possible complications from the treatment
of orientational motion in molecular systems.  All systems were
started with the same initial positions and velocities.  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) = \frac{\langle v_i(0)\cdot v_i(t)\rangle}{\langle v_i(0)\cdot v_i(0)\rangle}.
\label{eq:vCorr}
\end{equation}
Velocity autocorrelation functions require detailed short time data,
thus velocity information was saved every 2 fs over 10 ps
trajectories. Because the NaCl crystal is composed of two different
atom types, the average of the two resulting velocity autocorrelation
functions was used for comparisons.

\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods}

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
($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,
\label{eq:powerSpec}
\end{equation}
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the
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. 

\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 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:
\begin{enumerate}
\item liquid water (SPC/E),\cite{Berendsen87}
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E),
\item NaCl crystals,
\item NaCl melts,
\item a low ionic strength solution of NaCl in water (0.11 M),
\item a high ionic strength solution of NaCl in water (1.1 M), and
\item a 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 hope to discern under
which conditions it will be possible to use one of the alternative
summation methodologies instead of the Ewald sum.

For the solid and liquid water configurations, configurations were
taken at regular intervals from high temperature trajectories of 1000
SPC/E water molecules.  Each configuration was equilibrated
independently at a lower temperature (300~K for the liquid, 200~K for
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}).

These procedures guaranteed us a set of representative configurations
from chemically-relevant systems sampled from an appropriate
ensemble. 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 (SPME):
\begin{itemize}
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
and 0.3 \AA$^{-1}$,
\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
and 0.3 \AA$^{-1}$,
\item reaction field with an infinite dielectric constant, and 
\item an unmodified cutoff.
\end{itemize}
Group-based cutoffs with a fifth-order polynomial switching function
were utilized for the reaction field simulations.  Additionally, we
investigated the use of these cutoffs with the SP, SF, and pure
cutoff.  The SPME electrostatics 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} All other portions of the energy calculation
(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
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.
Typical molecular mechanics packages set 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 at the expense of increased time spent calculating
the reciprocal-space portion of the
summation.\cite{Perram88,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}

The most striking feature of this plot is how well the Shifted Force
({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy
differences.  For the undamped {\sc sf} method, and the
moderately-damped {\sc sp} methods, the results are nearly
indistinguishable from the Ewald results.  The other common methods do
significantly less well.  

The unmodified cutoff method is essentially unusable.  This is not
surprising since hard cutoffs give large energy fluctuations as atoms
or molecules move in and out of the cutoff
radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to
some degree by using group based cutoffs with a switching
function.\cite{Adams79,Steinbach94,Leach01} However, we do not see
significant improvement using the group-switched cutoff because the
salt and salt solution systems contain non-neutral groups.  Interested
readers can consult the accompanying supporting information for a
comparison where all groups are neutral.

For the {\sc sp} method, inclusion of potential damping improves the
agreement with Ewald, 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 {\sc sf} sets, increasing damping results in progressively
worse correlation with Ewald.  Overall, the undamped case is the best
performing set, as the correlation and quality of fits are
consistently superior regardless of the cutoff distance.  The undamped
case is also less computationally demanding (because no evaluation of
the complementary error function is required).

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 {\sc sp} 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
{\sc sf} 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 {\sc sp} and helpful
yet possibly unnecessary with the {\sc sf} 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 {\sc sf}, 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 {\sc sp}
and moderately for the {\sc sf} 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 {\sc sp} and {\sc sf} 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 {\sc sp} shows much narrower angular distributions for
both the force and torque vectors when using an $\alpha$ of 0.2
\AA$^{-1}$ or less, while {\sc sf} 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 {\sc sp} and
{\sc sf} 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 {\sc sf} method.

\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}

In the previous studies using a {\sc sf} 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 {\sc sf} 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]{./vCorrPlot.pdf}
\caption{Velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset is a magnification of the first trough. The times to first collision are nearly identical, but the differences can be seen in the peaks and troughs, where the undamped to weakly damped methods are stiffer than the moderately damped and SPME methods.}
\label{fig:vCorrPlot}
\end{figure}

The short-time decays through the first collision are nearly identical
in figure \ref{fig:vCorrPlot}, but the peaks and troughs of the
functions show how the methods differ.  The undamped {\sc sf} method
has deeper troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher
peaks than any of the other methods.  As the damping function is
increased, these peaks are smoothed out, and approach the SPME
curve. The damping acts as a distance dependent Gaussian screening of
the point charges for the pairwise summation methods; thus, the
collisions are more elastic in the undamped {\sc sf} potential, and the
stiffness of the potential is diminished as the electrostatic
interactions are softened by the damping function.  With $\alpha$
values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are
nearly identical and track the SPME features quite well.  This is not
too surprising in that the differences between the {\sc sf} and {\sc
sp} potentials are mitigated with increased damping.  However, this
appears to indicate that once damping is utilized, the form of the
potential seems to play a lesser role in the crystal dynamics.

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

The short time dynamics were extended to evaluate how the differences
between the methods affect the collective long-time motion.  The same
electrostatic summation methods were used as in the short time
velocity autocorrelation function evaluation, but the trajectories
were sampled over a much longer time. The power spectra of the
resulting velocity autocorrelation functions were calculated and are
displayed in figure \ref{fig:methodPS}.

\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, {\sc sf} ($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\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}

While high frequency peaks of the spectra in this figure 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 {\sc sf}.  When
using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the {\sc sf}
and {\sc sp} methods give near identical correlated motion behavior as
the Ewald method (which has a damping value of 0.3119).  This
weakening of the electrostatic interaction with increased damping
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 alone 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})}{r}\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 {\sc sf} method
producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite
respectable and shows that the shifted force procedure accounts for
most of the effect afforded through use of the Ewald summation.
However, it appears as though moderate damping is required for
accurate reproduction of crystal dynamics.
\begin{figure}
\centering
\includegraphics[width = \linewidth]{./comboSquare.pdf}
\caption{Regions of spectra showing the low-frequency correlated motions for NaCl crystals at 1000 K using various electrostatic summation methods.  The upper plot is a zoomed inset from figure \ref{fig:methodPS}.  As the damping value for the {\sc sf} potential increases, the low-frequency peaks red-shift.  The lower plot is of spectra 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
{\sc sf} method, reformulated above as eq. (\ref{eq:DSFPot}),
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 {\sc sf} and moderately damped {\sc sp} methods
produced results nearly identical to SPME.  Similarly for the dynamic
features, the undamped or moderately damped {\sc sf} and
moderately damped {\sc sp} 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 researchers have come to expect

\section{Acknowledgments}
\newpage 

\bibliographystyle{jcp2}
\bibliography{electrostaticMethods}


\end{document}
Revision:	2649
Committed:	Tue Mar 21 14:58:12 2006 UTC (18 years, 5 months ago) by chrisfen
Content type:	application/x-tex
File size:	60681 byte(s)
Log Message:	recip. space convergence citations