trunk/chrisDissertation/Electrostatics.tex

\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES}

In molecular simulations, proper accumulation of 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 repulsive Lennard-Jones 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 chapter, we focus on a new set of pairwise methods devised by
Wolf {\it et al.},\cite{Wolf99} which we further extend.  These
methods along with a few other mixed methods (i.e. reaction field) are
compared with the smooth particle mesh Ewald
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.

\section{The Ewald Sum}

The complete accumulation of the electrostatic interactions in a system with
periodic boundary conditions (PBC) requires the consideration of the
effect of all charges within a (cubic) simulation box as well as those
in the periodic replicas,
\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,
equation (\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
equation (\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$. 

\begin{figure}
\includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf}
\caption{The change in the need for the Ewald sum with
increasing computational power.  A:~Initially, only small systems
could be studied, and the Ewald sum replicated the simulation box to
convergence.  B:~Now, radial cutoff methods should be able to reach
convergence for the larger systems of charges that are common today.}
\label{fig:ewaldTime}
\end{figure}
Figure \ref{fig:ewaldTime} shows how the Ewald sum has been applied
over time.  Initially, due to the small system sizes that could be
simulated feasibly, the entire simulation box was replicated to
convergence.  In more modern simulations, the systems have grown large
enough that a real-space cutoff could potentially give convergent
behavior.  Indeed, it has been observed that with the choice of a
small $\alpha$, the reciprocal-space portion of the Ewald sum can be
rapidly convergent and small relative to the real-space
portion.\cite{Karasawa89,Kolafa92}

The original Ewald summation is an $\mathcal{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 $\mathcal{O}(N^2)$ calculation (with the other being
$\mathcal{O}(N)$).\cite{Sagui99} With the appropriate choice of
$\alpha$ and thoughtful algorithm development, this cost can be
reduced to $\mathcal{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 $\mathcal{O}(N^2)$ down to $\mathcal{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
2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} but these
methods are computationally expensive.\cite{Spohr97,Yeh99} More
recently, there have been several successful efforts toward reducing
the computational cost of 2-D lattice
summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04}
bringing them more in line with the cost of the full 3-D summation.

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.


\section{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 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}
Equation (\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}
\begin{split}
F_{\textrm{Wolf}}(r_{ij}) = q_i q_j \Biggr\{& 
\Biggr[\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}}
\Biggr]\\ 
&-\Biggr[
\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}}}
\Biggr]\Biggr\},
\end{split}
\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
equations (\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 equation (\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}
\begin{split}
V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\Biggr\{&
\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)\Biggr\},
\end{split}
\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.

\section{Simple Forms for Pairwise Electrostatics}\label{sec:PairwiseDerivation}

The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et
al.} are constructed using two different (and separable) computational
tricks: 

