| 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} |
| 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. |
| 80 |
> |
increases proportionally with the cutoff sphere, it quickly becomes |
| 81 |
> |
very time-consuming to perform these calculations. |
| 82 |
|
|
| 83 |
+ |
There have been many efforts to address this issue of both proper and |
| 84 |
+ |
practical handling of electrostatic interactions, and these have |
| 85 |
+ |
resulted in the availability of a variety of |
| 86 |
+ |
techniques.\cite{Roux99,Sagui99,Tobias01} These are typically |
| 87 |
+ |
classified as implicit methods (i.e., continuum dielectrics, static |
| 88 |
+ |
dipolar fields),\cite{Born20,Grossfield00} explicit methods (i.e., |
| 89 |
+ |
Ewald summations, interaction shifting or |
| 90 |
+ |
truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., |
| 91 |
+ |
reaction field type methods, fast multipole |
| 92 |
+ |
methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are |
| 93 |
+ |
often preferred because they incorporate dynamic solvent molecules in |
| 94 |
+ |
the system of interest, but these methods are sometimes difficult to |
| 95 |
+ |
utilize because of their high computational cost.\cite{Roux99} In |
| 96 |
+ |
addition to this cost, there has been some question of the inherent |
| 97 |
+ |
periodicity of the explicit Ewald summation artificially influencing |
| 98 |
+ |
systems dynamics.\cite{Tobias01} |
| 99 |
+ |
|
| 100 |
+ |
In this paper, we focus on the common mixed and explicit methods of |
| 101 |
+ |
reaction filed and smooth particle mesh |
| 102 |
+ |
Ewald\cite{Onsager36,Essmann99} and a new set of shifted methods |
| 103 |
+ |
devised by Wolf {\it et al.} which we further extend.\cite{Wolf99} |
| 104 |
+ |
These new methods for handling electrostatics are quite |
| 105 |
+ |
computationally efficient, since they involve only a simple |
| 106 |
+ |
modification to the direct pairwise sum, and they lack the added |
| 107 |
+ |
periodicity of the Ewald sum. Below, these methods are evaluated using |
| 108 |
+ |
a variety of model systems and comparison methodologies to establish |
| 109 |
+ |
their usability in molecular simulations. |
| 110 |
+ |
|
| 111 |
|
\subsection{The Ewald Sum} |
| 112 |
< |
blah blah blah Ewald Sum Important blah blah blah |
| 112 |
> |
The complete accumulation electrostatic interactions in a system with |
| 113 |
> |
periodic boundary conditions (PBC) requires the consideration of the |
| 114 |
> |
effect of all charges within a simulation box, as well as those in the |
| 115 |
> |
periodic replicas, |
| 116 |
> |
\begin{equation} |
| 117 |
> |
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], |
| 118 |
> |
\label{eq:PBCSum} |
| 119 |
> |
\end{equation} |
| 120 |
> |
where the sum over $\mathbf{n}$ is a sum over all periodic box |
| 121 |
> |
replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the |
| 122 |
> |
prime indicates $i = j$ are neglected for $\mathbf{n} = |
| 123 |
> |
0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic |
| 124 |
> |
particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is |
| 125 |
> |
the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and |
| 126 |
> |
$j$, and $\phi$ is Poisson's equation ($\phi(\mathbf{r}_{ij}) = q_i |
| 127 |
> |
q_j |\mathbf{r}_{ij}|^{-1}$ for charge-charge interactions). In the |
| 128 |
> |
case of monopole electrostatics, eq. (\ref{eq:PBCSum}) is |
| 129 |
> |
conditionally convergent and is discontinuous for non-neutral systems. |
| 130 |
|
|
| 131 |
+ |
This electrostatic summation problem was originally studied by Ewald |
| 132 |
+ |
for the case of an infinite crystal.\cite{Ewald21}. The approach he |
| 133 |
+ |
took was to convert this conditionally convergent sum into two |
| 134 |
+ |
absolutely convergent summations: a short-ranged real-space summation |
| 135 |
+ |
and a long-ranged reciprocal-space summation, |
| 136 |
+ |
\begin{equation} |
| 137 |
+ |
\begin{split} |
| 138 |
+ |
V_\textrm{elec} = \frac{1}{2}& \sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime \frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)}{|\mathbf{r}_{ij}+\mathbf{n}|} \\ &+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} \exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) \cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ &- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 + \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3}\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2, |
| 139 |
+ |
\end{split} |
| 140 |
+ |
\label{eq:EwaldSum} |
| 141 |
+ |
\end{equation} |
| 142 |
+ |
where $\alpha$ is a damping parameter, or separation constant, with |
| 143 |
+ |
units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and equal |
| 144 |
+ |
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 145 |
+ |
constant of the encompassing medium. The final two terms of |
| 146 |
+ |
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 147 |
+ |
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 148 |
+ |
dipolar term was neglected in early applications in molecular |
| 149 |
+ |
simulations,\cite{Brush66,Woodcock71} until it was introduced by de |
| 150 |
+ |
Leeuw {\it et al.} to address situations where the unit cell has a |
| 151 |
+ |
dipole moment and this dipole moment gets magnified through |
| 152 |
+ |
replication of the periodic images.\cite{deLeeuw80,Smith81} If this |
| 153 |
+ |
term is taken to be zero, the system is using conducting boundary |
| 154 |
+ |
conditions, $\epsilon_{\rm S} = \infty$. Figure |
| 155 |
+ |
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 156 |
+ |
time. Initially, due to the small size of systems, the entire |
| 157 |
+ |
simulation box was replicated to convergence. Currently, we balance a |
| 158 |
+ |
spherical real-space cutoff with the reciprocal sum and consider the |
| 159 |
+ |
surrounding dielectric. |
| 160 |
|
\begin{figure} |
| 161 |
|
\centering |
| 162 |
|
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 170 |
|
\label{fig:ewaldTime} |
| 171 |
|
\end{figure} |
| 172 |
|
|
| 173 |
+ |
The Ewald summation in the straight-forward form is an |
| 174 |
+ |
$\mathscr{O}(N^2)$ algorithm. The separation constant $(\alpha)$ |
| 175 |
+ |
plays an important role in the computational cost balance between the |
| 176 |
+ |
direct and reciprocal-space portions of the summation. The choice of |
| 177 |
+ |
the magnitude of this value allows one to select whether the |
| 178 |
+ |
real-space or reciprocal space portion of the summation is an |
| 179 |
+ |
$\mathscr{O}(N^2)$ calculation (with the other being |
| 180 |
+ |
$\mathscr{O}(N)$).\cite{Sagui99} With appropriate choice of $\alpha$ |
| 181 |
+ |
and thoughtful algorithm development, this cost can be brought down to |
| 182 |
+ |
$\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to |
| 183 |
+ |
reduce the cost of the Ewald summation further is to set $\alpha$ such |
| 184 |
+ |
that the real-space interactions decay rapidly, allowing for a short |
| 185 |
+ |
spherical cutoff, and then optimize the reciprocal space summation. |
| 186 |
+ |
These optimizations usually involve the utilization of the fast |
| 187 |
+ |
Fourier transform (FFT),\cite{Hockney81} leading to the |
| 188 |
+ |
particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) |
| 189 |
+ |
methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these |
| 190 |
+ |
methods, the cost of the reciprocal-space portion of the Ewald |
| 191 |
+ |
summation is from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N \log N)$. |
| 192 |
+ |
|
| 193 |
+ |
These developments and optimizations have led the use of the Ewald |
| 194 |
+ |
summation to become routine in simulations with periodic boundary |
| 195 |
+ |
conditions. However, in certain systems the intrinsic three |
| 196 |
+ |
dimensional periodicity can prove to be problematic, such as two |
| 197 |
+ |
dimensional surfaces and membranes. The Ewald sum has been |
| 198 |
+ |
reformulated to handle 2D |
| 199 |
+ |
systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the new |
| 200 |
+ |
methods have been found to be computationally |
| 201 |
+ |
expensive.\cite{Spohr97,Yeh99} Inclusion of a correction term in the |
| 202 |
+ |
full Ewald summation is a possible direction for enabling the handling |
| 203 |
+ |
of 2D systems and the inclusion of the optimizations described |
| 204 |
+ |
previously.\cite{Yeh99} |
| 205 |
+ |
|
| 206 |
+ |
Several studies have recognized that the inherent periodicity in the |
| 207 |
+ |
Ewald sum can also have an effect on systems that have the same |
| 208 |
+ |
dimensionality.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 209 |
+ |
Good examples are solvated proteins kept at high relative |
| 210 |
+ |
concentration due to the periodicity of the electrostatics. In these |
| 211 |
+ |
systems, the more compact folded states of a protein can be |
| 212 |
+ |
artificially stabilized by the periodic replicas introduced by the |
| 213 |
+ |
Ewald summation.\cite{Weber00} Thus, care ought to be taken when |
| 214 |
+ |
considering the use of the Ewald summation where the intrinsic |
| 215 |
+ |
periodicity may negatively affect the system dynamics. |
| 216 |
+ |
|
| 217 |
+ |
|
| 218 |
|
\subsection{The Wolf and Zahn Methods} |
| 219 |
|
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
| 220 |
|
for the accurate accumulation of electrostatic interactions in an |
| 221 |
< |
efficient pairwise fashion.\cite{Wolf99} Wolf \textit{et al.} observed |
| 222 |
< |
that the electrostatic interaction is effectively short-ranged in |
| 223 |
< |
condensed phase systems and that neutralization of the charge |
| 224 |
< |
contained within the cutoff radius is crucial for potential |
| 225 |
< |
stability. They devised a pairwise summation method that ensures |
| 226 |
< |
charge neutrality and gives results similar to those obtained with |
| 227 |
< |
the Ewald summation. The resulting shifted Coulomb potential |
| 228 |
< |
(Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through |
| 229 |
< |
placement on the cutoff sphere and a distance-dependent damping |
| 230 |
< |
function (identical to that seen in the real-space portion of the |
| 231 |
< |
Ewald sum) to aid convergence |
| 221 |
> |
efficient pairwise fashion and lacks the inherent periodicity of the |
| 222 |
> |
Ewald summation.\cite{Wolf99} Wolf \textit{et al.} observed that the |
| 223 |
> |
electrostatic interaction is effectively short-ranged in condensed |
| 224 |
> |
phase systems and that neutralization of the charge contained within |
| 225 |
> |
the cutoff radius is crucial for potential stability. They devised a |
| 226 |
> |
pairwise summation method that ensures charge neutrality and gives |
| 227 |
> |
results similar to those obtained with the Ewald summation. The |
| 228 |
> |
resulting shifted Coulomb potential (Eq. \ref{eq:WolfPot}) includes |
| 229 |
> |
image-charges subtracted out through placement on the cutoff sphere |
| 230 |
> |
and a distance-dependent damping function (identical to that seen in |
| 231 |
> |
the real-space portion of the Ewald sum) to aid convergence |
| 232 |
|
\begin{equation} |
| 233 |
< |
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\}. |
| 233 |
> |
V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. |
| 234 |
|
\label{eq:WolfPot} |
| 235 |
|
\end{equation} |
| 236 |
|
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
| 245 |
|
derivative of this potential prior to evaluation of the limit. This |
| 246 |
|
procedure gives an expression for the forces, |
| 247 |
|
\begin{equation} |
| 248 |
< |
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\}, |
| 248 |
> |
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\}, |
| 249 |
|
\label{eq:WolfForces} |
| 250 |
|
\end{equation} |
| 251 |
|
that incorporates both image charges and damping of the electrostatic |
| 325 |
|
The forces associated with the shifted potential are simply the forces |
| 326 |
|
of the unshifted potential itself (when inside the cutoff sphere), |
| 327 |
|
\begin{equation} |
| 328 |
< |
F_{\textrm{SP}} = \left( \frac{d v(r)}{dr} \right), |
| 328 |
> |
f_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), |
| 329 |
|
\end{equation} |
| 330 |
|
and are zero outside. Inside the cutoff sphere, the forces associated |
| 331 |
|
with the shifted force form can be written, |
| 332 |
|
\begin{equation} |
| 333 |
< |
F_{\textrm{SF}} = \left( \frac{d v(r)}{dr} \right) - \left(\frac{d |
| 333 |
> |
f_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d |
| 334 |
|
v(r)}{dr} \right)_{r=R_\textrm{c}}. |
| 335 |
|
\end{equation} |
| 336 |
|
|
| 342 |
|
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
| 343 |
|
al.}'s undamped prescription: |
| 344 |
|
\begin{equation} |
| 345 |
< |
V_\textrm{SP}(r) = |
| 345 |
> |
v_\textrm{SP}(r) = |
| 346 |
|
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
| 347 |
|
r\leqslant R_\textrm{c}, |
| 348 |
< |
\label{eq:WolfSP} |
| 348 |
> |
\label{eq:SPPot} |
| 349 |
|
\end{equation} |
| 350 |
|
with associated forces, |
| 351 |
|
\begin{equation} |
| 352 |
< |
F_\textrm{SP}(r) = q_iq_j\left(-\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
| 353 |
< |
\label{eq:FWolfSP} |
| 352 |
> |
f_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
| 353 |
> |
\label{eq:SPForces} |
| 354 |
|
\end{equation} |
| 355 |
|
These forces are identical to the forces of the standard Coulomb |
| 356 |
|
interaction, and cutting these off at $R_c$ was addressed by Wolf |
| 364 |
|
The shifted force ({\sc sf}) form using the normal Coulomb potential |
| 365 |
|
will give, |
| 366 |
|
\begin{equation} |
| 367 |
< |
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}. |
| 367 |
> |
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}. |
| 368 |
|
\label{eq:SFPot} |
| 369 |
|
\end{equation} |
| 370 |
|
with associated forces, |
| 371 |
|
\begin{equation} |
| 372 |
< |
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}. |
| 372 |
> |
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}. |
| 373 |
|
\label{eq:SFForces} |
| 374 |
|
\end{equation} |
| 375 |
|
This formulation has the benefits that there are no discontinuities at |
| 383 |
|
to gain functionality in dynamics simulations. |
| 384 |
|
|
| 385 |
|
Wolf \textit{et al.} originally discussed the energetics of the |
| 386 |
< |
shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that |
| 386 |
> |
shifted Coulomb potential (Eq. \ref{eq:SPPot}), and they found that |
| 387 |
|
it was still insufficient for accurate determination of the energy |
| 388 |
|
with reasonable cutoff distances. The calculated Madelung energies |
| 389 |
|
fluctuate around the expected value with increasing cutoff radius, but |
| 398 |
|
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
| 399 |
|
\label{eq:dampCoulomb} |
| 400 |
|
\end{equation} |
| 401 |
< |
the shifted potential (Eq. (\ref{eq:WolfSP})) can be recovered |
| 402 |
< |
using eq. (\ref{eq:shiftingForm}), |
| 401 |
> |
the shifted potential (Eq. (\ref{eq:SPPot})) can be reacquired using |
| 402 |
> |
eq. (\ref{eq:shiftingForm}), |
| 403 |
|
\begin{equation} |
| 404 |
< |
v_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}(\alpha r)}{r}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c}, |
| 404 |
> |
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}, |
| 405 |
|
\label{eq:DSPPot} |
| 406 |
|
\end{equation} |
| 407 |
|
with associated forces, |
| 408 |
|
\begin{equation} |
| 409 |
< |
f_{\textrm{DSP}}(r) = q_iq_j\left(-\frac{\textrm{erfc}(\alpha r)}{r^2}-\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r^2)}}{r}\right) \quad r\leqslant R_\textrm{c}. |
| 409 |
> |
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}. |
| 410 |
|
\label{eq:DSPForces} |
| 411 |
|
\end{equation} |
| 412 |
|
Again, this damped shifted potential suffers from a discontinuity and |
| 419 |
|
\label{eq:DSFPot} |
| 420 |
|
\end{split} |
| 421 |
|
\end{equation} |
| 422 |
< |
The derivative of the above potential gives the following forces, |
| 422 |
> |
The derivative of the above potential will lead to the following forces, |
| 423 |
|
\begin{equation} |
| 424 |
|
\begin{split} |
| 425 |
< |
f_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&-\left(\frac{\textrm{erfc}(\alpha r)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r^2)}}{r}\right) \\ &\left.+\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] \quad r\leqslant R_\textrm{c}. |
| 425 |
> |
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}. |
| 426 |
|
\label{eq:DSFForces} |
| 427 |
|
\end{split} |
| 428 |
|
\end{equation} |
| 429 |
+ |
If the damping parameter $(\alpha)$ is chosen to be zero, the undamped |
| 430 |
+ |
case, eqs. (\ref{eq:SPPot}-\ref{eq:SFForces}) are correctly recovered |
| 431 |
+ |
from eqs. (\ref{eq:DSPPot}-\ref{eq:DSFForces}). |
| 432 |
|
|
| 433 |
< |
This new {\sc sf} potential is similar to equation |
| 434 |
< |
\ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are |
| 435 |
< |
two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term |
| 436 |
< |
from eq. (\ref{eq:shiftingForm}) is equal to |
| 437 |
< |
eq. (\ref{eq:dampCoulomb}) with $r$ replaced by $R_\textrm{c}$. This |
| 438 |
< |
term is {\it not} present in the Zahn potential, resulting in a |
| 439 |
< |
potential discontinuity as particles cross $R_\textrm{c}$. Second, |
| 440 |
< |
the sign of the derivative portion is different. The missing |
| 441 |
< |
$v_\textrm{c}$ term would not affect molecular dynamics simulations |
| 442 |
< |
(although the computed energy would be expected to have sudden jumps |
| 443 |
< |
as particle distances crossed $R_c$). The sign problem would be a |
| 444 |
< |
potential source of errors, however. In fact, it introduces a |
| 445 |
< |
discontinuity in the forces at the cutoff, because the force function |
| 446 |
< |
is shifted in the wrong direction and doesn't cross zero at |
| 325 |
< |
$R_\textrm{c}$. |
| 433 |
> |
This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} |
| 434 |
> |
derived by Zahn \textit{et al.}; however, there are two important |
| 435 |
> |
differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from |
| 436 |
> |
eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb}) |
| 437 |
> |
with $r$ replaced by $R_\textrm{c}$. This term is {\it not} present |
| 438 |
> |
in the Zahn potential, resulting in a potential discontinuity as |
| 439 |
> |
particles cross $R_\textrm{c}$. Second, the sign of the derivative |
| 440 |
> |
portion is different. The missing $v_\textrm{c}$ term would not |
| 441 |
> |
affect molecular dynamics simulations (although the computed energy |
| 442 |
> |
would be expected to have sudden jumps as particle distances crossed |
| 443 |
> |
$R_c$). The sign problem would be a potential source of errors, |
| 444 |
> |
however. In fact, it introduces a discontinuity in the forces at the |
| 445 |
> |
cutoff, because the force function is shifted in the wrong direction |
| 446 |
> |
and doesn't cross zero at $R_\textrm{c}$. |
| 447 |
|
|
| 448 |
|
Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
| 449 |
|
electrostatic summation method that is continuous in both the |
| 480 |
|
|
| 481 |
|
Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$, |
| 482 |
|
and to incorporate their effect, a method like Reaction Field ({\sc |
| 483 |
< |
rf}) can be used. The orignal theory for {\sc rf} was originally |
| 483 |
> |
rf}) can be used. The original theory for {\sc rf} was originally |
| 484 |
|
developed by Onsager,\cite{Onsager36} and it was applied in |
| 485 |
|
simulations for the study of water by Barker and Watts.\cite{Barker73} |
| 486 |
|
In application, it is simply an extension of the group-based cutoff |
| 487 |
|
method where the net dipole within the cutoff sphere polarizes an |
| 488 |
|
external dielectric, which reacts back on the central dipole. The |
| 489 |
|
same switching function considerations for group-based cutoffs need to |
| 490 |
< |
made for {\sc rf}, with the additional prespecification of a |
| 490 |
> |
made for {\sc rf}, with the additional pre-specification of a |
| 491 |
|
dielectric constant. |
| 492 |
|
|
| 493 |
|
\section{Methods} |
| 573 |
|
investigated through measurement of the angle ($\theta$) formed |
| 574 |
|
between those computed from the particular method and those from SPME, |
| 575 |
|
\begin{equation} |
| 576 |
< |
\theta_f = \cos^{-1} \hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}, |
| 576 |
> |
\theta_f = \cos^{-1} \left(\hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}\right), |
| 577 |
|
\end{equation} |
| 578 |
|
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the |
| 579 |
|
force vector computed using method $M$. |
| 603 |
|
when using the reference method (SPME). |
| 604 |
|
|
| 605 |
|
\subsection{Short-time Dynamics} |
| 606 |
< |
|
| 607 |
< |
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 487 |
< |
Evaluation of the long-time dynamics of charged systems was performed |
| 488 |
< |
by considering the NaCl crystal system while using a subset of the |
| 606 |
> |
Evaluation of the short-time dynamics of charged systems was performed |
| 607 |
> |
by considering the 1000 K NaCl crystal system while using a subset of the |
| 608 |
|
best performing pairwise methods. The NaCl crystal was chosen to |
| 609 |
|
avoid possible complications involving the propagation techniques of |
| 610 |
< |
orientational motion in molecular systems. To enhance the atomic |
| 611 |
< |
motion, these crystals were equilibrated at 1000 K, near the |
| 612 |
< |
experimental $T_m$ for NaCl. Simulations were performed under the |
| 613 |
< |
microcanonical ensemble, and velocity autocorrelation functions |
| 614 |
< |
(Eq. \ref{eq:vCorr}) were computed for each of the trajectories, |
| 610 |
> |
orientational motion in molecular systems. All systems were started |
| 611 |
> |
with the same initial positions and velocities. Simulations were |
| 612 |
> |
performed under the microcanonical ensemble, and velocity |
| 613 |
> |
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
| 614 |
> |
of the trajectories, |
| 615 |
|
\begin{equation} |
| 616 |
< |
C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle. |
| 616 |
> |
C_v(t) = \frac{\langle v_i(0)\cdot v_i(t)\rangle}{\langle v_i(0)\cdot v_i(0)\rangle}. |
| 617 |
|
\label{eq:vCorr} |
| 618 |
|
\end{equation} |
| 619 |
< |
Velocity autocorrelation functions require detailed short time data |
| 620 |
< |
and long trajectories for good statistics, thus velocity information |
| 621 |
< |
was saved every 5 fs over 100 ps trajectories. The power spectrum |
| 622 |
< |
($I(\omega)$) is obtained via Fourier transform of the autocorrelation |
| 623 |
< |
function |
| 619 |
> |
Velocity autocorrelation functions require detailed short time data, |
| 620 |
> |
thus velocity information was saved every 2 fs over 10 ps |
| 621 |
> |
trajectories. Because the NaCl crystal is composed of two different |
| 622 |
> |
atom types, the average of the two resulting velocity autocorrelation |
| 623 |
> |
functions was used for comparisons. |
| 624 |
> |
|
| 625 |
> |
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 626 |
> |
Evaluation of the long-time dynamics of charged systems was performed |
| 627 |
> |
by considering the NaCl crystal system, again while using a subset of |
| 628 |
> |
the best performing pairwise methods. To enhance the atomic motion, |
| 629 |
> |
these crystals were equilibrated at 1000 K, near the experimental |
| 630 |
> |
$T_m$ for NaCl. Simulations were performed under the microcanonical |
| 631 |
> |
ensemble, and velocity information was saved every 5 fs over 100 ps |
| 632 |
> |
trajectories. The power spectrum ($I(\omega)$) was obtained via |
| 633 |
> |
Fourier transform of the velocity autocorrelation function |
| 634 |
|
\begin{equation} |
| 635 |
|
I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 636 |
|
\label{eq:powerSpec} |
| 637 |
|
\end{equation} |
| 638 |
< |
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. |
| 638 |
> |
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the |
| 639 |
> |
NaCl crystal is composed of two different atom types, the average of |
| 640 |
> |
the two resulting power spectra was used for comparisons. |
| 641 |
|
|
| 642 |
|
\subsection{Representative Simulations}\label{sec:RepSims} |
| 643 |
|
A variety of common and representative simulations were analyzed to |
| 661 |
|
Generation of the system configurations was dependent on the system |
| 662 |
|
type. For the solid and liquid water configurations, configuration |
| 663 |
|
snapshots were taken at regular intervals from higher temperature 1000 |
| 664 |
< |
SPC/E water molecule trajectories and each equilibrated individually. |
| 665 |
< |
The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- |
| 666 |
< |
ions and were selected and equilibrated in the same fashion as the |
| 667 |
< |
water systems. For the low and high ionic strength NaCl solutions, 4 |
| 668 |
< |
and 40 ions were first solvated in a 1000 water molecule boxes |
| 669 |
< |
respectively. Ion and water positions were then randomly swapped, and |
| 670 |
< |
the resulting configurations were again equilibrated individually. |
| 671 |
< |
Finally, for the Argon/Water "charge void" systems, the identities of |
| 672 |
< |
all the SPC/E waters within 6 \AA\ of the center of the equilibrated |
| 673 |
< |
water configurations were converted to argon |
| 674 |
< |
(Fig. \ref{fig:argonSlice}). |
| 664 |
> |
SPC/E water molecule trajectories and each equilibrated |
| 665 |
> |
individually.\cite{Berendsen87} The solid and liquid NaCl systems |
| 666 |
> |
consisted of 500 Na+ and 500 Cl- ions and were selected and |
| 667 |
> |
equilibrated in the same fashion as the water systems. For the low |
| 668 |
> |
and high ionic strength NaCl solutions, 4 and 40 ions were first |
| 669 |
> |
solvated in a 1000 water molecule boxes respectively. Ion and water |
| 670 |
> |
positions were then randomly swapped, and the resulting configurations |
| 671 |
> |
were again equilibrated individually. Finally, for the Argon/Water |
| 672 |
> |
"charge void" systems, the identities of all the SPC/E waters within 6 |
| 673 |
> |
\AA\ of the center of the equilibrated water configurations were |
| 674 |
> |
converted to argon (Fig. \ref{fig:argonSlice}). |
| 675 |
|
|
| 676 |
|
\begin{figure} |
| 677 |
|
\centering |
| 701 |
|
the energies and forces calculated. Typical molecular mechanics |
| 702 |
|
packages default this to a value dependent on the cutoff radius and a |
| 703 |
|
tolerance (typically less than $1 \times 10^{-4}$ kcal/mol). Smaller |
| 704 |
< |
tolerances are typically associated with increased accuracy in the |
| 705 |
< |
real-space portion of the summation.\cite{Essmann95} The default |
| 706 |
< |
TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME |
| 704 |
> |
tolerances are typically associated with increased accuracy, but this |
| 705 |
> |
usually means more time spent calculating the reciprocal-space portion |
| 706 |
> |
of the summation.\cite{Perram88,Essmann95} The default TINKER |
| 707 |
> |
tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME |
| 708 |
|
calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and |
| 709 |
|
0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 710 |
|
|
| 726 |
|
|
| 727 |
|
In this figure, it is apparent that it is unreasonable to expect |
| 728 |
|
realistic results using an unmodified cutoff. This is not all that |
| 729 |
< |
surprising since this results in large energy fluctuations as atoms |
| 730 |
< |
move in and out of the cutoff radius. These fluctuations can be |
| 731 |
< |
alleviated to some degree by using group based cutoffs with a |
| 732 |
< |
switching function.\cite{Steinbach94} The Group Switch Cutoff row |
| 733 |
< |
doesn't show a significant improvement in this plot because the salt |
| 734 |
< |
and salt solution systems contain non-neutral groups, see the |
| 729 |
> |
surprising since this results in large energy fluctuations as atoms or |
| 730 |
> |
molecules move in and out of the cutoff radius.\cite{Rahman71,Adams79} |
| 731 |
> |
These fluctuations can be alleviated to some degree by using group |
| 732 |
> |
based cutoffs with a switching |
| 733 |
> |
function.\cite{Adams79,Steinbach94,Leach01} The Group Switch Cutoff |
| 734 |
> |
row doesn't show a significant improvement in this plot because the |
| 735 |
> |
salt and salt solution systems contain non-neutral groups, see the |
| 736 |
|
accompanying supporting information for a comparison where all groups |
| 737 |
|
are neutral. |
| 738 |
|
|
| 739 |
|
Correcting the resulting charged cutoff sphere is one of the purposes |
| 740 |
|
of the damped Coulomb summation proposed by Wolf \textit{et |
| 741 |
|
al.},\cite{Wolf99} and this correction indeed improves the results as |
| 742 |
< |
seen in the Shifted-Potental rows. While the undamped case of this |
| 742 |
> |
seen in the {\sc sp} rows. While the undamped case of this |
| 743 |
|
method is a significant improvement over the pure cutoff, it still |
| 744 |
|
doesn't correlate that well with SPME. Inclusion of potential damping |
| 745 |
|
improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 924 |
|
up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably |
| 925 |
|
unnecessary when using the {\sc sf} method. |
| 926 |
|
|
| 927 |
< |
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 927 |
> |
\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals} |
| 928 |
|
|
| 929 |
|
In the previous studies using a {\sc sf} variant of the damped |
| 930 |
|
Wolf coulomb potential, the structure and dynamics of water were |
| 936 |
|
systems and simply recapitulate their results, we decided to look at |
| 937 |
|
the solid state dynamical behavior obtained using the best performing |
| 938 |
|
summation methods from the above results. |
| 939 |
+ |
|
| 940 |
+ |
\begin{figure} |
| 941 |
+ |
\centering |
| 942 |
+ |
\includegraphics[width = \linewidth]{./vCorrPlot.pdf} |
| 943 |
+ |
\caption{Velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset is a magnification of the first trough. The times to first collision are nearly identical, but the differences can be seen in the peaks and troughs, where the undamped to weakly damped methods are stiffer than the moderately damped and SPME methods.} |
| 944 |
+ |
\label{fig:vCorrPlot} |
| 945 |
+ |
\end{figure} |
| 946 |
+ |
|
| 947 |
+ |
The short-time decays through the first collision are nearly identical |
| 948 |
+ |
in figure \ref{fig:vCorrPlot}, but the peaks and troughs of the |
| 949 |
+ |
functions show how the methods differ. The undamped {\sc sf} method |
| 950 |
+ |
has deeper troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher |
| 951 |
+ |
peaks than any of the other methods. As the damping function is |
| 952 |
+ |
increased, these peaks are smoothed out, and approach the SPME |
| 953 |
+ |
curve. The damping acts as a distance dependent Gaussian screening of |
| 954 |
+ |
the point charges for the pairwise summation methods; thus, the |
| 955 |
+ |
collisions are more elastic in the undamped {\sc sf} potential, and the |
| 956 |
+ |
stiffness of the potential is diminished as the electrostatic |
| 957 |
+ |
interactions are softened by the damping function. With $\alpha$ |
| 958 |
+ |
values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are |
| 959 |
+ |
nearly identical and track the SPME features quite well. This is not |
| 960 |
+ |
too surprising in that the differences between the {\sc sf} and {\sc |
| 961 |
+ |
sp} potentials are mitigated with increased damping. However, this |
| 962 |
+ |
appears to indicate that once damping is utilized, the form of the |
| 963 |
+ |
potential seems to play a lesser role in the crystal dynamics. |
| 964 |
+ |
|
| 965 |
+ |
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 966 |
+ |
|
| 967 |
+ |
The short time dynamics were extended to evaluate how the differences |
| 968 |
+ |
between the methods affect the collective long-time motion. The same |
| 969 |
+ |
electrostatic summation methods were used as in the short time |
| 970 |
+ |
velocity autocorrelation function evaluation, but the trajectories |
| 971 |
+ |
were sampled over a much longer time. The power spectra of the |
| 972 |
+ |
resulting velocity autocorrelation functions were calculated and are |
| 973 |
+ |
displayed in figure \ref{fig:methodPS}. |
| 974 |
|
|
| 975 |
|
\begin{figure} |
| 976 |
|
\centering |
| 979 |
|
\label{fig:methodPS} |
| 980 |
|
\end{figure} |
| 981 |
|
|
| 982 |
< |
Figure \ref{fig:methodPS} shows the power spectra for the NaCl |
| 983 |
< |
crystals (from averaged Na and Cl ion velocity autocorrelation |
| 984 |
< |
functions) using the stated electrostatic summation methods. While |
| 985 |
< |
high frequency peaks of all the spectra overlap, showing the same |
| 986 |
< |
general features, the low frequency region shows how the summation |
| 987 |
< |
methods differ. Considering the low-frequency inset (expanded in the |
| 988 |
< |
upper frame of figure \ref{fig:dampInc}), at frequencies below 100 |
| 989 |
< |
cm$^{-1}$, the correlated motions are blue-shifted when using undamped |
| 990 |
< |
or weakly damped {\sc sf}. When using moderate damping ($\alpha |
| 991 |
< |
= 0.2$ \AA$^{-1}$) both the {\sc sf} and {\sc sp} |
| 992 |
< |
methods give near identical correlated motion behavior as the Ewald |
| 993 |
< |
method (which has a damping value of 0.3119). The damping acts as a |
| 994 |
< |
distance dependent Gaussian screening of the point charges for the |
| 995 |
< |
pairwise summation methods. This weakening of the electrostatic |
| 996 |
< |
interaction with distance explains why the long-ranged correlated |
| 829 |
< |
motions are at lower frequencies for the moderately damped methods |
| 830 |
< |
than for undamped or weakly damped methods. To see this effect more |
| 831 |
< |
clearly, we show how damping strength affects a simple real-space |
| 832 |
< |
electrostatic potential, |
| 982 |
> |
While high frequency peaks of the spectra in this figure overlap, |
| 983 |
> |
showing the same general features, the low frequency region shows how |
| 984 |
> |
the summation methods differ. Considering the low-frequency inset |
| 985 |
> |
(expanded in the upper frame of figure \ref{fig:dampInc}), at |
| 986 |
> |
frequencies below 100 cm$^{-1}$, the correlated motions are |
| 987 |
> |
blue-shifted when using undamped or weakly damped {\sc sf}. When |
| 988 |
> |
using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the {\sc sf} |
| 989 |
> |
and {\sc sp} methods give near identical correlated motion behavior as |
| 990 |
> |
the Ewald method (which has a damping value of 0.3119). This |
| 991 |
> |
weakening of the electrostatic interaction with increased damping |
| 992 |
> |
explains why the long-ranged correlated motions are at lower |
| 993 |
> |
frequencies for the moderately damped methods than for undamped or |
| 994 |
> |
weakly damped methods. To see this effect more clearly, we show how |
| 995 |
> |
damping strength alone affects a simple real-space electrostatic |
| 996 |
> |
potential, |
| 997 |
|
\begin{equation} |
| 998 |
|
V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r), |
| 999 |
|
\end{equation} |
| 1008 |
|
shift to higher frequency in exponential fashion. Though not shown, |
| 1009 |
|
the spectrum for the simple undamped electrostatic potential is |
| 1010 |
|
blue-shifted such that the lowest frequency peak resides near 325 |
| 1011 |
< |
cm$^{-1}$. In light of these results, the undamped {\sc sf} |
| 1012 |
< |
method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is |
| 1013 |
< |
quite respectable; however, it appears as though moderate damping is |
| 1014 |
< |
required for accurate reproduction of crystal dynamics. |
| 1011 |
> |
cm$^{-1}$. In light of these results, the undamped {\sc sf} method |
| 1012 |
> |
producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite |
| 1013 |
> |
respectable and shows that the shifted force procedure accounts for |
| 1014 |
> |
most of the effect afforded through use of the Ewald summation. |
| 1015 |
> |
However, it appears as though moderate damping is required for |
| 1016 |
> |
accurate reproduction of crystal dynamics. |
| 1017 |
|
\begin{figure} |
| 1018 |
|
\centering |
| 1019 |
|
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 1020 |
< |
\caption{Upper: Zoomed inset from figure \ref{fig:methodPS}. As the damping value for the {\sc sf} 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.} |
| 1020 |
> |
\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.} |
| 1021 |
|
\label{fig:dampInc} |
| 1022 |
|
\end{figure} |
| 1023 |
|
|
| 1057 |
|
standard by which these simple pairwise sums are judged. However, |
| 1058 |
|
these results do suggest that in the typical simulations performed |
| 1059 |
|
today, the Ewald summation may no longer be required to obtain the |
| 1060 |
< |
level of accuracy most researcher have come to expect |
| 1060 |
> |
level of accuracy most researchers have come to expect |
| 1061 |
|
|
| 1062 |
|
\section{Acknowledgments} |
| 1063 |
|
\newpage |
