| 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} |
| 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 |
< |
using eq. (\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} |
| 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 |
| 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 {\sc sf} 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 |
| 449 |
|
|
| 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 orignal theory for {\sc rf} was originally |
| 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 prespecification of a |
| 459 |
> |
made for {\sc rf}, with the additional pre-specification of a |
| 460 |
|
dielectric constant. |
| 461 |
|
|
| 462 |
|
\section{Methods} |
| 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 |
|
|
| 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 {\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.} |
| 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 |
|
|