\begin{enumerate}[itemsep=0pt]
\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 {\it et al.} and
Zahn {\it 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
equation (\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(\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)\ \Biggr{]} 
\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(\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)\Biggr{]} 
\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,
equations (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly
recovered from equations (\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 equation
(\ref{eq:shiftingForm}) is equal to equation (\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}$.

Equations (\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.


\section{Evaluating Pairwise Summation Techniques}

As mentioned in the introduction, there are two primary techniques
utilized to obtain information about the system of interest in
classical molecular mechanics simulations: 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}

\begin{figure}
\centering
\includegraphics[width = 3.5in]{./figures/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}
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
{\sc 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 {\sc spme}.
Sample correlation plots for two alternate methods are shown in
Fig. \ref{fig:linearFit}.

Each of the seven system types (detailed in section \ref{sec:RepSims})
were represented using 500 independent configurations.  Thus, each of
the alternative (non-Ewald) electrostatic summation methods was
evaluated using an accumulated 873,250 configurational energy
differences.

Results and discussion for the individual analysis of each of the
system types appear in appendix \ref{app:IndividualResults}, while the
cumulative results over all the investigated systems appear below in
sections \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 ({\sc 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 {\sc 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.  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. Thus, gaussian fits were used to measure the width of the
resulting distributions. The variance ($\sigma^2$) was extracted from
each of these fits and was used to compare distribution widths.
Values of $\sigma^2$ near zero indicate vector directions
indistinguishable from those calculated when using the reference
method ({\sc spme}).

\subsection{Short-time Dynamics}

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(0)\cdot v(t)\rangle}{\langle v^2\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 molecular simulations were analyzed to
determine the relative effectiveness of the pairwise summation
techniques in reproducing the energetics and dynamics exhibited by
{\sc spme}.  We wanted to span the space of typical molecular
simulations (i.e. from liquids of neutral molecules to ionic
crystals), so the systems studied were:

\begin{enumerate}[itemsep=0pt]
\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. In order to introduce measurable fluctuations in the
configuration energy differences, the crystalline simulations were
equilibrated at 1000~K, near the $T_\textrm{m}$ for NaCl. The liquid
NaCl configurations needed to represent a fully disordered array of
point charges, so the high temperature of 7000~K was selected for
equilibration. The ionic solutions were made by solvating 4 (or 40)
ions in a periodic box containing 1000 SPC/E water molecules.  Ion and
water positions were then randomly swapped, and the resulting
configurations were again equilibrated individually.  Finally, for the
Argon / Water ``charge void'' systems, the identities of all the SPC/E
waters within 6~\AA\ of the center of the equilibrated water
configurations were converted to argon.

These procedures guaranteed us a set of representative configurations
from chemically-relevant systems sampled from appropriate
ensembles. Force field parameters for the ions and Argon were taken
from the force field utilized by {\sc oopse}.\cite{Meineke05}

\subsection{Comparison of Summation Methods}\label{sec:ESMethods}
We compared the following alternative summation methods with results
from the reference method ({\sc spme}):

\begin{enumerate}[itemsep=0pt]
\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{enumerate}

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 {\sc sp}, {\sc sf}, and pure
cutoff.  The {\sc spme} electrostatics were performed using the {\sc tinker}
implementation of {\sc spme},\cite{Ponder87} while all other calculations
were performed using the {\sc 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 alternative methods were also evaluated with three different
cutoff radii (9, 12, and 15~\AA).  As noted previously, the
convergence parameter ($\alpha$) plays a role in the balance of the
real-space and reciprocal-space portions of the Ewald calculation.
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
increasing accuracy at the expense of computational time spent on the
reciprocal-space portion of the summation.\cite{Perram88,Essmann95}
The default {\sc tinker} tolerance of $1 \times 10^{-8}$~kcal/mol was used
in all {\sc 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{Configuration Energy Difference Results}\label{sec:EnergyResults}
In order to evaluate the performance of the pairwise electrostatic
summation methods for Monte Carlo (MC) simulations, the energy
differences between configurations were compared to the values
obtained when using {\sc spme}.  The results for the combined
regression analysis of all of the systems are shown in figure
\ref{fig:delE}.

\begin{figure}
\centering
\includegraphics[width=4.75in]{./figures/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
{\sc 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.  Appendix
\ref{app:IndividualResults} includes results for systems comprised
entirely of neutral groups.

For the {\sc sp} method, inclusion of electrostatic 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 {\sc 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 {\it
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 homogeneous
systems; although it does provide results that are an improvement over
those from an unmodified cutoff.

\section{Magnitude of the Force and Torque Vector Results}\label{sec:FTMagResults}

Evaluation of pairwise methods for use in Molecular Dynamics
simulations requires consideration of effects on the forces and
torques.  Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the
regression results for the force and torque vector magnitudes,
respectively.  The data in these figures was generated from an
accumulation of the statistics from all of the system types.

\begin{figure}
\centering
\includegraphics[width=4.75in]{./figures/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 {\sc 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}
Again, it is striking how well the Shifted Potential and Shifted Force
methods are doing at reproducing the {\sc spme} forces.  The undamped and
weakly-damped {\sc sf} method gives the best agreement with Ewald.
This is perhaps expected because this method explicitly incorporates a
smooth transition in the forces at the cutoff radius as well as the
neutralizing image charges.

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
an improvement much more significant than what was seen with $\Delta
E$.

With moderate damping and a large enough cutoff radius, the {\sc sp}
method is generating usable forces.  Further increases in damping,
while beneficial for simulations with a cutoff radius of 9~\AA\ , is
detrimental to simulations with larger cutoff radii.

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 the forces correlate
well with the Ewald forces in general.  We note that the reaction
field calculations do not include the pure NaCl systems, so these
results are partly biased towards conditions in which the method
performs more favorably.

\begin{figure}
\centering
\includegraphics[width=4.75in]{./figures/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 {\sc 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}
Molecular torques were only available from the systems which contained
rigid molecules (i.e. the systems containing water).  The data in
fig. \ref{fig:trqMag} is taken from this smaller sampling pool.

Torques appear to be much more sensitive to charges at a longer
distance.   The striking feature in comparing the new electrostatic
methods with {\sc spme} is how much the agreement improves with increasing
cutoff radius.  Again, the weakly damped and undamped {\sc sf} method
appears to reproduce the {\sc spme} torques most accurately.  

Water molecules are dipolar, and the reaction field method reproduces
the effect of the surrounding polarized medium on each of the
molecular bodies. Therefore it is not surprising that reaction field
performs best of all of the methods on molecular torques.

\section{Directionality of the Force and Torque Vector Results}\label{sec:FTDirResults}

It is clearly important that a new electrostatic method should be able
to reproduce the magnitudes of the force and torque vectors obtained
via the Ewald sum. However, the {\it directionality} of these vectors
will also be vital in calculating dynamical quantities accurately.
Force and torque directionalities were investigated by measuring the
angles formed between these vectors and the same vectors calculated
using {\sc 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=4.75in]{./figures/frcTrqAngplot.pdf}
\caption{Statistical analysis of the width of the angular distribution
that the force and torque vectors from a given electrostatic method
make with their counterparts obtained using the reference Ewald sum.
Results with a variance ($\sigma^2$) equal to zero (dashed line)
indicate force and torque directions indistinguishable from those
obtained using {\sc 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}. Here it is clear that the Shifted Potential ({\sc
sp}) method would be essentially unusable for molecular dynamics
unless the damping function is added.  The Shifted Force ({\sc sf})
method, however, is generating force and torque vectors which are
within a few degrees of the Ewald results even with weak (or no)
damping.

All of the sets (aside from the over-damped case) show the improvement
afforded by choosing a larger cutoff radius.  Increasing the cutoff
from 9 to 12~\AA\ typically results in a halving of the width of the
distribution, with a similar improvement when 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 the electrostatic damping improves the
angular behavior significantly for the {\sc sp} and moderately for the
{\sc sf} methods.  Over-damping is detrimental to both methods.  Again
it is important to recognize that the force vectors cover all
particles in all seven systems, while torque vectors are only
available for neutral molecular groups.  Damping is more beneficial to
charged bodies, and this observation is investigated further in
appendix \ref{app:IndividualResults}.

Although not discussed previously, group based cutoffs can be applied
to both the {\sc sp} and {\sc sf} methods.  The group-based cutoffs
will reintroduce small discontinuities at the cutoff radius, but the
effects of these can be minimized by utilizing a switching function.
Though there are no significant benefits or drawbacks observed in
$\Delta E$ and the force and torque magnitudes when doing this, there
is a measurable improvement in the directionality of the forces and
torques. Table \ref{tab:groupAngle} shows the angular variances
obtained both without (N) and with (Y) group based cutoffs and a
switching function.  Note that the $\alpha$ values have units of
\AA$^{-1}$ and the variance values have units of degrees$^2$.  The
{\sc sp} (with an $\alpha$ of 0.2~\AA$^{-1}$ or smaller) shows much
narrower angular distributions when using group-based cutoffs. The
{\sc sf} method likewise shows improvement in the undamped and lightly
damped cases.

\begin{table}
\caption{STATISTICAL ANALYSIS OF THE ANGULAR DISTRIBUTIONS ($\sigma^2$) 
THAT THE FORCE ({\it upper}) AND TORQUE ({\it lower}) VECTORS FROM A
GIVEN ELECTROSTATIC METHOD MAKE WITH THEIR COUNTERPARTS OBTAINED USING
THE REFERENCE EWALD SUMMATION}

\footnotesize
\begin{center}
\begin{tabular}{@{} ccrrrrrrrr @{}}
\toprule 
\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} 
\end{center}
\label{tab:groupAngle}
\end{table}

One additional trend 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 more obvious with group-based cutoffs.
The complimentary error function inserted into the potential weakens
the electrostatic interaction as the value of $\alpha$ is increased.
However, at larger values of $\alpha$, it is possible to over-damp the
electrostatic interaction and to remove it completely.  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 our
observations, empirical damping up to 0.2~\AA$^{-1}$ is beneficial,
but damping may be unnecessary when using the {\sc sf} method.


\section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}\label{sec:ShortTimeDynamics}

Zahn {\it et al.} investigated the structure and dynamics of water
using equations (\ref{eq:ZahnPot}) and
(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated
that a method similar (but not identical with) the damped {\sc sf}
method resulted in properties very similar to those obtained when
using the Ewald summation.  The properties they studied (pair
distribution functions, diffusion constants, and velocity and
orientational correlation functions) may not be particularly sensitive
to the long-range and collective behavior that governs the
low-frequency behavior in crystalline systems.  Additionally, the
ionic crystals are the worst case scenario for the pairwise methods
because they lack the reciprocal space contribution contained in the
Ewald summation.

We are using two separate measures to probe the effects of these
alternative electrostatic methods on the dynamics in crystalline
materials.  For short- and intermediate-time dynamics, we are
computing the velocity autocorrelation function, and for long-time
and large length-scale collective motions, we are looking at the
low-frequency portion of the power spectrum.

\begin{figure}
\centering
\includegraphics[width = \linewidth]{./figures/vCorrPlot.pdf}
\caption{Velocity autocorrelation functions of NaCl crystals at 
1000~K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \&
0.2~\AA$^{-1}$), and {\sc sp} ($\alpha$ = 0.2~\AA$^{-1}$). The inset is
a magnification of the area around the first minimum.  The times to
first collision are nearly identical, but differences can be seen in
the peaks and troughs, where the undamped and weakly damped methods
are stiffer than the moderately damped and {\sc spme} methods.}
\label{fig:vCorrPlot}
\end{figure}
The short-time decay of the velocity autocorrelation function 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 parameter ($\alpha$) is
increased, these peaks are smoothed out, and the {\sc sf} method
approaches the {\sc spme} results.  With $\alpha$ values of 0.2~\AA$^{-1}$,
the {\sc sf} and {\sc sp} functions are nearly identical and track the
{\sc spme} features quite well.  This is not surprising because the {\sc sf}
and {\sc sp} potentials become nearly identical with increased
damping.  However, this appears to indicate that once damping is
utilized, the details of the form of the potential (and forces)
constructed out of the damped electrostatic interaction are less
important.

\section{Collective Motion: Power Spectra of NaCl Crystals}\label{sec:LongTimeDynamics}

\begin{figure}
\centering
\includegraphics[width = \linewidth]{./figures/spectraSquare.pdf}
\caption{Power spectra obtained from the velocity auto-correlation 
functions of NaCl crystals at 1000~K while using {\sc spme}, {\sc sf}
($\alpha$ = 0, 0.1, \& 0.2~\AA$^{-1}$), and {\sc sp} ($\alpha$ =
0.2~\AA$^{-1}$).  The inset shows the frequency region below
100~cm$^{-1}$ to highlight where the spectra differ.}
\label{fig:methodPS}
\end{figure}
To evaluate how the differences between the methods affect the
collective long-time motion, we computed power spectra from long-time
traces of the velocity autocorrelation function. The power spectra for
the best-performing alternative methods are shown in
fig. \ref{fig:methodPS}.  Apodization of the correlation functions via
a cubic switching function between 40 and 50~ps was used to reduce the
ringing resulting from data truncation.  This procedure had no
noticeable effect on peak location or magnitude.

While the high frequency regions of the power spectra for the
alternative methods are quantitatively identical with Ewald spectrum,
the low frequency region shows how the summation methods differ.
Considering the low-frequency inset (expanded in the upper frame of
figure \ref{fig:methodPS}), 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 nearly
identical correlated motion to the Ewald method (which has a
convergence parameter of 0.3119~\AA$^{-1}$).  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.

\begin{figure}
\centering
\includegraphics[width = \linewidth]{./figures/increasedDamping.pdf}
\caption{Effect of damping on the two lowest-frequency phonon modes in
the NaCl crystal at 1000~K.  The undamped shifted force ({\sc sf})
method is off by less than 10~cm$^{-1}$, and increasing the
electrostatic damping to 0.25~\AA$^{-1}$ gives quantitative agreement
with the power spectrum obtained using the Ewald sum.  Over-damping can
result in underestimates of frequencies of the long-wavelength
motions.}
\label{fig:dampInc}
\end{figure}
To isolate the role of the damping constant, we have computed the
spectra for a single method ({\sc sf}) with a range of damping
constants and compared this with the {\sc spme} spectrum.
Fig. \ref{fig:dampInc} shows more clearly that increasing the
electrostatic damping red-shifts the lowest frequency phonon modes.
However, even without any electrostatic damping, the {\sc sf} method
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode.
Without the {\sc sf} modifications, an undamped (pure cutoff) method
would predict the lowest frequency peak near 325~cm$^{-1}$.  {\it
Most} of the collective behavior in the crystal is accurately captured
using the {\sc sf} method.  Quantitative agreement with Ewald can be
obtained using moderate damping in addition to the shifting at the
cutoff distance.

\section{An Application: TIP5P-E Water}\label{sec:t5peApplied}

The above sections focused on the energetics and dynamics of a variety
of systems when utilizing the {\sc sp} and {\sc sf} pairwise
techniques.  A unitary correlation with results obtained using the
Ewald summation should result in a successful reproduction of both the
static and dynamic properties of any selected system.  To test this,
we decided to calculate a series of properties for the TIP5P-E water
model when using the {\sc sf} technique.

The TIP5P-E water model is a variant of Mahoney and Jorgensen's
five-point transferable intermolecular potential (TIP5P) model for
water.\cite{Mahoney00} TIP5P was developed to reproduce the density
maximum anomaly present in liquid water near 4$^\circ$C. As with many
previous point charge water models (such as ST2, TIP3P, TIP4P, SPC,
and SPC/E), TIP5P was parametrized using a simple cutoff with no
long-range electrostatic
correction.\cite{Stillinger74,Jorgensen83,Berendsen81,Berendsen87}
Without this correction, the pressure term on the central particle
from the surroundings is missing. Because they expand to compensate
for this added pressure term when this correction is included, systems
composed of these particles tend to under-predict the density of water
under standard conditions. When using any form of long-range
electrostatic correction, it has become common practice to develop or
utilize a reparametrized water model that corrects for this
effect.\cite{vanderSpoel98,Fennell04,Horn04} The TIP5P-E model follows
this practice and was optimized specifically for use with the Ewald
summation.\cite{Rick04} In his publication, Rick preserved the
geometry and point charge magnitudes in TIP5P and focused on altering
the Lennard-Jones parameters to correct the density at
298K.\cite{Rick04} With the density corrected, he compared common
water properties for TIP5P-E using the Ewald sum with TIP5P using a
9~\AA\ cutoff.

In the following sections, we compared these same water properties
calculated from TIP5P-E using the Ewald sum with TIP5P-E using the
{\sc sf} technique.  In the above evaluation of the pairwise
techniques, we observed some flexibility in the choice of parameters.
Because of this, the following comparisons include the {\sc sf}
technique with a 12~\AA\ cutoff and an $\alpha$ = 0.0, 0.1, and
0.2~\AA$^{-1}$, as well as a 9~\AA\ cutoff with an $\alpha$ =
0.2~\AA$^{-1}$.

\subsection{Density}\label{sec:t5peDensity}

As stated previously, the property that prompted the development of
TIP5P-E was the density at 1 atm.  The density depends upon the
internal pressure of the system in the $NPT$ ensemble, and the
calculation of the pressure includes a components from both the
kinetic energy and the virial. More specifically, the instantaneous
molecular pressure ($p(t)$) is given by
\begin{equation}
p(t) =  \frac{1}{\textrm{d}V}\sum_\mu
        \left[\frac{\mathbf{P}_{\mu}^2}{M_{\mu}} 
        + \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right],
\label{eq:MolecularPressure}
\end{equation}
where d is the dimensionality of the system, $V$ is the volume,
$\mathbf{P}_{\mu}$ is the momentum of molecule $\mu$, $\mathbf{R}_\mu$
is the position of the center of mass ($M_\mu$) of molecule $\mu$, and
$\mathbf{F}_{\mu i}$ is the force on atom $i$ of molecule
$\mu$.\cite{Melchionna93} The virial term (the right term in the
brackets of equation
\ref{eq:MolecularPressure}) is directly dependent on the interatomic
forces.  Since the {\sc sp} method does not modify the forces (see
section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp}
will be identical to that obtained without an electrostatic
correction.  The {\sc sf} method does alter the virial component and,
by way of the modified pressures, should provide densities more in
line with those obtained using the Ewald summation.

To compare densities, $NPT$ simulations were performed with the same
temperatures as those selected by Rick in his Ewald summation
simulations.\cite{Rick04} In order to improve statistics around the
density maximum, 3~ns trajectories were accumulated at 0, 12.5, and
25$^\circ$C, while 2~ns trajectories were obtained at all other
temperatures. The average densities were calculated from the later
three-fourths of each trajectory. Similar to Mahoney and Jorgensen's
method for accumulating statistics, these sequences were spliced into
200 segments to calculate the average density and standard deviation
at each temperature.\cite{Mahoney00}

\begin{figure}
\includegraphics[width=\linewidth]{./figures/tip5peDensities.pdf}
\caption{Density versus temperature for the TIP5P-E water model when 
using the Ewald summation (Ref. \cite{Rick04}) and the {\sc sf} method
with various parameters. The pressure term from the image-charge shell
is larger than that provided by the reciprocal-space portion of the
Ewald summation, leading to slightly lower densities. This effect is
more visible with the 9~\AA\ cutoff, where the image charges exert a
greater force on the central particle. The error bars for the {\sc sf}
methods show the average one-sigma uncertainty of the density
measurement, and this uncertainty is the same for all the {\sc sf}
curves.}
\label{fig:t5peDensities}
\end{figure}
Figure \ref{fig:t5peDensities} shows the densities calculated for
TIP5P-E using differing electrostatic corrections overlaid on the
experimental values.\cite{CRC80} The densities when using the {\sc sf}
technique are close to, though typically lower than, those calculated
while using the Ewald summation. These slightly reduced densities
indicate that the pressure component from the image charges at
R$_\textrm{c}$ is larger than that exerted by the reciprocal-space
portion of the Ewald summation. Bringing the image charges closer to
the central particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than
the preferred 12~\AA\ R$_\textrm{c}$) increases the strength of their
interactions, resulting in a further reduction of the densities.

Because the strength of the image charge interactions has a noticeable
effect on the density, we would expect the use of electrostatic
damping to also play a role in these calculations. Larger values of
$\alpha$ weaken the pair-interactions; and since electrostatic damping
is distance-dependent, force components from the image charges will be
reduced more than those from particles close the the central
charge. This effect is visible in figure \ref{fig:t5peDensities} with
the damped {\sc sf} sums showing slightly higher densities; however,
it is apparent that the choice of cutoff radius plays a much more
important role in the resulting densities.

As a final note, all of the above density calculations were performed
with systems of 512 water molecules. Rick observed a system size
dependence of the computed densities when using the Ewald summation,
most likely due to his tying of the convergence parameter to the box
dimensions.\cite{Rick04} For systems of 256 water molecules, the
calculated densities were on average 0.002 to 0.003~g/cm$^3$ lower. A
system size of 256 molecules would force the use of a shorter
R$_\textrm{c}$ when using the {\sc sf} technique, and this would also
lower the densities. Moving to larger systems, as long as the
R$_\textrm{c}$ remains at a fixed value, we would expect the densities
to remain constant.

\subsection{Liquid Structure}\label{sec:t5peLiqStructure}

A common function considered when developing and comparing water
models is the oxygen-oxygen radial distribution function
($g_\textrm{OO}(r)$). The $g_\textrm{OO}(r)$ is the probability of
finding a pair of oxygen atoms some distance ($r$) apart relative to a
random distribution at the same density.\cite{Allen87} It is
calculated via
\begin{equation}
g_\textrm{OO}(r) = \frac{V}{N^2}\left\langle\sum_i\sum_{j\ne i}
        \delta(\mathbf{r}-\mathbf{r}_{ij})\right\rangle,
\label{eq:GOOofR}
\end{equation}
where the double sum is over all $i$ and $j$ pairs of $N$ oxygen
atoms. The $g_\textrm{OO}(r)$ can be directly compared to X-ray or
neutron scattering experiments through the oxygen-oxygen structure
factor ($S_\textrm{OO}(k)$) by the following relationship:
\begin{equation}
S_\textrm{OO}(k) = 1 + 4\pi\rho
        \int_0^\infty r^2\frac{\sin kr}{kr}g_\textrm{OO}(r)\textrm{d}r.
\label{eq:SOOofK}
\end{equation}
Thus, $S_\textrm{OO}(k)$ is related to the Fourier-Laplace transform
of $g_\textrm{OO}(r)$.

The experimentally determined $g_\textrm{OO}(r)$ for liquid water has
been compared in great detail with the various common water models,
and TIP5P was found to be in better agreement than other rigid,
non-polarizable models.\cite{Sorenson00} This excellent agreement with
experiment was maintained when Rick developed TIP5P-E.\cite{Rick04} To
check whether the choice of using the Ewald summation or the {\sc sf}
technique alters the liquid structure, the $g_\textrm{OO}(r)$s at 298K
and 1atm were determined for the systems compared in the previous
section.

\begin{figure}
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf}
\caption{The $g_\textrm{OO}(r)$s calculated for TIP5P-E at 298~K and 
1atm while using the Ewald summation (Ref. \cite{Rick04}) and the {\sc
sf} technique with varying parameters. Even with the reduced densities
using the {\sc sf} technique, the $g_\textrm{OO}(r)$s are essentially
identical.}
\label{fig:t5peGofRs}
\end{figure}
The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc
sf} technique with a various parameters are overlaid on the
$g_\textrm{OO}(r)$ while using the Ewald summation in figure
\ref{fig:t5peGofRs}. The differences in density do not appear to have
any effect on the liquid structure as the $g_\textrm{OO}(r)$s are
indistinguishable. These results indicate that the $g_\textrm{OO}(r)$
is insensitive to the choice of electrostatic correction.

\subsection{Thermodynamic Properties}\label{sec:t5peThermo}

In addition to the density, there are a variety of thermodynamic
quantities that can be calculated for water and compared directly to
experimental values. Some of these additional quantities include the
latent heat of vaporization ($\Delta H_\textrm{vap}$), the constant
pressure heat capacity ($C_p$), the isothermal compressibility
($\kappa_T$), the thermal expansivity ($\alpha_p$), and the static
dielectric constant ($\epsilon$). All of these properties were
calculated for TIP5P-E with the Ewald summation, so they provide a
good set for comparisons involving the {\sc sf} technique.

The $\Delta H_\textrm{vap}$ is the enthalpy change required to
transform one mol of substance from the liquid phase to the gas
phase.\cite{Berry00} In molecular simulations, this quantity can be
determined via
\begin{equation}
\begin{split}
\Delta H_\textrm{vap} &= H_\textrm{gas} - H_\textrm{liq.} \\
                      &= E_\textrm{gas} - E_\textrm{liq.} 
                        + p(V_\textrm{gas} - V_\textrm{liq.}) \\
                      &\approx -\frac{\langle U_\textrm{liq.}\rangle}{N}+ RT,
\end{split}
\label{eq:DeltaHVap}
\end{equation}
where $E$ is the total energy, $U$ is the potential energy, $p$ is the
pressure, $V$ is the volume, $N$ is the number of molecules, $R$ is
the gas constant, and $T$ is the temperature.\cite{Jorgensen98b} As
seen in the last line of equation (\ref{eq:DeltaHVap}), we can
approximate $\Delta H_\textrm{vap}$ by using the ideal gas for the gas
state. This allows us to cancel the kinetic energy terms, leaving only
the $U_\textrm{liq.}$ term. Additionally, the $pV$ term for the gas is
several orders of magnitude larger than that of the liquid, so we can
neglect the liquid $pV$ term.

The remaining thermodynamic properties can all be calculated from
fluctuations of the enthalpy, volume, and system dipole
moment.\cite{Allen87} $C_p$ can be calculated from fluctuations in the
enthalpy in constant pressure simulations via
\begin{equation}
\begin{split}
C_p     = \left(\frac{\partial H}{\partial T}\right)_{N,p} 
        = \frac{(\langle H^2\rangle - \langle H\rangle^2)}{Nk_BT^2},
\end{split}
\label{eq:Cp}
\end{equation}
where $k_B$ is Boltzmann's constant.  $\kappa_T$ can be calculated via
\begin{equation}
\begin{split}
\kappa_T = \frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{N,T}
         = \frac{(\langle V^2\rangle_{N,P,T} - \langle V\rangle^2_{N,P,T})}
                {k_BT\langle V\rangle_{N,P,T}},
\end{split}
\label{eq:kappa}
\end{equation}
and $\alpha_p$ can be calculated via
\begin{equation}
\begin{split}
\alpha_p        = \frac{1}{V}\left(\frac{\partial V}{\partial T}\right)_{N,P}
                = \frac{(\langle VH\rangle_{N,P,T} 
                        - \langle V\rangle_{N,P,T}\langle H\rangle_{N,P,T})}
                        {k_BT^2\langle V\rangle_{N,P,T}}.
\end{split}
\label{eq:alpha}
\end{equation}
Using the Ewald sum under tin-foil boundary conditions, $\epsilon$ can
be calculated for systems of non-polarizable substances via
\begin{equation}
\epsilon = 1 + \frac{\langle M^2\rangle}{3\epsilon_0k_\textrm{B}TV},
\label{eq:staticDielectric}
\end{equation}
where $\epsilon_0$ is the permittivity of free space and $\langle
M^2\rangle$ is the fluctuation of the system dipole
moment.\cite{Allen87} The numerator in the fractional term in equation
(\ref{eq:staticDielectric}) is the fluctuation of the simulation-box
dipole moment, identical to the quantity calculated in the
finite-system Kirkwood $g$ factor ($G_k$):
\begin{equation}
G_k = \frac{\langle M^2\rangle}{N\mu^2},
\label{eq:KirkwoodFactor}
\end{equation}
where $\mu$ is the dipole moment of a single molecule of the
homogeneous system.\cite{Neumann83,Neumann84,Neumann85} The
fluctuation term in both equation (\ref{eq:staticDielectric}) and
\ref{eq:KirkwoodFactor} is calculated as follows,
\begin{equation}
\begin{split}
\langle M^2\rangle &= \langle\bm{M}\cdot\bm{M}\rangle 
                        - \langle\bm{M}\rangle\cdot\langle\bm{M}\rangle \\
                   &= \langle M_x^2+M_y^2+M_z^2\rangle
                        - (\langle M_x\rangle^2 + \langle M_x\rangle^2 
                                + \langle M_x\rangle^2).        
\end{split}
\label{eq:fluctBoxDipole}
\end{equation}
This fluctuation term can be accumulated during the simulation;
however, it converges rather slowly, thus requiring multi-nanosecond
simulation times.\cite{Horn04} In the case of tin-foil boundary
conditions, the dielectric/surface term of equation (\ref{eq:EwaldSum})
is equal to zero. Since the {\sc sf} method also lacks this
dielectric/surface term, equation (\ref{eq:staticDielectric}) is still
valid for determining static dielectric constants.

All of the above properties were calculated from the same trajectories
used to determine the densities in section \ref{sec:t5peDensity}
except for the static dielectric constants. The $\epsilon$ values were
accumulated from 2~ns $NVE$ ensemble trajectories with system densities
fixed at the average values from the $NPT$ simulations at each of the
temperatures. The resulting values are displayed in figure
\ref{fig:t5peThermo}.
\begin{figure}
\centering
\includegraphics[width=4.5in]{./figures/t5peThermo.pdf}
\caption{Thermodynamic properties for TIP5P-E using the Ewald summation 
and the {\sc sf} techniques along with the experimental values. Units
for the properties are kcal mol$^{-1}$ for $\Delta H_\textrm{vap}$,
cal mol$^{-1}$ K$^{-1}$ for $C_p$, 10$^6$ atm$^{-1}$ for $\kappa_T$,
and 10$^5$ K$^{-1}$ for $\alpha_p$. Ewald summation results are from
reference \cite{Rick04}. Experimental values for $\Delta
H_\textrm{vap}$, $\kappa_T$, and $\alpha_p$ are from reference
\cite{Kell75}. Experimental values for $C_p$ are from reference
\cite{Wagner02}. Experimental values for $\epsilon$ are from reference
\cite{Malmberg56}.}
\label{fig:t5peThermo}
\end{figure}

As observed for the density in section \ref{sec:t5peDensity}, the
property trends with temperature seen when using the Ewald summation
are reproduced with the {\sc sf} technique. One noticable difference
between the properties calculated using the two methods are the lower
$\Delta H_\textrm{vap}$ values when using {\sc sf}. This is to be
expected due to the direct weakening of the electrostatic interaction
through forced neutralization. This results in an increase of the
intermolecular potential producing lower values from equation
(\ref{eq:DeltaHVap}). The slopes of these values with temperature are
similar to that seen using the Ewald summation; however, they are both
steeper than the experimental trend, indirectly resulting in the
inflated $C_p$ values at all temperatures.

Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$
values all overlap within error. As indicated for the $\Delta
H_\textrm{vap}$ and $C_p$ results discussed in the previous paragraph,
the deviations between experiment and simulation in this region are
not the fault of the electrostatic summation methods but are due to
the TIP5P class model itself. Like most rigid, non-polarizable,
point-charge water models, the density decreases with temperature at a
much faster rate than experiment (see figure
\ref{fig:t5peDensities}). The reduced density leads to the inflated
compressibility and expansivity values at higher temperatures seen
here in figure \ref{fig:t5peThermo}. Incorporation of polarizability
and many-body effects are required in order for simulation to overcome
these differences with experiment.\cite{Laasonen93,Donchev06}

At temperatures below the freezing point for experimental water, the
differences between {\sc sf} and the Ewald summation results are more
apparent. The larger $C_p$ and lower $\alpha_p$ values in this region
indicate a more pronounced transition in the supercooled regime,
particularly in the case of {\sc sf} without damping. This points to
the onset of a more frustrated or glassy behavior for TIP5P-E at
temperatures below 250~K in the {\sc sf} simulations, indicating that
disorder in the reciprical-space term of the Ewald summation might act
to loosen up the local structure more than the image-charges in {\sc
sf}. Because the systems are locked in different regions of
phase-space, comparisons between properties at these temperatures are
not exactly fair. This observation is explored in more detail in
section \ref{sec:t5peDynamics}.

The final thermodynamic property displayed in figure
\ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy
between the Ewald summation and the {\sc sf} technique (and experiment
for that matter). It is known that the dielectric constant is
dependent upon and quite sensitive to the imposed boundary
conditions.\cite{Neumann80,Neumann83} This is readily apparent in the
converged $\epsilon$ values accumulated for the {\sc sf}
simulations. Lack of a damping function results in dielectric
constants significantly smaller than that obtained using the Ewald
sum. Increasing the damping coefficient to 0.2~\AA$^{-1}$ improves the
agreement considerably. It should be noted that the choice of the
``Ewald coefficient'' value also has a significant effect on the
calculated value when using the Ewald summation. In the simulations of
TIP5P-E with the Ewald sum, this screening parameter was tethered to
the simulation box size (as was the $R_\textrm{c}$).\cite{Rick04} In
general, systems with larger screening parameters reported larger
dielectric constant values, the same behavior we see here with {\sc
sf}; however, the choice of cutoff radius also plays an important
role. In section \ref{sec:dampingDielectric}, this connection is
further explored as optimal damping coefficients for different choices
of $R_\textrm{c}$ are determined for {\sc sf} for capturing the
dielectric behavior.

\subsection{Dynamic Properties}\label{sec:t5peDynamics}

To look at the dynamic properties of TIP5P-E when using the {\sc sf}
method, 200~ps $NVE$ simulations were performed for each temperature at
the average density reported by the $NPT$ simulations. The
self-diffusion constants ($D$) were calculated with the Einstein
relation using the mean square displacement (MSD),
\begin{equation}
D = \lim_{t\rightarrow\infty} 
    \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
\label{eq:MSD}
\end{equation}
where $t$ is time, and $\mathbf{r}_i$ is the position of particle
$i$.\cite{Allen87} Figure \ref{fig:ExampleMSD} shows an example MSD
plot. As labeled in the figure, MSD plots consist of three distinct
regions:

\begin{enumerate}[itemsep=0pt]
\item parabolic short-time ballistic motion, 
\item linear diffusive regime, and
\item a region with poor statistics. 
\end{enumerate}
The slope from the linear region (region 2) is used to calculate $D$.
\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/ExampleMSD.pdf}
\caption{Example plot of mean square displacement verses time. The 
left red region is the ballistic motion regime, the middle green
region is the linear diffusive regime, and the right blue region is
the region with poor statistics.}
\label{fig:ExampleMSD}
\end{figure}

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/waterFrame.pdf}
\caption{Body-fixed coordinate frame for a water molecule. The 
respective molecular principle axes point in the direction of the
labeled frame axes.}
\label{fig:waterFrame}
\end{figure}
In addition to translational diffusion, reorientational time constants
were calculated for comparisons with the Ewald simulations and with
experiments. These values were determined from 25~ps $NVE$ trajectories
through calculation of the orientational time correlation function,
\begin{equation}
C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\alpha(t)
                \cdot\hat{\mathbf{u}}_i^\alpha(0)\right]\right\rangle,
\label{eq:OrientCorr}
\end{equation}
where $P_l$ is the Legendre polynomial of order $l$ and
$\hat{\mathbf{u}}_i^\alpha$ is the unit vector of molecule $i$ along
principle axis $\alpha$. The principle axis frame for these water
molecules is shown in figure \ref{fig:waterFrame}. As an example,
$C_l^y$ is calculated from the time evolution of the unit vector
connecting the two hydrogen atoms.

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/ExampleOrientCorr.pdf}
\caption{Example plots of the orientational autocorrelation functions 
for the first and second Legendre polynomials. These curves show the
time decay of the unit vector along the $y$ principle axis.}
\label{fig:OrientCorr}
\end{figure}
From the orientation autocorrelation functions, we can obtain time
constants for rotational relaxation. Figure \ref{fig:OrientCorr} shows
some example plots of orientational autocorrelation functions for the
first and second Legendre polynomials. The relatively short time
portions (between 1 and 3~ps for water) of these curves can be fit to
an exponential decay to obtain these constants, and they are directly
comparable to water orientational relaxation times from nuclear
magnetic resonance (NMR). The relaxation constant obtained from
$C_2^y(t)$ is of particular interest because it describes the
relaxation of the principle axis connecting the hydrogen atoms. Thus,
$C_2^y(t)$ can be compared to the intermolecular portion of the
dipole-dipole relaxation from a proton NMR signal and should provide
the best estimate of the NMR relaxation time constant.\cite{Impey82}

\begin{figure}
\centering
\includegraphics[width=3.5in]{./figures/t5peDynamics.pdf}
\caption{Diffusion constants ({\it upper}) and reorientational time 
constants ({\it lower}) for TIP5P-E using the Ewald sum and {\sc sf}
technique compared with experiment. Data at temperatures less that
0$^\circ$C were not included in the $\tau_2^y$ plot to allow for
easier comparisons in the more relevant temperature regime.}
\label{fig:t5peDynamics}
\end{figure}
Results for the diffusion constants and orientational relaxation times
are shown in figure \ref{fig:t5peDynamics}. From this figure, it is
apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using
the Ewald sum are reproduced with the {\sc sf} technique. The enhanced
diffusion at high temperatures are again the product of the lower
densities in comparison with experiment and do not provide any special
insight into differences between the electrostatic summation
techniques. With the undamped {\sc sf} technique, TIP5P-E tends to
diffuse a little faster than with the Ewald sum; however, use of light
to moderate damping results in indistinguishable $D$ values. Though
not apparent in this figure, {\sc sf} values at the lowest temperature
are approximately an order of magnitude lower than with Ewald. These
values support the observation from section \ref{sec:t5peThermo} that
there appeared to be a change to a more glassy-like phase with the
{\sc sf} technique at these lower temperatures.

The $\tau_2^y$ results in the lower frame of figure
\ref{fig:t5peDynamics} show a much greater difference between the {\sc
sf} results and the Ewald results. At all temperatures shown, TIP5P-E
relaxes faster than experiment with the Ewald sum while tracking
experiment fairly well when using the {\sc sf} technique, independent
of the choice of damping constant. Their are several possible reasons
for this deviation between techniques. The Ewald results were taken
shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick
calculation from a 10~ps trajectory with {\sc sf} with an $\alpha$ of
0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$,
placing the result more in line with that obtained using the Ewald
sum. These results support this explanation; however, recomputing the
results to meet a poorer statistical standard is
counter-productive. Assuming the Ewald results are not the product of
poor statistics, differences in techniques to integrate the
orientational motion could also play a role. {\sc shake} is the most
commonly used technique for approximating rigid-body orientational
motion,\cite{Ryckaert77} where as in {\sc oopse}, we maintain and
integrate the entire rotation matrix using the {\sc dlm}
method.\cite{Meineke05} Since {\sc shake} is an iterative constraint
technique, if the convergence tolerances are raised for increased
performance, error will accumulate in the orientational
motion. Finally, the Ewald results were calculated using the $NVT$
ensemble, while the $NVE$ ensemble was used for {\sc sf}
calculations. The additional mode of motion due to the thermostat will
alter the dynamics, resulting in differences between $NVT$ and $NVE$
results. These differences are increasingly noticeable as the
thermostat time constant decreases.

\section{Damping of Point Multipoles}\label{sec:dampingMultipoles}

As discussed above, the {\sc sp} and {\sc sf} methods operate by
neutralizing the cutoff sphere with charge-charge interaction shifting
and by damping the electrostatic interactions. Now we would like to
consider an extension of these techniques to include point multipole
interactions. How will the shifting and damping need to be modified in
order to accommodate point multipoles?

Of the two techniques, the easiest to adapt is shifting. Shifting is
employed to neutralize the cutoff sphere; however, in a system
composed purely of point multipoles, the cutoff sphere is already
neutralized. This means that shifting is not necessary between point
multipoles. In a mixed system of monopoles and multipoles, the
undamped {\sc sf} potential needs only to shift the force terms of the
monopole (and use the monopole potential of equation (\ref{eq:SFPot}))
and smoothly cutoff the multipole interactions with a switching
function. The switching function is required in order to conserve
energy, because a discontinuity will exist at $R_\textrm{c}$ in the
absence of shifting terms.

If we consider damping the {\sc sf} potential (Eq. (\ref{eq:DSFPot})),
then we need to incorporate the complimentary error function term into
the multipole potentials. The most direct way to do this is by
replacing $r^{-1}$ with erfc$(\alpha r)\cdot r^{-1}$ in the multipole
expansion.\cite{Hirschfelder67} In the multipole expansion, rather
than considering only the interactions between single point charges,
the electrostatic interaction is reformulated such that it describes
the interaction between charge distributions about central sites of
the respective sets of charges. This procedure is what leads to the
familiar charge-dipole,
\begin{equation}
V_\textrm{cd} = -q_i\frac{\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij}}{r^3_{ij}}
                = -q_i\mu_j\frac{\cos\theta}{r^2_{ij}}, 
\label{eq:chargeDipole}
\end{equation}
and dipole-dipole,
\begin{equation}
V_\textrm{dd} = 3\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})
                        (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r^5_{ij}} -
                \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r^3_{ij}},
\label{eq:dipoleDipole}
\end{equation}
interaction potentials.

Using the charge-dipole interaction as an example, if we insert
erfc$(\alpha r)\cdot r^{-1}$ in place of $r^{-1}$, a damped
charge-dipole results,
\begin{equation}
V_\textrm{Dcd} = -q_i\mu_j\frac{\cos\theta}{r^2_{ij}} c_1(r_{ij}),
\label{eq:dChargeDipole}
\end{equation}
where $c_1(r_{ij})$ is
\begin{equation}
c_1(r_{ij}) = \frac{2\alpha r_{ij}e^{-\alpha^2r^2_{ij}}}{\sqrt{\pi}}
                + \textrm{erfc}(\alpha r_{ij}).
\label{eq:c1Func}
\end{equation}
Thus, $c_1(r_{ij})$ is the resulting damping term that modifies the
standard charge-dipole potential (Eq. (\ref{eq:chargeDipole})).  Note
that this damping term is dependent upon distance and not upon
orientation, and that it is acting on what was originally an
$r^{-3}$ function. By writing the damped form in this manner, we
can collect the damping into one function and apply it to the original
potential when damping is desired. This works well for potentials that
have only one $r^{-n}$ term (where $n$ is an odd positive integer);
but in the case of the dipole-dipole potential, there is one part
dependent on $r^{-3}$ and another dependent on $r^{-5}$. In order to
properly damping this potential, each of these parts is dampened with
separate damping functions. We can determine the necessary damping
functions by continuing with the multipole expansion; however, it
quickly becomes more complex with ``two-center'' systems, like the
dipole-dipole potential, and is typically approached with a spherical
harmonic formalism.\cite{Hirschfelder67} A simpler method for
determining these functions arises from adopting the tensor formalism
for expressing the electrostatic interactions.\cite{Stone02}

The tensor formalism for electrostatic interactions involves obtaining
the multipole interactions from successive gradients of the monopole
potential. Thus, tensors of rank one through three are
\begin{equation}
T = \frac{1}{4\pi\epsilon_0r_{ij}},
\label{eq:tensorRank1}
\end{equation}
\begin{equation}
T_\alpha = \frac{1}{4\pi\epsilon_0}\nabla_\alpha \frac{1}{r_{ij}},
\label{eq:tensorRank2}
\end{equation}
\begin{equation}
T_{\alpha\beta} = \frac{1}{4\pi\epsilon_0}
                        \nabla_\alpha\nabla_\beta \frac{1}{r_{ij}},
\label{eq:tensorRank3}
\end{equation}
where the form of the first tensor gives the monopole-monopole
potential, the second gives the monopole-dipole potential, and the
third gives the monopole-quadrupole and dipole-dipole
potentials.\cite{Stone02} Since the force is $-\nabla V$, the forces
for each potential come from the next higher tensor.

To obtain the damped electrostatic forms, we replace $r^{-1}$ with
erfc$(\alpha r)\cdot r^{-1}$. Equation \ref{eq:tensorRank2} generates
$c_1(r_{ij})$, just like the multipole expansion, while equation
\ref{eq:tensorRank3} generates $c_2(r_{ij})$, where
\begin{equation}
c_2(r_{ij}) = \frac{4\alpha^3r^3_{ij}e^{-\alpha^2r^2_{ij}}}{3\sqrt{\pi}}
                + \frac{2\alpha r_{ij}e^{-\alpha^2r^2_{ij}}}{\sqrt{\pi}}
                + \textrm{erfc}(\alpha r_{ij}).
\end{equation} 
Note that $c_2(r_{ij})$ is equal to $c_1(r_{ij})$ plus an additional
term. Continuing with higher rank tensors, we can obtain the damping
functions for higher multipoles as well as the forces. Each subsequent
damping function includes one additional term, and we can simplify the
procedure for obtaining these terms by writing out the following
generating function,
\begin{equation}
c_n(r_{ij}) = \frac{2^n(\alpha r_{ij})^{2n-1}e^{-\alpha^2r^2_{ij}}}
                {(2n-1)!!\sqrt{\pi}} + c_{n-1}(r_{ij}),
\label{eq:dampingGeneratingFunc}
\end{equation}
where,
\begin{equation}
m!! \equiv \left\{ \begin{array}{l@{\quad\quad}l}
                m\cdot(m-2)\ldots 5\cdot3\cdot1 & m > 0\textrm{ odd} \\
                m\cdot(m-2)\ldots 6\cdot4\cdot2 & m > 0\textrm{ even} \\
                1 & m = -1\textrm{ or }0,
                \end{array}\right.
\label{eq:doubleFactorial}              
\end{equation}
and $c_0(r_{ij})$ is erfc$(\alpha r_{ij})$. This generating function
is similar in form to those obtained by researchers for the
application of the Ewald sum to
multipoles.\cite{Smith82,Smith98,Aguado03}

Returning to the dipole-dipole example, the potential consists of a
portion dependent upon $r^{-5}$ and another dependent upon
$r^{-3}$. In the damped dipole-dipole potential,
\begin{equation}
V_\textrm{Ddd} = 3\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})
                (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})}{r^5_{ij}}
                c_2(r_{ij}) -
                \frac{\boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j}{r^3_{ij}}
                c_1(r_{ij}),
