| 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} |
| 5 |
> |
\usepackage{amsmath,bm} |
| 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{berkeley} |
| 13 |
> |
\usepackage{booktabs} |
| 14 |
> |
\usepackage{bibentry} |
| 15 |
> |
\usepackage{mathrsfs} |
| 16 |
|
\usepackage[ref]{overcite} |
| 17 |
|
\pagestyle{plain} |
| 18 |
|
\pagenumbering{arabic} |
| 25 |
|
|
| 26 |
|
\begin{document} |
| 27 |
|
|
| 28 |
< |
\title{On the necessity of the Ewald Summation in molecular simulations.} |
| 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 |
| 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 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 61 |
> |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 62 |
|
% BODY OF TEXT |
| 63 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 63 |
> |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 64 |
|
|
| 65 |
|
\section{Introduction} |
| 66 |
|
|
| 67 |
< |
In this paper, a variety of simulation situations were analyzed to determine the relative effectiveness of the adapted Wolf spherical truncation schemes at reproducing the results obtained using a smooth particle mesh Ewald (SPME) summation technique. In addition to the Shifted-Potential (SP) and Shifted-Force (SF) adapted Wolf methods, both reaction field and uncorrected cutoff methods were included for comparison purposes. The general usability of these methods in both Monte Carlo (MC) and Molecular Dynamics (MD) calculations was assessed through statistical analysis over the combined results from all of the following studied systems: |
| 68 |
< |
\begin{list}{-}{} |
| 69 |
< |
\item Liquid Water |
| 70 |
< |
\item Crystalline Water (Ice I$_\textrm{c}$) |
| 71 |
< |
\item 1 M Solution of NaCl in Water |
| 72 |
< |
\item 10 M Solution of NaCl in Water |
| 73 |
< |
\item 6 \AA\ Radius Sphere of Argon in Water |
| 74 |
< |
\item NaCl Crystal |
| 75 |
< |
\item NaCl Melt |
| 76 |
< |
\end{list} |
| 77 |
< |
Additional discussion on the results from the individual systems was also performed to identify limitations of the considered methods in specific systems. |
| 67 |
> |
In molecular simulations, proper accumulation of the electrostatic |
| 68 |
> |
interactions is considered one of the most essential and |
| 69 |
> |
computationally demanding tasks. The common molecular mechanics force |
| 70 |
> |
fields are founded on representation of the atomic sites centered on |
| 71 |
> |
full or partial charges shielded by Lennard-Jones type interactions. |
| 72 |
> |
This means that nearly every pair interaction involves an |
| 73 |
> |
charge-charge calculation. Coupled with $r^{-1}$ decay, the monopole |
| 74 |
> |
interactions quickly become a burden for molecular systems of all |
| 75 |
> |
sizes. For example, in small systems, the electrostatic pair |
| 76 |
> |
interaction may not have decayed appreciably within the box length |
| 77 |
> |
leading to an effect excluded from the pair interactions within a unit |
| 78 |
> |
box. In large systems, excessively large cutoffs need to be used to |
| 79 |
> |
accurately incorporate their effect, and since the computational cost |
| 80 |
> |
increases proportionally with the cutoff sphere, it quickly becomes an |
| 81 |
> |
impractical task to perform these calculations. |
| 82 |
|
|
| 83 |
< |
\section{Methods} |
| 83 |
> |
\subsection{The Ewald Sum} |
| 84 |
> |
The complete accumulation electrostatic interactions in a system with periodic boundary conditions (PBC) requires the consideration of the effect of all charges within a simulation box, as well as those in the periodic replicas, |
| 85 |
> |
\begin{equation} |
| 86 |
> |
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], |
| 87 |
> |
\label{eq:PBCSum} |
| 88 |
> |
\end{equation} |
| 89 |
> |
where the sum over $\mathbf{n}$ is a sum over all periodic box replicas |
| 90 |
> |
with integer coordinates $\mathbf{n} = (l,m,n)$, and the prime indicates |
| 91 |
> |
$i = j$ are neglected for $\mathbf{n} = 0$.\cite{deLeeuw80} Within the |
| 92 |
> |
sum, $N$ is the number of electrostatic particles, $\mathbf{r}_{ij}$ is |
| 93 |
> |
$\mathbf{r}_j - \mathbf{r}_i$, $L$ is the cell length, $\bm{\Omega}_{i,j}$ are |
| 94 |
> |
the Euler angles for $i$ and $j$, and $\phi$ is Poisson's equation |
| 95 |
> |
($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge |
| 96 |
> |
interactions). In the case of monopole electrostatics, |
| 97 |
> |
eq. (\ref{eq:PBCSum}) is conditionally convergent and is discontiuous |
| 98 |
> |
for non-neutral systems. |
| 99 |
|
|
| 100 |
< |
In each of the simulated systems, 500 distinct configurations were generated, and the electrostatic summation methods were compared via sequential application on each of these fixed configurations. The methods compared include SPME, the aforementioned SP and SF methods - both with damping parameters ($\alpha$) of 0, 0.1, 0.2, and 0.3 \AA$^{-1}$, reaction field with an infinite dielectric constant, and an unmodified cutoff. Group-based cutoffs with a fifth-order polynomial switching function were necessary for the reaction field simulations and were utilized in the SP, SF, and pure cutoff methods for comparison to the standard lack of group-based cutoffs with a hard truncation. |
| 100 |
> |
This electrostatic summation problem was originally studied by Ewald |
| 101 |
> |
for the case of an infinite crystal.\cite{Ewald21}. The approach he |
| 102 |
> |
took was to convert this conditionally convergent sum into two |
| 103 |
> |
absolutely convergent summations: a short-ranged real-space summation |
| 104 |
> |
and a long-ranged reciprocal-space summation, |
| 105 |
> |
\begin{equation} |
| 106 |
> |
\begin{split} |
| 107 |
> |
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}{3L^3}\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2, |
| 108 |
> |
\end{split} |
| 109 |
> |
\label{eq:EwaldSum} |
| 110 |
> |
\end{equation} |
| 111 |
> |
where $\alpha$ is a damping parameter, or separation constant, with |
| 112 |
> |
units of \AA$^{-1}$, and $\mathbf{k}$ are the reciprocal vectors and |
| 113 |
> |
equal $2\pi\mathbf{n}/L^2$. The final two terms of |
| 114 |
> |
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 115 |
> |
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 116 |
> |
dipolar term was neglected in early applications in molecular |
| 117 |
> |
simulations,\cite{Brush66,Woodcock71} until it was introduced by de |
| 118 |
> |
Leeuw {\it et al.} to address situations where the unit cell has a |
| 119 |
> |
dipole moment and this dipole moment gets magnified through |
| 120 |
> |
replication of the periodic images.\cite{deLeeuw80} This term is zero |
| 121 |
> |
for systems where $\epsilon_{\rm S} = \infty$. Figure |
| 122 |
> |
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 123 |
> |
time. Initially, due to the small size of systems, the entire |
| 124 |
> |
simulation box was replicated to convergence. Currently, we balance a |
| 125 |
> |
spherical real-space cutoff with the reciprocal sum and consider the |
| 126 |
> |
surrounding dielectric. |
| 127 |
> |
\begin{figure} |
| 128 |
> |
\centering |
| 129 |
> |
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 130 |
> |
\caption{How the application of the Ewald summation has changed with |
| 131 |
> |
the increase in computer power. Initially, only small numbers of |
| 132 |
> |
particles could be studied, and the Ewald sum acted to replicate the |
| 133 |
> |
unit cell charge distribution out to convergence. Now, much larger |
| 134 |
> |
systems of charges are investigated with fixed distance cutoffs. The |
| 135 |
> |
calculated structure factor is used to sum out to great distance, and |
| 136 |
> |
a surrounding dielectric term is included.} |
| 137 |
> |
\label{fig:ewaldTime} |
| 138 |
> |
\end{figure} |
| 139 |
|
|
| 140 |
< |
Generation of the system configurations was dependent on the system type. For the solid and liquid water configurations, configuration snapshots were taken at regular intervals from higher temperature 1000 SPC/E water molecule trajectories and individually equilibrated. The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- ions and were selected and equilibrated in the same fashion as the water systems. For the 1 and 10 M NaCl solutions, 4 and 40 ions, respectively, were first solvated in a 1000 water molecule boxes. Ion and water positions were then randomly swapped, and the resulting configurations were again individually equilibrated. Finally, for the Argon/Water "charge void" systems, the identities of all the SPC/E waters within 6 \AA\ of the center of the equilibrated water configurations were converted to argon (Fig. \ref{argonSlice}). |
| 140 |
> |
The Ewald summation in the straight-forward form is an |
| 141 |
> |
$\mathscr{O}(N^2)$ algorithm. The separation constant $(\alpha)$ |
| 142 |
> |
plays an important role in the computational cost balance between the |
| 143 |
> |
direct and reciprocal-space portions of the summation. The choice of |
| 144 |
> |
the magnitude of this value allows one to whether the real-space or |
| 145 |
> |
reciprocal space portion of the summation is an $\mathscr{O}(N^2)$ |
| 146 |
> |
calcualtion, with the other being $\mathscr{O}(N)$.\cite{Sagui99} With |
| 147 |
> |
appropriate choice of $\alpha$ and thoughtful algorithm development, |
| 148 |
> |
this cost can be brought down to |
| 149 |
> |
$\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to |
| 150 |
> |
accelerate the Ewald summation is to se |
| 151 |
|
|
| 152 |
+ |
\subsection{The Wolf and Zahn Methods} |
| 153 |
+ |
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
| 154 |
+ |
for the accurate accumulation of electrostatic interactions in an |
| 155 |
+ |
efficient pairwise fashion.\cite{Wolf99} Wolf \textit{et al.} observed |
| 156 |
+ |
that the electrostatic interaction is effectively short-ranged in |
| 157 |
+ |
condensed phase systems and that neutralization of the charge |
| 158 |
+ |
contained within the cutoff radius is crucial for potential |
| 159 |
+ |
stability. They devised a pairwise summation method that ensures |
| 160 |
+ |
charge neutrality and gives results similar to those obtained with |
| 161 |
+ |
the Ewald summation. The resulting shifted Coulomb potential |
| 162 |
+ |
(Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through |
| 163 |
+ |
placement on the cutoff sphere and a distance-dependent damping |
| 164 |
+ |
function (identical to that seen in the real-space portion of the |
| 165 |
+ |
Ewald sum) to aid convergence |
| 166 |
+ |
\begin{equation} |
| 167 |
+ |
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\}. |
| 168 |
+ |
\label{eq:WolfPot} |
| 169 |
+ |
\end{equation} |
| 170 |
+ |
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
| 171 |
+ |
potential. However, neutralizing the charge contained within each |
| 172 |
+ |
cutoff sphere requires the placement of a self-image charge on the |
| 173 |
+ |
surface of the cutoff sphere. This additional self-term in the total |
| 174 |
+ |
potential enabled Wolf {\it et al.} to obtain excellent estimates of |
| 175 |
+ |
Madelung energies for many crystals. |
| 176 |
+ |
|
| 177 |
+ |
In order to use their charge-neutralized potential in molecular |
| 178 |
+ |
dynamics simulations, Wolf \textit{et al.} suggested taking the |
| 179 |
+ |
derivative of this potential prior to evaluation of the limit. This |
| 180 |
+ |
procedure gives an expression for the forces, |
| 181 |
+ |
\begin{equation} |
| 182 |
+ |
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\}, |
| 183 |
+ |
\label{eq:WolfForces} |
| 184 |
+ |
\end{equation} |
| 185 |
+ |
that incorporates both image charges and damping of the electrostatic |
| 186 |
+ |
interaction. |
| 187 |
+ |
|
| 188 |
+ |
More recently, Zahn \textit{et al.} investigated these potential and |
| 189 |
+ |
force expressions for use in simulations involving water.\cite{Zahn02} |
| 190 |
+ |
In their work, they pointed out that the forces and derivative of |
| 191 |
+ |
the potential are not commensurate. Attempts to use both |
| 192 |
+ |
Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead |
| 193 |
+ |
to poor energy conservation. They correctly observed that taking the |
| 194 |
+ |
limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the |
| 195 |
+ |
derivatives gives forces for a different potential energy function |
| 196 |
+ |
than the one shown in Eq. (\ref{eq:WolfPot}). |
| 197 |
+ |
|
| 198 |
+ |
Zahn \textit{et al.} proposed a modified form of this ``Wolf summation |
| 199 |
+ |
method'' as a way to use this technique in Molecular Dynamics |
| 200 |
+ |
simulations. Taking the integral of the forces shown in equation |
| 201 |
+ |
\ref{eq:WolfForces}, they proposed a new damped Coulomb |
| 202 |
+ |
potential, |
| 203 |
+ |
\begin{equation} |
| 204 |
+ |
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\}. |
| 205 |
+ |
\label{eq:ZahnPot} |
| 206 |
+ |
\end{equation} |
| 207 |
+ |
They showed that this potential does fairly well at capturing the |
| 208 |
+ |
structural and dynamic properties of water compared the same |
| 209 |
+ |
properties obtained using the Ewald sum. |
| 210 |
+ |
|
| 211 |
+ |
\subsection{Simple Forms for Pairwise Electrostatics} |
| 212 |
+ |
|
| 213 |
+ |
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et |
| 214 |
+ |
al.} are constructed using two different (and separable) computational |
| 215 |
+ |
tricks: \begin{enumerate} |
| 216 |
+ |
\item shifting through the use of image charges, and |
| 217 |
+ |
\item damping the electrostatic interaction. |
| 218 |
+ |
\end{enumerate} Wolf \textit{et al.} treated the |
| 219 |
+ |
development of their summation method as a progressive application of |
| 220 |
+ |
these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded |
| 221 |
+ |
their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the |
| 222 |
+ |
post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using |
| 223 |
+ |
both techniques. It is possible, however, to separate these |
| 224 |
+ |
tricks and study their effects independently. |
| 225 |
+ |
|
| 226 |
+ |
Starting with the original observation that the effective range of the |
| 227 |
+ |
electrostatic interaction in condensed phases is considerably less |
| 228 |
+ |
than $r^{-1}$, either the cutoff sphere neutralization or the |
| 229 |
+ |
distance-dependent damping technique could be used as a foundation for |
| 230 |
+ |
a new pairwise summation method. Wolf \textit{et al.} made the |
| 231 |
+ |
observation that charge neutralization within the cutoff sphere plays |
| 232 |
+ |
a significant role in energy convergence; therefore we will begin our |
| 233 |
+ |
analysis with the various shifted forms that maintain this charge |
| 234 |
+ |
neutralization. We can evaluate the methods of Wolf |
| 235 |
+ |
\textit{et al.} and Zahn \textit{et al.} by considering the standard |
| 236 |
+ |
shifted potential, |
| 237 |
+ |
\begin{equation} |
| 238 |
+ |
v_\textrm{SP}(r) = \begin{cases} |
| 239 |
+ |
v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > |
| 240 |
+ |
R_\textrm{c} |
| 241 |
+ |
\end{cases}, |
| 242 |
+ |
\label{eq:shiftingPotForm} |
| 243 |
+ |
\end{equation} |
| 244 |
+ |
and shifted force, |
| 245 |
+ |
\begin{equation} |
| 246 |
+ |
v_\textrm{SF}(r) = \begin{cases} |
| 247 |
+ |
v(r)-v_\textrm{c}-\left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) |
| 248 |
+ |
&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} |
| 249 |
+ |
\end{cases}, |
| 250 |
+ |
\label{eq:shiftingForm} |
| 251 |
+ |
\end{equation} |
| 252 |
+ |
functions where $v(r)$ is the unshifted form of the potential, and |
| 253 |
+ |
$v_c$ is $v(R_\textrm{c})$. The Shifted Force ({\sc sf}) form ensures |
| 254 |
+ |
that both the potential and the forces goes to zero at the cutoff |
| 255 |
+ |
radius, while the Shifted Potential ({\sc sp}) form only ensures the |
| 256 |
+ |
potential is smooth at the cutoff radius |
| 257 |
+ |
($R_\textrm{c}$).\cite{Allen87} |
| 258 |
+ |
|
| 259 |
+ |
The forces associated with the shifted potential are simply the forces |
| 260 |
+ |
of the unshifted potential itself (when inside the cutoff sphere), |
| 261 |
+ |
\begin{equation} |
| 262 |
+ |
f_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), |
| 263 |
+ |
\end{equation} |
| 264 |
+ |
and are zero outside. Inside the cutoff sphere, the forces associated |
| 265 |
+ |
with the shifted force form can be written, |
| 266 |
+ |
\begin{equation} |
| 267 |
+ |
f_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d |
| 268 |
+ |
v(r)}{dr} \right)_{r=R_\textrm{c}}. |
| 269 |
+ |
\end{equation} |
| 270 |
+ |
|
| 271 |
+ |
If the potential ($v(r)$) is taken to be the normal Coulomb potential, |
| 272 |
+ |
\begin{equation} |
| 273 |
+ |
v(r) = \frac{q_i q_j}{r}, |
| 274 |
+ |
\label{eq:Coulomb} |
| 275 |
+ |
\end{equation} |
| 276 |
+ |
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
| 277 |
+ |
al.}'s undamped prescription: |
| 278 |
+ |
\begin{equation} |
| 279 |
+ |
v_\textrm{SP}(r) = |
| 280 |
+ |
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
| 281 |
+ |
r\leqslant R_\textrm{c}, |
| 282 |
+ |
\label{eq:SPPot} |
| 283 |
+ |
\end{equation} |
| 284 |
+ |
with associated forces, |
| 285 |
+ |
\begin{equation} |
| 286 |
+ |
f_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
| 287 |
+ |
\label{eq:SPForces} |
| 288 |
+ |
\end{equation} |
| 289 |
+ |
These forces are identical to the forces of the standard Coulomb |
| 290 |
+ |
interaction, and cutting these off at $R_c$ was addressed by Wolf |
| 291 |
+ |
\textit{et al.} as undesirable. They pointed out that the effect of |
| 292 |
+ |
the image charges is neglected in the forces when this form is |
| 293 |
+ |
used,\cite{Wolf99} thereby eliminating any benefit from the method in |
| 294 |
+ |
molecular dynamics. Additionally, there is a discontinuity in the |
| 295 |
+ |
forces at the cutoff radius which results in energy drift during MD |
| 296 |
+ |
simulations. |
| 297 |
+ |
|
| 298 |
+ |
The shifted force ({\sc sf}) form using the normal Coulomb potential |
| 299 |
+ |
will give, |
| 300 |
+ |
\begin{equation} |
| 301 |
+ |
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}. |
| 302 |
+ |
\label{eq:SFPot} |
| 303 |
+ |
\end{equation} |
| 304 |
+ |
with associated forces, |
| 305 |
+ |
\begin{equation} |
| 306 |
+ |
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}. |
| 307 |
+ |
\label{eq:SFForces} |
| 308 |
+ |
\end{equation} |
| 309 |
+ |
This formulation has the benefits that there are no discontinuities at |
| 310 |
+ |
the cutoff distance, while the neutralizing image charges are present |
| 311 |
+ |
in both the energy and force expressions. It would be simple to add |
| 312 |
+ |
the self-neutralizing term back when computing the total energy of the |
| 313 |
+ |
system, thereby maintaining the agreement with the Madelung energies. |
| 314 |
+ |
A side effect of this treatment is the alteration in the shape of the |
| 315 |
+ |
potential that comes from the derivative term. Thus, a degree of |
| 316 |
+ |
clarity about agreement with the empirical potential is lost in order |
| 317 |
+ |
to gain functionality in dynamics simulations. |
| 318 |
+ |
|
| 319 |
+ |
Wolf \textit{et al.} originally discussed the energetics of the |
| 320 |
+ |
shifted Coulomb potential (Eq. \ref{eq:SPPot}), and they found that |
| 321 |
+ |
it was still insufficient for accurate determination of the energy |
| 322 |
+ |
with reasonable cutoff distances. The calculated Madelung energies |
| 323 |
+ |
fluctuate around the expected value with increasing cutoff radius, but |
| 324 |
+ |
the oscillations converge toward the correct value.\cite{Wolf99} A |
| 325 |
+ |
damping function was incorporated to accelerate the convergence; and |
| 326 |
+ |
though alternative functional forms could be |
| 327 |
+ |
used,\cite{Jones56,Heyes81} the complimentary error function was |
| 328 |
+ |
chosen to mirror the effective screening used in the Ewald summation. |
| 329 |
+ |
Incorporating this error function damping into the simple Coulomb |
| 330 |
+ |
potential, |
| 331 |
+ |
\begin{equation} |
| 332 |
+ |
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
| 333 |
+ |
\label{eq:dampCoulomb} |
| 334 |
+ |
\end{equation} |
| 335 |
+ |
the shifted potential (Eq. (\ref{eq:SPPot})) can be reacquired using |
| 336 |
+ |
eq. (\ref{eq:shiftingForm}), |
| 337 |
+ |
\begin{equation} |
| 338 |
+ |
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}, |
| 339 |
+ |
\label{eq:DSPPot} |
| 340 |
+ |
\end{equation} |
| 341 |
+ |
with associated forces, |
| 342 |
+ |
\begin{equation} |
| 343 |
+ |
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}. |
| 344 |
+ |
\label{eq:DSPForces} |
| 345 |
+ |
\end{equation} |
| 346 |
+ |
Again, this damped shifted potential suffers from a discontinuity and |
| 347 |
+ |
a lack of the image charges in the forces. To remedy these concerns, |
| 348 |
+ |
one may derive a {\sc sf} variant by including the derivative |
| 349 |
+ |
term in eq. (\ref{eq:shiftingForm}), |
| 350 |
+ |
\begin{equation} |
| 351 |
+ |
\begin{split} |
| 352 |
+ |
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}. |
| 353 |
+ |
\label{eq:DSFPot} |
| 354 |
+ |
\end{split} |
| 355 |
+ |
\end{equation} |
| 356 |
+ |
The derivative of the above potential will lead to the following forces, |
| 357 |
+ |
\begin{equation} |
| 358 |
+ |
\begin{split} |
| 359 |
+ |
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}. |
| 360 |
+ |
\label{eq:DSFForces} |
| 361 |
+ |
\end{split} |
| 362 |
+ |
\end{equation} |
| 363 |
+ |
If the damping parameter $(\alpha)$ is chosen to be zero, the undamped |
| 364 |
+ |
case, eqs. (\ref{eq:SPPot}-\ref{eq:SFForces}) are correctly recovered |
| 365 |
+ |
from eqs. (\ref{eq:DSPPot}-\ref{eq:DSFForces}). |
| 366 |
+ |
|
| 367 |
+ |
This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} |
| 368 |
+ |
derived by Zahn \textit{et al.}; however, there are two important |
| 369 |
+ |
differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from |
| 370 |
+ |
eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb}) |
| 371 |
+ |
with $r$ replaced by $R_\textrm{c}$. This term is {\it not} present |
| 372 |
+ |
in the Zahn potential, resulting in a potential discontinuity as |
| 373 |
+ |
particles cross $R_\textrm{c}$. Second, the sign of the derivative |
| 374 |
+ |
portion is different. The missing $v_\textrm{c}$ term would not |
| 375 |
+ |
affect molecular dynamics simulations (although the computed energy |
| 376 |
+ |
would be expected to have sudden jumps as particle distances crossed |
| 377 |
+ |
$R_c$). The sign problem would be a potential source of errors, |
| 378 |
+ |
however. In fact, it introduces a discontinuity in the forces at the |
| 379 |
+ |
cutoff, because the force function is shifted in the wrong direction |
| 380 |
+ |
and doesn't cross zero at $R_\textrm{c}$. |
| 381 |
+ |
|
| 382 |
+ |
Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
| 383 |
+ |
electrostatic summation method that is continuous in both the |
| 384 |
+ |
potential and forces and which incorporates the damping function |
| 385 |
+ |
proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of this |
| 386 |
+ |
paper, we will evaluate exactly how good these methods ({\sc sp}, {\sc |
| 387 |
+ |
sf}, damping) are at reproducing the correct electrostatic summation |
| 388 |
+ |
performed by the Ewald sum. |
| 389 |
+ |
|
| 390 |
+ |
\subsection{Other alternatives} |
| 391 |
+ |
In addition to the methods described above, we will consider some |
| 392 |
+ |
other techniques that commonly get used in molecular simulations. The |
| 393 |
+ |
simplest of these is group-based cutoffs. Though of little use for |
| 394 |
+ |
non-neutral molecules, collecting atoms into neutral groups takes |
| 395 |
+ |
advantage of the observation that the electrostatic interactions decay |
| 396 |
+ |
faster than those for monopolar pairs.\cite{Steinbach94} When |
| 397 |
+ |
considering these molecules as groups, an orientational aspect is |
| 398 |
+ |
introduced to the interactions. Consequently, as these molecular |
| 399 |
+ |
particles move through $R_\textrm{c}$, the energy will drift upward |
| 400 |
+ |
due to the anisotropy of the net molecular dipole |
| 401 |
+ |
interactions.\cite{Rahman71} To maintain good energy conservation, |
| 402 |
+ |
both the potential and derivative need to be smoothly switched to zero |
| 403 |
+ |
at $R_\textrm{c}$.\cite{Adams79} This is accomplished using a |
| 404 |
+ |
switching function, |
| 405 |
+ |
\begin{equation} |
| 406 |
+ |
S(r) = \begin{cases} 1 &\quad r\leqslant r_\textrm{sw} \\ |
| 407 |
+ |
\frac{(R_\textrm{c}+2r-3r_\textrm{sw})(R_\textrm{c}-r)^2}{(R_\textrm{c}-r_\textrm{sw})^3} &\quad r_\textrm{sw}<r\leqslant R_\textrm{c} \\ |
| 408 |
+ |
0 &\quad r>R_\textrm{c} |
| 409 |
+ |
\end{cases}, |
| 410 |
+ |
\end{equation} |
| 411 |
+ |
where the above form is for a cubic function. If a smooth second |
| 412 |
+ |
derivative is desired, a fifth (or higher) order polynomial can be |
| 413 |
+ |
used.\cite{Andrea83} |
| 414 |
+ |
|
| 415 |
+ |
Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$, |
| 416 |
+ |
and to incorporate their effect, a method like Reaction Field ({\sc |
| 417 |
+ |
rf}) can be used. The original theory for {\sc rf} was originally |
| 418 |
+ |
developed by Onsager,\cite{Onsager36} and it was applied in |
| 419 |
+ |
simulations for the study of water by Barker and Watts.\cite{Barker73} |
| 420 |
+ |
In application, it is simply an extension of the group-based cutoff |
| 421 |
+ |
method where the net dipole within the cutoff sphere polarizes an |
| 422 |
+ |
external dielectric, which reacts back on the central dipole. The |
| 423 |
+ |
same switching function considerations for group-based cutoffs need to |
| 424 |
+ |
made for {\sc rf}, with the additional pre-specification of a |
| 425 |
+ |
dielectric constant. |
| 426 |
+ |
|
| 427 |
+ |
\section{Methods} |
| 428 |
+ |
|
| 429 |
+ |
In classical molecular mechanics simulations, there are two primary |
| 430 |
+ |
techniques utilized to obtain information about the system of |
| 431 |
+ |
interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these |
| 432 |
+ |
techniques utilize pairwise summations of interactions between |
| 433 |
+ |
particle sites, but they use these summations in different ways. |
| 434 |
+ |
|
| 435 |
+ |
In MC, the potential energy difference between two subsequent |
| 436 |
+ |
configurations dictates the progression of MC sampling. Going back to |
| 437 |
+ |
the origins of this method, the acceptance criterion for the canonical |
| 438 |
+ |
ensemble laid out by Metropolis \textit{et al.} states that a |
| 439 |
+ |
subsequent configuration is accepted if $\Delta E < 0$ or if $\xi < |
| 440 |
+ |
\exp(-\Delta E/kT)$, where $\xi$ is a random number between 0 and |
| 441 |
+ |
1.\cite{Metropolis53} Maintaining the correct $\Delta E$ when using an |
| 442 |
+ |
alternate method for handling the long-range electrostatics will |
| 443 |
+ |
ensure proper sampling from the ensemble. |
| 444 |
+ |
|
| 445 |
+ |
In MD, the derivative of the potential governs how the system will |
| 446 |
+ |
progress in time. Consequently, the force and torque vectors on each |
| 447 |
+ |
body in the system dictate how the system evolves. If the magnitude |
| 448 |
+ |
and direction of these vectors are similar when using alternate |
| 449 |
+ |
electrostatic summation techniques, the dynamics in the short term |
| 450 |
+ |
will be indistinguishable. Because error in MD calculations is |
| 451 |
+ |
cumulative, one should expect greater deviation at longer times, |
| 452 |
+ |
although methods which have large differences in the force and torque |
| 453 |
+ |
vectors will diverge from each other more rapidly. |
| 454 |
+ |
|
| 455 |
+ |
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 456 |
+ |
The pairwise summation techniques (outlined in section |
| 457 |
+ |
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 458 |
+ |
studying the energy differences between conformations. We took the |
| 459 |
+ |
SPME-computed energy difference between two conformations to be the |
| 460 |
+ |
correct behavior. An ideal performance by an alternative method would |
| 461 |
+ |
reproduce these energy differences exactly. Since none of the methods |
| 462 |
+ |
provide exact energy differences, we used linear least squares |
| 463 |
+ |
regressions of the $\Delta E$ values between configurations using SPME |
| 464 |
+ |
against $\Delta E$ values using tested methods provides a quantitative |
| 465 |
+ |
comparison of this agreement. Unitary results for both the |
| 466 |
+ |
correlation and correlation coefficient for these regressions indicate |
| 467 |
+ |
equivalent energetic results between the method under consideration |
| 468 |
+ |
and electrostatics handled using SPME. Sample correlation plots for |
| 469 |
+ |
two alternate methods are shown in Fig. \ref{fig:linearFit}. |
| 470 |
+ |
|
| 471 |
|
\begin{figure} |
| 472 |
|
\centering |
| 473 |
< |
\includegraphics[width=3.25in]{./slice.pdf} |
| 473 |
> |
\includegraphics[width = \linewidth]{./dualLinear.pdf} |
| 474 |
> |
\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.} |
| 475 |
> |
\label{fig:linearFit} |
| 476 |
> |
\end{figure} |
| 477 |
> |
|
| 478 |
> |
Each system type (detailed in section \ref{sec:RepSims}) was |
| 479 |
> |
represented using 500 independent configurations. Additionally, we |
| 480 |
> |
used seven different system types, so each of the alternate |
| 481 |
> |
(non-Ewald) electrostatic summation methods was evaluated using |
| 482 |
> |
873,250 configurational energy differences. |
| 483 |
> |
|
| 484 |
> |
Results and discussion for the individual analysis of each of the |
| 485 |
> |
system types appear in the supporting information, while the |
| 486 |
> |
cumulative results over all the investigated systems appears below in |
| 487 |
> |
section \ref{sec:EnergyResults}. |
| 488 |
> |
|
| 489 |
> |
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
| 490 |
> |
We evaluated the pairwise methods (outlined in section |
| 491 |
> |
\ref{sec:ESMethods}) for use in MD simulations by |
| 492 |
> |
comparing the force and torque vectors with those obtained using the |
| 493 |
> |
reference Ewald summation (SPME). Both the magnitude and the |
| 494 |
> |
direction of these vectors on each of the bodies in the system were |
| 495 |
> |
analyzed. For the magnitude of these vectors, linear least squares |
| 496 |
> |
regression analyses were performed as described previously for |
| 497 |
> |
comparing $\Delta E$ values. Instead of a single energy difference |
| 498 |
> |
between two system configurations, we compared the magnitudes of the |
| 499 |
> |
forces (and torques) on each molecule in each configuration. For a |
| 500 |
> |
system of 1000 water molecules and 40 ions, there are 1040 force |
| 501 |
> |
vectors and 1000 torque vectors. With 500 configurations, this |
| 502 |
> |
results in 520,000 force and 500,000 torque vector comparisons. |
| 503 |
> |
Additionally, data from seven different system types was aggregated |
| 504 |
> |
before the comparison was made. |
| 505 |
> |
|
| 506 |
> |
The {\it directionality} of the force and torque vectors was |
| 507 |
> |
investigated through measurement of the angle ($\theta$) formed |
| 508 |
> |
between those computed from the particular method and those from SPME, |
| 509 |
> |
\begin{equation} |
| 510 |
> |
\theta_f = \cos^{-1} \hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}, |
| 511 |
> |
\end{equation} |
| 512 |
> |
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the |
| 513 |
> |
force vector computed using method $M$. |
| 514 |
> |
|
| 515 |
> |
Each of these $\theta$ values was accumulated in a distribution |
| 516 |
> |
function, weighted by the area on the unit sphere. Non-linear |
| 517 |
> |
Gaussian fits were used to measure the width of the resulting |
| 518 |
> |
distributions. |
| 519 |
> |
|
| 520 |
> |
\begin{figure} |
| 521 |
> |
\centering |
| 522 |
> |
\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 523 |
> |
\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.} |
| 524 |
> |
\label{fig:gaussian} |
| 525 |
> |
\end{figure} |
| 526 |
> |
|
| 527 |
> |
Figure \ref{fig:gaussian} shows an example distribution with applied |
| 528 |
> |
non-linear fits. The solid line is a Gaussian profile, while the |
| 529 |
> |
dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 530 |
> |
Lorentzian. Since this distribution is a measure of angular error |
| 531 |
> |
between two different electrostatic summation methods, there is no |
| 532 |
> |
{\it a priori} reason for the profile to adhere to any specific shape. |
| 533 |
> |
Gaussian fits was used to compare all the tested methods. The |
| 534 |
> |
variance ($\sigma^2$) was extracted from each of these fits and was |
| 535 |
> |
used to compare distribution widths. Values of $\sigma^2$ near zero |
| 536 |
> |
indicate vector directions indistinguishable from those calculated |
| 537 |
> |
when using the reference method (SPME). |
| 538 |
> |
|
| 539 |
> |
\subsection{Short-time Dynamics} |
| 540 |
> |
|
| 541 |
> |
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 542 |
> |
Evaluation of the long-time dynamics of charged systems was performed |
| 543 |
> |
by considering the NaCl crystal system while using a subset of the |
| 544 |
> |
best performing pairwise methods. The NaCl crystal was chosen to |
| 545 |
> |
avoid possible complications involving the propagation techniques of |
| 546 |
> |
orientational motion in molecular systems. To enhance the atomic |
| 547 |
> |
motion, these crystals were equilibrated at 1000 K, near the |
| 548 |
> |
experimental $T_m$ for NaCl. Simulations were performed under the |
| 549 |
> |
microcanonical ensemble, and velocity autocorrelation functions |
| 550 |
> |
(Eq. \ref{eq:vCorr}) were computed for each of the trajectories, |
| 551 |
> |
\begin{equation} |
| 552 |
> |
C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle. |
| 553 |
> |
\label{eq:vCorr} |
| 554 |
> |
\end{equation} |
| 555 |
> |
Velocity autocorrelation functions require detailed short time data |
| 556 |
> |
and long trajectories for good statistics, thus velocity information |
| 557 |
> |
was saved every 5 fs over 100 ps trajectories. The power spectrum |
| 558 |
> |
($I(\omega)$) is obtained via Fourier transform of the autocorrelation |
| 559 |
> |
function |
| 560 |
> |
\begin{equation} |
| 561 |
> |
I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 562 |
> |
\label{eq:powerSpec} |
| 563 |
> |
\end{equation} |
| 564 |
> |
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. |
| 565 |
> |
|
| 566 |
> |
\subsection{Representative Simulations}\label{sec:RepSims} |
| 567 |
> |
A variety of common and representative simulations were analyzed to |
| 568 |
> |
determine the relative effectiveness of the pairwise summation |
| 569 |
> |
techniques in reproducing the energetics and dynamics exhibited by |
| 570 |
> |
SPME. The studied systems were as follows: |
| 571 |
> |
\begin{enumerate} |
| 572 |
> |
\item Liquid Water |
| 573 |
> |
\item Crystalline Water (Ice I$_\textrm{c}$) |
| 574 |
> |
\item NaCl Crystal |
| 575 |
> |
\item NaCl Melt |
| 576 |
> |
\item Low Ionic Strength Solution of NaCl in Water |
| 577 |
> |
\item High Ionic Strength Solution of NaCl in Water |
| 578 |
> |
\item 6 \AA\ Radius Sphere of Argon in Water |
| 579 |
> |
\end{enumerate} |
| 580 |
> |
By utilizing the pairwise techniques (outlined in section |
| 581 |
> |
\ref{sec:ESMethods}) in systems composed entirely of neutral groups, |
| 582 |
> |
charged particles, and mixtures of the two, we can comment on possible |
| 583 |
> |
system dependence and/or universal applicability of the techniques. |
| 584 |
> |
|
| 585 |
> |
Generation of the system configurations was dependent on the system |
| 586 |
> |
type. For the solid and liquid water configurations, configuration |
| 587 |
> |
snapshots were taken at regular intervals from higher temperature 1000 |
| 588 |
> |
SPC/E water molecule trajectories and each equilibrated individually. |
| 589 |
> |
The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- |
| 590 |
> |
ions and were selected and equilibrated in the same fashion as the |
| 591 |
> |
water systems. For the low and high ionic strength NaCl solutions, 4 |
| 592 |
> |
and 40 ions were first solvated in a 1000 water molecule boxes |
| 593 |
> |
respectively. Ion and water positions were then randomly swapped, and |
| 594 |
> |
the resulting configurations were again equilibrated individually. |
| 595 |
> |
Finally, for the Argon/Water "charge void" systems, the identities of |
| 596 |
> |
all the SPC/E waters within 6 \AA\ of the center of the equilibrated |
| 597 |
> |
water configurations were converted to argon |
| 598 |
> |
(Fig. \ref{fig:argonSlice}). |
| 599 |
> |
|
| 600 |
> |
\begin{figure} |
| 601 |
> |
\centering |
| 602 |
> |
\includegraphics[width = \linewidth]{./slice.pdf} |
| 603 |
|
\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.} |
| 604 |
< |
\label{argonSlice} |
| 604 |
> |
\label{fig:argonSlice} |
| 605 |
|
\end{figure} |
| 606 |
|
|
| 607 |
< |
All of these comparisons were performed with three different cutoff radii (9, 12, and 15 \AA) to investigate the cutoff radius dependence of the various techniques. It should be noted that the damping parameter chosen in SPME, or so called ``Ewald Coefficient", has a significant effect on the energies and forces calculated. Typical molecular mechanics packages default this to a value dependent on the cutoff radius and a tolerance (typically less than $1 \times 10^{-5}$ kcal/mol). We chose a tolerance of $1 \times 10^{-8}$, resulting in Ewald Coefficients of 0.4200, 0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 607 |
> |
\subsection{Electrostatic Summation Methods}\label{sec:ESMethods} |
| 608 |
> |
Electrostatic summation method comparisons were performed using SPME, |
| 609 |
> |
the {\sc sp} and {\sc sf} methods - both with damping |
| 610 |
> |
parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak, |
| 611 |
> |
moderate, and strong damping respectively), reaction field with an |
| 612 |
> |
infinite dielectric constant, and an unmodified cutoff. Group-based |
| 613 |
> |
cutoffs with a fifth-order polynomial switching function were |
| 614 |
> |
necessary for the reaction field simulations and were utilized in the |
| 615 |
> |
SP, SF, and pure cutoff methods for comparison to the standard lack of |
| 616 |
> |
group-based cutoffs with a hard truncation. The SPME calculations |
| 617 |
> |
were performed using the TINKER implementation of SPME,\cite{Ponder87} |
| 618 |
> |
while all other method calculations were performed using the OOPSE |
| 619 |
> |
molecular mechanics package.\cite{Meineke05} |
| 620 |
|
|
| 621 |
+ |
These methods were additionally evaluated with three different cutoff |
| 622 |
+ |
radii (9, 12, and 15 \AA) to investigate possible cutoff radius |
| 623 |
+ |
dependence. It should be noted that the damping parameter chosen in |
| 624 |
+ |
SPME, or so called ``Ewald Coefficient", has a significant effect on |
| 625 |
+ |
the energies and forces calculated. Typical molecular mechanics |
| 626 |
+ |
packages default this to a value dependent on the cutoff radius and a |
| 627 |
+ |
tolerance (typically less than $1 \times 10^{-4}$ kcal/mol). Smaller |
| 628 |
+ |
tolerances are typically associated with increased accuracy, but this |
| 629 |
+ |
usually means more time spent calculating the reciprocal-space portion |
| 630 |
+ |
of the summation.\cite{Perram88,Essmann95} The default TINKER |
| 631 |
+ |
tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME |
| 632 |
+ |
calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and |
| 633 |
+ |
0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 634 |
+ |
|
| 635 |
|
\section{Results and Discussion} |
| 636 |
|
|
| 637 |
< |
In order to evaluate the performance of the adapted Wolf SP and SF electrostatic summation methods for Monte Carlo simulations, the energy differences between configurations need to be compared to the results using SPME. Considering the SPME results to be the correct or desired behavior, ideal performance of a tested method is taken to be agreement between the energy differences calculated. Linear least squares regression of the $\Delta$E values between configurations using SPME against $\Delta$E values using tested methods provides a quantitative comparison of this agreement. Unitary results for both the correlation and correlation coefficient for these regressions indicate equivalent energetic results between the methods. The correlation is the slope of the plotted data while the correlation coefficient ($R^2$) is a measure of the of the data scatter around the fitted line and gives an idea of the quality of the fit (Fig. \ref{linearFit}). |
| 637 |
> |
\subsection{Configuration Energy Differences}\label{sec:EnergyResults} |
| 638 |
> |
In order to evaluate the performance of the pairwise electrostatic |
| 639 |
> |
summation methods for Monte Carlo simulations, the energy differences |
| 640 |
> |
between configurations were compared to the values obtained when using |
| 641 |
> |
SPME. The results for the subsequent regression analysis are shown in |
| 642 |
> |
figure \ref{fig:delE}. |
| 643 |
|
|
| 644 |
|
\begin{figure} |
| 645 |
|
\centering |
| 646 |
< |
\includegraphics[width=3.25in]{./linearFit.pdf} |
| 647 |
< |
\caption{Example least squares regression of the $\Delta$E between configurations for the SF method against SPME in the pure water system. } |
| 648 |
< |
\label{linearFit} |
| 646 |
> |
\includegraphics[width=5.5in]{./delEplot.pdf} |
| 647 |
> |
\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).} |
| 648 |
> |
\label{fig:delE} |
| 649 |
|
\end{figure} |
| 650 |
|
|
| 651 |
< |
With 500 independent configurations, 124,750 $\Delta$E data points are used in a regression of a single system. A table with the results for analysis of To gauge the applicability of each method in the general case, all the different system types were included in a separate Figure \ref{delEplot} shows the results for analysis of all the simulation types |
| 651 |
> |
In this figure, it is apparent that it is unreasonable to expect |
| 652 |
> |
realistic results using an unmodified cutoff. This is not all that |
| 653 |
> |
surprising since this results in large energy fluctuations as atoms |
| 654 |
> |
move in and out of the cutoff radius. These fluctuations can be |
| 655 |
> |
alleviated to some degree by using group based cutoffs with a |
| 656 |
> |
switching function.\cite{Steinbach94} The Group Switch Cutoff row |
| 657 |
> |
doesn't show a significant improvement in this plot because the salt |
| 658 |
> |
and salt solution systems contain non-neutral groups, see the |
| 659 |
> |
accompanying supporting information for a comparison where all groups |
| 660 |
> |
are neutral. |
| 661 |
|
|
| 662 |
+ |
Correcting the resulting charged cutoff sphere is one of the purposes |
| 663 |
+ |
of the damped Coulomb summation proposed by Wolf \textit{et |
| 664 |
+ |
al.},\cite{Wolf99} and this correction indeed improves the results as |
| 665 |
+ |
seen in the Shifted-Potental rows. While the undamped case of this |
| 666 |
+ |
method is a significant improvement over the pure cutoff, it still |
| 667 |
+ |
doesn't correlate that well with SPME. Inclusion of potential damping |
| 668 |
+ |
improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 669 |
+ |
an excellent correlation and quality of fit with the SPME results, |
| 670 |
+ |
particularly with a cutoff radius greater than 12 \AA . Use of a |
| 671 |
+ |
larger damping parameter is more helpful for the shortest cutoff |
| 672 |
+ |
shown, but it has a detrimental effect on simulations with larger |
| 673 |
+ |
cutoffs. In the {\sc sf} sets, increasing damping results in |
| 674 |
+ |
progressively poorer correlation. Overall, the undamped case is the |
| 675 |
+ |
best performing set, as the correlation and quality of fits are |
| 676 |
+ |
consistently superior regardless of the cutoff distance. This result |
| 677 |
+ |
is beneficial in that the undamped case is less computationally |
| 678 |
+ |
prohibitive do to the lack of complimentary error function calculation |
| 679 |
+ |
when performing the electrostatic pair interaction. The reaction |
| 680 |
+ |
field results illustrates some of that method's limitations, primarily |
| 681 |
+ |
that it was developed for use in homogenous systems; although it does |
| 682 |
+ |
provide results that are an improvement over those from an unmodified |
| 683 |
+ |
cutoff. |
| 684 |
+ |
|
| 685 |
+ |
\subsection{Magnitudes of the Force and Torque Vectors} |
| 686 |
+ |
|
| 687 |
+ |
Evaluation of pairwise methods for use in Molecular Dynamics |
| 688 |
+ |
simulations requires consideration of effects on the forces and |
| 689 |
+ |
torques. Investigation of the force and torque vector magnitudes |
| 690 |
+ |
provides a measure of the strength of these values relative to SPME. |
| 691 |
+ |
Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the |
| 692 |
+ |
force and torque vector magnitude regression results for the |
| 693 |
+ |
accumulated analysis over all the system types. |
| 694 |
+ |
|
| 695 |
+ |
\begin{figure} |
| 696 |
+ |
\centering |
| 697 |
+ |
\includegraphics[width=5.5in]{./frcMagplot.pdf} |
| 698 |
+ |
\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).} |
| 699 |
+ |
\label{fig:frcMag} |
| 700 |
+ |
\end{figure} |
| 701 |
+ |
|
| 702 |
+ |
Figure \ref{fig:frcMag}, for the most part, parallels the results seen |
| 703 |
+ |
in the previous $\Delta E$ section. The unmodified cutoff results are |
| 704 |
+ |
poor, but using group based cutoffs and a switching function provides |
| 705 |
+ |
a improvement much more significant than what was seen with $\Delta |
| 706 |
+ |
E$. Looking at the {\sc sp} sets, the slope and $R^2$ |
| 707 |
+ |
improve with the use of damping to an optimal result of 0.2 \AA |
| 708 |
+ |
$^{-1}$ for the 12 and 15 \AA\ cutoffs. Further increases in damping, |
| 709 |
+ |
while beneficial for simulations with a cutoff radius of 9 \AA\ , is |
| 710 |
+ |
detrimental to simulations with larger cutoff radii. The undamped |
| 711 |
+ |
{\sc sf} method gives forces in line with those obtained using |
| 712 |
+ |
SPME, and use of a damping function results in minor improvement. The |
| 713 |
+ |
reaction field results are surprisingly good, considering the poor |
| 714 |
+ |
quality of the fits for the $\Delta E$ results. There is still a |
| 715 |
+ |
considerable degree of scatter in the data, but it correlates well in |
| 716 |
+ |
general. To be fair, we again note that the reaction field |
| 717 |
+ |
calculations do not encompass NaCl crystal and melt systems, so these |
| 718 |
+ |
results are partly biased towards conditions in which the method |
| 719 |
+ |
performs more favorably. |
| 720 |
+ |
|
| 721 |
+ |
\begin{figure} |
| 722 |
+ |
\centering |
| 723 |
+ |
\includegraphics[width=5.5in]{./trqMagplot.pdf} |
| 724 |
+ |
\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).} |
| 725 |
+ |
\label{fig:trqMag} |
| 726 |
+ |
\end{figure} |
| 727 |
+ |
|
| 728 |
+ |
To evaluate the torque vector magnitudes, the data set from which |
| 729 |
+ |
values are drawn is limited to rigid molecules in the systems |
| 730 |
+ |
(i.e. water molecules). In spite of this smaller sampling pool, the |
| 731 |
+ |
torque vector magnitude results in figure \ref{fig:trqMag} are still |
| 732 |
+ |
similar to those seen for the forces; however, they more clearly show |
| 733 |
+ |
the improved behavior that comes with increasing the cutoff radius. |
| 734 |
+ |
Moderate damping is beneficial to the {\sc sp} and helpful |
| 735 |
+ |
yet possibly unnecessary with the {\sc sf} method, and they also |
| 736 |
+ |
show that over-damping adversely effects all cutoff radii rather than |
| 737 |
+ |
showing an improvement for systems with short cutoffs. The reaction |
| 738 |
+ |
field method performs well when calculating the torques, better than |
| 739 |
+ |
the Shifted Force method over this limited data set. |
| 740 |
+ |
|
| 741 |
+ |
\subsection{Directionality of the Force and Torque Vectors} |
| 742 |
+ |
|
| 743 |
+ |
Having force and torque vectors with magnitudes that are well |
| 744 |
+ |
correlated to SPME is good, but if they are not pointing in the proper |
| 745 |
+ |
direction the results will be incorrect. These vector directions were |
| 746 |
+ |
investigated through measurement of the angle formed between them and |
| 747 |
+ |
those from SPME. The results (Fig. \ref{fig:frcTrqAng}) are compared |
| 748 |
+ |
through the variance ($\sigma^2$) of the Gaussian fits of the angle |
| 749 |
+ |
error distributions of the combined set over all system types. |
| 750 |
+ |
|
| 751 |
+ |
\begin{figure} |
| 752 |
+ |
\centering |
| 753 |
+ |
\includegraphics[width=5.5in]{./frcTrqAngplot.pdf} |
| 754 |
+ |
\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).} |
| 755 |
+ |
\label{fig:frcTrqAng} |
| 756 |
+ |
\end{figure} |
| 757 |
+ |
|
| 758 |
+ |
Both the force and torque $\sigma^2$ results from the analysis of the |
| 759 |
+ |
total accumulated system data are tabulated in figure |
| 760 |
+ |
\ref{fig:frcTrqAng}. All of the sets, aside from the over-damped case |
| 761 |
+ |
show the improvement afforded by choosing a longer simulation cutoff. |
| 762 |
+ |
Increasing the cutoff from 9 to 12 \AA\ typically results in a halving |
| 763 |
+ |
of the distribution widths, with a similar improvement going from 12 |
| 764 |
+ |
to 15 \AA . The undamped {\sc sf}, Group Based Cutoff, and |
| 765 |
+ |
Reaction Field methods all do equivalently well at capturing the |
| 766 |
+ |
direction of both the force and torque vectors. Using damping |
| 767 |
+ |
improves the angular behavior significantly for the {\sc sp} |
| 768 |
+ |
and moderately for the {\sc sf} methods. Increasing the damping |
| 769 |
+ |
too far is destructive for both methods, particularly to the torque |
| 770 |
+ |
vectors. Again it is important to recognize that the force vectors |
| 771 |
+ |
cover all particles in the systems, while torque vectors are only |
| 772 |
+ |
available for neutral molecular groups. Damping appears to have a |
| 773 |
+ |
more beneficial effect on non-neutral bodies, and this observation is |
| 774 |
+ |
investigated further in the accompanying supporting information. |
| 775 |
+ |
|
| 776 |
+ |
\begin{table}[htbp] |
| 777 |
+ |
\centering |
| 778 |
+ |
\caption{Variance ($\sigma^2$) of the force (top set) and torque (bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods. Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function. The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.} |
| 779 |
+ |
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 780 |
+ |
\\ |
| 781 |
+ |
\toprule |
| 782 |
+ |
& & \multicolumn{4}{c}{Shifted Potential} & \multicolumn{4}{c}{Shifted Force} \\ |
| 783 |
+ |
\cmidrule(lr){3-6} |
| 784 |
+ |
\cmidrule(l){7-10} |
| 785 |
+ |
$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$\\ |
| 786 |
+ |
\midrule |
| 787 |
+ |
|
| 788 |
+ |
9 \AA & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\ |
| 789 |
+ |
& \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} \\ |
| 790 |
+ |
12 \AA & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\ |
| 791 |
+ |
& \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} \\ |
| 792 |
+ |
15 \AA & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\ |
| 793 |
+ |
& \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} \\ |
| 794 |
+ |
|
| 795 |
+ |
\midrule |
| 796 |
+ |
|
| 797 |
+ |
9 \AA & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\ |
| 798 |
+ |
& \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} \\ |
| 799 |
+ |
12 \AA & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\ |
| 800 |
+ |
& \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} \\ |
| 801 |
+ |
15 \AA & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\ |
| 802 |
+ |
& \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} \\ |
| 803 |
+ |
|
| 804 |
+ |
\bottomrule |
| 805 |
+ |
\end{tabular} |
| 806 |
+ |
\label{tab:groupAngle} |
| 807 |
+ |
\end{table} |
| 808 |
+ |
|
| 809 |
+ |
Although not discussed previously, group based cutoffs can be applied |
| 810 |
+ |
to both the {\sc sp} and {\sc sf} methods. Use off a |
| 811 |
+ |
switching function corrects for the discontinuities that arise when |
| 812 |
+ |
atoms of a group exit the cutoff before the group's center of mass. |
| 813 |
+ |
Though there are no significant benefit or drawbacks observed in |
| 814 |
+ |
$\Delta E$ and vector magnitude results when doing this, there is a |
| 815 |
+ |
measurable improvement in the vector angle results. Table |
| 816 |
+ |
\ref{tab:groupAngle} shows the angular variance values obtained using |
| 817 |
+ |
group based cutoffs and a switching function alongside the standard |
| 818 |
+ |
results seen in figure \ref{fig:frcTrqAng} for comparison purposes. |
| 819 |
+ |
The {\sc sp} shows much narrower angular distributions for |
| 820 |
+ |
both the force and torque vectors when using an $\alpha$ of 0.2 |
| 821 |
+ |
\AA$^{-1}$ or less, while {\sc sf} shows improvements in the |
| 822 |
+ |
undamped and lightly damped cases. Thus, by calculating the |
| 823 |
+ |
electrostatic interactions in terms of molecular pairs rather than |
| 824 |
+ |
atomic pairs, the direction of the force and torque vectors are |
| 825 |
+ |
determined more accurately. |
| 826 |
+ |
|
| 827 |
+ |
One additional trend to recognize in table \ref{tab:groupAngle} is |
| 828 |
+ |
that the $\sigma^2$ values for both {\sc sp} and |
| 829 |
+ |
{\sc sf} converge as $\alpha$ increases, something that is easier |
| 830 |
+ |
to see when using group based cutoffs. Looking back on figures |
| 831 |
+ |
\ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this |
| 832 |
+ |
behavior clearly at large $\alpha$ and cutoff values. The reason for |
| 833 |
+ |
this is that the complimentary error function inserted into the |
| 834 |
+ |
potential weakens the electrostatic interaction as $\alpha$ increases. |
| 835 |
+ |
Thus, at larger values of $\alpha$, both the summation method types |
| 836 |
+ |
progress toward non-interacting functions, so care is required in |
| 837 |
+ |
choosing large damping functions lest one generate an undesirable loss |
| 838 |
+ |
in the pair interaction. Kast \textit{et al.} developed a method for |
| 839 |
+ |
choosing appropriate $\alpha$ values for these types of electrostatic |
| 840 |
+ |
summation methods by fitting to $g(r)$ data, and their methods |
| 841 |
+ |
indicate optimal values of 0.34, 0.25, and 0.16 \AA$^{-1}$ for cutoff |
| 842 |
+ |
values of 9, 12, and 15 \AA\ respectively.\cite{Kast03} These appear |
| 843 |
+ |
to be reasonable choices to obtain proper MC behavior |
| 844 |
+ |
(Fig. \ref{fig:delE}); however, based on these findings, choices this |
| 845 |
+ |
high would introduce error in the molecular torques, particularly for |
| 846 |
+ |
the shorter cutoffs. Based on the above findings, empirical damping |
| 847 |
+ |
up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably |
| 848 |
+ |
unnecessary when using the {\sc sf} method. |
| 849 |
+ |
|
| 850 |
+ |
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 851 |
+ |
|
| 852 |
+ |
In the previous studies using a {\sc sf} variant of the damped |
| 853 |
+ |
Wolf coulomb potential, the structure and dynamics of water were |
| 854 |
+ |
investigated rather extensively.\cite{Zahn02,Kast03} Their results |
| 855 |
+ |
indicated that the damped {\sc sf} method results in properties |
| 856 |
+ |
very similar to those obtained when using the Ewald summation. |
| 857 |
+ |
Considering the statistical results shown above, the good performance |
| 858 |
+ |
of this method is not that surprising. Rather than consider the same |
| 859 |
+ |
systems and simply recapitulate their results, we decided to look at |
| 860 |
+ |
the solid state dynamical behavior obtained using the best performing |
| 861 |
+ |
summation methods from the above results. |
| 862 |
+ |
|
| 863 |
+ |
\begin{figure} |
| 864 |
+ |
\centering |
| 865 |
+ |
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 866 |
+ |
\caption{Power spectra obtained from the velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, {\sc sf} ($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). Apodization of the correlation functions via a cubic switching function between 40 and 50 ps was used to clear up the spectral noise resulting from data truncation, and had no noticeable effect on peak location or magnitude. The inset shows the frequency region below 100 cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 867 |
+ |
\label{fig:methodPS} |
| 868 |
+ |
\end{figure} |
| 869 |
+ |
|
| 870 |
+ |
Figure \ref{fig:methodPS} shows the power spectra for the NaCl |
| 871 |
+ |
crystals (from averaged Na and Cl ion velocity autocorrelation |
| 872 |
+ |
functions) using the stated electrostatic summation methods. While |
| 873 |
+ |
high frequency peaks of all the spectra overlap, showing the same |
| 874 |
+ |
general features, the low frequency region shows how the summation |
| 875 |
+ |
methods differ. Considering the low-frequency inset (expanded in the |
| 876 |
+ |
upper frame of figure \ref{fig:dampInc}), at frequencies below 100 |
| 877 |
+ |
cm$^{-1}$, the correlated motions are blue-shifted when using undamped |
| 878 |
+ |
or weakly damped {\sc sf}. When using moderate damping ($\alpha |
| 879 |
+ |
= 0.2$ \AA$^{-1}$) both the {\sc sf} and {\sc sp} |
| 880 |
+ |
methods give near identical correlated motion behavior as the Ewald |
| 881 |
+ |
method (which has a damping value of 0.3119). The damping acts as a |
| 882 |
+ |
distance dependent Gaussian screening of the point charges for the |
| 883 |
+ |
pairwise summation methods. This weakening of the electrostatic |
| 884 |
+ |
interaction with distance explains why the long-ranged correlated |
| 885 |
+ |
motions are at lower frequencies for the moderately damped methods |
| 886 |
+ |
than for undamped or weakly damped methods. To see this effect more |
| 887 |
+ |
clearly, we show how damping strength affects a simple real-space |
| 888 |
+ |
electrostatic potential, |
| 889 |
+ |
\begin{equation} |
| 890 |
+ |
V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r), |
| 891 |
+ |
\end{equation} |
| 892 |
+ |
where $S(r)$ is a switching function that smoothly zeroes the |
| 893 |
+ |
potential at the cutoff radius. Figure \ref{fig:dampInc} shows how |
| 894 |
+ |
the low frequency motions are dependent on the damping used in the |
| 895 |
+ |
direct electrostatic sum. As the damping increases, the peaks drop to |
| 896 |
+ |
lower frequencies. Incidentally, use of an $\alpha$ of 0.25 |
| 897 |
+ |
\AA$^{-1}$ on a simple electrostatic summation results in low |
| 898 |
+ |
frequency correlated dynamics equivalent to a simulation using SPME. |
| 899 |
+ |
When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks |
| 900 |
+ |
shift to higher frequency in exponential fashion. Though not shown, |
| 901 |
+ |
the spectrum for the simple undamped electrostatic potential is |
| 902 |
+ |
blue-shifted such that the lowest frequency peak resides near 325 |
| 903 |
+ |
cm$^{-1}$. In light of these results, the undamped {\sc sf} |
| 904 |
+ |
method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is |
| 905 |
+ |
quite respectable; however, it appears as though moderate damping is |
| 906 |
+ |
required for accurate reproduction of crystal dynamics. |
| 907 |
+ |
\begin{figure} |
| 908 |
+ |
\centering |
| 909 |
+ |
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 910 |
+ |
\caption{Regions of spectra showing the low-frequency correlated motions for NaCl crystals at 1000 K using various electrostatic summation methods. The upper plot is a zoomed inset from figure \ref{fig:methodPS}. As the damping value for the {\sc sf} potential increases, the low-frequency peaks red-shift. The lower plot is of spectra when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$. As $\alpha$ increases, the peaks are red-shifted toward and eventually beyond the values given by SPME. The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.} |
| 911 |
+ |
\label{fig:dampInc} |
| 912 |
+ |
\end{figure} |
| 913 |
+ |
|
| 914 |
|
\section{Conclusions} |
| 915 |
|
|
| 916 |
< |
\section{Acknowledgments} |
| 916 |
> |
This investigation of pairwise electrostatic summation techniques |
| 917 |
> |
shows that there are viable and more computationally efficient |
| 918 |
> |
electrostatic summation techniques than the Ewald summation, chiefly |
| 919 |
> |
methods derived from the damped Coulombic sum originally proposed by |
| 920 |
> |
Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the |
| 921 |
> |
{\sc sf} method, reformulated above as eq. (\ref{eq:DSFPot}), |
| 922 |
> |
shows a remarkable ability to reproduce the energetic and dynamic |
| 923 |
> |
characteristics exhibited by simulations employing lattice summation |
| 924 |
> |
techniques. The cumulative energy difference results showed the |
| 925 |
> |
undamped {\sc sf} and moderately damped {\sc sp} methods |
| 926 |
> |
produced results nearly identical to SPME. Similarly for the dynamic |
| 927 |
> |
features, the undamped or moderately damped {\sc sf} and |
| 928 |
> |
moderately damped {\sc sp} methods produce force and torque |
| 929 |
> |
vector magnitude and directions very similar to the expected values. |
| 930 |
> |
These results translate into long-time dynamic behavior equivalent to |
| 931 |
> |
that produced in simulations using SPME. |
| 932 |
|
|
| 933 |
< |
\newpage |
| 933 |
> |
Aside from the computational cost benefit, these techniques have |
| 934 |
> |
applicability in situations where the use of the Ewald sum can prove |
| 935 |
> |
problematic. Primary among them is their use in interfacial systems, |
| 936 |
> |
where the unmodified lattice sum techniques artificially accentuate |
| 937 |
> |
the periodicity of the system in an undesirable manner. There have |
| 938 |
> |
been alterations to the standard Ewald techniques, via corrections and |
| 939 |
> |
reformulations, to compensate for these systems; but the pairwise |
| 940 |
> |
techniques discussed here require no modifications, making them |
| 941 |
> |
natural tools to tackle these problems. Additionally, this |
| 942 |
> |
transferability gives them benefits over other pairwise methods, like |
| 943 |
> |
reaction field, because estimations of physical properties (e.g. the |
| 944 |
> |
dielectric constant) are unnecessary. |
| 945 |
|
|
| 946 |
< |
\bibliographystyle{achemso} |
| 946 |
> |
We are not suggesting any flaw with the Ewald sum; in fact, it is the |
| 947 |
> |
standard by which these simple pairwise sums are judged. However, |
| 948 |
> |
these results do suggest that in the typical simulations performed |
| 949 |
> |
today, the Ewald summation may no longer be required to obtain the |
| 950 |
> |
level of accuracy most researcher have come to expect |
| 951 |
> |
|
| 952 |
> |
\section{Acknowledgments} |
| 953 |
> |
\newpage |
| 954 |
> |
|
| 955 |
> |
\bibliographystyle{jcp2} |
| 956 |
|
\bibliography{electrostaticMethods} |
| 957 |
|
|
| 958 |
|
|
