| 1 |
|
%\documentclass[prb,aps,twocolumn,tabularx]{revtex4} |
| 2 |
< |
\documentclass[12pt]{article} |
| 2 |
> |
%\documentclass[aps,prb,preprint]{revtex4} |
| 3 |
> |
\documentclass[11pt]{article} |
| 4 |
|
\usepackage{endfloat} |
| 5 |
|
\usepackage{amsmath} |
| 6 |
|
\usepackage{amssymb} |
| 7 |
|
\usepackage{epsf} |
| 8 |
|
\usepackage{times} |
| 9 |
< |
\usepackage{mathptm} |
| 9 |
> |
\usepackage{mathptmx} |
| 10 |
|
\usepackage{setspace} |
| 11 |
|
\usepackage{tabularx} |
| 12 |
|
\usepackage{graphicx} |
| 13 |
|
\usepackage{booktabs} |
| 14 |
|
\usepackage{bibentry} |
| 15 |
|
\usepackage{mathrsfs} |
| 15 |
– |
%\usepackage{berkeley} |
| 16 |
|
\usepackage[ref]{overcite} |
| 17 |
|
\pagestyle{plain} |
| 18 |
|
\pagenumbering{arabic} |
| 25 |
|
|
| 26 |
|
\begin{document} |
| 27 |
|
|
| 28 |
< |
\title{Is the Ewald Summation necessary? : Pairwise alternatives to the accepted standard for long-range electrostatics} |
| 28 |
> |
\title{Is the Ewald Summation necessary? Pairwise alternatives to the accepted standard for long-range electrostatics} |
| 29 |
|
|
| 30 |
< |
\author{Christopher J. Fennell and J. Daniel Gezelter \\ |
| 30 |
> |
\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
| 31 |
> |
gezelter@nd.edu} \\ |
| 32 |
|
Department of Chemistry and Biochemistry\\ |
| 33 |
|
University of Notre Dame\\ |
| 34 |
|
Notre Dame, Indiana 46556} |
| 36 |
|
\date{\today} |
| 37 |
|
|
| 38 |
|
\maketitle |
| 39 |
< |
%\doublespacing |
| 39 |
> |
\doublespacing |
| 40 |
> |
|
| 41 |
|
\nobibliography{} |
| 42 |
|
\begin{abstract} |
| 43 |
< |
A new method for accumulating electrostatic interactions was derived from the previous efforts described in \bibentry{Wolf99} and \bibentry{Zahn02} as a possible replacement for lattice sum methods in molecular simulations. Comparisons were performed with this and other pairwise electrostatic summation techniques against the smooth particle mesh Ewald (SPME) summation to see how well they reproduce the energetics and dynamics of a variety of simulation types. The newly derived Shifted-Force technique shows a remarkable ability to reproduce the behavior exhibited in simulations using SPME with an $\mathscr{O}(N)$ computational cost, equivalent to merely the real-space portion of the lattice summation. |
| 43 |
> |
A new method for accumulating electrostatic interactions was derived |
| 44 |
> |
from the previous efforts described in \bibentry{Wolf99} and |
| 45 |
> |
\bibentry{Zahn02} as a possible replacement for lattice sum methods in |
| 46 |
> |
molecular simulations. Comparisons were performed with this and other |
| 47 |
> |
pairwise electrostatic summation techniques against the smooth |
| 48 |
> |
particle mesh Ewald (SPME) summation to see how well they reproduce |
| 49 |
> |
the energetics and dynamics of a variety of simulation types. The |
| 50 |
> |
newly derived Shifted-Force technique shows a remarkable ability to |
| 51 |
> |
reproduce the behavior exhibited in simulations using SPME with an |
| 52 |
> |
$\mathscr{O}(N)$ computational cost, equivalent to merely the |
| 53 |
> |
real-space portion of the lattice summation. |
| 54 |
> |
|
| 55 |
|
\end{abstract} |
| 56 |
|
|
| 57 |
+ |
\newpage |
| 58 |
+ |
|
| 59 |
|
%\narrowtext |
| 60 |
|
|
| 61 |
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 64 |
|
|
| 65 |
|
\section{Introduction} |
| 66 |
|
|
| 67 |
< |
In molecular simulations, proper accumulation of the electrostatic interactions is considered one of the most essential and computationally demanding tasks. |
| 67 |
> |
In molecular simulations, proper accumulation of the electrostatic |
| 68 |
> |
interactions is considered one of the most essential and |
| 69 |
> |
computationally demanding tasks. |
| 70 |
|
|
| 71 |
|
\subsection{The Ewald Sum} |
| 72 |
|
blah blah blah Ewald Sum Important blah blah blah |
| 73 |
|
|
| 74 |
|
\begin{figure} |
| 75 |
|
\centering |
| 76 |
< |
\includegraphics[width = 3.25in]{./ewaldProgression.pdf} |
| 77 |
< |
\caption{How the application of the Ewald summation has changed with the increase in computer power. Initially, only small numbers of particles could be studied, and the Ewald sum acted to replicate the unit cell charge distribution out to infinity. Now, much larger systems of charges are investigated with fixed distance cutoffs. The calculated structure factor is used to sum out to great distance, and a surrounding dielectric term is included.} |
| 76 |
> |
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 77 |
> |
\caption{How the application of the Ewald summation has changed with |
| 78 |
> |
the increase in computer power. Initially, only small numbers of |
| 79 |
> |
particles could be studied, and the Ewald sum acted to replicate the |
| 80 |
> |
unit cell charge distribution out to convergence. Now, much larger |
| 81 |
> |
systems of charges are investigated with fixed distance cutoffs. The |
| 82 |
> |
calculated structure factor is used to sum out to great distance, and |
| 83 |
> |
a surrounding dielectric term is included.} |
| 84 |
|
\label{fig:ewaldTime} |
| 85 |
|
\end{figure} |
| 86 |
|
|
| 87 |
|
\subsection{The Wolf and Zahn Methods} |
| 88 |
< |
In a recent paper by Wolf \textit{et al.}, a procedure was outlined for accumulation of electrostatic interactions in a simple pairwise fashion.\cite{Wolf99} They took the observation that the effective electrostatic interaction is short-ranged in systems of charges and that charge neutralization within the cutoff spheres is crucial for potential stability. They devised a pairwise summation method that ensures charge neutrality and gives results similar to those obtained using the Ewald summation. The resulting shifted Coulomb potential (Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through placement on the cutoff sphere and a distance-dependent damping function (identical to that seen in the real-space portion of the Ewald sum) to aid energetic convergence |
| 88 |
> |
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
| 89 |
> |
for an accurate accumulation of electrostatic interactions in an |
| 90 |
> |
efficient pairwise fashion.\cite{Wolf99} Wolf \textit{et al.} observed |
| 91 |
> |
that the electrostatic interaction is effectively short-ranged in |
| 92 |
> |
condensed phase systems and that neutralization of the charge |
| 93 |
> |
contained within the cutoff radius is crucial for potential |
| 94 |
> |
stability. They devised a pairwise summation method that ensures |
| 95 |
> |
charge neutrality and gives results similar to those obtained with |
| 96 |
> |
the Ewald summation. The resulting shifted Coulomb potential |
| 97 |
> |
(Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through |
| 98 |
> |
placement on the cutoff sphere and a distance-dependent damping |
| 99 |
> |
function (identical to that seen in the real-space portion of the |
| 100 |
> |
Ewald sum) to aid convergence |
| 101 |
|
\begin{equation} |
| 102 |
|
V^{\textrm{Wolf}}(r_{ij})= \frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. |
| 103 |
|
\label{eq:WolfPot} |
| 104 |
|
\end{equation} |
| 105 |
< |
In order to use this potential in molecular dynamics simulations, Wolf \textit{et al.} suggested taking the derivative of this potential, followed by evaluation of the limit to give the following forces, |
| 105 |
> |
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
| 106 |
> |
potential. However, neutralizing the charge contained within each |
| 107 |
> |
cutoff sphere requires the placement of a self-image charge on the |
| 108 |
> |
surface of the cutoff sphere. This additional self-term in the total |
| 109 |
> |
potential enables Wolf {\it et al.} to obtain excellent estimates of |
| 110 |
> |
Madelung energies for many crystals. |
| 111 |
> |
|
| 112 |
> |
In order to use their charge-neutralized potential in molecular |
| 113 |
> |
dynamics simulations, Wolf \textit{et al.} suggested taking the |
| 114 |
> |
derivative of this potential prior to evaluation of the limit. This |
| 115 |
> |
procedure gives an expression for the forces, |
| 116 |
|
\begin{equation} |
| 117 |
< |
F^{\textrm{Wolf}}(r_{ij}) = q_iq_j\left[\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right)-\left(\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right]. |
| 117 |
> |
F^{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{-\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right]+\left[\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right]\right\}, |
| 118 |
|
\label{eq:WolfForces} |
| 119 |
|
\end{equation} |
| 120 |
< |
More recently, Zahn \textit{et al.} investigated this electrostatic summation method for use in simulations involving water.\cite{Zahn02} In their work, they point out that the method as proposed is problematic for use in Molecular Dynamics simulations, because the forces and derivative of the potential are not equivalent. This comes about from the procedure of taking the limit shown in equation \ref{eq:WolfPot} after calculating the derivatives.\cite{Wolf99} Zahn \textit{et al.} proposed a shifted force adaptation of this ``Wolf summation method" as a way to use this technique in Molecular Dynamics simulations. Taking the integral of the forces shown in equation \ref{eq:WolfForces}, they obtained a new shifted damped Coulomb potential |
| 120 |
> |
that incorporates both image charges and damping of the electrostatic |
| 121 |
> |
interaction. |
| 122 |
> |
|
| 123 |
> |
More recently, Zahn \textit{et al.} investigated these potential and |
| 124 |
> |
force expressions for use in simulations involving water.\cite{Zahn02} |
| 125 |
> |
In their work, they pointed out that the method that the forces and |
| 126 |
> |
derivative of the potential are not commensurate. Attempts to use |
| 127 |
> |
both Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will |
| 128 |
> |
lead to poor energy conservation. They correctly observed that taking |
| 129 |
> |
the limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating |
| 130 |
> |
the derivatives is mathematically invalid. |
| 131 |
> |
|
| 132 |
> |
Zahn \textit{et al.} proposed a modified form of this ``Wolf summation |
| 133 |
> |
method'' as a way to use this technique in Molecular Dynamics |
| 134 |
> |
simulations. Taking the integral of the forces shown in equation |
| 135 |
> |
\ref{eq:WolfForces}, they proposed a new damped Coulomb |
| 136 |
> |
potential, |
| 137 |
|
\begin{equation} |
| 138 |
|
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\}. |
| 139 |
|
\label{eq:ZahnPot} |
| 140 |
|
\end{equation} |
| 141 |
< |
They showed that this new potential does well in capturing the structural and dynamic properties present when using the Ewald sum with the models of water used in their simulations. |
| 141 |
> |
They showed that this potential does fairly well at capturing the |
| 142 |
> |
structural and dynamic properties of water compared the same |
| 143 |
> |
properties obtained using the Ewald sum. |
| 144 |
|
|
| 145 |
|
\subsection{Simple Forms for Pairwise Electrostatics} |
| 83 |
– |
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et al.} are constructed using two different (and separable) computational tricks: shifting through use of image charges and damping of the electrostatic interaction. |
| 146 |
|
|
| 147 |
< |
While implementing these methods for use in our own work, we discovered the potential presented in equation \ref{eq:ZahnPot} is still not entirely correct. The derivative of this equation leads to a sign error in the forces, resulting in erroneous dynamics. We can apply the standard shifted force potential, |
| 147 |
> |
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et |
| 148 |
> |
al.} are constructed using two different (and separable) computational |
| 149 |
> |
tricks: \begin{itemize} |
| 150 |
> |
\item shifting through the use of image charges, and |
| 151 |
> |
\item damping the electrostatic interaction. |
| 152 |
> |
\end{itemize} Wolf \textit{et al.} treated the |
| 153 |
> |
development of their summation method as a progressive application of |
| 154 |
> |
these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded |
| 155 |
> |
their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the |
| 156 |
> |
post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using |
| 157 |
> |
both techniques. It is possible, however, to separate these |
| 158 |
> |
tricks and study their effects independently. |
| 159 |
> |
|
| 160 |
> |
Starting with the original observation that the effective range of the |
| 161 |
> |
electrostatic interaction in condensed phases is considerably less |
| 162 |
> |
than $r^{-1}$, either the cutoff sphere neutralization or the |
| 163 |
> |
distance-dependent damping technique could be used as a foundation for |
| 164 |
> |
a new pairwise summation method. Wolf \textit{et al.} made the |
| 165 |
> |
observation that charge neutralization within the cutoff sphere plays |
| 166 |
> |
a significant role in energy convergence; therefore we will begin our |
| 167 |
> |
analysis with the various shifted forms that maintain this charge |
| 168 |
> |
neutralization. We can evaluate the methods of Wolf |
| 169 |
> |
\textit{et al.} and Zahn \textit{et al.} by considering the standard |
| 170 |
> |
shifted potential, |
| 171 |
|
\begin{equation} |
| 172 |
< |
V^\textrm{SF}(r_{ij}) = \begin{cases} v(r_{ij})-v_\textrm{c}-\left(\frac{\textrm{d}v(r_{ij})}{\textrm{d}r_{ij}}\right)_{r_{ij}=R_\textrm{c}}(r_{ij}-R_\textrm{c}) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 172 |
> |
v^\textrm{SP}(r) = \begin{cases} |
| 173 |
> |
v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > |
| 174 |
> |
R_\textrm{c} |
| 175 |
> |
\end{cases}, |
| 176 |
> |
\label{eq:shiftingPotForm} |
| 177 |
> |
\end{equation} |
| 178 |
> |
and shifted force, |
| 179 |
> |
\begin{equation} |
| 180 |
> |
v^\textrm{SF}(r) = \begin{cases} |
| 181 |
> |
v(r)-v_\textrm{c}-\left(\frac{\textrm{d}v(r)}{\textrm{d}r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) |
| 182 |
> |
&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} |
| 183 |
|
\end{cases}, |
| 184 |
+ |
\label{eq:shiftingForm} |
| 185 |
|
\end{equation} |
| 186 |
< |
where $v(r_{ij})$ is the unshifted form of the potential, and $v_c$ is $v(R_\textrm{c})$ and insures the potential goes to zero at the cutoff radius.\cite{Allen87} Using the simple damped Coulomb potential as the starting point, |
| 186 |
> |
functions where $v(r)$ is the unshifted form of the potential, and |
| 187 |
> |
$v_c$ is $v(R_\textrm{c})$. The Shifted Force ({\sc sf}) form ensures |
| 188 |
> |
that both the potential and the forces goes to zero at the cutoff |
| 189 |
> |
radius, while the Shifted Potential ({\sc sp}) form only ensures the |
| 190 |
> |
potential is smooth at the cutoff radius |
| 191 |
> |
($R_\textrm{c}$).\cite{Allen87} |
| 192 |
> |
|
| 193 |
> |
|
| 194 |
> |
|
| 195 |
> |
|
| 196 |
> |
If the derivative term is taken to be zero, we are left with the shifted Coulomb potential devised by Wolf \textit{et al.},\cite{Wolf99} |
| 197 |
|
\begin{equation} |
| 198 |
< |
v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}, |
| 199 |
< |
\label{eq:dampCoulomb} |
| 198 |
> |
V^\textrm{WSP}(r_{ij}) = \begin{cases} q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 199 |
> |
\end{cases}. |
| 200 |
> |
\label{eq:WolfSP} |
| 201 |
|
\end{equation} |
| 202 |
< |
the resulting shifted force potential is |
| 202 |
> |
The forces associated with this potential are obtained by taking the derivative, resulting in the following, |
| 203 |
|
\begin{equation} |
| 204 |
< |
V^\mathrm{SF}\left(r_{ij}\right)=q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\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_{ij}-R_\mathrm{c}\right)\right\}. |
| 204 |
> |
F^\textrm{WSP}(r_{ij}) = \begin{cases} q_iq_j\left(-\frac{1}{r_{ij}^2}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 205 |
> |
\end{cases}. |
| 206 |
> |
\label{eq:FWolfSP} |
| 207 |
> |
\end{equation} |
| 208 |
> |
These forces are identical to the forces of the standard electrostatic interaction, and this was addressed by Wolf \textit{et al.} as undesirable. They pointed out that the effect of the image charges is neglected in the forces when they would expect there to be some pressure exerted due to their presence.\cite{Wolf99} As a consequence the forces, though mathematically valid, may not be physically correct due to this missing component. Additionally, there is a discontinuity in the forces. This can be remedied with the use of a switching function to zero the potential and forces smoothly as particles near $R_\textrm{c}$. |
| 209 |
> |
|
| 210 |
> |
If the derivative term in equation \ref{eq:shiftingForm} is evaluated, we obtain an hitherto undiscussed shifted force Coulomb potential, |
| 211 |
> |
\begin{equation} |
| 212 |
> |
V^\textrm{SF}(r_{ij}) = \begin{cases} q_iq_j\left\{\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}+\left[\frac{1}{R_\textrm{c}^2}\right](r_{ij}-R_\textrm{c})\right\} &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 213 |
> |
\end{cases}. |
| 214 |
|
\label{eq:SFPot} |
| 215 |
|
\end{equation} |
| 216 |
< |
Equation \ref{eq:SFPot} 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 is simply equation \ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This term is not present in equation \ref{eq:ZahnPot}, resulting in a discontinuity in the potential as particles cross $R_\textrm{c}$. Second, the sign of the derivative portion is different. The constant $v_\textrm{c}$ term is not going to have a presence in the forces after performing the derivative, but the negative sign does effect the derivative. In fact, it introduces a discontinuity in the forces at the cutoff, because the force function is shifted in the wrong direction and doesn't cross zero at $R_\textrm{c}$. Thus, these alterations make for an electrostatic summation method that is continuous in both the potential and forces and incorporates the pairwise sum considerations stressed by Wolf \textit{et al.}\cite{Wolf99} |
| 216 |
> |
Taking the derivative of this shifted force potential gives the following forces, |
| 217 |
> |
\begin{equation} |
| 218 |
> |
F^\textrm{SF}(r_{ij}) = \begin{cases} q_iq_j\left(-\frac{1}{r_{ij}^2}+\frac{1}{R_\textrm{c}^2}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 219 |
> |
\end{cases}. |
| 220 |
> |
\label{eq:SFForces} |
| 221 |
> |
\end{equation} |
| 222 |
> |
Using this formulation rather than the simple shifted potential (Eq. \ref{eq:WolfSP}) means that there are no discontinuities in the forces in addition to the potential. This form also has the benefit that the image charges have a force presence, addressing the concerns about a missing physical component. One side effect of this treatment is a slight alteration in the shape of the potential that comes about from the derivative term. Thus, a degree of clarity about the original formulation of the potential is lost in order to gain functionality in dynamics simulations. |
| 223 |
|
|
| 224 |
< |
It is important to note that shifted force techniques have a drawback in that they alter the shape of the original potential. We thereby lose a degree of clarity about the original formulation of the potential in order to gain functionality in dynamics simulations. An alternative direction would be use the derivatives of the original potential for the forces. This was addressed by Wolf \textit{et al.} as undesirable, because the effect of the image charges is neglected in the forces when they would expect there to be some pressure exerted due to their presence.\cite{Wolf99} As a consequence the forces, though mathematically valid, may not be physically correct due to this missing component. In Monte Carlo simulations, this argument is moot, because forces are not evaluated. We decided to consider both the Shifted-Force technique described above and this Shifted-Potential technique to determine their usability in the evaluation of both energetic and dynamic results in simulations with electrostatics. |
| 224 |
> |
Wolf \textit{et al.} originally discussed the energetics of the shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that it was still insufficient for accurate determination of the energy. The energy would fluctuate around the expected value with increasing cutoff radius, but the oscillations appeared to be converging toward the correct value.\cite{Wolf99} A damping function was incorporated to accelerate convergence; and though alternative functional forms could be used,\cite{Jones56,Heyes81} the complimentary error function was chosen to draw parallels to the Ewald summation. Incorporating damping into the simple Coulomb potential, |
| 225 |
> |
\begin{equation} |
| 226 |
> |
v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}, |
| 227 |
> |
\label{eq:dampCoulomb} |
| 228 |
> |
\end{equation} |
| 229 |
> |
the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived \textit{via} equation \ref{eq:shiftingForm}, |
| 230 |
> |
\begin{equation} |
| 231 |
> |
V^{\textrm{DSP}}(r_{ij}) = \begin{cases} q_iq_j\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right] &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 232 |
> |
\end{cases}. |
| 233 |
> |
\label{eq:DSPPot} |
| 234 |
> |
\end{equation} |
| 235 |
> |
The derivative of this Shifted-Potential can be taken to obtain forces for use in MD, |
| 236 |
> |
\begin{equation} |
| 237 |
> |
F^{\textrm{DSP}}(r_{ij}) = \begin{cases} q_iq_j\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right] &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 238 |
> |
\end{cases}. |
| 239 |
> |
\label{eq:DSPForces} |
| 240 |
> |
\end{equation} |
| 241 |
> |
Again, this Shifted-Potential suffers from a discontinuity in the forces, and a lack of an image-charge component in the forces. To remedy these concerns, a Shifted-Force variant is obtained by inclusion of the derivative term in equation \ref{eq:shiftingForm} to give, |
| 242 |
> |
\begin{equation} |
| 243 |
> |
V^\mathrm{DSF}(r_{ij}) = \begin{cases} q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\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_{ij}-R_\mathrm{c}\right)\right\} &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 244 |
> |
\end{cases}. |
| 245 |
> |
\label{eq:DSFPot} |
| 246 |
> |
\end{equation} |
| 247 |
> |
The derivative of the above potential gives the following forces, |
| 248 |
> |
\begin{equation} |
| 249 |
> |
F^\mathrm{DSF}(r_{ij}) = \begin{cases} q_iq_j\left\{-\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right]+\left[\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2R_{\textrm{c}}^2)}}{R_{\textrm{c}}}\right]\right\} &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 250 |
> |
\end{cases}. |
| 251 |
> |
\label{eq:DSFForces} |
| 252 |
> |
\end{equation} |
| 253 |
|
|
| 254 |
+ |
This new Shifted-Force potential is similar to equation \ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from equation \ref{eq:shiftingForm} is equal to equation \ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This term is not present in the Zahn potential, resulting in a discontinuity as particles cross $R_\textrm{c}$. Second, the sign of the derivative portion is different. The constant $v_\textrm{c}$ term is not going to have a presence in the forces after performing the derivative, but the negative sign does effect the derivative. In fact, it introduces a discontinuity in the forces at the cutoff, because the force function is shifted in the wrong direction and doesn't cross zero at $R_\textrm{c}$. Thus, these alterations make for an electrostatic summation method that is continuous in both the potential and forces and incorporates the pairwise sum considerations stressed by Wolf \textit{et al.}\cite{Wolf99} |
| 255 |
+ |
|
| 256 |
|
\section{Methods} |
| 257 |
|
|
| 258 |
|
\subsection{What Qualities are Important?}\label{sec:Qualities} |
| 267 |
|
|
| 268 |
|
\begin{figure} |
| 269 |
|
\centering |
| 270 |
< |
\includegraphics[width=3.25in]{./linearFit.pdf} |
| 271 |
< |
\caption{Example least squares regression of the $\Delta E$ between configurations for the SF method against SPME in the pure water system. } |
| 270 |
> |
\includegraphics[width = \linewidth]{./dualLinear.pdf} |
| 271 |
> |
\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.} |
| 272 |
|
\label{fig:linearFit} |
| 273 |
|
\end{figure} |
| 274 |
|
|
| 275 |
< |
Each system type (detailed in section \ref{sec:RepSims}) studied consisted of 500 independent configurations, each equilibrated from higher temperature trajectories. Thus, 124,750 $\Delta E$ data points are used in a regression of a single system type. Results and discussion for the individual analysis of each of the system types appear in the supporting information, while the cumulative results over all the investigated systems appears below in section \ref{sec:EnergyResults}. |
| 275 |
> |
Each system type (detailed in section \ref{sec:RepSims}) studied |
| 276 |
> |
consisted of 500 independent configurations, each equilibrated from |
| 277 |
> |
higher temperature trajectories. Thus, 124,750 $\Delta E$ data points |
| 278 |
> |
are used in a regression of a single system type. Results and |
| 279 |
> |
discussion for the individual analysis of each of the system types |
| 280 |
> |
appear in the supporting information, while the cumulative results |
| 281 |
> |
over all the investigated systems appears below in section |
| 282 |
> |
\ref{sec:EnergyResults}. |
| 283 |
|
|
| 284 |
|
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
| 285 |
|
Evaluation of the pairwise methods (outlined in section \ref{sec:ESMethods}) for use in MD simulations was performed through comparison of the force and torque vectors obtained with those from SPME. Both the magnitude and the direction of these vectors on each of the bodies in the system were analyzed. For the magnitude of these vectors, linear least squares regression analysis can be performed as described previously for comparing $\Delta E$ values. Instead of a single value between two system configurations, there is a value for each particle in each configuration. For a system of 1000 water molecules and 40 ions, there are 1040 force vectors and 1000 torque vectors. With 500 configurations, this results in 520,000 force and 500,000 torque vector comparisons samples for each system type. |
| 292 |
|
|
| 293 |
|
\begin{figure} |
| 294 |
|
\centering |
| 295 |
< |
\includegraphics[width=3.25in]{./gaussFit.pdf} |
| 295 |
> |
\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 296 |
|
\caption{Sample fit of the angular distribution of the force vectors over all of the studied systems. Gaussian fits were used to obtain values for the variance in force and torque vectors used in the following figure.} |
| 297 |
|
\label{fig:gaussian} |
| 298 |
|
\end{figure} |
| 329 |
|
|
| 330 |
|
\begin{figure} |
| 331 |
|
\centering |
| 332 |
< |
\includegraphics[width=3.25in]{./slice.pdf} |
| 332 |
> |
\includegraphics[width = \linewidth]{./slice.pdf} |
| 333 |
|
\caption{A slice from the center of a water box used in a charge void simulation. The darkened region represents the boundary sphere within which the water molecules were converted to argon atoms.} |
| 334 |
|
\label{fig:argonSlice} |
| 335 |
|
\end{figure} |
| 346 |
|
|
| 347 |
|
\begin{figure} |
| 348 |
|
\centering |
| 349 |
< |
\includegraphics[width=3.25in]{./delEplot.pdf} |
| 349 |
> |
\includegraphics[width=5.5in]{./delEplot.pdf} |
| 350 |
|
\caption{Statistical analysis of the quality of configurational energy differences for a given electrostatic method compared with the reference Ewald sum. Results with a value equal to 1 (dashed line) indicate $\Delta E$ values indistinguishable from those obtained using SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 351 |
|
\label{fig:delE} |
| 352 |
|
\end{figure} |
| 361 |
|
|
| 362 |
|
\begin{figure} |
| 363 |
|
\centering |
| 364 |
< |
\includegraphics[width=3.25in]{./frcMagplot.pdf} |
| 364 |
> |
\includegraphics[width=5.5in]{./frcMagplot.pdf} |
| 365 |
|
\caption{Statistical analysis of the quality of the force vector magnitudes for a given electrostatic method compared with the reference Ewald sum. Results with a value equal to 1 (dashed line) indicate force magnitude values indistinguishable from those obtained using SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 366 |
|
\label{fig:frcMag} |
| 367 |
|
\end{figure} |
| 370 |
|
|
| 371 |
|
\begin{figure} |
| 372 |
|
\centering |
| 373 |
< |
\includegraphics[width=3.25in]{./trqMagplot.pdf} |
| 373 |
> |
\includegraphics[width=5.5in]{./trqMagplot.pdf} |
| 374 |
|
\caption{Statistical analysis of the quality of the torque vector magnitudes for a given electrostatic method compared with the reference Ewald sum. Results with a value equal to 1 (dashed line) indicate torque magnitude values indistinguishable from those obtained using SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 375 |
|
\label{fig:trqMag} |
| 376 |
|
\end{figure} |
| 383 |
|
|
| 384 |
|
\begin{figure} |
| 385 |
|
\centering |
| 386 |
< |
\includegraphics[width=3.25in]{./frcTrqAngplot.pdf} |
| 386 |
> |
\includegraphics[width=5.5in]{./frcTrqAngplot.pdf} |
| 387 |
|
\caption{Statistical analysis of the quality of the Gaussian fit of the force and torque vector angular distributions for a given electrostatic method compared with the reference Ewald sum. Results with a variance ($\sigma^2$) equal to zero (dashed line) indicate force and torque directions indistinguishable from those obtained using SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 388 |
|
\label{fig:frcTrqAng} |
| 389 |
|
\end{figure} |
| 433 |
|
|
| 434 |
|
\begin{figure} |
| 435 |
|
\centering |
| 436 |
< |
\includegraphics[width = 3.25in]{./nModeFTPlotDot.pdf} |
| 437 |
< |
\caption{Power spectra obtained from the velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, Shifted-Force ($\alpha$ = 0, 0.1, \& 0.2), and Shifted-Potential ($\alpha$ = 0.2). Apodization of the correlation functions via a cubic switching function between 40 and 50 ps was used to clear up the spectral noise resulting from data truncation, and had no noticeable effect on peak location or magnitude. The inset shows the frequency region below 100 cm$^{-1}$ to highlight where the spectra begin to differentiate.} |
| 436 |
> |
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 437 |
> |
\caption{Power spectra obtained from the velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, Shifted-Force ($\alpha$ = 0, 0.1, \& 0.2), and Shifted-Potential ($\alpha$ = 0.2). Apodization of the correlation functions via a cubic switching function between 40 and 50 ps was used to clear up the spectral noise resulting from data truncation, and had no noticeable effect on peak location or magnitude. The inset shows the frequency region below 100 cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 438 |
|
\label{fig:methodPS} |
| 439 |
|
\end{figure} |
| 440 |
|
|
| 445 |
|
where $S(r)$ is a switching function that smoothly zeroes the potential at the cutoff radius. Figure \ref{fig:dampInc} shows how the low frequency motions are dependent on the damping used in the direct electrostatic sum. As the damping increases, the peaks drop to lower frequencies. Incidentally, use of an $\alpha$ of 0.25 \AA$^{-1}$ on a simple electrostatic summation results in low frequency correlated dynamics equivalent to a simulation using SPME. When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks shift to higher frequency in exponential fashion. Though not shown, the spectrum for the simple undamped electrostatic potential is blue-shifted such that the lowest frequency peak resides near 325 cm$^{-1}$. In light of these results, the undamped Shifted-Force method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite respectable; however, it appears as though moderate damping is required for accurate reproduction of crystal dynamics. |
| 446 |
|
\begin{figure} |
| 447 |
|
\centering |
| 448 |
< |
\includegraphics[width = 3.25in]{./alphaCompare.pdf} |
| 449 |
< |
\caption{Normal modes for an NaCl crystal at 1000 K when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$)ranging from 0.15 to 0.3 \AA$^{-1}$. As $\alpha$ increases, the peaks are red-shifted toward and eventually beyond the values given by SPME. The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.} |
| 448 |
> |
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 449 |
> |
\caption{Upper: Zoomed inset from figure \ref{fig:methodPS}. As the damping value for the Shifted-Force potential increases, the low-frequency peaks red-shift. Lower: Low-frequency correlated motions for NaCl crystals at 1000 K when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$. As $\alpha$ increases, the peaks are red-shifted toward and eventually beyond the values given by SPME. The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.} |
| 450 |
|
\label{fig:dampInc} |
| 451 |
|
\end{figure} |
| 452 |
|
|
| 459 |
|
We are not suggesting any flaw with the Ewald sum; in fact, it is the standard by which these simple pairwise sums are judged. However, these results do suggest that in the typical simulations performed today, the Ewald summation may no longer be required to obtain the level of accuracy most researcher have come to expect |
| 460 |
|
|
| 461 |
|
\section{Acknowledgments} |
| 303 |
– |
|
| 462 |
|
\newpage |
| 463 |
|
|
| 464 |
< |
\bibliographystyle{achemso} |
| 464 |
> |
\bibliographystyle{jcp2} |
| 465 |
|
\bibliography{electrostaticMethods} |
| 466 |
|
|
| 467 |
|
|