\label{eq:dampDipoleDipole}
\end{equation}
$c_2(r_{ij})$ and $c_1(r_{ij})$ respectively dampen these two
parts. The forces for the damped dipole-dipole interaction,
\begin{equation}
\begin{split}
F_\textrm{Ddd} = &15\frac{(\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij})
        (\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij})\mathbf{r}_{ij}}{r^7_{ij}} 
        c_3(r_{ij})\\ &- 
3\frac{\boldsymbol{\mu}_i\cdot\mathbf{r}_{ij}\cdot\hat{\boldsymbol{\mu}}_j +
        \boldsymbol{\mu}_j\cdot\mathbf{r}_{ij}\cdot\hat{\boldsymbol{\mu}}_i +
        \boldsymbol{\mu}_i\cdot\boldsymbol{\mu}_j\cdot\mathbf{r}_{ij}}
        {r^5_{ij}} c_2(r_{ij}),
\end{split}
\label{eq:dampDipoleDipoleForces}
\end{equation}
rely on higher order damping functions because we perform another
gradient operation. In this manner, we can dampen higher order
multipolar interactions along with the monopole interactions, allowing
us to include multipoles in simulations involving damped electrostatic
interactions.


\section{Damping and Dielectric Constants}\label{sec:dampingDielectric}

