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 |