| 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}{(2\epsilon_\textrm{S}+1)L^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}$, $\mathbf{k}$ are the reciprocal vectors and equal |
| 113 |
+ |
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 114 |
+ |
constant of the encompassing medium. The final two terms of |
| 115 |
+ |
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 116 |
+ |
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 117 |
+ |
dipolar term was neglected in early applications in molecular |
| 118 |
+ |
simulations,\cite{Brush66,Woodcock71} until it was introduced by de |
| 119 |
+ |
Leeuw {\it et al.} to address situations where the unit cell has a |
| 120 |
+ |
dipole moment and this dipole moment gets magnified through |
| 121 |
+ |
replication of the periodic images.\cite{deLeeuw80,Smith81} If this |
| 122 |
+ |
term is taken to be zero, the system is using conducting boundary |
| 123 |
+ |
conditions, $\epsilon_{\rm S} = \infty$. Figure |
| 124 |
+ |
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 125 |
+ |
time. Initially, due to the small size of systems, the entire |
| 126 |
+ |
simulation box was replicated to convergence. Currently, we balance a |
| 127 |
+ |
spherical real-space cutoff with the reciprocal sum and consider the |
| 128 |
+ |
surrounding dielectric. |
| 129 |
|
\begin{figure} |
| 130 |
|
\centering |
| 131 |
|
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 138 |
|
a surrounding dielectric term is included.} |
| 139 |
|
\label{fig:ewaldTime} |
| 140 |
|
\end{figure} |
| 141 |
+ |
|
| 142 |
+ |
The Ewald summation in the straight-forward form is an |
| 143 |
+ |
$\mathscr{O}(N^2)$ algorithm. The separation constant $(\alpha)$ |
| 144 |
+ |
plays an important role in the computational cost balance between the |
| 145 |
+ |
direct and reciprocal-space portions of the summation. The choice of |
| 146 |
+ |
the magnitude of this value allows one to select whether the |
| 147 |
+ |
real-space or reciprocal space portion of the summation is an |
| 148 |
+ |
$\mathscr{O}(N^2)$ calcualtion (with the other being |
| 149 |
+ |
$\mathscr{O}(N)$).\cite{Sagui99} With appropriate choice of $\alpha$ |
| 150 |
+ |
and thoughtful algorithm development, this cost can be brought down to |
| 151 |
+ |
$\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route taken to |
| 152 |
+ |
reduce the cost of the Ewald summation further is to set $\alpha$ such |
| 153 |
+ |
that the real-space interactions decay rapidly, allowing for a short |
| 154 |
+ |
spherical cutoff, and then optimize the reciprocal space summation. |
| 155 |
+ |
These optimizations usually involve the utilization of the fast |
| 156 |
+ |
Fourier transform (FFT),\cite{Hockney81} leading to the |
| 157 |
+ |
particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) |
| 158 |
+ |
methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these |
| 159 |
+ |
methods, the cost of the reciprocal-space portion of the Ewald |
| 160 |
+ |
summation is from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N \log N)$. |
| 161 |
+ |
|
| 162 |
+ |
These developments and optimizations have led the use of the Ewald |
| 163 |
+ |
summation to become routine in simulations with periodic boundary |
| 164 |
+ |
conditions. However, in certain systems the intrinsic three |
| 165 |
+ |
dimensional periodicity can prove to be problematic, such as two |
| 166 |
+ |
dimensional surfaces and membranes. The Ewald sum has been |
| 167 |
+ |
reformulated to handle 2D |
| 168 |
+ |
systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the new |
| 169 |
+ |
methods have been found to be computationally |
| 170 |
+ |
expensive.\cite{Spohr97,Yeh99} Inclusion of a correction term in the |
| 171 |
+ |
full Ewald summation is a possible direction for enabling the handling |
| 172 |
+ |
of 2D systems and the inclusion of the optimizations described |
| 173 |
+ |
previously.\cite{Yeh99} |
| 174 |
|
|
| 175 |
+ |
Several studies have recognized that the inherent periodicity in the |
| 176 |
+ |
Ewald sum can also have an effect on systems that have the same |
| 177 |
+ |
dimensionality.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
| 178 |
+ |
Good examples are solvated proteins kept at high relative |
| 179 |
+ |
concentration due to the periodicity of the electrostatics. In these |
| 180 |
+ |
systems, the more compact folded states of a protein can be |
| 181 |
+ |
artificially stabilized by the periodic replicas introduced by the |
| 182 |
+ |
Ewald summation.\cite{Weber00} Thus, care ought to be taken when |
| 183 |
+ |
considering the use of the Ewald summation where the intrinsic |
| 184 |
+ |
perodicity may negatively affect the system dynamics. |
| 185 |
+ |
|
| 186 |
+ |
|
| 187 |
|
\subsection{The Wolf and Zahn Methods} |
| 188 |
|
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
| 189 |
|
for the accurate accumulation of electrostatic interactions in an |
| 190 |
< |
efficient pairwise fashion.\cite{Wolf99} Wolf \textit{et al.} observed |
| 191 |
< |
that the electrostatic interaction is effectively short-ranged in |
| 192 |
< |
condensed phase systems and that neutralization of the charge |
| 193 |
< |
contained within the cutoff radius is crucial for potential |
| 194 |
< |
stability. They devised a pairwise summation method that ensures |
| 195 |
< |
charge neutrality and gives results similar to those obtained with |
| 196 |
< |
the Ewald summation. The resulting shifted Coulomb potential |
| 197 |
< |
(Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through |
| 198 |
< |
placement on the cutoff sphere and a distance-dependent damping |
| 199 |
< |
function (identical to that seen in the real-space portion of the |
| 200 |
< |
Ewald sum) to aid convergence |
| 190 |
> |
efficient pairwise fashion and lacks the inherent periodicity of the |
| 191 |
> |
Ewald summation.\cite{Wolf99} Wolf \textit{et al.} observed that the |
| 192 |
> |
electrostatic interaction is effectively short-ranged in condensed |
| 193 |
> |
phase systems and that neutralization of the charge contained within |
| 194 |
> |
the cutoff radius is crucial for potential stability. They devised a |
| 195 |
> |
pairwise summation method that ensures charge neutrality and gives |
| 196 |
> |
results similar to those obtained with the Ewald summation. The |
| 197 |
> |
resulting shifted Coulomb potential (Eq. \ref{eq:WolfPot}) includes |
| 198 |
> |
image-charges subtracted out through placement on the cutoff sphere |
| 199 |
> |
and a distance-dependent damping function (identical to that seen in |
| 200 |
> |
the real-space portion of the Ewald sum) to aid convergence |
| 201 |
|
\begin{equation} |
| 202 |
|
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\}. |
| 203 |
|
\label{eq:WolfPot} |
| 214 |
|
derivative of this potential prior to evaluation of the limit. This |
| 215 |
|
procedure gives an expression for the forces, |
| 216 |
|
\begin{equation} |
| 217 |
< |
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\}, |
| 217 |
> |
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\}, |
| 218 |
|
\label{eq:WolfForces} |
| 219 |
|
\end{equation} |
| 220 |
|
that incorporates both image charges and damping of the electrostatic |
| 294 |
|
The forces associated with the shifted potential are simply the forces |
| 295 |
|
of the unshifted potential itself (when inside the cutoff sphere), |
| 296 |
|
\begin{equation} |
| 297 |
< |
F_{\textrm{SP}} = \left( \frac{d v(r)}{dr} \right), |
| 297 |
> |
f_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), |
| 298 |
|
\end{equation} |
| 299 |
|
and are zero outside. Inside the cutoff sphere, the forces associated |
| 300 |
|
with the shifted force form can be written, |
| 301 |
|
\begin{equation} |
| 302 |
< |
F_{\textrm{SF}} = \left( \frac{d v(r)}{dr} \right) - \left(\frac{d |
| 302 |
> |
f_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d |
| 303 |
|
v(r)}{dr} \right)_{r=R_\textrm{c}}. |
| 304 |
|
\end{equation} |
| 305 |
|
|
| 311 |
|
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
| 312 |
|
al.}'s undamped prescription: |
| 313 |
|
\begin{equation} |
| 314 |
< |
V_\textrm{SP}(r) = |
| 314 |
> |
v_\textrm{SP}(r) = |
| 315 |
|
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
| 316 |
|
r\leqslant R_\textrm{c}, |
| 317 |
< |
\label{eq:WolfSP} |
| 317 |
> |
\label{eq:SPPot} |
| 318 |
|
\end{equation} |
| 319 |
|
with associated forces, |
| 320 |
|
\begin{equation} |
| 321 |
< |
F_\textrm{SP}(r) = q_iq_j\left(-\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
| 322 |
< |
\label{eq:FWolfSP} |
| 321 |
> |
f_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
| 322 |
> |
\label{eq:SPForces} |
| 323 |
|
\end{equation} |
| 324 |
|
These forces are identical to the forces of the standard Coulomb |
| 325 |
|
interaction, and cutting these off at $R_c$ was addressed by Wolf |
| 333 |
|
The shifted force ({\sc sf}) form using the normal Coulomb potential |
| 334 |
|
will give, |
| 335 |
|
\begin{equation} |
| 336 |
< |
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}. |
| 336 |
> |
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}. |
| 337 |
|
\label{eq:SFPot} |
| 338 |
|
\end{equation} |
| 339 |
|
with associated forces, |
| 340 |
|
\begin{equation} |
| 341 |
< |
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}. |
| 341 |
> |
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}. |
| 342 |
|
\label{eq:SFForces} |
| 343 |
|
\end{equation} |
| 344 |
|
This formulation has the benefits that there are no discontinuities at |
| 352 |
|
to gain functionality in dynamics simulations. |
| 353 |
|
|
| 354 |
|
Wolf \textit{et al.} originally discussed the energetics of the |
| 355 |
< |
shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that |
| 355 |
> |
shifted Coulomb potential (Eq. \ref{eq:SPPot}), and they found that |
| 356 |
|
it was still insufficient for accurate determination of the energy |
| 357 |
|
with reasonable cutoff distances. The calculated Madelung energies |
| 358 |
|
fluctuate around the expected value with increasing cutoff radius, but |
| 367 |
|
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
| 368 |
|
\label{eq:dampCoulomb} |
| 369 |
|
\end{equation} |
| 370 |
< |
the shifted potential (Eq. \ref{eq:WolfSP}) can be recovered |
| 371 |
< |
\textit{via} equation \ref{eq:shiftingForm}, |
| 370 |
> |
the shifted potential (Eq. (\ref{eq:SPPot})) can be reacquired using |
| 371 |
> |
eq. (\ref{eq:shiftingForm}), |
| 372 |
|
\begin{equation} |
| 373 |
< |
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}. |
| 373 |
> |
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}, |
| 374 |
|
\label{eq:DSPPot} |
| 375 |
< |
\end{equation}, |
| 375 |
> |
\end{equation} |
| 376 |
|
with associated forces, |
| 377 |
|
\begin{equation} |
| 378 |
< |
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}. |
| 378 |
> |
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}. |
| 379 |
|
\label{eq:DSPForces} |
| 380 |
|
\end{equation} |
| 381 |
|
Again, this damped shifted potential suffers from a discontinuity and |
| 382 |
|
a lack of the image charges in the forces. To remedy these concerns, |
| 383 |
< |
one may derive a Shifted-Force variant by including the derivative |
| 384 |
< |
term in equation \ref{eq:shiftingForm}, |
| 383 |
> |
one may derive a {\sc sf} variant by including the derivative |
| 384 |
> |
term in eq. (\ref{eq:shiftingForm}), |
| 385 |
|
\begin{equation} |
| 386 |
|
\begin{split} |
| 387 |
|
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}. |
| 388 |
|
\label{eq:DSFPot} |
| 389 |
|
\end{split} |
| 390 |
|
\end{equation} |
| 391 |
< |
The derivative of the above potential gives the following forces, |
| 391 |
> |
The derivative of the above potential will lead to the following forces, |
| 392 |
|
\begin{equation} |
| 393 |
|
\begin{split} |
| 394 |
< |
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}. |
| 394 |
> |
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}. |
| 395 |
|
\label{eq:DSFForces} |
| 396 |
|
\end{split} |
| 397 |
|
\end{equation} |
| 398 |
+ |
If the damping parameter $(\alpha)$ is chosen to be zero, the undamped |
| 399 |
+ |
case, eqs. (\ref{eq:SPPot}-\ref{eq:SFForces}) are correctly recovered |
| 400 |
+ |
from eqs. (\ref{eq:DSPPot}-\ref{eq:DSFForces}). |
| 401 |
|
|
| 402 |
< |
This new Shifted-Force potential is similar to equation |
| 403 |
< |
\ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are |
| 404 |
< |
two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term |
| 405 |
< |
from eq. (\ref{eq:shiftingForm}) is equal to |
| 406 |
< |
eq. (\ref{eq:dampCoulomb}) with $r$ replaced by $R_\textrm{c}$. This |
| 407 |
< |
term is {\it not} present in the Zahn potential, resulting in a |
| 408 |
< |
potential discontinuity as particles cross $R_\textrm{c}$. Second, |
| 409 |
< |
the sign of the derivative portion is different. The missing |
| 410 |
< |
$v_\textrm{c}$ term would not affect molecular dynamics simulations |
| 411 |
< |
(although the computed energy would be expected to have sudden jumps |
| 412 |
< |
as particle distances crossed $R_c$). The sign problem would be a |
| 413 |
< |
potential source of errors, however. In fact, it introduces a |
| 414 |
< |
discontinuity in the forces at the cutoff, because the force function |
| 415 |
< |
is shifted in the wrong direction and doesn't cross zero at |
| 325 |
< |
$R_\textrm{c}$. |
| 402 |
> |
This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} |
| 403 |
> |
derived by Zahn \textit{et al.}; however, there are two important |
| 404 |
> |
differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from |
| 405 |
> |
eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb}) |
| 406 |
> |
with $r$ replaced by $R_\textrm{c}$. This term is {\it not} present |
| 407 |
> |
in the Zahn potential, resulting in a potential discontinuity as |
| 408 |
> |
particles cross $R_\textrm{c}$. Second, the sign of the derivative |
| 409 |
> |
portion is different. The missing $v_\textrm{c}$ term would not |
| 410 |
> |
affect molecular dynamics simulations (although the computed energy |
| 411 |
> |
would be expected to have sudden jumps as particle distances crossed |
| 412 |
> |
$R_c$). The sign problem would be a potential source of errors, |
| 413 |
> |
however. In fact, it introduces a discontinuity in the forces at the |
| 414 |
> |
cutoff, because the force function is shifted in the wrong direction |
| 415 |
> |
and doesn't cross zero at $R_\textrm{c}$. |
| 416 |
|
|
| 417 |
|
Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
| 418 |
|
electrostatic summation method that is continuous in both the |
| 423 |
|
performed by the Ewald sum. |
| 424 |
|
|
| 425 |
|
\subsection{Other alternatives} |
| 426 |
+ |
In addition to the methods described above, we will consider some |
| 427 |
+ |
other techniques that commonly get used in molecular simulations. The |
| 428 |
+ |
simplest of these is group-based cutoffs. Though of little use for |
| 429 |
+ |
non-neutral molecules, collecting atoms into neutral groups takes |
| 430 |
+ |
advantage of the observation that the electrostatic interactions decay |
| 431 |
+ |
faster than those for monopolar pairs.\cite{Steinbach94} When |
| 432 |
+ |
considering these molecules as groups, an orientational aspect is |
| 433 |
+ |
introduced to the interactions. Consequently, as these molecular |
| 434 |
+ |
particles move through $R_\textrm{c}$, the energy will drift upward |
| 435 |
+ |
due to the anisotropy of the net molecular dipole |
| 436 |
+ |
interactions.\cite{Rahman71} To maintain good energy conservation, |
| 437 |
+ |
both the potential and derivative need to be smoothly switched to zero |
| 438 |
+ |
at $R_\textrm{c}$.\cite{Adams79} This is accomplished using a |
| 439 |
+ |
switching function, |
| 440 |
+ |
\begin{equation} |
| 441 |
+ |
S(r) = \begin{cases} 1 &\quad r\leqslant r_\textrm{sw} \\ |
| 442 |
+ |
\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} \\ |
| 443 |
+ |
0 &\quad r>R_\textrm{c} |
| 444 |
+ |
\end{cases}, |
| 445 |
+ |
\end{equation} |
| 446 |
+ |
where the above form is for a cubic function. If a smooth second |
| 447 |
+ |
derivative is desired, a fifth (or higher) order polynomial can be |
| 448 |
+ |
used.\cite{Andrea83} |
| 449 |
|
|
| 450 |
< |
Reaction Field blah |
| 451 |
< |
|
| 452 |
< |
Group-based cutoff blah |
| 453 |
< |
|
| 450 |
> |
Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$, |
| 451 |
> |
and to incorporate their effect, a method like Reaction Field ({\sc |
| 452 |
> |
rf}) can be used. The original theory for {\sc rf} was originally |
| 453 |
> |
developed by Onsager,\cite{Onsager36} and it was applied in |
| 454 |
> |
simulations for the study of water by Barker and Watts.\cite{Barker73} |
| 455 |
> |
In application, it is simply an extension of the group-based cutoff |
| 456 |
> |
method where the net dipole within the cutoff sphere polarizes an |
| 457 |
> |
external dielectric, which reacts back on the central dipole. The |
| 458 |
> |
same switching function considerations for group-based cutoffs need to |
| 459 |
> |
made for {\sc rf}, with the additional pre-specification of a |
| 460 |
> |
dielectric constant. |
| 461 |
|
|
| 462 |
|
\section{Methods} |
| 463 |
|
|
| 641 |
|
|
| 642 |
|
\subsection{Electrostatic Summation Methods}\label{sec:ESMethods} |
| 643 |
|
Electrostatic summation method comparisons were performed using SPME, |
| 644 |
< |
the Shifted-Potential and Shifted-Force methods - both with damping |
| 644 |
> |
the {\sc sp} and {\sc sf} methods - both with damping |
| 645 |
|
parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak, |
| 646 |
|
moderate, and strong damping respectively), reaction field with an |
| 647 |
|
infinite dielectric constant, and an unmodified cutoff. Group-based |
| 660 |
|
the energies and forces calculated. Typical molecular mechanics |
| 661 |
|
packages default this to a value dependent on the cutoff radius and a |
| 662 |
|
tolerance (typically less than $1 \times 10^{-4}$ kcal/mol). Smaller |
| 663 |
< |
tolerances are typically associated with increased accuracy in the |
| 664 |
< |
real-space portion of the summation.\cite{Essmann95} The default |
| 665 |
< |
TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME |
| 663 |
> |
tolerances are typically associated with increased accuracy, but this |
| 664 |
> |
usually means more time spent calculating the reciprocal-space portion |
| 665 |
> |
of the summation.\cite{Perram88,Essmann95} The default TINKER |
| 666 |
> |
tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME |
| 667 |
|
calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and |
| 668 |
|
0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 669 |
|
|
| 705 |
|
particularly with a cutoff radius greater than 12 \AA . Use of a |
| 706 |
|
larger damping parameter is more helpful for the shortest cutoff |
| 707 |
|
shown, but it has a detrimental effect on simulations with larger |
| 708 |
< |
cutoffs. In the Shifted-Force sets, increasing damping results in |
| 708 |
> |
cutoffs. In the {\sc sf} sets, increasing damping results in |
| 709 |
|
progressively poorer correlation. Overall, the undamped case is the |
| 710 |
|
best performing set, as the correlation and quality of fits are |
| 711 |
|
consistently superior regardless of the cutoff distance. This result |
| 738 |
|
in the previous $\Delta E$ section. The unmodified cutoff results are |
| 739 |
|
poor, but using group based cutoffs and a switching function provides |
| 740 |
|
a improvement much more significant than what was seen with $\Delta |
| 741 |
< |
E$. Looking at the Shifted-Potential sets, the slope and $R^2$ |
| 741 |
> |
E$. Looking at the {\sc sp} sets, the slope and $R^2$ |
| 742 |
|
improve with the use of damping to an optimal result of 0.2 \AA |
| 743 |
|
$^{-1}$ for the 12 and 15 \AA\ cutoffs. Further increases in damping, |
| 744 |
|
while beneficial for simulations with a cutoff radius of 9 \AA\ , is |
| 745 |
|
detrimental to simulations with larger cutoff radii. The undamped |
| 746 |
< |
Shifted-Force method gives forces in line with those obtained using |
| 746 |
> |
{\sc sf} method gives forces in line with those obtained using |
| 747 |
|
SPME, and use of a damping function results in minor improvement. The |
| 748 |
|
reaction field results are surprisingly good, considering the poor |
| 749 |
|
quality of the fits for the $\Delta E$ results. There is still a |
| 766 |
|
torque vector magnitude results in figure \ref{fig:trqMag} are still |
| 767 |
|
similar to those seen for the forces; however, they more clearly show |
| 768 |
|
the improved behavior that comes with increasing the cutoff radius. |
| 769 |
< |
Moderate damping is beneficial to the Shifted-Potential and helpful |
| 770 |
< |
yet possibly unnecessary with the Shifted-Force method, and they also |
| 769 |
> |
Moderate damping is beneficial to the {\sc sp} and helpful |
| 770 |
> |
yet possibly unnecessary with the {\sc sf} method, and they also |
| 771 |
|
show that over-damping adversely effects all cutoff radii rather than |
| 772 |
|
showing an improvement for systems with short cutoffs. The reaction |
| 773 |
|
field method performs well when calculating the torques, better than |
| 796 |
|
show the improvement afforded by choosing a longer simulation cutoff. |
| 797 |
|
Increasing the cutoff from 9 to 12 \AA\ typically results in a halving |
| 798 |
|
of the distribution widths, with a similar improvement going from 12 |
| 799 |
< |
to 15 \AA . The undamped Shifted-Force, Group Based Cutoff, and |
| 799 |
> |
to 15 \AA . The undamped {\sc sf}, Group Based Cutoff, and |
| 800 |
|
Reaction Field methods all do equivalently well at capturing the |
| 801 |
|
direction of both the force and torque vectors. Using damping |
| 802 |
< |
improves the angular behavior significantly for the Shifted-Potential |
| 803 |
< |
and moderately for the Shifted-Force methods. Increasing the damping |
| 802 |
> |
improves the angular behavior significantly for the {\sc sp} |
| 803 |
> |
and moderately for the {\sc sf} methods. Increasing the damping |
| 804 |
|
too far is destructive for both methods, particularly to the torque |
| 805 |
|
vectors. Again it is important to recognize that the force vectors |
| 806 |
|
cover all particles in the systems, while torque vectors are only |
| 842 |
|
\end{table} |
| 843 |
|
|
| 844 |
|
Although not discussed previously, group based cutoffs can be applied |
| 845 |
< |
to both the Shifted-Potential and Shifted-Force methods. Use off a |
| 845 |
> |
to both the {\sc sp} and {\sc sf} methods. Use off a |
| 846 |
|
switching function corrects for the discontinuities that arise when |
| 847 |
|
atoms of a group exit the cutoff before the group's center of mass. |
| 848 |
|
Though there are no significant benefit or drawbacks observed in |
| 851 |
|
\ref{tab:groupAngle} shows the angular variance values obtained using |
| 852 |
|
group based cutoffs and a switching function alongside the standard |
| 853 |
|
results seen in figure \ref{fig:frcTrqAng} for comparison purposes. |
| 854 |
< |
The Shifted-Potential shows much narrower angular distributions for |
| 854 |
> |
The {\sc sp} shows much narrower angular distributions for |
| 855 |
|
both the force and torque vectors when using an $\alpha$ of 0.2 |
| 856 |
< |
\AA$^{-1}$ or less, while Shifted-Force shows improvements in the |
| 856 |
> |
\AA$^{-1}$ or less, while {\sc sf} shows improvements in the |
| 857 |
|
undamped and lightly damped cases. Thus, by calculating the |
| 858 |
|
electrostatic interactions in terms of molecular pairs rather than |
| 859 |
|
atomic pairs, the direction of the force and torque vectors are |
| 860 |
|
determined more accurately. |
| 861 |
|
|
| 862 |
|
One additional trend to recognize in table \ref{tab:groupAngle} is |
| 863 |
< |
that the $\sigma^2$ values for both Shifted-Potential and |
| 864 |
< |
Shifted-Force converge as $\alpha$ increases, something that is easier |
| 863 |
> |
that the $\sigma^2$ values for both {\sc sp} and |
| 864 |
> |
{\sc sf} converge as $\alpha$ increases, something that is easier |
| 865 |
|
to see when using group based cutoffs. Looking back on figures |
| 866 |
|
\ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this |
| 867 |
|
behavior clearly at large $\alpha$ and cutoff values. The reason for |
| 880 |
|
high would introduce error in the molecular torques, particularly for |
| 881 |
|
the shorter cutoffs. Based on the above findings, empirical damping |
| 882 |
|
up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably |
| 883 |
< |
unnecessary when using the Shifted-Force method. |
| 883 |
> |
unnecessary when using the {\sc sf} method. |
| 884 |
|
|
| 885 |
|
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 886 |
|
|
| 887 |
< |
In the previous studies using a Shifted-Force variant of the damped |
| 887 |
> |
In the previous studies using a {\sc sf} variant of the damped |
| 888 |
|
Wolf coulomb potential, the structure and dynamics of water were |
| 889 |
|
investigated rather extensively.\cite{Zahn02,Kast03} Their results |
| 890 |
< |
indicated that the damped Shifted-Force method results in properties |
| 890 |
> |
indicated that the damped {\sc sf} method results in properties |
| 891 |
|
very similar to those obtained when using the Ewald summation. |
| 892 |
|
Considering the statistical results shown above, the good performance |
| 893 |
|
of this method is not that surprising. Rather than consider the same |
| 898 |
|
\begin{figure} |
| 899 |
|
\centering |
| 900 |
|
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 901 |
< |
\caption{Power spectra obtained from the velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, Shifted-Force ($\alpha$ = 0, 0.1, \& 0.2), and Shifted-Potential ($\alpha$ = 0.2). Apodization of the correlation functions via a cubic switching function between 40 and 50 ps was used to clear up the spectral noise resulting from data truncation, and had no noticeable effect on peak location or magnitude. The inset shows the frequency region below 100 cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 901 |
> |
\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.} |
| 902 |
|
\label{fig:methodPS} |
| 903 |
|
\end{figure} |
| 904 |
|
|
| 910 |
|
methods differ. Considering the low-frequency inset (expanded in the |
| 911 |
|
upper frame of figure \ref{fig:dampInc}), at frequencies below 100 |
| 912 |
|
cm$^{-1}$, the correlated motions are blue-shifted when using undamped |
| 913 |
< |
or weakly damped Shifted-Force. When using moderate damping ($\alpha |
| 914 |
< |
= 0.2$ \AA$^{-1}$) both the Shifted-Force and Shifted-Potential |
| 913 |
> |
or weakly damped {\sc sf}. When using moderate damping ($\alpha |
| 914 |
> |
= 0.2$ \AA$^{-1}$) both the {\sc sf} and {\sc sp} |
| 915 |
|
methods give near identical correlated motion behavior as the Ewald |
| 916 |
|
method (which has a damping value of 0.3119). The damping acts as a |
| 917 |
|
distance dependent Gaussian screening of the point charges for the |
| 935 |
|
shift to higher frequency in exponential fashion. Though not shown, |
| 936 |
|
the spectrum for the simple undamped electrostatic potential is |
| 937 |
|
blue-shifted such that the lowest frequency peak resides near 325 |
| 938 |
< |
cm$^{-1}$. In light of these results, the undamped Shifted-Force |
| 938 |
> |
cm$^{-1}$. In light of these results, the undamped {\sc sf} |
| 939 |
|
method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is |
| 940 |
|
quite respectable; however, it appears as though moderate damping is |
| 941 |
|
required for accurate reproduction of crystal dynamics. |
| 942 |
|
\begin{figure} |
| 943 |
|
\centering |
| 944 |
|
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 945 |
< |
\caption{Upper: Zoomed inset from figure \ref{fig:methodPS}. As the damping value for the Shifted-Force potential increases, the low-frequency peaks red-shift. Lower: Low-frequency correlated motions for NaCl crystals at 1000 K when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$. As $\alpha$ increases, the peaks are red-shifted toward and eventually beyond the values given by SPME. The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.} |
| 945 |
> |
\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.} |
| 946 |
|
\label{fig:dampInc} |
| 947 |
|
\end{figure} |
| 948 |
|
|
| 953 |
|
electrostatic summation techniques than the Ewald summation, chiefly |
| 954 |
|
methods derived from the damped Coulombic sum originally proposed by |
| 955 |
|
Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the |
| 956 |
< |
Shifted-Force method, reformulated above as equation \ref{eq:SFPot}, |
| 956 |
> |
{\sc sf} method, reformulated above as eq. (\ref{eq:DSFPot}), |
| 957 |
|
shows a remarkable ability to reproduce the energetic and dynamic |
| 958 |
|
characteristics exhibited by simulations employing lattice summation |
| 959 |
|
techniques. The cumulative energy difference results showed the |
| 960 |
< |
undamped Shifted-Force and moderately damped Shifted-Potential methods |
| 960 |
> |
undamped {\sc sf} and moderately damped {\sc sp} methods |
| 961 |
|
produced results nearly identical to SPME. Similarly for the dynamic |
| 962 |
< |
features, the undamped or moderately damped Shifted-Force and |
| 963 |
< |
moderately damped Shifted-Potential methods produce force and torque |
| 962 |
> |
features, the undamped or moderately damped {\sc sf} and |
| 963 |
> |
moderately damped {\sc sp} methods produce force and torque |
| 964 |
|
vector magnitude and directions very similar to the expected values. |
| 965 |
|
These results translate into long-time dynamic behavior equivalent to |
| 966 |
|
that produced in simulations using SPME. |