In section \ref{sec:t5peThermo}, we observed that the choice of
damping coefficient plays a major role in the calculated dielectric
constant. This is not too surprising given the results for damping
parameter influence on the long-time correlated motions of the NaCl
crystal in section \ref{sec:LongTimeDynamics}. The static dielectric
constant is calculated from the long-time fluctuations of the system's
accumulated dipole moment (Eq. (\ref{eq:staticDielectric})), so it is
going to be quite sensitive to the choice of damping parameter. We
would like to choose an optimal damping constant for any particular
cutoff radius choice that would properly capture the dielectric
behavior of the liquid.

In order to find these optimal values, we mapped out the static
dielectric constant as a function of both the damping parameter and
cutoff radius for several different water models. To calculate the
static dielectric constant, we performed 5~ns $NPT$ calculations on
systems of 512 water molecules, using the TIP5P-E, TIP4P-Ew, SPC/E,
and SSD/RF water models. TIP4P-Ew is a reparametrized version of the
four-point transferable intermolecular potential (TIP4P) for water
targeted for use with the Ewald summation.\cite{Horn04} SSD/RF is the
reaction field modified variant of the soft sticky dipole (SSD) model
for water\cite{Fennell04} This model is discussed in more detail in
the next chapter. One thing to note about it, electrostatic
interactions are handled via dipole-dipole interactions rather than
charge-charge interactions like the other three models. Damping of the
dipole-dipole interaction was handled as described in section
\ref{sec:dampingMultipoles}. Each of these systems were studied with
cutoff radii of 9, 10, 11, and 12~\AA\ and with damping parameter values
ranging from 0 to 0.35~\AA$^{-1}$.

