| 98 |
|
|
| 99 |
|
\subsection{The Wolf and Zahn Methods} |
| 100 |
|
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
| 101 |
< |
for an accurate accumulation of electrostatic interactions in an |
| 101 |
> |
for the accurate accumulation of electrostatic interactions in an |
| 102 |
|
efficient pairwise fashion.\cite{Wolf99} Wolf \textit{et al.} observed |
| 103 |
|
that the electrostatic interaction is effectively short-ranged in |
| 104 |
|
condensed phase systems and that neutralization of the charge |
| 111 |
|
function (identical to that seen in the real-space portion of the |
| 112 |
|
Ewald sum) to aid convergence |
| 113 |
|
\begin{equation} |
| 114 |
< |
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\}. |
| 114 |
> |
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\}. |
| 115 |
|
\label{eq:WolfPot} |
| 116 |
|
\end{equation} |
| 117 |
|
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
| 118 |
|
potential. However, neutralizing the charge contained within each |
| 119 |
|
cutoff sphere requires the placement of a self-image charge on the |
| 120 |
|
surface of the cutoff sphere. This additional self-term in the total |
| 121 |
< |
potential enables Wolf {\it et al.} to obtain excellent estimates of |
| 121 |
> |
potential enabled Wolf {\it et al.} to obtain excellent estimates of |
| 122 |
|
Madelung energies for many crystals. |
| 123 |
|
|
| 124 |
|
In order to use their charge-neutralized potential in molecular |
| 126 |
|
derivative of this potential prior to evaluation of the limit. This |
| 127 |
|
procedure gives an expression for the forces, |
| 128 |
|
\begin{equation} |
| 129 |
< |
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\}, |
| 129 |
> |
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\}, |
| 130 |
|
\label{eq:WolfForces} |
| 131 |
|
\end{equation} |
| 132 |
|
that incorporates both image charges and damping of the electrostatic |
| 134 |
|
|
| 135 |
|
More recently, Zahn \textit{et al.} investigated these potential and |
| 136 |
|
force expressions for use in simulations involving water.\cite{Zahn02} |
| 137 |
< |
In their work, they pointed out that the method that the forces and |
| 138 |
< |
derivative of the potential are not commensurate. Attempts to use |
| 139 |
< |
both Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will |
| 140 |
< |
lead to poor energy conservation. They correctly observed that taking |
| 141 |
< |
the limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating |
| 142 |
< |
the derivatives is mathematically invalid. |
| 137 |
> |
In their work, they pointed out that the forces and derivative of |
| 138 |
> |
the potential are not commensurate. Attempts to use both |
| 139 |
> |
Eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead |
| 140 |
> |
to poor energy conservation. They correctly observed that taking the |
| 141 |
> |
limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the |
| 142 |
> |
derivatives gives forces for a different potential energy function |
| 143 |
> |
than the one shown in Eq. (\ref{eq:WolfPot}). |
| 144 |
|
|
| 145 |
|
Zahn \textit{et al.} proposed a modified form of this ``Wolf summation |
| 146 |
|
method'' as a way to use this technique in Molecular Dynamics |
| 148 |
|
\ref{eq:WolfForces}, they proposed a new damped Coulomb |
| 149 |
|
potential, |
| 150 |
|
\begin{equation} |
| 151 |
< |
V^{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\}. |
| 151 |
> |
V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\}. |
| 152 |
|
\label{eq:ZahnPot} |
| 153 |
|
\end{equation} |
| 154 |
|
They showed that this potential does fairly well at capturing the |
| 159 |
|
|
| 160 |
|
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et |
| 161 |
|
al.} are constructed using two different (and separable) computational |
| 162 |
< |
tricks: \begin{itemize} |
| 162 |
> |
tricks: \begin{enumerate} |
| 163 |
|
\item shifting through the use of image charges, and |
| 164 |
|
\item damping the electrostatic interaction. |
| 165 |
< |
\end{itemize} Wolf \textit{et al.} treated the |
| 165 |
> |
\end{enumerate} Wolf \textit{et al.} treated the |
| 166 |
|
development of their summation method as a progressive application of |
| 167 |
|
these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded |
| 168 |
|
their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the |
| 182 |
|
\textit{et al.} and Zahn \textit{et al.} by considering the standard |
| 183 |
|
shifted potential, |
| 184 |
|
\begin{equation} |
| 185 |
< |
v^\textrm{SP}(r) = \begin{cases} |
| 185 |
> |
v_\textrm{SP}(r) = \begin{cases} |
| 186 |
|
v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > |
| 187 |
|
R_\textrm{c} |
| 188 |
|
\end{cases}, |
| 190 |
|
\end{equation} |
| 191 |
|
and shifted force, |
| 192 |
|
\begin{equation} |
| 193 |
< |
v^\textrm{SF}(r) = \begin{cases} |
| 194 |
< |
v(r)-v_\textrm{c}-\left(\frac{\textrm{d}v(r)}{\textrm{d}r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) |
| 193 |
> |
v_\textrm{SF}(r) = \begin{cases} |
| 194 |
> |
v(r)-v_\textrm{c}-\left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) |
| 195 |
|
&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} |
| 196 |
|
\end{cases}, |
| 197 |
|
\label{eq:shiftingForm} |
| 203 |
|
potential is smooth at the cutoff radius |
| 204 |
|
($R_\textrm{c}$).\cite{Allen87} |
| 205 |
|
|
| 206 |
< |
|
| 207 |
< |
|
| 207 |
< |
|
| 208 |
< |
If the derivative term is taken to be zero, we are left with the shifted Coulomb potential devised by Wolf \textit{et al.},\cite{Wolf99} |
| 206 |
> |
The forces associated with the shifted potential are simply the forces |
| 207 |
> |
of the unshifted potential itself (when inside the cutoff sphere), |
| 208 |
|
\begin{equation} |
| 209 |
< |
V^\textrm{SP}(r_{ij}) = q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}. \label{eq:WolfSP} |
| 209 |
> |
F_{\textrm{SP}} = \left( \frac{d v(r)}{dr} \right), |
| 210 |
|
\end{equation} |
| 211 |
< |
The forces associated with this potential are obtained by taking the derivative, resulting in the following, |
| 211 |
> |
and are zero outside. Inside the cutoff sphere, the forces associated |
| 212 |
> |
with the shifted force form can be written, |
| 213 |
|
\begin{equation} |
| 214 |
< |
F^\textrm{SP}(r_{ij}) = q_iq_j\left(-\frac{1}{r_{ij}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}. |
| 214 |
> |
F_{\textrm{SF}} = \left( \frac{d v(r)}{dr} \right) - \left(\frac{d |
| 215 |
> |
v(r)}{dr} \right)_{r=R_\textrm{c}}. |
| 216 |
> |
\end{equation} |
| 217 |
> |
|
| 218 |
> |
If the potential ($v(r)$) is taken to be the normal Coulomb potential, |
| 219 |
> |
\begin{equation} |
| 220 |
> |
v(r) = \frac{q_i q_j}{r}, |
| 221 |
> |
\label{eq:Coulomb} |
| 222 |
> |
\end{equation} |
| 223 |
> |
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
| 224 |
> |
al.}'s undamped prescription: |
| 225 |
> |
\begin{equation} |
| 226 |
> |
V_\textrm{SP}(r) = |
| 227 |
> |
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
| 228 |
> |
r\leqslant R_\textrm{c}, |
| 229 |
> |
\label{eq:WolfSP} |
| 230 |
> |
\end{equation} |
| 231 |
> |
with associated forces, |
| 232 |
> |
\begin{equation} |
| 233 |
> |
F_\textrm{SP}(r) = q_iq_j\left(-\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
| 234 |
|
\label{eq:FWolfSP} |
| 235 |
|
\end{equation} |
| 236 |
< |
These forces are identical to the forces of the standard electrostatic interaction, and this was addressed by Wolf \textit{et al.} as undesirable. They pointed out that the effect of the image charges is neglected in the forces when they would expect there to be some pressure exerted due to their presence.\cite{Wolf99} As a consequence the forces, though mathematically valid, may not be physically correct due to this missing component. Additionally, there is a discontinuity in the forces. This can be remedied with the use of a switching function to zero the potential and forces smoothly as particles near $R_\textrm{c}$. |
| 236 |
> |
These forces are identical to the forces of the standard Coulomb |
| 237 |
> |
interaction, and cutting these off at $R_c$ was addressed by Wolf |
| 238 |
> |
\textit{et al.} as undesirable. They pointed out that the effect of |
| 239 |
> |
the image charges is neglected in the forces when this form is |
| 240 |
> |
used,\cite{Wolf99} thereby eliminating any benefit from the method in |
| 241 |
> |
molecular dynamics. Additionally, there is a discontinuity in the |
| 242 |
> |
forces at the cutoff radius which results in energy drift during MD |
| 243 |
> |
simulations. |
| 244 |
|
|
| 245 |
< |
If the derivative term in equation \ref{eq:shiftingForm} is evaluated, we obtain an hitherto undiscussed shifted force Coulomb potential, |
| 245 |
> |
The shifted force ({\sc sf}) form using the normal Coulomb potential |
| 246 |
> |
will give, |
| 247 |
|
\begin{equation} |
| 248 |
< |
V^\textrm{SF}(r_{ij}) = q_iq_j\left[\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r_{ij}-R_\textrm{c})\right] \quad r_{ij}\leqslant R_\textrm{c}. |
| 248 |
> |
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}. |
| 249 |
|
\label{eq:SFPot} |
| 250 |
|
\end{equation} |
| 251 |
< |
Taking the derivative of this shifted force potential gives the |
| 225 |
< |
following forces, |
| 251 |
> |
with associated forces, |
| 252 |
|
\begin{equation} |
| 253 |
< |
F^\textrm{SF}(r_{ij} = q_iq_j\left(-\frac{1}{r_{ij}^2}+\frac{1}{R_\textrm{c}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}. |
| 253 |
> |
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}. |
| 254 |
|
\label{eq:SFForces} |
| 255 |
|
\end{equation} |
| 256 |
< |
Using this formulation rather than the simple shifted potential |
| 257 |
< |
(Eq. \ref{eq:WolfSP}) means that there are no discontinuities in the |
| 258 |
< |
forces in addition to the potential. This form also has the benefit |
| 259 |
< |
that the image charges have a force presence, addressing the concerns |
| 260 |
< |
about a missing physical component. One side effect of this treatment |
| 261 |
< |
is a slight alteration in the shape of the potential that comes about |
| 262 |
< |
from the derivative term. Thus, a degree of clarity about the |
| 263 |
< |
original formulation of the potential is lost in order to gain |
| 264 |
< |
functionality in dynamics simulations. |
| 256 |
> |
This formulation has the benefits that there are no discontinuities at |
| 257 |
> |
the cutoff distance, while the neutralizing image charges are present |
| 258 |
> |
in both the energy and force expressions. It would be simple to add |
| 259 |
> |
the self-neutralizing term back when computing the total energy of the |
| 260 |
> |
system, thereby maintaining the agreement with the Madelung energies. |
| 261 |
> |
A side effect of this treatment is the alteration in the shape of the |
| 262 |
> |
potential that comes from the derivative term. Thus, a degree of |
| 263 |
> |
clarity about agreement with the empirical potential is lost in order |
| 264 |
> |
to gain functionality in dynamics simulations. |
| 265 |
|
|
| 266 |
|
Wolf \textit{et al.} originally discussed the energetics of the |
| 267 |
|
shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that |
| 268 |
< |
it was still insufficient for accurate determination of the energy. |
| 269 |
< |
The energy would fluctuate around the expected value with increasing |
| 270 |
< |
cutoff radius, but the oscillations appeared to be converging toward |
| 271 |
< |
the correct value.\cite{Wolf99} A damping function was incorporated to |
| 272 |
< |
accelerate convergence; and though alternative functional forms could |
| 273 |
< |
be used,\cite{Jones56,Heyes81} the complimentary error function was |
| 274 |
< |
chosen to draw parallels to the Ewald summation. Incorporating |
| 275 |
< |
damping into the simple Coulomb potential, |
| 268 |
> |
it was still insufficient for accurate determination of the energy |
| 269 |
> |
with reasonable cutoff distances. The calculated Madelung energies |
| 270 |
> |
fluctuate around the expected value with increasing cutoff radius, but |
| 271 |
> |
the oscillations converge toward the correct value.\cite{Wolf99} A |
| 272 |
> |
damping function was incorporated to accelerate the convergence; and |
| 273 |
> |
though alternative functional forms could be |
| 274 |
> |
used,\cite{Jones56,Heyes81} the complimentary error function was |
| 275 |
> |
chosen to mirror the effective screening used in the Ewald summation. |
| 276 |
> |
Incorporating this error function damping into the simple Coulomb |
| 277 |
> |
potential, |
| 278 |
|
\begin{equation} |
| 279 |
< |
v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}, |
| 279 |
> |
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
| 280 |
|
\label{eq:dampCoulomb} |
| 281 |
|
\end{equation} |
| 282 |
< |
the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived |
| 282 |
> |
the shifted potential (Eq. \ref{eq:WolfSP}) can be recovered |
| 283 |
|
\textit{via} equation \ref{eq:shiftingForm}, |
| 284 |
|
\begin{equation} |
| 285 |
< |
V^{\textrm{DSP}}(r_{ij}) = q_iq_j\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}. |
| 285 |
> |
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}. |
| 286 |
|
\label{eq:DSPPot} |
| 287 |
< |
\end{equation} |
| 288 |
< |
The derivative of this Shifted-Potential can be taken to obtain forces |
| 261 |
< |
for use in MD, |
| 287 |
> |
\end{equation}, |
| 288 |
> |
with associated forces, |
| 289 |
|
\begin{equation} |
| 290 |
< |
F^{\textrm{DSP}}(r_{ij}) = q_iq_j\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) \quad r_{ij}\leqslant R_\textrm{c}. |
| 290 |
> |
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}. |
| 291 |
|
\label{eq:DSPForces} |
| 292 |
|
\end{equation} |
| 293 |
< |
Again, this Shifted-Potential suffers from a discontinuity in the |
| 294 |
< |
forces, and a lack of an image-charge component in the forces. To |
| 295 |
< |
remedy these concerns, a Shifted-Force variant is obtained by |
| 296 |
< |
inclusion of the derivative term in equation \ref{eq:shiftingForm} to |
| 270 |
< |
give, |
| 293 |
> |
Again, this damped shifted potential suffers from a discontinuity and |
| 294 |
> |
a lack of the image charges in the forces. To remedy these concerns, |
| 295 |
> |
one may derive a Shifted-Force variant by including the derivative |
| 296 |
> |
term in equation \ref{eq:shiftingForm}, |
| 297 |
|
\begin{equation} |
| 298 |
|
\begin{split} |
| 299 |
< |
V^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\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_{ij}-R_\mathrm{c}\right)\ \right] \quad r_{ij}\leqslant R_\textrm{c}. |
| 299 |
> |
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}. |
| 300 |
|
\label{eq:DSFPot} |
| 301 |
|
\end{split} |
| 302 |
|
\end{equation} |
| 303 |
|
The derivative of the above potential gives the following forces, |
| 304 |
|
\begin{equation} |
| 305 |
|
\begin{split} |
| 306 |
< |
F^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&-\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.+\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_{ij}\leqslant R_\textrm{c}. |
| 306 |
> |
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}. |
| 307 |
|
\label{eq:DSFForces} |
| 308 |
|
\end{split} |
| 309 |
|
\end{equation} |
| 311 |
|
This new Shifted-Force potential is similar to equation |
| 312 |
|
\ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are |
| 313 |
|
two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term |
| 314 |
< |
from equation \ref{eq:shiftingForm} is equal to equation |
| 315 |
< |
\ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This |
| 316 |
< |
term is not present in the Zahn potential, resulting in a |
| 317 |
< |
discontinuity as particles cross $R_\textrm{c}$. Second, the sign of |
| 318 |
< |
the derivative portion is different. The constant $v_\textrm{c}$ term |
| 319 |
< |
is not going to have a presence in the forces after performing the |
| 320 |
< |
derivative, but the negative sign does effect the derivative. In |
| 321 |
< |
fact, it introduces a discontinuity in the forces at the cutoff, |
| 322 |
< |
because the force function is shifted in the wrong direction and |
| 323 |
< |
doesn't cross zero at $R_\textrm{c}$. Thus, these alterations make |
| 324 |
< |
for an electrostatic summation method that is continuous in both the |
| 325 |
< |
potential and forces and incorporates the pairwise sum considerations |
| 300 |
< |
stressed by Wolf \textit{et al.}\cite{Wolf99} |
| 314 |
> |
from eq. (\ref{eq:shiftingForm}) is equal to |
| 315 |
> |
eq. (\ref{eq:dampCoulomb}) with $r$ replaced by $R_\textrm{c}$. This |
| 316 |
> |
term is {\it not} present in the Zahn potential, resulting in a |
| 317 |
> |
potential discontinuity as particles cross $R_\textrm{c}$. Second, |
| 318 |
> |
the sign of the derivative portion is different. The missing |
| 319 |
> |
$v_\textrm{c}$ term would not affect molecular dynamics simulations |
| 320 |
> |
(although the computed energy would be expected to have sudden jumps |
| 321 |
> |
as particle distances crossed $R_c$). The sign problem would be a |
| 322 |
> |
potential source of errors, however. In fact, it introduces a |
| 323 |
> |
discontinuity in the forces at the cutoff, because the force function |
| 324 |
> |
is shifted in the wrong direction and doesn't cross zero at |
| 325 |
> |
$R_\textrm{c}$. |
| 326 |
|
|
| 327 |
+ |
Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
| 328 |
+ |
electrostatic summation method that is continuous in both the |
| 329 |
+ |
potential and forces and which incorporates the damping function |
| 330 |
+ |
proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of this |
| 331 |
+ |
paper, we will evaluate exactly how good these methods ({\sc sp}, {\sc |
| 332 |
+ |
sf}, damping) are at reproducing the correct electrostatic summation |
| 333 |
+ |
performed by the Ewald sum. |
| 334 |
+ |
|
| 335 |
+ |
\subsection{Other alternatives} |
| 336 |
+ |
|
| 337 |
+ |
Reaction Field blah |
| 338 |
+ |
|
| 339 |
+ |
Group-based cutoff blah |
| 340 |
+ |
|
| 341 |
+ |
|
| 342 |
|
\section{Methods} |
| 343 |
|
|
| 304 |
– |
\subsection{What Qualities are Important?}\label{sec:Qualities} |
| 344 |
|
In classical molecular mechanics simulations, there are two primary |
| 345 |
|
techniques utilized to obtain information about the system of |
| 346 |
|
interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these |
| 349 |
|
|
| 350 |
|
In MC, the potential energy difference between two subsequent |
| 351 |
|
configurations dictates the progression of MC sampling. Going back to |
| 352 |
< |
the origins of this method, the Canonical ensemble acceptance criteria |
| 353 |
< |
laid out by Metropolis \textit{et al.} states that a subsequent |
| 354 |
< |
configuration is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta |
| 355 |
< |
E/kT)$, where $\xi$ is a random number between 0 and |
| 356 |
< |
1.\cite{Metropolis53} Maintaining a consistent $\Delta E$ when using |
| 357 |
< |
an alternate method for handling the long-range electrostatics ensures |
| 358 |
< |
proper sampling within the ensemble. |
| 352 |
> |
the origins of this method, the acceptance criterion for the canonical |
| 353 |
> |
ensemble laid out by Metropolis \textit{et al.} states that a |
| 354 |
> |
subsequent configuration is accepted if $\Delta E < 0$ or if $\xi < |
| 355 |
> |
\exp(-\Delta E/kT)$, where $\xi$ is a random number between 0 and |
| 356 |
> |
1.\cite{Metropolis53} Maintaining the correct $\Delta E$ when using an |
| 357 |
> |
alternate method for handling the long-range electrostatics will |
| 358 |
> |
ensure proper sampling from the ensemble. |
| 359 |
|
|
| 360 |
< |
In MD, the derivative of the potential directs how the system will |
| 360 |
> |
In MD, the derivative of the potential governs how the system will |
| 361 |
|
progress in time. Consequently, the force and torque vectors on each |
| 362 |
< |
body in the system dictate how it develops as a whole. If the |
| 363 |
< |
magnitude and direction of these vectors are similar when using |
| 364 |
< |
alternate electrostatic summation techniques, the dynamics in the near |
| 365 |
< |
term will be indistinguishable. Because error in MD calculations is |
| 366 |
< |
cumulative, one should expect greater deviation in the long term |
| 367 |
< |
trajectories with greater differences in these vectors between |
| 368 |
< |
configurations using different long-range electrostatics. |
| 362 |
> |
body in the system dictate how the system evolves. If the magnitude |
| 363 |
> |
and direction of these vectors are similar when using alternate |
| 364 |
> |
electrostatic summation techniques, the dynamics in the short term |
| 365 |
> |
will be indistinguishable. Because error in MD calculations is |
| 366 |
> |
cumulative, one should expect greater deviation at longer times, |
| 367 |
> |
although methods which have large differences in the force and torque |
| 368 |
> |
vectors will diverge from each other more rapidly. |
| 369 |
|
|
| 370 |
|
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 371 |
< |
Evaluation of the pairwise summation techniques (outlined in section |
| 372 |
< |
\ref{sec:ESMethods}) for use in MC simulations was performed through |
| 373 |
< |
study of the energy differences between conformations. Considering |
| 374 |
< |
the SPME results to be the correct or desired behavior, ideal |
| 375 |
< |
performance of a tested method was taken to be agreement between the |
| 376 |
< |
energy differences calculated. Linear least squares regression of the |
| 377 |
< |
$\Delta E$ values between configurations using SPME against $\Delta E$ |
| 378 |
< |
values using tested methods provides a quantitative comparison of this |
| 379 |
< |
agreement. Unitary results for both the correlation and correlation |
| 380 |
< |
coefficient for these regressions indicate equivalent energetic |
| 381 |
< |
results between the methods. The correlation is the slope of the |
| 382 |
< |
plotted data while the correlation coefficient ($R^2$) is a measure of |
| 383 |
< |
the of the data scatter around the fitted line and tells about the |
| 384 |
< |
quality of the fit (Fig. \ref{fig:linearFit}). |
| 371 |
> |
The pairwise summation techniques (outlined in section |
| 372 |
> |
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 373 |
> |
studying the energy differences between conformations. We took the |
| 374 |
> |
SPME-computed energy difference between two conformations to be the |
| 375 |
> |
correct behavior. An ideal performance by an alternative method would |
| 376 |
> |
reproduce these energy differences exactly. Since none of the methods |
| 377 |
> |
provide exact energy differences, we used linear least squares |
| 378 |
> |
regressions of the $\Delta E$ values between configurations using SPME |
| 379 |
> |
against $\Delta E$ values using tested methods provides a quantitative |
| 380 |
> |
comparison of this agreement. Unitary results for both the |
| 381 |
> |
correlation and correlation coefficient for these regressions indicate |
| 382 |
> |
equivalent energetic results between the method under consideration |
| 383 |
> |
and electrostatics handled using SPME. Sample correlation plots for |
| 384 |
> |
two alternate methods are shown in Fig. \ref{fig:linearFit}. |
| 385 |
|
|
| 386 |
|
\begin{figure} |
| 387 |
|
\centering |
| 390 |
|
\label{fig:linearFit} |
| 391 |
|
\end{figure} |
| 392 |
|
|
| 393 |
< |
Each system type (detailed in section \ref{sec:RepSims}) studied |
| 394 |
< |
consisted of 500 independent configurations, each equilibrated from |
| 395 |
< |
higher temperature trajectories. Thus, 124,750 $\Delta E$ data points |
| 396 |
< |
are used in a regression of a single system type. Results and |
| 397 |
< |
discussion for the individual analysis of each of the system types |
| 359 |
< |
appear in the supporting information, while the cumulative results |
| 360 |
< |
over all the investigated systems appears below in section |
| 361 |
< |
\ref{sec:EnergyResults}. |
| 393 |
> |
Each system type (detailed in section \ref{sec:RepSims}) was |
| 394 |
> |
represented using 500 independent configurations. Additionally, we |
| 395 |
> |
used seven different system types, so each of the alternate |
| 396 |
> |
(non-Ewald) electrostatic summation methods was evaluated using |
| 397 |
> |
873,250 configurational energy differences. |
| 398 |
|
|
| 399 |
< |
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
| 400 |
< |
Evaluation of the pairwise methods (outlined in section |
| 401 |
< |
\ref{sec:ESMethods}) for use in MD simulations was performed through |
| 402 |
< |
comparison of the force and torque vectors obtained with those from |
| 367 |
< |
SPME. Both the magnitude and the direction of these vectors on each |
| 368 |
< |
of the bodies in the system were analyzed. For the magnitude of these |
| 369 |
< |
vectors, linear least squares regression analysis can be performed as |
| 370 |
< |
described previously for comparing $\Delta E$ values. Instead of a |
| 371 |
< |
single value between two system configurations, there is a value for |
| 372 |
< |
each particle in each configuration. For a system of 1000 water |
| 373 |
< |
molecules and 40 ions, there are 1040 force vectors and 1000 torque |
| 374 |
< |
vectors. With 500 configurations, this results in 520,000 force and |
| 375 |
< |
500,000 torque vector comparisons samples for each system type. |
| 399 |
> |
Results and discussion for the individual analysis of each of the |
| 400 |
> |
system types appear in the supporting information, while the |
| 401 |
> |
cumulative results over all the investigated systems appears below in |
| 402 |
> |
section \ref{sec:EnergyResults}. |
| 403 |
|
|
| 404 |
< |
The force and torque vector directions were investigated through |
| 405 |
< |
measurement of the angle ($\theta$) formed between those from the |
| 406 |
< |
particular method and those from SPME |
| 404 |
> |
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
| 405 |
> |
We evaluated the pairwise methods (outlined in section |
| 406 |
> |
\ref{sec:ESMethods}) for use in MD simulations by |
| 407 |
> |
comparing the force and torque vectors with those obtained using the |
| 408 |
> |
reference Ewald summation (SPME). Both the magnitude and the |
| 409 |
> |
direction of these vectors on each of the bodies in the system were |
| 410 |
> |
analyzed. For the magnitude of these vectors, linear least squares |
| 411 |
> |
regression analyses were performed as described previously for |
| 412 |
> |
comparing $\Delta E$ values. Instead of a single energy difference |
| 413 |
> |
between two system configurations, we compared the magnitudes of the |
| 414 |
> |
forces (and torques) on each molecule in each configuration. For a |
| 415 |
> |
system of 1000 water molecules and 40 ions, there are 1040 force |
| 416 |
> |
vectors and 1000 torque vectors. With 500 configurations, this |
| 417 |
> |
results in 520,000 force and 500,000 torque vector comparisons. |
| 418 |
> |
Additionally, data from seven different system types was aggregated |
| 419 |
> |
before the comparison was made. |
| 420 |
> |
|
| 421 |
> |
The {\it directionality} of the force and torque vectors was |
| 422 |
> |
investigated through measurement of the angle ($\theta$) formed |
| 423 |
> |
between those computed from the particular method and those from SPME, |
| 424 |
|
\begin{equation} |
| 425 |
< |
\theta_F = \frac{\vec{F}_\textrm{SPME}}{|\vec{F}_\textrm{SPME}|}\cdot\frac{\vec{F}_\textrm{Method}}{|\vec{F}_\textrm{Method}|}. |
| 425 |
> |
\theta_f = \cos^{-1} \hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}, |
| 426 |
|
\end{equation} |
| 427 |
+ |
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the |
| 428 |
+ |
force vector computed using method $M$. |
| 429 |
+ |
|
| 430 |
|
Each of these $\theta$ values was accumulated in a distribution |
| 431 |
< |
function, weighted by the area on the unit sphere. Non-linear fits |
| 432 |
< |
were used to measure the shape of the resulting distributions. |
| 431 |
> |
function, weighted by the area on the unit sphere. Non-linear |
| 432 |
> |
Gaussian fits were used to measure the width of the resulting |
| 433 |
> |
distributions. |
| 434 |
|
|
| 435 |
|
\begin{figure} |
| 436 |
|
\centering |
| 443 |
|
non-linear fits. The solid line is a Gaussian profile, while the |
| 444 |
|
dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 445 |
|
Lorentzian. Since this distribution is a measure of angular error |
| 446 |
< |
between two different electrostatic summation methods, there is |
| 447 |
< |
particular reason for the profile to adhere to a specific shape. |
| 448 |
< |
Because of this and the Gaussian profile's more statistically |
| 449 |
< |
meaningful properties, Gaussian fits was used to compare all the |
| 450 |
< |
tested methods. The variance ($\sigma^2$) was extracted from each of |
| 451 |
< |
these fits and was used to compare distribution widths. Values of |
| 452 |
< |
$\sigma^2$ near zero indicate vector directions indistinguishable from |
| 453 |
< |
those calculated when using SPME. |
| 446 |
> |
between two different electrostatic summation methods, there is no |
| 447 |
> |
{\it a priori} reason for the profile to adhere to any specific shape. |
| 448 |
> |
Gaussian fits was used to compare all the tested methods. The |
| 449 |
> |
variance ($\sigma^2$) was extracted from each of these fits and was |
| 450 |
> |
used to compare distribution widths. Values of $\sigma^2$ near zero |
| 451 |
> |
indicate vector directions indistinguishable from those calculated |
| 452 |
> |
when using the reference method (SPME). |
| 453 |
> |
|
| 454 |
> |
\subsection{Short-time Dynamics} |
| 455 |
|
|
| 456 |
|
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 457 |
|
Evaluation of the long-time dynamics of charged systems was performed |
| 801 |
|
clearly, we show how damping strength affects a simple real-space |
| 802 |
|
electrostatic potential, |
| 803 |
|
\begin{equation} |
| 804 |
< |
V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r_{ij}})}{r_{ij}}\right]S(r), |
| 804 |
> |
V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r), |
| 805 |
|
\end{equation} |
| 806 |
|
where $S(r)$ is a switching function that smoothly zeroes the |
| 807 |
|
potential at the cutoff radius. Figure \ref{fig:dampInc} shows how |
