[ViewVC] Diff of: group/trunk/electrostaticMethodsPaper/electrostaticMethods.tex

Comparing trunk/electrostaticMethodsPaper/electrostaticMethods.tex (file contents):
Revision 2575 by chrisfen, Sat Jan 28 22:42:46 2006 UTC vs.
Revision 2742 by chrisfen, Wed Apr 26 17:41:25 2006 UTC

+%\documentclass[prb,aps,twocolumn,tabularx]{revtex4}
-<
+\documentclass[12pt]{article}
->
+%\documentclass[aps,prb,preprint]{revtex4}
->
+\documentclass[11pt]{article}
+\usepackage{endfloat}
-<
+\usepackage{amsmath}
->
+\usepackage{amsmath,bm}
->
+\usepackage{amssymb}
+\usepackage{epsf}
+\usepackage{times}
-<
+\usepackage{mathptm}
->
+\usepackage{mathptmx}
+\usepackage{setspace}
+\usepackage{tabularx}
+\usepackage{graphicx}
-+
+\usepackage{booktabs}
-+
+\usepackage{bibentry}
-+
+\usepackage{mathrsfs}
+\usepackage[ref]{overcite}
+\pagestyle{plain}
+\pagenumbering{arabic}
+\topmargin -21pt \headsep 10pt
+\textheight 9.0in \textwidth 6.5in
+\brokenpenalty=10000
-<
+\renewcommand{\baselinestretch}{1.2}
->
+%\renewcommand{\baselinestretch}{1.2}
->
+\renewcommand{\baselinestretch}{2}
+\renewcommand\citemid{\ } % no comma in optional reference note
-+
+\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{2}} %doublespace captions
-+
+\let\Caption\caption
-+
+\renewcommand\caption[1]{%
-+
+        \Caption[#1]{}%
-+
+}
-+
+\begin{document}
-<
+\title{On the necessity of the Ewald Summation in molecular simulations.}
->
+\title{Is the Ewald summation still necessary? \\
->
+Pairwise alternatives to the accepted standard for
->
+long-range electrostatics in molecular simulations}
-<
+\author{Christopher J. Fennell and J. Daniel Gezelter \\
->
+\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail:
->
+gezelter@nd.edu} \\
+Department of Chemistry and Biochemistry\\
+University of Notre Dame\\
+Notre Dame, Indiana 46556}
+%\doublespacing
+\begin{abstract}
-+
+We investigate pairwise electrostatic interaction methods and show
-+
+that there are viable and computationally efficient $(\mathscr{O}(N))$
-+
+alternatives to the Ewald summation for typical modern molecular
-+
+simulations.  These methods are extended from the damped and
-+
+cutoff-neutralized Coulombic sum originally proposed by
-+
+[D. Wolf, P. Keblinski, S.~R. Phillpot, and J. Eggebrecht, {\it J. Chem. Phys.} {\bf 110}, 8255 (1999)] One of these, the damped shifted force method, shows
-+
+a remarkable ability to reproduce the energetic and dynamic
-+
+characteristics exhibited by simulations employing lattice summation
-+
+techniques.  Comparisons were performed with this and other pairwise
-+
+methods against the smooth particle mesh Ewald ({\sc spme}) summation
-+
+to see how well they reproduce the energetics and dynamics of a
-+
+variety of simulation types.
+\end{abstract}
-+
+\newpage
-+
+%\narrowtext
-<
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
->
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%                              BODY OF TEXT
-<
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
->
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+\section{Introduction}
-+
+In molecular simulations, proper accumulation of the electrostatic
-+
+interactions is essential and is one of the most
-+
+computationally-demanding tasks.  The common molecular mechanics force
-+
+fields represent atomic sites with full or partial charges protected
-+
+by Lennard-Jones (short range) interactions.  This means that nearly
-+
+every pair interaction involves a calculation of charge-charge forces.
-+
+Coupled with relatively long-ranged $r^{-1}$ decay, the monopole
-+
+interactions quickly become the most expensive part of molecular
-+
+simulations.  Historically, the electrostatic pair interaction would
-+
+not have decayed appreciably within the typical box lengths that could
-+
+be feasibly simulated.  In the larger systems that are more typical of
-+
+modern simulations, large cutoffs should be used to incorporate
-+
+electrostatics correctly.
-+
-+
+There have been many efforts to address the proper and practical
-+
+handling of electrostatic interactions, and these have resulted in a
-+
+variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are
-+
+typically classified as implicit methods (i.e., continuum dielectrics,
-+
+static dipolar fields),\cite{Born20,Grossfield00} explicit methods
-+
+(i.e., Ewald summations, interaction shifting or
-+
+truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e.,
-+
+reaction field type methods, fast multipole
-+
+methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are
-+
+often preferred because they physically incorporate solvent molecules
-+
+in the system of interest, but these methods are sometimes difficult
-+
+to utilize because of their high computational cost.\cite{Roux99} In
-+
+addition to the computational cost, there have been some questions
-+
+regarding possible artifacts caused by the inherent periodicity of the
-+
+explicit Ewald summation.\cite{Tobias01}
-+
-+
+In this paper, we focus on a new set of 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.
-+
-+
+\subsection{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} =
-+
+$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic
-+
+particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is
-+
+the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and
-+
+$j$, and $\phi$ is the solution to Poisson's equation
-+
+($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for
-+
+charge-charge interactions). In the case of monopole electrostatics,
-+
+eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for
-+
+non-neutral systems.
-+
-+
+The electrostatic summation problem was originally studied by Ewald
-+
+for the case of an infinite crystal.\cite{Ewald21}. The approach he
-+
+took was to convert this conditionally convergent sum into two
-+
+absolutely convergent summations: a short-ranged real-space summation
-+
+and a long-ranged reciprocal-space summation,
-+
+\begin{equation}
-+
+\begin{split}
-+
+V_\textrm{elec} = \frac{1}{2}& \sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime \frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)}{|\mathbf{r}_{ij}+\mathbf{n}|} \\ &+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} \exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) \cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ &- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 + \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3}\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2,
-+
+\end{split}
-+
+\label{eq:EwaldSum}
-+
+\end{equation}
-+
+where $\alpha$ is the damping or convergence parameter with units of
-+
+\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to
-+
+$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric
-+
+constant of the surrounding medium. The final two terms of
-+
+eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term
-+
+for interacting with a surrounding dielectric.\cite{Allen87} This
-+
+dipolar term was neglected in early applications in molecular
-+
+simulations,\cite{Brush66,Woodcock71} until it was introduced by de
-+
+Leeuw {\it et al.} to address situations where the unit cell has a
-+
+dipole moment which is magnified through replication of the periodic
-+
+images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the
-+
+system is said to be using conducting (or ``tin-foil'') boundary
-+
+conditions, $\epsilon_{\rm S} = \infty$. Figure
-+
+\ref{fig:ewaldTime} shows how the Ewald sum has been applied over
-+
+time.  Initially, due to the small 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}
-+
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width = \linewidth]{./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}
-+
-+
+The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm.  The
-+
+convergence parameter $(\alpha)$ plays an important role in balancing
-+
+the computational cost between the direct and reciprocal-space
-+
+portions of the summation.  The choice of this value allows one to
-+
+select whether the real-space or reciprocal space portion of the
-+
+summation is an $\mathscr{O}(N^2)$ calculation (with the other being
-+
+$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of
-+
+$\alpha$ and thoughtful algorithm development, this cost can be
-+
+reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route
-+
+taken to reduce the cost of the Ewald summation even further is to set
-+
+$\alpha$ such that the real-space interactions decay rapidly, allowing
-+
+for a short spherical cutoff. Then the reciprocal space summation is
-+
+optimized.  These optimizations usually involve utilization of the
-+
+fast Fourier transform (FFT),\cite{Hockney81} leading to the
-+
+particle-particle particle-mesh (P3M) and particle mesh Ewald (PME)
-+
+methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these
-+
+methods, the cost of the reciprocal-space portion of the Ewald
-+
+summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N
-+
+\log N)$.
-+
-+
+These developments and optimizations have made the use of the Ewald
-+
+summation routine in simulations with periodic boundary
-+
+conditions. However, in certain systems, such as vapor-liquid
-+
+interfaces and membranes, the intrinsic three-dimensional periodicity
-+
+can prove problematic.  The Ewald sum has been reformulated to handle
-+
+D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the
-+
+new methods are computationally expensive.\cite{Spohr97,Yeh99} More
-+
+recently, there have been several successful efforts toward reducing
-+
+the computational cost of 2D lattice summations, often enabling the
-+
+use of the mentioned
-+
+optimizations.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04}
-+
-+
+Several studies have recognized that the inherent periodicity in the
-+
+Ewald sum can also have an effect on three-dimensional
-+
+systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00}
-+
+Solvated proteins are essentially kept at high concentration due to
-+
+the periodicity of the electrostatic summation method.  In these
-+
+systems, the more compact folded states of a protein can be
-+
+artificially stabilized by the periodic replicas introduced by the
-+
+Ewald summation.\cite{Weber00} Thus, care must be taken when
-+
+considering the use of the Ewald summation where the assumed
-+
+periodicity would introduce spurious effects in the system dynamics.
-+
-+
+\subsection{The Wolf and Zahn Methods}
-+
+In a recent paper by Wolf \textit{et al.}, a procedure was outlined
-+
+for the accurate accumulation of electrostatic interactions in an
-+
+efficient pairwise fashion.  This procedure lacks the inherent
-+
+periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.}
-+
+observed that the electrostatic interaction is effectively
-+
+short-ranged in condensed phase systems and that neutralization of the
-+
+charge contained within the cutoff radius is crucial for potential
-+
+stability. They devised a pairwise summation method that ensures
-+
+charge neutrality and gives results similar to those obtained with the
-+
+Ewald summation.  The resulting shifted Coulomb potential includes
-+
+image-charges subtracted out through placement on the cutoff sphere
-+
+and a distance-dependent damping function (identical to that seen in
-+
+the real-space portion of the Ewald sum) to aid convergence
-+
+\begin{equation}
-+
+V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}.
-+
+\label{eq:WolfPot}
-+
+\end{equation}
-+
+Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted
-+
+potential.  However, neutralizing the charge contained within each
-+
+cutoff sphere requires the placement of a self-image charge on the
-+
+surface of the cutoff sphere.  This additional self-term in the total
-+
+potential enabled Wolf {\it et al.}  to obtain excellent estimates of
-+
+Madelung energies for many crystals.
-+
-+
+In order to use their charge-neutralized potential in molecular
-+
+dynamics simulations, Wolf \textit{et al.} suggested taking the
-+
+derivative of this potential prior to evaluation of the limit.  This
-+
+procedure gives an expression for the forces,
-+
+\begin{equation}
-+
+F_{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{\left[\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{r_{ij}}\right]-\left[\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right]\right\},
-+
+\label{eq:WolfForces}
-+
+\end{equation}
-+
+that incorporates both image charges and damping of the electrostatic
-+
+interaction.
-+
-+
+More recently, Zahn \textit{et al.} investigated these potential and
-+
+force expressions for use in simulations involving water.\cite{Zahn02}
-+
+In their work, they pointed out that the forces and derivative of
-+
+the potential are not commensurate.  Attempts to use both
-+
+eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead
-+
+to poor energy conservation.  They correctly observed that taking the
-+
+limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the
-+
+derivatives gives forces for a different potential energy function
-+
+than the one shown in eq. (\ref{eq:WolfPot}).
-+
-+
+Zahn \textit{et al.} introduced a modified form of this summation
-+
+method as a way to use the technique in Molecular Dynamics
-+
+simulations.  They proposed a new damped Coulomb potential,
-+
+\begin{equation}
-+
+V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\},
-+
+\label{eq:ZahnPot}
-+
+\end{equation}
-+
+and showed that this potential does fairly well at capturing the
-+
+structural and dynamic properties of water compared the same
-+
+properties obtained using the Ewald sum.
-+
-+
+\subsection{Simple Forms for Pairwise Electrostatics}
-+
-+
+The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et
-+
+al.} are constructed using two different (and separable) computational
-+
+tricks: \begin{enumerate}
-+
+\item shifting through the use of image charges, and
-+
+\item damping the electrostatic interaction.
-+
+\end{enumerate}  Wolf \textit{et al.} treated the
-+
+development of their summation method as a progressive application of
-+
+these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded
-+
+their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the
-+
+post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using
-+
+both techniques.  It is possible, however, to separate these
-+
+tricks and study their effects independently.
-+
-+
+Starting with the original observation that the effective range of the
-+
+electrostatic interaction in condensed phases is considerably less
-+
+than $r^{-1}$, either the cutoff sphere neutralization or the
-+
+distance-dependent damping technique could be used as a foundation for
-+
+a new pairwise summation method.  Wolf \textit{et al.} made the
-+
+observation that charge neutralization within the cutoff sphere plays
-+
+a significant role in energy convergence; therefore we will begin our
-+
+analysis with the various shifted forms that maintain this charge
-+
+neutralization.  We can evaluate the methods of Wolf
-+
+\textit{et al.}  and Zahn \textit{et al.} by considering the standard
-+
+shifted potential,
-+
+\begin{equation}
-+
+V_\textrm{SP}(r) =      \begin{cases}
-+
+v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r >
-+
+R_\textrm{c}
-+
+\end{cases},
-+
+\label{eq:shiftingPotForm}
-+
+\end{equation}
-+
+and shifted force,
-+
+\begin{equation}
-+
+V_\textrm{SF}(r) =      \begin{cases}
-+
+v(r)-v_\textrm{c}-\left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c})
-+
+&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c}
-+
+                                                \end{cases},
-+
+\label{eq:shiftingForm}
-+
+\end{equation}
-+
+functions where $v(r)$ is the unshifted form of the potential, and
-+
+$v_c$ is $v(R_\textrm{c})$.  The Shifted Force ({\sc sf}) form ensures
-+
+that both the potential and the forces goes to zero at the cutoff
-+
+radius, while the Shifted Potential ({\sc sp}) form only ensures the
-+
+potential is smooth at the cutoff radius
-+
+($R_\textrm{c}$).\cite{Allen87}
-+
-+
+The forces associated with the shifted potential are simply the forces
-+
+of the unshifted potential itself (when inside the cutoff sphere),
-+
+\begin{equation}
-+
+F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right),
-+
+\end{equation}
-+
+and are zero outside.  Inside the cutoff sphere, the forces associated
-+
+with the shifted force form can be written,
-+
+\begin{equation}
-+
+F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d
-+
+v(r)}{dr} \right)_{r=R_\textrm{c}}.
-+
+\end{equation}
-+
-+
+If the potential, $v(r)$, is taken to be the normal Coulomb potential,
-+
+\begin{equation}
-+
+v(r) = \frac{q_i q_j}{r},
-+
+\label{eq:Coulomb}
-+
+\end{equation}
-+
+then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et
-+
+al.}'s undamped prescription:
-+
+\begin{equation}
-+
+V_\textrm{SP}(r) =
-+
+q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad
-+
+r\leqslant R_\textrm{c},
-+
+\label{eq:SPPot}
-+
+\end{equation}
-+
+with associated forces,
-+
+\begin{equation}
-+
+F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}.
-+
+\label{eq:SPForces}
-+
+\end{equation}
-+
+These forces are identical to the forces of the standard Coulomb
-+
+interaction, and cutting these off at $R_c$ was addressed by Wolf
-+
+\textit{et al.} as undesirable.  They pointed out that the effect of
-+
+the image charges is neglected in the forces when this form is
-+
+used,\cite{Wolf99} thereby eliminating any benefit from the method in
-+
+molecular dynamics.  Additionally, there is a discontinuity in the
-+
+forces at the cutoff radius which results in energy drift during MD
-+
+simulations.
-+
-+
+The shifted force ({\sc sf}) form using the normal Coulomb potential
-+
+will give,
-+
+\begin{equation}
-+
+V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] \quad r\leqslant R_\textrm{c}.
-+
+\label{eq:SFPot}
-+
+\end{equation}
-+
+with associated forces,
-+
+\begin{equation}
-+
+F_\textrm{SF}(r) =  q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) \quad r\leqslant R_\textrm{c}.
-+
+\label{eq:SFForces}
-+
+\end{equation}
-+
+This formulation has the benefits that there are no discontinuities at
-+
+the cutoff radius, while the neutralizing image charges are present in
-+
+both the energy and force expressions.  It would be simple to add the
-+
+self-neutralizing term back when computing the total energy of the
-+
+system, thereby maintaining the agreement with the Madelung energies.
-+
+A side effect of this treatment is the alteration in the shape of the
-+
+potential that comes from the derivative term.  Thus, a degree of
-+
+clarity about agreement with the empirical potential is lost in order
-+
+to gain functionality in dynamics simulations.
-+
-+
+Wolf \textit{et al.} originally discussed the energetics of the
-+
+shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was
-+
+insufficient for accurate determination of the energy with reasonable
-+
+cutoff distances.  The calculated Madelung energies fluctuated around
-+
+the expected value as the cutoff radius was increased, but the
-+
+oscillations converged toward the correct value.\cite{Wolf99} A
-+
+damping function was incorporated to accelerate the convergence; and
-+
+though alternative forms for the damping function could be
-+
+used,\cite{Jones56,Heyes81} the complimentary error function was
-+
+chosen to mirror the effective screening used in the Ewald summation.
-+
+Incorporating this error function damping into the simple Coulomb
-+
+potential,
-+
+\begin{equation}
-+
+v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r},
-+
+\label{eq:dampCoulomb}
-+
+\end{equation}
-+
+the shifted potential (eq. (\ref{eq:SPPot})) becomes
-+
+\begin{equation}
-+
+V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}-\frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c},
-+
+\label{eq:DSPPot}
-+
+\end{equation}
-+
+with associated forces,
-+
+\begin{equation}
-+
+F_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \quad r\leqslant R_\textrm{c}.
-+
+\label{eq:DSPForces}
-+
+\end{equation}
-+
+Again, this damped shifted potential suffers from a
-+
+force-discontinuity at the cutoff radius, and the image charges play
-+
+no role in the forces.  To remedy these concerns, one may derive a
-+
+{\sc sf} variant by including the derivative term in
-+
+eq. (\ref{eq:shiftingForm}),
-+
+\begin{equation}
-+
+\begin{split}
-+
+V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r\right)}{r}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ &\left.+\left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right)\left(r-R_\mathrm{c}\right)\ \right] \quad r\leqslant R_\textrm{c}.
-+
+\label{eq:DSFPot}
-+
+\end{split}
-+
+\end{equation}
-+
+The derivative of the above potential will lead to the following forces,
-+
+\begin{equation}
-+
+\begin{split}
-+
+F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ &\left.-\left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right] \quad r\leqslant R_\textrm{c}.
-+
+\label{eq:DSFForces}
-+
+\end{split}
-+
+\end{equation}
-+
+If the damping parameter $(\alpha)$ is set to zero, the undamped case,
-+
+eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly
-+
+recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}).
-+
-+
+This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot}
-+
+derived by Zahn \textit{et al.}; however, there are two important
-+
+differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from
-+
+eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb})
-+
+with $r$ replaced by $R_\textrm{c}$.  This term is {\it not} present
-+
+in the Zahn potential, resulting in a potential discontinuity as
-+
+particles cross $R_\textrm{c}$.  Second, the sign of the derivative
-+
+portion is different.  The missing $v_\textrm{c}$ term would not
-+
+affect molecular dynamics simulations (although the computed energy
-+
+would be expected to have sudden jumps as particle distances crossed
-+
+$R_c$).  The sign problem is a potential source of errors, however.
-+
+In fact, it introduces a discontinuity in the forces at the cutoff,
-+
+because the force function is shifted in the wrong direction and
-+
+doesn't cross zero at $R_\textrm{c}$.
-+
-+
+Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an
-+
+electrostatic summation method in which the potential and forces are
-+
+continuous at the cutoff radius and which incorporates the damping
-+
+function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of
-+
+this paper, we will evaluate exactly how good these methods ({\sc sp},
-+
+{\sc sf}, damping) are at reproducing the correct electrostatic
-+
+summation performed by the Ewald sum.
-+
-+
+\subsection{Other alternatives}
-+
+In addition to the methods described above, we considered some other
-+
+techniques that are commonly used in molecular simulations.  The
-+
+simplest of these is group-based cutoffs.  Though of little use for
-+
+charged molecules, collecting atoms into neutral groups takes
-+
+advantage of the observation that the electrostatic interactions decay
-+
+faster than those for monopolar pairs.\cite{Steinbach94} When
-+
+considering these molecules as neutral groups, the relative
-+
+orientations of the molecules control the strength of the interactions
-+
+at the cutoff radius.  Consequently, as these molecular particles move
-+
+through $R_\textrm{c}$, the energy will drift upward due to the
-+
+anisotropy of the net molecular dipole interactions.\cite{Rahman71} To
-+
+maintain good energy conservation, both the potential and derivative
-+
+need to be smoothly switched to zero at $R_\textrm{c}$.\cite{Adams79}
-+
+This is accomplished using a standard switching function.  If a smooth
-+
+second derivative is desired, a fifth (or higher) order polynomial can
-+
+be used.\cite{Andrea83}
-+
-+
+Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$,
-+
+and to incorporate the effects of the surroundings, a method like
-+
+Reaction Field ({\sc rf}) can be used.  The original theory for {\sc
-+
+rf} was originally developed by Onsager,\cite{Onsager36} and it was
-+
+applied in simulations for the study of water by Barker and
-+
+Watts.\cite{Barker73} In modern simulation codes, {\sc rf} is simply
-+
+an extension of the group-based cutoff method where the net dipole
-+
+within the cutoff sphere polarizes an external dielectric, which
-+
+reacts back on the central dipole.  The same switching function
-+
+considerations for group-based cutoffs need to made for {\sc rf}, with
-+
+the additional pre-specification of a dielectric constant.
-+
+\section{Methods}
-+
+In classical molecular mechanics simulations, there are two primary
-+
+techniques utilized to obtain information about the system of
-+
+interest: Monte Carlo (MC) and Molecular Dynamics (MD).  Both of these
-+
+techniques utilize pairwise summations of interactions between
-+
+particle sites, but they use these summations in different ways.
-+
-+
+In MC, the potential energy difference between configurations dictates
-+
+the progression of MC sampling.  Going back to the origins of this
-+
+method, the acceptance criterion for the canonical ensemble laid out
-+
+by Metropolis \textit{et al.} states that a subsequent configuration
-+
+is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where
-+
+$\xi$ is a random number between 0 and 1.\cite{Metropolis53}
-+
+Maintaining the correct $\Delta E$ when using an alternate method for
-+
+handling the long-range electrostatics will ensure proper sampling
-+
+from the ensemble.
-+
-+
+In MD, the derivative of the potential governs how the system will
-+
+progress in time.  Consequently, the force and torque vectors on each
-+
+body in the system dictate how the system evolves.  If the magnitude
-+
+and direction of these vectors are similar when using alternate
-+
+electrostatic summation techniques, the dynamics in the short term
-+
+will be indistinguishable.  Because error in MD calculations is
-+
+cumulative, one should expect greater deviation at longer times,
-+
+although methods which have large differences in the force and torque
-+
+vectors will diverge from each other more rapidly.
-+
-+
+\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
-+
-+
+The pairwise summation techniques (outlined in section
-+
+\ref{sec:ESMethods}) were evaluated for use in MC simulations by
-+
+studying the energy differences between conformations.  We took the
-+
+{\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}.
-+
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width = \linewidth]{./dualLinear.pdf}
-+
+\caption{Example least squares regressions of the configuration energy
-+
+differences for SPC/E water systems. The upper plot shows a data set
-+
+with a poor correlation coefficient ($R^2$), while the lower plot
-+
+shows a data set with a good correlation coefficient.}
-+
+\label{fig:linearFit}
-+
+\end{figure}
-+
-+
+Each of the seven system types (detailed in section \ref{sec:RepSims})
-+
+were represented using 500 independent configurations.  Thus, each of
-+
+the alternative (non-Ewald) electrostatic summation methods was
-+
+evaluated using an accumulated 873,250 configurational energy
-+
+differences.
-+
-+
+Results and discussion for the individual analysis of each of the
-+
+system types appear in the supporting information,\cite{EPAPSdeposit}
-+
+while the cumulative results over all the investigated systems appears
-+
+below in section \ref{sec:EnergyResults}.
-+
-+
+\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods}
-+
+We evaluated the pairwise methods (outlined in section
-+
+\ref{sec:ESMethods}) for use in MD simulations by
-+
+comparing the force and torque vectors with those obtained using the
-+
+reference Ewald summation ({\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 1000~K).  The power spectrum
-+
+($I(\omega)$) was obtained via Fourier transform of the velocity
-+
+autocorrelation function, \begin{equation} I(\omega) =
-+
+\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt,
-+
+\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}
-+
+\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{itemize}
-+
+\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
-+
+and 0.3 \AA$^{-1}$,
-+
+\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2,
-+
+and 0.3 \AA$^{-1}$,
-+
+\item reaction field with an infinite dielectric constant, and
-+
+\item an unmodified cutoff.
-+
+\end{itemize}
-+
+Group-based cutoffs with a fifth-order polynomial switching function
-+
+were utilized for the reaction field simulations.  Additionally, we
-+
+investigated the use of these cutoffs with the {\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
-+
+^{-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,
-+
+.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\
-+
+respectively.
-+
+\section{Results and Discussion}
-+
+\subsection{Configuration Energy Differences}\label{sec:EnergyResults}
-+
+In order to evaluate the performance of the pairwise electrostatic
-+
+summation methods for Monte Carlo simulations, the energy differences
-+
+between configurations were compared to the values obtained when using
-+
+{\sc spme}.  The results for the subsequent regression analysis are shown in
-+
+figure \ref{fig:delE}.
-+
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width=5.5in]{./delEplot.pdf}
-+
+\caption{Statistical analysis of the quality of configurational energy
-+
+differences for a given electrostatic method compared with the
-+
+reference Ewald sum.  Results with a value equal to 1 (dashed line)
-+
+indicate $\Delta E$ values indistinguishable from those obtained using
-+
+{\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.  Interested
-+
+readers can consult the accompanying supporting information for a
-+
+comparison where all groups are neutral.\cite{EPAPSdeposit}
-+
-+
+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 homogenous
-+
+systems; although it does provide results that are an improvement over
-+
+those from an unmodified cutoff.
-+
-+
+\subsection{Magnitudes of the Force and Torque Vectors}
-+
-+
+Evaluation of pairwise methods for use in Molecular Dynamics
-+
+simulations requires consideration of effects on the forces and
-+
+torques.  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=5.5in]{./frcMagplot.pdf}
-+
+\caption{Statistical analysis of the quality of the force vector
-+
+magnitudes for a given electrostatic method compared with the
-+
+reference Ewald sum.  Results with a value equal to 1 (dashed line)
-+
+indicate force magnitude values indistinguishable from those obtained
-+
+using {\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=5.5in]{./trqMagplot.pdf}
-+
+\caption{Statistical analysis of the quality of the torque vector
-+
+magnitudes for a given electrostatic method compared with the
-+
+reference Ewald sum.  Results with a value equal to 1 (dashed line)
-+
+indicate torque magnitude values indistinguishable from those obtained
-+
+using {\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 be reproducing 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.
-+
-+
+\subsection{Directionality of the Force and Torque Vectors}
-+
-+
+It is clearly important that a new electrostatic method can 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=5.5in]{./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.  Overdamping 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 the
-+
+accompanying supporting information.\cite{EPAPSdeposit}
-+
-+
+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 using group based cutoffs along with the results seen in
-+
+figure \ref{fig:frcTrqAng}.  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}[htbp]
-+
+   \centering
-+
+   \caption{Statistical analysis of the angular
-+
+   distributions that the force (upper) and torque (lower) vectors
-+
+   from a given electrostatic method make with their counterparts
-+
+   obtained using the reference Ewald sum.  Calculations were
-+
+   performed both with (Y) and without (N) group based cutoffs and a
-+
+   switching function.  The $\alpha$ values have units of \AA$^{-1}$
-+
+   and the variance values have units of degrees$^2$.}
-+
-+
+   \begin{tabular}{@{} ccrrrrrrrr @{}}
-+
+      \\
-+
+      \toprule
-+
+      & & \multicolumn{4}{c}{Shifted Potential} & \multicolumn{4}{c}{Shifted Force} \\
-+
+      \cmidrule(lr){3-6}
-+
+      \cmidrule(l){7-10}
-+
+            $R_\textrm{c}$ & Groups & $\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$ & $\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$\\
-+
+      \midrule
-+
-+
+\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} \\
-+
+\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} \\
-+
+\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
-+
-+
+\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} \\
-+
+\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} \\
-+
+\AA  & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\
-+
+        & \textbf{Y} & \textbf{0.282} & \textbf{0.294} & \textbf{1.272} & \textbf{4.999} & \textbf{0.324} & \textbf{0.318} & \textbf{1.272} & \textbf{4.999} \\
-+
-+
+      \bottomrule
-+
+   \end{tabular}
-+
+   \label{tab:groupAngle}
-+
+\end{table}
-+
-+
+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 overdamp 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,
-+
+.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.
-+
-+
+\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}
-+
-+
+Zahn {\it et al.} investigated the structure and dynamics of water
-+
+using eqs. (\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]{./vCorrPlot.pdf}
-+
+\caption{Velocity autocorrelation functions of NaCl crystals at
-+
+K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc
-+
+sp} ($\alpha$ = 0.2). The inset is a magnification of the area around
-+
+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.
-+
-+
+\subsection{Collective Motion: Power Spectra of NaCl Crystals}
-+
-+
+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.
-+
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width = \linewidth]{./spectraSquare.pdf}
-+
+\caption{Power spectra obtained from the velocity auto-correlation
-+
+functions of NaCl crystals at 1000 K while using {\sc spme}, {\sc sf}
-+
+($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2).  The inset
-+
+shows the frequency region below 100 cm$^{-1}$ to highlight where the
-+
+spectra differ.}
-+
+\label{fig:methodPS}
-+
+\end{figure}
-+
-+
+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:dampInc}), at frequencies below 100 cm$^{-1}$, the
-+
+correlated motions are blue-shifted when using undamped or weakly
-+
+damped {\sc sf}.  When using moderate damping ($\alpha = 0.2$
-+
+\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give 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.
-+
-+
+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.
-+
-+
+\begin{figure}
-+
+\centering
-+
+\includegraphics[width = \linewidth]{./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.  Overdamping can
-+
+result in underestimates of frequencies of the long-wavelength
-+
+motions.}
-+
+\label{fig:dampInc}
-+
+\end{figure}
-+
+\section{Conclusions}
-+
+This 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 eqs. (\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 {\sc spme}.  Similarly
-+
+for the dynamic features, the undamped or moderately damped {\sc sf}
-+
+and moderately damped {\sc sp} methods produce force and torque vector
-+
+magnitude and directions very similar to the expected values.  These
-+
+results translate into long-time dynamic behavior equivalent to that
-+
+produced in simulations using {\sc spme}.
-+
-+
+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}$.
-+
-+
+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.
-+
+\section{Acknowledgments}
-+
+Support for this project was provided by the National Science
-+
+Foundation under grant CHE-0134881.  The authors would like to thank
-+
+Steve Corcelli and Ed Maginn for helpful discussions and comments.
-<
+\newpage
->
+\newpage
-<
+\bibliographystyle{achemso}
->
+\bibliographystyle{jcp2}
+\bibliography{electrostaticMethods}

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing trunk/electrostaticMethodsPaper/electrostaticMethods.tex (file contents): Revision 2575 by chrisfen, Sat Jan 28 22:42:46 2006 UTC vs. Revision 2742 by chrisfen, Wed Apr 26 17:41:25 2006 UTC

Diff Legend

Comparing trunk/electrostaticMethodsPaper/electrostaticMethods.tex (file contents):
Revision 2575 by chrisfen, Sat Jan 28 22:42:46 2006 UTC vs.
Revision 2742 by chrisfen, Wed Apr 26 17:41:25 2006 UTC