\begin{figure}
\centering
\includegraphics[width=\linewidth]{./figures/dielectricMap.pdf}
\caption{The static dielectric constant for the TIP5P-E (A), TIP4P-Ew 
(B), SPC/E (C), and SSD/RF (D) water models as a function of cutoff
radius ($R_\textrm{c}$) and damping coefficient ($\alpha$).}
\label{fig:dielectricMap}
\end{figure}
The results of these calculations are displayed in figure
\ref{fig:dielectricMap} in the form of shaded contour plots. An
interesting aspect of all four contour plots is that the dielectric
constant is effectively linear with respect to $\alpha$ and
$R_\textrm{c}$ in the low to moderate damping regions, and the slope
is 0.025~\AA$^{-1}\cdot R_\textrm{c}^{-1}$. Another point to note is
that choosing $\alpha$ and $R_\textrm{c}$ identical to those used in
studies with the Ewald summation results in the same calculated
dielectric constant. As an example, in the paper outlining the
development of TIP5P-E, the real-space cutoff and Ewald coefficient
were tethered to the system size, and for a 512 molecule system are
approximately 12~\AA\ and 0.25~\AA$^{-1}$ respectively.\cite{Rick04}
These parameters resulted in a dielectric constant of 92$\pm$14, while
with {\sc sf} these parameters give a dielectric constant of
90.8$\pm$0.9. Another example comes from the TIP4P-Ew paper where
$\alpha$ and $R_\textrm{c}$ were chosen to be 9.5~\AA\ and
0.35~\AA$^{-1}$, and these parameters resulted in a $\epsilon_0$ equal
to 63$\pm$1.\cite{Horn04} We did not perform calculations with these
exact parameters, but interpolating between surrounding values gives a
$\epsilon_0$ of 61$\pm$1. Seeing a dependence of the dielectric
constant on $\alpha$ and $R_\textrm{c}$ with the {\sc sf} technique,
it might be interesting to investigate the dielectric dependence of
the real-space Ewald parameters.

