| 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} |
| 81 |
|
impractical task to perform these calculations. |
| 82 |
|
|
| 83 |
|
\subsection{The Ewald Sum} |
| 84 |
< |
blah blah blah Ewald Sum Important blah blah blah |
| 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 |
+ |
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} |
| 137 |
|
\label{fig:ewaldTime} |
| 138 |
|
\end{figure} |
| 139 |
|
|
| 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 |
| 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}(\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\}, |
| 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 |
| 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), |
| 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 |
| 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 |
|
|
| 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) = |
| 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:WolfSP} |
| 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:FWolfSP} |
| 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 |
| 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}. |
| 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}. |
| 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 |
| 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:WolfSP}), and they found that |
| 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 |
| 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:WolfSP})) can be recovered |
| 336 |
< |
using eq. (\ref{eq:shiftingForm}), |
| 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}(\alpha r)}{r}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c}, |
| 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}(\alpha r)}{r^2}-\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r^2)}}{r}\right) \quad r\leqslant R_\textrm{c}. |
| 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 |
| 353 |
|
\label{eq:DSFPot} |
| 354 |
|
\end{split} |
| 355 |
|
\end{equation} |
| 356 |
< |
The derivative of the above potential gives the following forces, |
| 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}(\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}. |
| 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 |
| 368 |
< |
\ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are |
| 369 |
< |
two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term |
| 370 |
< |
from eq. (\ref{eq:shiftingForm}) is equal to |
| 371 |
< |
eq. (\ref{eq:dampCoulomb}) with $r$ replaced by $R_\textrm{c}$. This |
| 372 |
< |
term is {\it not} present in the Zahn potential, resulting in a |
| 373 |
< |
potential discontinuity as particles cross $R_\textrm{c}$. Second, |
| 374 |
< |
the sign of the derivative portion is different. The missing |
| 375 |
< |
$v_\textrm{c}$ term would not affect molecular dynamics simulations |
| 376 |
< |
(although the computed energy would be expected to have sudden jumps |
| 377 |
< |
as particle distances crossed $R_c$). The sign problem would be a |
| 378 |
< |
potential source of errors, however. In fact, it introduces a |
| 379 |
< |
discontinuity in the forces at the cutoff, because the force function |
| 380 |
< |
is shifted in the wrong direction and doesn't cross zero at |
| 325 |
< |
$R_\textrm{c}$. |
| 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 |
| 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 in the |
| 629 |
< |
real-space portion of the summation.\cite{Essmann95} The default |
| 630 |
< |
TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME |
| 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 |
|
|
| 907 |
|
\begin{figure} |
| 908 |
|
\centering |
| 909 |
|
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 910 |
< |
\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.} |
| 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 |
|
|