Although it is tempting to choose damping parameters equivalent to
these Ewald examples, the results discussed in sections
\ref{sec:EnergyResults} through \ref{sec:FTDirResults} and appendix
\ref{app:IndividualResults} indicate that values this high are
destructive to both the energetics and dynamics. Ideally, $\alpha$
should not exceed 0.3~\AA$^{-1}$ for any of the cutoff values in this
range. If the optimal damping parameter is chosen to be midway between
0.275 and 0.3~\AA$^{-1}$ (0.2875~\AA$^{-1}$) at the 9~\AA\ cutoff,
then the linear trend with $R_\textrm{c}$ will always keep $\alpha$
below 0.3~\AA$^{-1}$. This linear progression would give values of
0.2875, 0.2625, 0.2375, and 0.2125~\AA$^{-1}$ for cutoff radii of 9,
10, 11, and 12~\AA. Setting this to be the default behavior for the
damped {\sc sf} technique will result in consistent dielectric
behavior for these and other condensed molecular systems, regardless
of the chosen cutoff radius. The static dielectric constants for
TIP5P-E, TIP4P-Ew, SPC/E, and SSD/RF will be fixed at approximately
74, 52, 58, and 89 respectively. These values are generally lower than
the values reported in the literature; however, the relative
dielectric behavior scales as expected when comparing the models to
one another.

\section{Conclusions}\label{sec:PairwiseConclusions}

The above investigation of pairwise electrostatic summation techniques
shows that there are viable and computationally efficient alternatives
to the Ewald summation.  These methods are derived from the damped and
cutoff-neutralized Coulombic sum originally proposed by Wolf
\textit{et al.}\cite{Wolf99} In particular, the {\sc sf}
method, reformulated above as equations (\ref{eq:DSFPot}) and
(\ref{eq:DSFForces}), 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 the Ewald
summation.  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 the
Ewald summation. A detailed study of water simulations showed that
liquid properties calculated when using {\sc sf} will also be
equivalent to those obtained using the Ewald summation.

As in all purely-pairwise cutoff methods, these methods are expected
to scale approximately {\it linearly} with system size, and they are
easily parallelizable.  This should result in substantial reductions
in the computational cost of performing large simulations.

Aside from the computational cost benefit, these techniques have
applicability in situations where the use of the Ewald sum can prove
problematic.  Of greatest interest is their potential 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.

If a researcher is using Monte Carlo simulations of large chemical
systems containing point charges, most structural features will be
accurately captured using the undamped {\sc sf} method or the {\sc sp}
method with an electrostatic damping of 0.2~\AA$^{-1}$.  These methods
would also be appropriate for molecular dynamics simulations where the
data of interest is either structural or short-time dynamical
quantities.  For long-time dynamics and collective motions, the safest
pairwise method we have evaluated is the {\sc sf} method with an
electrostatic damping between 0.2 and 0.25~\AA$^{-1}$. It is also
important to note that the static dielectric constant in water
simulations is highly dependent on both $\alpha$ and
$R_\textrm{c}$. For consistent dielectric behavior, the damped {\sc
sf} method should use an $\alpha$ of 0.2175~\AA$^{-1}$ for an
$R_\textrm{c}$ of 12~\AA, and $\alpha$ should decrease by
0.025~\AA$^{-1}$ for every 1~\AA\ increase in cutoff radius.

We are not suggesting that there is any flaw with the Ewald sum; in
fact, it is the standard by which these simple pairwise sums have been
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.
Revision:	3023
Committed:	Tue Sep 26 03:07:59 2006 UTC (18 years, 7 months ago) by chrisfen
Content type:	application/x-tex
File size:	96315 byte(s)
Log Message:	Refinements. Just a little proof checking left...