| 58 |
|
|
| 59 |
|
%\narrowtext |
| 60 |
|
|
| 61 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 61 |
> |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 62 |
|
% BODY OF TEXT |
| 63 |
< |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 63 |
> |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
| 64 |
|
|
| 65 |
|
\section{Introduction} |
| 66 |
|
|
| 67 |
|
In molecular simulations, proper accumulation of the electrostatic |
| 68 |
|
interactions is considered one of the most essential and |
| 69 |
< |
computationally demanding tasks. |
| 69 |
> |
computationally demanding tasks. The common molecular mechanics force |
| 70 |
> |
fields are founded on representation of the atomic sites centered on |
| 71 |
> |
full or partial charges shielded by Lennard-Jones type interactions. |
| 72 |
> |
This means that nearly every pair interaction involves an |
| 73 |
> |
charge-charge calculation. Coupled with $r^{-1}$ decay, the monopole |
| 74 |
> |
interactions quickly become a burden for molecular systems of all |
| 75 |
> |
sizes. For example, in small systems, the electrostatic pair |
| 76 |
> |
interaction may not have decayed appreciably within the box length |
| 77 |
> |
leading to an effect excluded from the pair interactions within a unit |
| 78 |
> |
box. In large systems, excessively large cutoffs need to be used to |
| 79 |
> |
accurately incorporate their effect, and since the computational cost |
| 80 |
> |
increases proportionally with the cutoff sphere, it quickly becomes an |
| 81 |
> |
impractical task to perform these calculations. |
| 82 |
|
|
| 83 |
|
\subsection{The Ewald Sum} |
| 84 |
|
blah blah blah Ewald Sum Important blah blah blah |
| 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} |
| 209 |
|
\begin{equation} |
| 210 |
< |
V^\textrm{WSP}(r_{ij}) = \begin{cases} q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 199 |
< |
\end{cases}. |
| 200 |
< |
\label{eq:WolfSP} |
| 210 |
> |
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} |
| 211 |
|
\end{equation} |
| 212 |
|
The forces associated with this potential are obtained by taking the derivative, resulting in the following, |
| 213 |
|
\begin{equation} |
| 214 |
< |
F^\textrm{WSP}(r_{ij}) = \begin{cases} q_iq_j\left(-\frac{1}{r_{ij}^2}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 205 |
< |
\end{cases}. |
| 214 |
> |
F^\textrm{SP}(r_{ij}) = q_iq_j\left(-\frac{1}{r_{ij}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}. |
| 215 |
|
\label{eq:FWolfSP} |
| 216 |
|
\end{equation} |
| 217 |
|
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}$. |
| 218 |
|
|
| 219 |
|
If the derivative term in equation \ref{eq:shiftingForm} is evaluated, we obtain an hitherto undiscussed shifted force Coulomb potential, |
| 220 |
|
\begin{equation} |
| 221 |
< |
V^\textrm{SF}(r_{ij}) = \begin{cases} 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} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 213 |
< |
\end{cases}. |
| 221 |
> |
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}. |
| 222 |
|
\label{eq:SFPot} |
| 223 |
|
\end{equation} |
| 224 |
< |
Taking the derivative of this shifted force potential gives the following forces, |
| 224 |
> |
Taking the derivative of this shifted force potential gives the |
| 225 |
> |
following forces, |
| 226 |
|
\begin{equation} |
| 227 |
< |
F^\textrm{SF}(r_{ij}) = \begin{cases} q_iq_j\left(-\frac{1}{r_{ij}^2}+\frac{1}{R_\textrm{c}^2}\right) &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 219 |
< |
\end{cases}. |
| 227 |
> |
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}. |
| 228 |
|
\label{eq:SFForces} |
| 229 |
|
\end{equation} |
| 230 |
< |
Using this formulation rather than the simple shifted potential (Eq. \ref{eq:WolfSP}) means that there are no discontinuities in the forces in addition to the potential. This form also has the benefit that the image charges have a force presence, addressing the concerns about a missing physical component. One side effect of this treatment is a slight alteration in the shape of the potential that comes about from the derivative term. Thus, a degree of clarity about the original formulation of the potential is lost in order to gain functionality in dynamics simulations. |
| 230 |
> |
Using this formulation rather than the simple shifted potential |
| 231 |
> |
(Eq. \ref{eq:WolfSP}) means that there are no discontinuities in the |
| 232 |
> |
forces in addition to the potential. This form also has the benefit |
| 233 |
> |
that the image charges have a force presence, addressing the concerns |
| 234 |
> |
about a missing physical component. One side effect of this treatment |
| 235 |
> |
is a slight alteration in the shape of the potential that comes about |
| 236 |
> |
from the derivative term. Thus, a degree of clarity about the |
| 237 |
> |
original formulation of the potential is lost in order to gain |
| 238 |
> |
functionality in dynamics simulations. |
| 239 |
|
|
| 240 |
< |
Wolf \textit{et al.} originally discussed the energetics of the shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that it was still insufficient for accurate determination of the energy. The energy would fluctuate around the expected value with increasing cutoff radius, but the oscillations appeared to be converging toward the correct value.\cite{Wolf99} A damping function was incorporated to accelerate convergence; and though alternative functional forms could be used,\cite{Jones56,Heyes81} the complimentary error function was chosen to draw parallels to the Ewald summation. Incorporating damping into the simple Coulomb potential, |
| 240 |
> |
Wolf \textit{et al.} originally discussed the energetics of the |
| 241 |
> |
shifted Coulomb potential (Eq. \ref{eq:WolfSP}), and they found that |
| 242 |
> |
it was still insufficient for accurate determination of the energy. |
| 243 |
> |
The energy would fluctuate around the expected value with increasing |
| 244 |
> |
cutoff radius, but the oscillations appeared to be converging toward |
| 245 |
> |
the correct value.\cite{Wolf99} A damping function was incorporated to |
| 246 |
> |
accelerate convergence; and though alternative functional forms could |
| 247 |
> |
be used,\cite{Jones56,Heyes81} the complimentary error function was |
| 248 |
> |
chosen to draw parallels to the Ewald summation. Incorporating |
| 249 |
> |
damping into the simple Coulomb potential, |
| 250 |
|
\begin{equation} |
| 251 |
|
v(r_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}, |
| 252 |
|
\label{eq:dampCoulomb} |
| 253 |
|
\end{equation} |
| 254 |
< |
the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived \textit{via} equation \ref{eq:shiftingForm}, |
| 254 |
> |
the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived |
| 255 |
> |
\textit{via} equation \ref{eq:shiftingForm}, |
| 256 |
|
\begin{equation} |
| 257 |
< |
V^{\textrm{DSP}}(r_{ij}) = \begin{cases} 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} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 232 |
< |
\end{cases}. |
| 257 |
> |
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}. |
| 258 |
|
\label{eq:DSPPot} |
| 259 |
|
\end{equation} |
| 260 |
< |
The derivative of this Shifted-Potential can be taken to obtain forces for use in MD, |
| 260 |
> |
The derivative of this Shifted-Potential can be taken to obtain forces |
| 261 |
> |
for use in MD, |
| 262 |
|
\begin{equation} |
| 263 |
< |
F^{\textrm{DSP}}(r_{ij}) = \begin{cases} 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} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 238 |
< |
\end{cases}. |
| 263 |
> |
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}. |
| 264 |
|
\label{eq:DSPForces} |
| 265 |
|
\end{equation} |
| 266 |
< |
Again, this Shifted-Potential suffers from a discontinuity in the forces, and a lack of an image-charge component in the forces. To remedy these concerns, a Shifted-Force variant is obtained by inclusion of the derivative term in equation \ref{eq:shiftingForm} to give, |
| 266 |
> |
Again, this Shifted-Potential suffers from a discontinuity in the |
| 267 |
> |
forces, and a lack of an image-charge component in the forces. To |
| 268 |
> |
remedy these concerns, a Shifted-Force variant is obtained by |
| 269 |
> |
inclusion of the derivative term in equation \ref{eq:shiftingForm} to |
| 270 |
> |
give, |
| 271 |
|
\begin{equation} |
| 272 |
< |
V^\mathrm{DSF}(r_{ij}) = \begin{cases} q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}}\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} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 273 |
< |
\end{cases}. |
| 272 |
> |
\begin{split} |
| 273 |
> |
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}. |
| 274 |
|
\label{eq:DSFPot} |
| 275 |
+ |
\end{split} |
| 276 |
|
\end{equation} |
| 277 |
|
The derivative of the above potential gives the following forces, |
| 278 |
|
\begin{equation} |
| 279 |
< |
F^\mathrm{DSF}(r_{ij}) = \begin{cases} q_iq_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}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2R_{\textrm{c}}^2)}}{R_{\textrm{c}}}\right]\right\} &\quad r_{ij}\leqslant R_\textrm{c} \\ 0 &\quad r_{ij}>R_\textrm{c} |
| 280 |
< |
\end{cases}. |
| 279 |
> |
\begin{split} |
| 280 |
> |
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}. |
| 281 |
|
\label{eq:DSFForces} |
| 282 |
+ |
\end{split} |
| 283 |
|
\end{equation} |
| 284 |
|
|
| 285 |
< |
This new Shifted-Force potential is similar to equation \ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from equation \ref{eq:shiftingForm} is equal to equation \ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This term is not present in the Zahn potential, resulting in a discontinuity as particles cross $R_\textrm{c}$. Second, the sign of the derivative portion is different. The constant $v_\textrm{c}$ term is not going to have a presence in the forces after performing the derivative, but the negative sign does effect the derivative. In fact, it introduces a discontinuity in the forces at the cutoff, because the force function is shifted in the wrong direction and doesn't cross zero at $R_\textrm{c}$. Thus, these alterations make for an electrostatic summation method that is continuous in both the potential and forces and incorporates the pairwise sum considerations stressed by Wolf \textit{et al.}\cite{Wolf99} |
| 285 |
> |
This new Shifted-Force potential is similar to equation |
| 286 |
> |
\ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are |
| 287 |
> |
two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term |
| 288 |
> |
from equation \ref{eq:shiftingForm} is equal to equation |
| 289 |
> |
\ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This |
| 290 |
> |
term is not present in the Zahn potential, resulting in a |
| 291 |
> |
discontinuity as particles cross $R_\textrm{c}$. Second, the sign of |
| 292 |
> |
the derivative portion is different. The constant $v_\textrm{c}$ term |
| 293 |
> |
is not going to have a presence in the forces after performing the |
| 294 |
> |
derivative, but the negative sign does effect the derivative. In |
| 295 |
> |
fact, it introduces a discontinuity in the forces at the cutoff, |
| 296 |
> |
because the force function is shifted in the wrong direction and |
| 297 |
> |
doesn't cross zero at $R_\textrm{c}$. Thus, these alterations make |
| 298 |
> |
for an electrostatic summation method that is continuous in both the |
| 299 |
> |
potential and forces and incorporates the pairwise sum considerations |
| 300 |
> |
stressed by Wolf \textit{et al.}\cite{Wolf99} |
| 301 |
|
|
| 302 |
|
\section{Methods} |
| 303 |
|
|
| 304 |
|
\subsection{What Qualities are Important?}\label{sec:Qualities} |
| 305 |
< |
In classical molecular mechanics simulations, there are two primary techniques utilized to obtain information about the system of interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these techniques utilize pairwise summations of interactions between particle sites, but they use these summations in different ways. |
| 305 |
> |
In classical molecular mechanics simulations, there are two primary |
| 306 |
> |
techniques utilized to obtain information about the system of |
| 307 |
> |
interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these |
| 308 |
> |
techniques utilize pairwise summations of interactions between |
| 309 |
> |
particle sites, but they use these summations in different ways. |
| 310 |
|
|
| 311 |
< |
In MC, the potential energy difference between two subsequent configurations dictates the progression of MC sampling. Going back to the origins of this method, the Canonical ensemble acceptance criteria laid out by Metropolis \textit{et al.} states that a subsequent configuration is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where $\xi$ is a random number between 0 and 1.\cite{Metropolis53} Maintaining a consistent $\Delta E$ when using an alternate method for handling the long-range electrostatics ensures proper sampling within the ensemble. |
| 311 |
> |
In MC, the potential energy difference between two subsequent |
| 312 |
> |
configurations dictates the progression of MC sampling. Going back to |
| 313 |
> |
the origins of this method, the Canonical ensemble acceptance criteria |
| 314 |
> |
laid out by Metropolis \textit{et al.} states that a subsequent |
| 315 |
> |
configuration is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta |
| 316 |
> |
E/kT)$, where $\xi$ is a random number between 0 and |
| 317 |
> |
1.\cite{Metropolis53} Maintaining a consistent $\Delta E$ when using |
| 318 |
> |
an alternate method for handling the long-range electrostatics ensures |
| 319 |
> |
proper sampling within the ensemble. |
| 320 |
|
|
| 321 |
< |
In MD, the derivative of the potential directs how the system will progress in time. Consequently, the force and torque vectors on each body in the system dictate how it develops as a whole. If the magnitude and direction of these vectors are similar when using alternate electrostatic summation techniques, the dynamics in the near term will be indistinguishable. Because error in MD calculations is cumulative, one should expect greater deviation in the long term trajectories with greater differences in these vectors between configurations using different long-range electrostatics. |
| 321 |
> |
In MD, the derivative of the potential directs how the system will |
| 322 |
> |
progress in time. Consequently, the force and torque vectors on each |
| 323 |
> |
body in the system dictate how it develops as a whole. If the |
| 324 |
> |
magnitude and direction of these vectors are similar when using |
| 325 |
> |
alternate electrostatic summation techniques, the dynamics in the near |
| 326 |
> |
term will be indistinguishable. Because error in MD calculations is |
| 327 |
> |
cumulative, one should expect greater deviation in the long term |
| 328 |
> |
trajectories with greater differences in these vectors between |
| 329 |
> |
configurations using different long-range electrostatics. |
| 330 |
|
|
| 331 |
|
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 332 |
< |
Evaluation of the pairwise summation techniques (outlined in section \ref{sec:ESMethods}) for use in MC simulations was performed through study of the energy differences between conformations. Considering the SPME results to be the correct or desired behavior, ideal performance of a tested method was taken to be agreement between the energy differences calculated. Linear least squares regression of the $\Delta E$ values between configurations using SPME against $\Delta E$ values using tested methods provides a quantitative comparison of this agreement. Unitary results for both the correlation and correlation coefficient for these regressions indicate equivalent energetic results between the methods. The correlation is the slope of the plotted data while the correlation coefficient ($R^2$) is a measure of the of the data scatter around the fitted line and tells about the quality of the fit (Fig. \ref{fig:linearFit}). |
| 332 |
> |
Evaluation of the pairwise summation techniques (outlined in section |
| 333 |
> |
\ref{sec:ESMethods}) for use in MC simulations was performed through |
| 334 |
> |
study of the energy differences between conformations. Considering |
| 335 |
> |
the SPME results to be the correct or desired behavior, ideal |
| 336 |
> |
performance of a tested method was taken to be agreement between the |
| 337 |
> |
energy differences calculated. Linear least squares regression of the |
| 338 |
> |
$\Delta E$ values between configurations using SPME against $\Delta E$ |
| 339 |
> |
values using tested methods provides a quantitative comparison of this |
| 340 |
> |
agreement. Unitary results for both the correlation and correlation |
| 341 |
> |
coefficient for these regressions indicate equivalent energetic |
| 342 |
> |
results between the methods. The correlation is the slope of the |
| 343 |
> |
plotted data while the correlation coefficient ($R^2$) is a measure of |
| 344 |
> |
the of the data scatter around the fitted line and tells about the |
| 345 |
> |
quality of the fit (Fig. \ref{fig:linearFit}). |
| 346 |
|
|
| 347 |
|
\begin{figure} |
| 348 |
|
\centering |
| 361 |
|
\ref{sec:EnergyResults}. |
| 362 |
|
|
| 363 |
|
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
| 364 |
< |
Evaluation of the pairwise methods (outlined in section \ref{sec:ESMethods}) for use in MD simulations was performed through comparison of the force and torque vectors obtained with those from SPME. Both the magnitude and the direction of these vectors on each of the bodies in the system were analyzed. For the magnitude of these vectors, linear least squares regression analysis can be performed as described previously for comparing $\Delta E$ values. Instead of a single value between two system configurations, there is a value for each particle in each configuration. For a system of 1000 water molecules and 40 ions, there are 1040 force vectors and 1000 torque vectors. With 500 configurations, this results in 520,000 force and 500,000 torque vector comparisons samples for each system type. |
| 364 |
> |
Evaluation of the pairwise methods (outlined in section |
| 365 |
> |
\ref{sec:ESMethods}) for use in MD simulations was performed through |
| 366 |
> |
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. |
| 376 |
|
|
| 377 |
< |
The force and torque vector directions were investigated through measurement of the angle ($\theta$) formed between those from the particular method and those from SPME |
| 377 |
> |
The force and torque vector directions were investigated through |
| 378 |
> |
measurement of the angle ($\theta$) formed between those from the |
| 379 |
> |
particular method and those from SPME |
| 380 |
|
\begin{equation} |
| 381 |
|
\theta_F = \frac{\vec{F}_\textrm{SPME}}{|\vec{F}_\textrm{SPME}|}\cdot\frac{\vec{F}_\textrm{Method}}{|\vec{F}_\textrm{Method}|}. |
| 382 |
|
\end{equation} |
| 383 |
< |
Each of these $\theta$ values was accumulated in a distribution function, weighted by the area on the unit sphere. Non-linear fits were used to measure the shape of the resulting distributions. |
| 383 |
> |
Each of these $\theta$ values was accumulated in a distribution |
| 384 |
> |
function, weighted by the area on the unit sphere. Non-linear fits |
| 385 |
> |
were used to measure the shape of the resulting distributions. |
| 386 |
|
|
| 387 |
|
\begin{figure} |
| 388 |
|
\centering |
| 391 |
|
\label{fig:gaussian} |
| 392 |
|
\end{figure} |
| 393 |
|
|
| 394 |
< |
Figure \ref{fig:gaussian} shows an example distribution with applied non-linear fits. The solid line is a Gaussian profile, while the dotted line is a Voigt profile, a convolution of a Gaussian and a Lorentzian. Since this distribution is a measure of angular error between two different electrostatic summation methods, there is particular reason for the profile to adhere to a specific shape. Because of this and the Gaussian profile's more statistically meaningful properties, Gaussian fits was used to compare all the tested methods. The variance ($\sigma^2$) was extracted from each of these fits and was used to compare distribution widths. Values of $\sigma^2$ near zero indicate vector directions indistinguishable from those calculated when using SPME. |
| 394 |
> |
Figure \ref{fig:gaussian} shows an example distribution with applied |
| 395 |
> |
non-linear fits. The solid line is a Gaussian profile, while the |
| 396 |
> |
dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 397 |
> |
Lorentzian. Since this distribution is a measure of angular error |
| 398 |
> |
between two different electrostatic summation methods, there is |
| 399 |
> |
particular reason for the profile to adhere to a specific shape. |
| 400 |
> |
Because of this and the Gaussian profile's more statistically |
| 401 |
> |
meaningful properties, Gaussian fits was used to compare all the |
| 402 |
> |
tested methods. The variance ($\sigma^2$) was extracted from each of |
| 403 |
> |
these fits and was used to compare distribution widths. Values of |
| 404 |
> |
$\sigma^2$ near zero indicate vector directions indistinguishable from |
| 405 |
> |
those calculated when using SPME. |
| 406 |
|
|
| 407 |
|
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 408 |
< |
Evaluation of the long-time dynamics of charged systems was performed by considering the NaCl crystal system while using a subset of the best performing pairwise methods. The NaCl crystal was chosen to avoid possible complications involving the propagation techniques of orientational motion in molecular systems. To enhance the atomic motion, these crystals were equilibrated at 1000 K, near the experimental $T_m$ for NaCl. Simulations were performed under the microcanonical ensemble, and velocity autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each of the trajectories, |
| 408 |
> |
Evaluation of the long-time dynamics of charged systems was performed |
| 409 |
> |
by considering the NaCl crystal system while using a subset of the |
| 410 |
> |
best performing pairwise methods. The NaCl crystal was chosen to |
| 411 |
> |
avoid possible complications involving the propagation techniques of |
| 412 |
> |
orientational motion in molecular systems. To enhance the atomic |
| 413 |
> |
motion, these crystals were equilibrated at 1000 K, near the |
| 414 |
> |
experimental $T_m$ for NaCl. Simulations were performed under the |
| 415 |
> |
microcanonical ensemble, and velocity autocorrelation functions |
| 416 |
> |
(Eq. \ref{eq:vCorr}) were computed for each of the trajectories, |
| 417 |
|
\begin{equation} |
| 418 |
|
C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle. |
| 419 |
|
\label{eq:vCorr} |
| 420 |
|
\end{equation} |
| 421 |
< |
Velocity autocorrelation functions require detailed short time data and long trajectories for good statistics, thus velocity information was saved every 5 fs over 100 ps trajectories. The power spectrum ($I(\omega)$) is obtained via Fourier transform of the autocorrelation function |
| 421 |
> |
Velocity autocorrelation functions require detailed short time data |
| 422 |
> |
and long trajectories for good statistics, thus velocity information |
| 423 |
> |
was saved every 5 fs over 100 ps trajectories. The power spectrum |
| 424 |
> |
($I(\omega)$) is obtained via Fourier transform of the autocorrelation |
| 425 |
> |
function |
| 426 |
|
\begin{equation} |
| 427 |
|
I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 428 |
|
\label{eq:powerSpec} |
| 430 |
|
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. |
| 431 |
|
|
| 432 |
|
\subsection{Representative Simulations}\label{sec:RepSims} |
| 433 |
< |
A variety of common and representative simulations were analyzed to determine the relative effectiveness of the pairwise summation techniques in reproducing the energetics and dynamics exhibited by SPME. The studied systems were as follows: |
| 433 |
> |
A variety of common and representative simulations were analyzed to |
| 434 |
> |
determine the relative effectiveness of the pairwise summation |
| 435 |
> |
techniques in reproducing the energetics and dynamics exhibited by |
| 436 |
> |
SPME. The studied systems were as follows: |
| 437 |
|
\begin{enumerate} |
| 438 |
|
\item Liquid Water |
| 439 |
|
\item Crystalline Water (Ice I$_\textrm{c}$) |
| 443 |
|
\item High Ionic Strength Solution of NaCl in Water |
| 444 |
|
\item 6 \AA\ Radius Sphere of Argon in Water |
| 445 |
|
\end{enumerate} |
| 446 |
< |
By utilizing the pairwise techniques (outlined in section \ref{sec:ESMethods}) in systems composed entirely of neutral groups, charged particles, and mixtures of the two, we can comment on possible system dependence and/or universal applicability of the techniques. |
| 446 |
> |
By utilizing the pairwise techniques (outlined in section |
| 447 |
> |
\ref{sec:ESMethods}) in systems composed entirely of neutral groups, |
| 448 |
> |
charged particles, and mixtures of the two, we can comment on possible |
| 449 |
> |
system dependence and/or universal applicability of the techniques. |
| 450 |
|
|
| 451 |
< |
Generation of the system configurations was dependent on the system type. For the solid and liquid water configurations, configuration snapshots were taken at regular intervals from higher temperature 1000 SPC/E water molecule trajectories and each equilibrated individually. The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- ions and were selected and equilibrated in the same fashion as the water systems. For the low and high ionic strength NaCl solutions, 4 and 40 ions were first solvated in a 1000 water molecule boxes respectively. Ion and water positions were then randomly swapped, and the resulting configurations were again equilibrated individually. Finally, for the Argon/Water "charge void" systems, the identities of all the SPC/E waters within 6 \AA\ of the center of the equilibrated water configurations were converted to argon (Fig. \ref{fig:argonSlice}). |
| 451 |
> |
Generation of the system configurations was dependent on the system |
| 452 |
> |
type. For the solid and liquid water configurations, configuration |
| 453 |
> |
snapshots were taken at regular intervals from higher temperature 1000 |
| 454 |
> |
SPC/E water molecule trajectories and each equilibrated individually. |
| 455 |
> |
The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- |
| 456 |
> |
ions and were selected and equilibrated in the same fashion as the |
| 457 |
> |
water systems. For the low and high ionic strength NaCl solutions, 4 |
| 458 |
> |
and 40 ions were first solvated in a 1000 water molecule boxes |
| 459 |
> |
respectively. Ion and water positions were then randomly swapped, and |
| 460 |
> |
the resulting configurations were again equilibrated individually. |
| 461 |
> |
Finally, for the Argon/Water "charge void" systems, the identities of |
| 462 |
> |
all the SPC/E waters within 6 \AA\ of the center of the equilibrated |
| 463 |
> |
water configurations were converted to argon |
| 464 |
> |
(Fig. \ref{fig:argonSlice}). |
| 465 |
|
|
| 466 |
|
\begin{figure} |
| 467 |
|
\centering |
| 471 |
|
\end{figure} |
| 472 |
|
|
| 473 |
|
\subsection{Electrostatic Summation Methods}\label{sec:ESMethods} |
| 474 |
< |
Electrostatic summation method comparisons were performed using SPME, the Shifted-Potential and Shifted-Force methods - both with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak, moderate, and strong damping respectively), reaction field with an infinite dielectric constant, and an unmodified cutoff. Group-based cutoffs with a fifth-order polynomial switching function were necessary for the reaction field simulations and were utilized in the SP, SF, and pure cutoff methods for comparison to the standard lack of group-based cutoffs with a hard truncation. The SPME calculations were performed using the TINKER implementation of SPME,\cite{Ponder87} while all other method calculations were performed using the OOPSE molecular mechanics package.\cite{Meineke05} |
| 474 |
> |
Electrostatic summation method comparisons were performed using SPME, |
| 475 |
> |
the Shifted-Potential and Shifted-Force methods - both with damping |
| 476 |
> |
parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak, |
| 477 |
> |
moderate, and strong damping respectively), reaction field with an |
| 478 |
> |
infinite dielectric constant, and an unmodified cutoff. Group-based |
| 479 |
> |
cutoffs with a fifth-order polynomial switching function were |
| 480 |
> |
necessary for the reaction field simulations and were utilized in the |
| 481 |
> |
SP, SF, and pure cutoff methods for comparison to the standard lack of |
| 482 |
> |
group-based cutoffs with a hard truncation. The SPME calculations |
| 483 |
> |
were performed using the TINKER implementation of SPME,\cite{Ponder87} |
| 484 |
> |
while all other method calculations were performed using the OOPSE |
| 485 |
> |
molecular mechanics package.\cite{Meineke05} |
| 486 |
|
|
| 487 |
< |
These methods were additionally evaluated with three different cutoff radii (9, 12, and 15 \AA) to investigate possible cutoff radius dependence. It should be noted that the damping parameter chosen in SPME, or so called ``Ewald Coefficient", has a significant effect on the energies and forces calculated. Typical molecular mechanics packages default this to a value dependent on the cutoff radius and a tolerance (typically less than $1 \times 10^{-4}$ kcal/mol). Smaller tolerances are typically associated with increased accuracy in the real-space portion of the summation.\cite{Essmann95} The default TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 487 |
> |
These methods were additionally evaluated with three different cutoff |
| 488 |
> |
radii (9, 12, and 15 \AA) to investigate possible cutoff radius |
| 489 |
> |
dependence. It should be noted that the damping parameter chosen in |
| 490 |
> |
SPME, or so called ``Ewald Coefficient", has a significant effect on |
| 491 |
> |
the energies and forces calculated. Typical molecular mechanics |
| 492 |
> |
packages default this to a value dependent on the cutoff radius and a |
| 493 |
> |
tolerance (typically less than $1 \times 10^{-4}$ kcal/mol). Smaller |
| 494 |
> |
tolerances are typically associated with increased accuracy in the |
| 495 |
> |
real-space portion of the summation.\cite{Essmann95} The default |
| 496 |
> |
TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME |
| 497 |
> |
calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and |
| 498 |
> |
0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 499 |
|
|
| 500 |
|
\section{Results and Discussion} |
| 501 |
|
|
| 502 |
|
\subsection{Configuration Energy Differences}\label{sec:EnergyResults} |
| 503 |
< |
In order to evaluate the performance of the pairwise electrostatic summation methods for Monte Carlo simulations, the energy differences between configurations were compared to the values obtained when using SPME. The results for the subsequent regression analysis are shown in figure \ref{fig:delE}. |
| 503 |
> |
In order to evaluate the performance of the pairwise electrostatic |
| 504 |
> |
summation methods for Monte Carlo simulations, the energy differences |
| 505 |
> |
between configurations were compared to the values obtained when using |
| 506 |
> |
SPME. The results for the subsequent regression analysis are shown in |
| 507 |
> |
figure \ref{fig:delE}. |
| 508 |
|
|
| 509 |
|
\begin{figure} |
| 510 |
|
\centering |
| 513 |
|
\label{fig:delE} |
| 514 |
|
\end{figure} |
| 515 |
|
|
| 516 |
< |
In this figure, it is apparent that it is unreasonable to expect realistic results using an unmodified cutoff. This is not all that surprising since this results in large energy fluctuations as atoms move in and out of the cutoff radius. These fluctuations can be alleviated to some degree by using group based cutoffs with a switching function.\cite{Steinbach94} The Group Switch Cutoff row doesn't show a significant improvement in this plot because the salt and salt solution systems contain non-neutral groups, see the accompanying supporting information for a comparison where all groups are neutral. |
| 516 |
> |
In this figure, it is apparent that it is unreasonable to expect |
| 517 |
> |
realistic results using an unmodified cutoff. This is not all that |
| 518 |
> |
surprising since this results in large energy fluctuations as atoms |
| 519 |
> |
move in and out of the cutoff radius. These fluctuations can be |
| 520 |
> |
alleviated to some degree by using group based cutoffs with a |
| 521 |
> |
switching function.\cite{Steinbach94} The Group Switch Cutoff row |
| 522 |
> |
doesn't show a significant improvement in this plot because the salt |
| 523 |
> |
and salt solution systems contain non-neutral groups, see the |
| 524 |
> |
accompanying supporting information for a comparison where all groups |
| 525 |
> |
are neutral. |
| 526 |
|
|
| 527 |
< |
Correcting the resulting charged cutoff sphere is one of the purposes of the damped Coulomb summation proposed by Wolf \textit{et al.},\cite{Wolf99} and this correction indeed improves the results as seen in the Shifted-Potental rows. While the undamped case of this method is a significant improvement over the pure cutoff, it still doesn't correlate that well with SPME. Inclusion of potential damping improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows an excellent correlation and quality of fit with the SPME results, particularly with a cutoff radius greater than 12 \AA . Use of a larger damping parameter is more helpful for the shortest cutoff shown, but it has a detrimental effect on simulations with larger cutoffs. In the Shifted-Force sets, increasing damping results in progressively poorer correlation. Overall, the undamped case is the best performing set, as the correlation and quality of fits are consistently superior regardless of the cutoff distance. This result is beneficial in that the undamped case is less computationally prohibitive do to the lack of complimentary error function calculation when performing the electrostatic pair interaction. The reaction field results illustrates some of that method's limitations, primarily that it was developed for use in homogenous systems; although it does provide results that are an improvement over those from an unmodified cutoff. |
| 527 |
> |
Correcting the resulting charged cutoff sphere is one of the purposes |
| 528 |
> |
of the damped Coulomb summation proposed by Wolf \textit{et |
| 529 |
> |
al.},\cite{Wolf99} and this correction indeed improves the results as |
| 530 |
> |
seen in the Shifted-Potental rows. While the undamped case of this |
| 531 |
> |
method is a significant improvement over the pure cutoff, it still |
| 532 |
> |
doesn't correlate that well with SPME. Inclusion of potential damping |
| 533 |
> |
improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 534 |
> |
an excellent correlation and quality of fit with the SPME results, |
| 535 |
> |
particularly with a cutoff radius greater than 12 \AA . Use of a |
| 536 |
> |
larger damping parameter is more helpful for the shortest cutoff |
| 537 |
> |
shown, but it has a detrimental effect on simulations with larger |
| 538 |
> |
cutoffs. In the Shifted-Force sets, increasing damping results in |
| 539 |
> |
progressively poorer correlation. Overall, the undamped case is the |
| 540 |
> |
best performing set, as the correlation and quality of fits are |
| 541 |
> |
consistently superior regardless of the cutoff distance. This result |
| 542 |
> |
is beneficial in that the undamped case is less computationally |
| 543 |
> |
prohibitive do to the lack of complimentary error function calculation |
| 544 |
> |
when performing the electrostatic pair interaction. The reaction |
| 545 |
> |
field results illustrates some of that method's limitations, primarily |
| 546 |
> |
that it was developed for use in homogenous systems; although it does |
| 547 |
> |
provide results that are an improvement over those from an unmodified |
| 548 |
> |
cutoff. |
| 549 |
|
|
| 550 |
|
\subsection{Magnitudes of the Force and Torque Vectors} |
| 551 |
|
|
| 552 |
< |
Evaluation of pairwise methods for use in Molecular Dynamics simulations requires consideration of effects on the forces and torques. Investigation of the force and torque vector magnitudes provides a measure of the strength of these values relative to SPME. Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the force and torque vector magnitude regression results for the accumulated analysis over all the system types. |
| 552 |
> |
Evaluation of pairwise methods for use in Molecular Dynamics |
| 553 |
> |
simulations requires consideration of effects on the forces and |
| 554 |
> |
torques. Investigation of the force and torque vector magnitudes |
| 555 |
> |
provides a measure of the strength of these values relative to SPME. |
| 556 |
> |
Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the |
| 557 |
> |
force and torque vector magnitude regression results for the |
| 558 |
> |
accumulated analysis over all the system types. |
| 559 |
|
|
| 560 |
|
\begin{figure} |
| 561 |
|
\centering |
| 564 |
|
\label{fig:frcMag} |
| 565 |
|
\end{figure} |
| 566 |
|
|
| 567 |
< |
Figure \ref{fig:frcMag}, for the most part, parallels the results seen in the previous $\Delta E$ section. The unmodified cutoff results are poor, but using group based cutoffs and a switching function provides a improvement much more significant than what was seen with $\Delta E$. Looking at the Shifted-Potential sets, the slope and $R^2$ improve with the use of damping to an optimal result of 0.2 \AA $^{-1}$ for the 12 and 15 \AA\ cutoffs. Further increases in damping, while beneficial for simulations with a cutoff radius of 9 \AA\ , is detrimental to simulations with larger cutoff radii. The undamped Shifted-Force method gives forces in line with those obtained using SPME, and use of a damping function results in minor improvement. The reaction field results are surprisingly good, considering the poor quality of the fits for the $\Delta E$ results. There is still a considerable degree of scatter in the data, but it correlates well in general. To be fair, we again note that the reaction field calculations do not encompass NaCl crystal and melt systems, so these results are partly biased towards conditions in which the method performs more favorably. |
| 567 |
> |
Figure \ref{fig:frcMag}, for the most part, parallels the results seen |
| 568 |
> |
in the previous $\Delta E$ section. The unmodified cutoff results are |
| 569 |
> |
poor, but using group based cutoffs and a switching function provides |
| 570 |
> |
a improvement much more significant than what was seen with $\Delta |
| 571 |
> |
E$. Looking at the Shifted-Potential sets, the slope and $R^2$ |
| 572 |
> |
improve with the use of damping to an optimal result of 0.2 \AA |
| 573 |
> |
$^{-1}$ for the 12 and 15 \AA\ cutoffs. Further increases in damping, |
| 574 |
> |
while beneficial for simulations with a cutoff radius of 9 \AA\ , is |
| 575 |
> |
detrimental to simulations with larger cutoff radii. The undamped |
| 576 |
> |
Shifted-Force method gives forces in line with those obtained using |
| 577 |
> |
SPME, and use of a damping function results in minor improvement. The |
| 578 |
> |
reaction field results are surprisingly good, considering the poor |
| 579 |
> |
quality of the fits for the $\Delta E$ results. There is still a |
| 580 |
> |
considerable degree of scatter in the data, but it correlates well in |
| 581 |
> |
general. To be fair, we again note that the reaction field |
| 582 |
> |
calculations do not encompass NaCl crystal and melt systems, so these |
| 583 |
> |
results are partly biased towards conditions in which the method |
| 584 |
> |
performs more favorably. |
| 585 |
|
|
| 586 |
|
\begin{figure} |
| 587 |
|
\centering |
| 590 |
|
\label{fig:trqMag} |
| 591 |
|
\end{figure} |
| 592 |
|
|
| 593 |
< |
To evaluate the torque vector magnitudes, the data set from which values are drawn is limited to rigid molecules in the systems (i.e. water molecules). In spite of this smaller sampling pool, the torque vector magnitude results in figure \ref{fig:trqMag} are still similar to those seen for the forces; however, they more clearly show the improved behavior that comes with increasing the cutoff radius. Moderate damping is beneficial to the Shifted-Potential and helpful yet possibly unnecessary with the Shifted-Force method, and they also show that over-damping adversely effects all cutoff radii rather than showing an improvement for systems with short cutoffs. The reaction field method performs well when calculating the torques, better than the Shifted Force method over this limited data set. |
| 593 |
> |
To evaluate the torque vector magnitudes, the data set from which |
| 594 |
> |
values are drawn is limited to rigid molecules in the systems |
| 595 |
> |
(i.e. water molecules). In spite of this smaller sampling pool, the |
| 596 |
> |
torque vector magnitude results in figure \ref{fig:trqMag} are still |
| 597 |
> |
similar to those seen for the forces; however, they more clearly show |
| 598 |
> |
the improved behavior that comes with increasing the cutoff radius. |
| 599 |
> |
Moderate damping is beneficial to the Shifted-Potential and helpful |
| 600 |
> |
yet possibly unnecessary with the Shifted-Force method, and they also |
| 601 |
> |
show that over-damping adversely effects all cutoff radii rather than |
| 602 |
> |
showing an improvement for systems with short cutoffs. The reaction |
| 603 |
> |
field method performs well when calculating the torques, better than |
| 604 |
> |
the Shifted Force method over this limited data set. |
| 605 |
|
|
| 606 |
|
\subsection{Directionality of the Force and Torque Vectors} |
| 607 |
|
|
| 608 |
< |
Having force and torque vectors with magnitudes that are well correlated to SPME is good, but if they are not pointing in the proper direction the results will be incorrect. These vector directions were investigated through measurement of the angle formed between them and those from SPME. The results (Fig. \ref{fig:frcTrqAng}) are compared through the variance ($\sigma^2$) of the Gaussian fits of the angle error distributions of the combined set over all system types. |
| 608 |
> |
Having force and torque vectors with magnitudes that are well |
| 609 |
> |
correlated to SPME is good, but if they are not pointing in the proper |
| 610 |
> |
direction the results will be incorrect. These vector directions were |
| 611 |
> |
investigated through measurement of the angle formed between them and |
| 612 |
> |
those from SPME. The results (Fig. \ref{fig:frcTrqAng}) are compared |
| 613 |
> |
through the variance ($\sigma^2$) of the Gaussian fits of the angle |
| 614 |
> |
error distributions of the combined set over all system types. |
| 615 |
|
|
| 616 |
|
\begin{figure} |
| 617 |
|
\centering |
| 620 |
|
\label{fig:frcTrqAng} |
| 621 |
|
\end{figure} |
| 622 |
|
|
| 623 |
< |
Both the force and torque $\sigma^2$ results from the analysis of the total accumulated system data are tabulated in figure \ref{fig:frcTrqAng}. All of the sets, aside from the over-damped case show the improvement afforded by choosing a longer simulation cutoff. Increasing the cutoff from 9 to 12 \AA\ typically results in a halving of the distribution widths, with a similar improvement going from 12 to 15 \AA . The undamped Shifted-Force, Group Based Cutoff, and Reaction Field methods all do equivalently well at capturing the direction of both the force and torque vectors. Using damping improves the angular behavior significantly for the Shifted-Potential and moderately for the Shifted-Force methods. Increasing the damping too far is destructive for both methods, particularly to the torque vectors. Again it is important to recognize that the force vectors cover all particles in the systems, while torque vectors are only available for neutral molecular groups. Damping appears to have a more beneficial effect on non-neutral bodies, and this observation is investigated further in the accompanying supporting information. |
| 623 |
> |
Both the force and torque $\sigma^2$ results from the analysis of the |
| 624 |
> |
total accumulated system data are tabulated in figure |
| 625 |
> |
\ref{fig:frcTrqAng}. All of the sets, aside from the over-damped case |
| 626 |
> |
show the improvement afforded by choosing a longer simulation cutoff. |
| 627 |
> |
Increasing the cutoff from 9 to 12 \AA\ typically results in a halving |
| 628 |
> |
of the distribution widths, with a similar improvement going from 12 |
| 629 |
> |
to 15 \AA . The undamped Shifted-Force, Group Based Cutoff, and |
| 630 |
> |
Reaction Field methods all do equivalently well at capturing the |
| 631 |
> |
direction of both the force and torque vectors. Using damping |
| 632 |
> |
improves the angular behavior significantly for the Shifted-Potential |
| 633 |
> |
and moderately for the Shifted-Force methods. Increasing the damping |
| 634 |
> |
too far is destructive for both methods, particularly to the torque |
| 635 |
> |
vectors. Again it is important to recognize that the force vectors |
| 636 |
> |
cover all particles in the systems, while torque vectors are only |
| 637 |
> |
available for neutral molecular groups. Damping appears to have a |
| 638 |
> |
more beneficial effect on non-neutral bodies, and this observation is |
| 639 |
> |
investigated further in the accompanying supporting information. |
| 640 |
|
|
| 641 |
|
\begin{table}[htbp] |
| 642 |
|
\centering |
| 671 |
|
\label{tab:groupAngle} |
| 672 |
|
\end{table} |
| 673 |
|
|
| 674 |
< |
Although not discussed previously, group based cutoffs can be applied to both the Shifted-Potential and Shifted-Force methods. Use off a switching function corrects for the discontinuities that arise when atoms of a group exit the cutoff before the group's center of mass. Though there are no significant benefit or drawbacks observed in $\Delta E$ and vector magnitude results when doing this, there is a measurable improvement in the vector angle results. Table \ref{tab:groupAngle} shows the angular variance values obtained using group based cutoffs and a switching function alongside the standard results seen in figure \ref{fig:frcTrqAng} for comparison purposes. The Shifted-Potential shows much narrower angular distributions for both the force and torque vectors when using an $\alpha$ of 0.2 \AA$^{-1}$ or less, while Shifted-Force shows improvements in the undamped and lightly damped cases. Thus, by calculating the electrostatic interactions in terms of molecular pairs rather than atomic pairs, the direction of the force and torque vectors are determined more accurately. |
| 674 |
> |
Although not discussed previously, group based cutoffs can be applied |
| 675 |
> |
to both the Shifted-Potential and Shifted-Force methods. Use off a |
| 676 |
> |
switching function corrects for the discontinuities that arise when |
| 677 |
> |
atoms of a group exit the cutoff before the group's center of mass. |
| 678 |
> |
Though there are no significant benefit or drawbacks observed in |
| 679 |
> |
$\Delta E$ and vector magnitude results when doing this, there is a |
| 680 |
> |
measurable improvement in the vector angle results. Table |
| 681 |
> |
\ref{tab:groupAngle} shows the angular variance values obtained using |
| 682 |
> |
group based cutoffs and a switching function alongside the standard |
| 683 |
> |
results seen in figure \ref{fig:frcTrqAng} for comparison purposes. |
| 684 |
> |
The Shifted-Potential shows much narrower angular distributions for |
| 685 |
> |
both the force and torque vectors when using an $\alpha$ of 0.2 |
| 686 |
> |
\AA$^{-1}$ or less, while Shifted-Force shows improvements in the |
| 687 |
> |
undamped and lightly damped cases. Thus, by calculating the |
| 688 |
> |
electrostatic interactions in terms of molecular pairs rather than |
| 689 |
> |
atomic pairs, the direction of the force and torque vectors are |
| 690 |
> |
determined more accurately. |
| 691 |
|
|
| 692 |
< |
One additional trend to recognize in table \ref{tab:groupAngle} is that the $\sigma^2$ values for both Shifted-Potential and Shifted-Force converge as $\alpha$ increases, something that is easier to see when using group based cutoffs. Looking back on figures \ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this behavior clearly at large $\alpha$ and cutoff values. The reason for this is that the complimentary error function inserted into the potential weakens the electrostatic interaction as $\alpha$ increases. Thus, at larger values of $\alpha$, both the summation method types progress toward non-interacting functions, so care is required in choosing large damping functions lest one generate an undesirable loss in the pair interaction. Kast \textit{et al.} developed a method for choosing appropriate $\alpha$ values for these types of electrostatic summation methods by fitting to $g(r)$ data, and their methods indicate optimal values of 0.34, 0.25, and 0.16 \AA$^{-1}$ for cutoff values of 9, 12, and 15 \AA\ respectively.\cite{Kast03} These appear to be reasonable choices to obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on these findings, choices this high would introduce error in the molecular torques, particularly for the shorter cutoffs. Based on the above findings, empirical damping up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably unnecessary when using the Shifted-Force method. |
| 692 |
> |
One additional trend to recognize in table \ref{tab:groupAngle} is |
| 693 |
> |
that the $\sigma^2$ values for both Shifted-Potential and |
| 694 |
> |
Shifted-Force converge as $\alpha$ increases, something that is easier |
| 695 |
> |
to see when using group based cutoffs. Looking back on figures |
| 696 |
> |
\ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this |
| 697 |
> |
behavior clearly at large $\alpha$ and cutoff values. The reason for |
| 698 |
> |
this is that the complimentary error function inserted into the |
| 699 |
> |
potential weakens the electrostatic interaction as $\alpha$ increases. |
| 700 |
> |
Thus, at larger values of $\alpha$, both the summation method types |
| 701 |
> |
progress toward non-interacting functions, so care is required in |
| 702 |
> |
choosing large damping functions lest one generate an undesirable loss |
| 703 |
> |
in the pair interaction. Kast \textit{et al.} developed a method for |
| 704 |
> |
choosing appropriate $\alpha$ values for these types of electrostatic |
| 705 |
> |
summation methods by fitting to $g(r)$ data, and their methods |
| 706 |
> |
indicate optimal values of 0.34, 0.25, and 0.16 \AA$^{-1}$ for cutoff |
| 707 |
> |
values of 9, 12, and 15 \AA\ respectively.\cite{Kast03} These appear |
| 708 |
> |
to be reasonable choices to obtain proper MC behavior |
| 709 |
> |
(Fig. \ref{fig:delE}); however, based on these findings, choices this |
| 710 |
> |
high would introduce error in the molecular torques, particularly for |
| 711 |
> |
the shorter cutoffs. Based on the above findings, empirical damping |
| 712 |
> |
up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably |
| 713 |
> |
unnecessary when using the Shifted-Force method. |
| 714 |
|
|
| 715 |
|
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 716 |
|
|
| 717 |
< |
In the previous studies using a Shifted-Force variant of the damped Wolf coulomb potential, the structure and dynamics of water were investigated rather extensively.\cite{Zahn02,Kast03} Their results indicated that the damped Shifted-Force method results in properties very similar to those obtained when using the Ewald summation. Considering the statistical results shown above, the good performance of this method is not that surprising. Rather than consider the same systems and simply recapitulate their results, we decided to look at the solid state dynamical behavior obtained using the best performing summation methods from the above results. |
| 717 |
> |
In the previous studies using a Shifted-Force variant of the damped |
| 718 |
> |
Wolf coulomb potential, the structure and dynamics of water were |
| 719 |
> |
investigated rather extensively.\cite{Zahn02,Kast03} Their results |
| 720 |
> |
indicated that the damped Shifted-Force method results in properties |
| 721 |
> |
very similar to those obtained when using the Ewald summation. |
| 722 |
> |
Considering the statistical results shown above, the good performance |
| 723 |
> |
of this method is not that surprising. Rather than consider the same |
| 724 |
> |
systems and simply recapitulate their results, we decided to look at |
| 725 |
> |
the solid state dynamical behavior obtained using the best performing |
| 726 |
> |
summation methods from the above results. |
| 727 |
|
|
| 728 |
|
\begin{figure} |
| 729 |
|
\centering |
| 732 |
|
\label{fig:methodPS} |
| 733 |
|
\end{figure} |
| 734 |
|
|
| 735 |
< |
Figure \ref{fig:methodPS} shows the power spectra for the NaCl crystals (from averaged Na and Cl ion velocity autocorrelation functions) using the stated electrostatic summation methods. While high frequency peaks of all the spectra overlap, showing the same general features, the low frequency region shows how the summation methods differ. Considering the low-frequency inset (expanded in the upper frame of figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the correlated motions are blue-shifted when using undamped or weakly damped Shifted-Force. When using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the Shifted-Force and Shifted-Potential methods give near identical correlated motion behavior as the Ewald method (which has a damping value of 0.3119). The damping acts as a distance dependent Gaussian screening of the point charges for the pairwise summation methods. This weakening of the electrostatic interaction with distance explains why the long-ranged correlated motions are at lower frequencies for the moderately damped methods than for undamped or weakly damped methods. To see this effect more clearly, we show how damping strength affects a simple real-space electrostatic potential, |
| 735 |
> |
Figure \ref{fig:methodPS} shows the power spectra for the NaCl |
| 736 |
> |
crystals (from averaged Na and Cl ion velocity autocorrelation |
| 737 |
> |
functions) using the stated electrostatic summation methods. While |
| 738 |
> |
high frequency peaks of all the spectra overlap, showing the same |
| 739 |
> |
general features, the low frequency region shows how the summation |
| 740 |
> |
methods differ. Considering the low-frequency inset (expanded in the |
| 741 |
> |
upper frame of figure \ref{fig:dampInc}), at frequencies below 100 |
| 742 |
> |
cm$^{-1}$, the correlated motions are blue-shifted when using undamped |
| 743 |
> |
or weakly damped Shifted-Force. When using moderate damping ($\alpha |
| 744 |
> |
= 0.2$ \AA$^{-1}$) both the Shifted-Force and Shifted-Potential |
| 745 |
> |
methods give near identical correlated motion behavior as the Ewald |
| 746 |
> |
method (which has a damping value of 0.3119). The damping acts as a |
| 747 |
> |
distance dependent Gaussian screening of the point charges for the |
| 748 |
> |
pairwise summation methods. This weakening of the electrostatic |
| 749 |
> |
interaction with distance explains why the long-ranged correlated |
| 750 |
> |
motions are at lower frequencies for the moderately damped methods |
| 751 |
> |
than for undamped or weakly damped methods. To see this effect more |
| 752 |
> |
clearly, we show how damping strength affects a simple real-space |
| 753 |
> |
electrostatic potential, |
| 754 |
|
\begin{equation} |
| 755 |
|
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), |
| 756 |
|
\end{equation} |
| 757 |
< |
where $S(r)$ is a switching function that smoothly zeroes the potential at the cutoff radius. Figure \ref{fig:dampInc} shows how the low frequency motions are dependent on the damping used in the direct electrostatic sum. As the damping increases, the peaks drop to lower frequencies. Incidentally, use of an $\alpha$ of 0.25 \AA$^{-1}$ on a simple electrostatic summation results in low frequency correlated dynamics equivalent to a simulation using SPME. When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks shift to higher frequency in exponential fashion. Though not shown, the spectrum for the simple undamped electrostatic potential is blue-shifted such that the lowest frequency peak resides near 325 cm$^{-1}$. In light of these results, the undamped Shifted-Force method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite respectable; however, it appears as though moderate damping is required for accurate reproduction of crystal dynamics. |
| 757 |
> |
where $S(r)$ is a switching function that smoothly zeroes the |
| 758 |
> |
potential at the cutoff radius. Figure \ref{fig:dampInc} shows how |
| 759 |
> |
the low frequency motions are dependent on the damping used in the |
| 760 |
> |
direct electrostatic sum. As the damping increases, the peaks drop to |
| 761 |
> |
lower frequencies. Incidentally, use of an $\alpha$ of 0.25 |
| 762 |
> |
\AA$^{-1}$ on a simple electrostatic summation results in low |
| 763 |
> |
frequency correlated dynamics equivalent to a simulation using SPME. |
| 764 |
> |
When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks |
| 765 |
> |
shift to higher frequency in exponential fashion. Though not shown, |
| 766 |
> |
the spectrum for the simple undamped electrostatic potential is |
| 767 |
> |
blue-shifted such that the lowest frequency peak resides near 325 |
| 768 |
> |
cm$^{-1}$. In light of these results, the undamped Shifted-Force |
| 769 |
> |
method producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is |
| 770 |
> |
quite respectable; however, it appears as though moderate damping is |
| 771 |
> |
required for accurate reproduction of crystal dynamics. |
| 772 |
|
\begin{figure} |
| 773 |
|
\centering |
| 774 |
|
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 778 |
|
|
| 779 |
|
\section{Conclusions} |
| 780 |
|
|
| 781 |
< |
This investigation of pairwise electrostatic summation techniques shows that there are viable and more computationally efficient electrostatic summation techniques than the Ewald summation, chiefly methods derived from the damped Coulombic sum originally proposed by Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the Shifted-Force method, reformulated above as equation \ref{eq:SFPot}, shows a remarkable ability to reproduce the energetic and dynamic characteristics exhibited by simulations employing lattice summation techniques. The cumulative energy difference results showed the undamped Shifted-Force and moderately damped Shifted-Potential methods produced results nearly identical to SPME. Similarly for the dynamic features, the undamped or moderately damped Shifted-Force and moderately damped Shifted-Potential methods produce force and torque vector magnitude and directions very similar to the expected values. These results translate into long-time dynamic behavior equivalent to that produced in simulations using SPME. |
| 781 |
> |
This investigation of pairwise electrostatic summation techniques |
| 782 |
> |
shows that there are viable and more computationally efficient |
| 783 |
> |
electrostatic summation techniques than the Ewald summation, chiefly |
| 784 |
> |
methods derived from the damped Coulombic sum originally proposed by |
| 785 |
> |
Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the |
| 786 |
> |
Shifted-Force method, reformulated above as equation \ref{eq:SFPot}, |
| 787 |
> |
shows a remarkable ability to reproduce the energetic and dynamic |
| 788 |
> |
characteristics exhibited by simulations employing lattice summation |
| 789 |
> |
techniques. The cumulative energy difference results showed the |
| 790 |
> |
undamped Shifted-Force and moderately damped Shifted-Potential methods |
| 791 |
> |
produced results nearly identical to SPME. Similarly for the dynamic |
| 792 |
> |
features, the undamped or moderately damped Shifted-Force and |
| 793 |
> |
moderately damped Shifted-Potential methods produce force and torque |
| 794 |
> |
vector magnitude and directions very similar to the expected values. |
| 795 |
> |
These results translate into long-time dynamic behavior equivalent to |
| 796 |
> |
that produced in simulations using SPME. |
| 797 |
|
|
| 798 |
< |
Aside from the computational cost benefit, these techniques have applicability in situations where the use of the Ewald sum can prove problematic. Primary among them is their use in interfacial systems, where the unmodified lattice sum techniques artificially accentuate the periodicity of the system in an undesirable manner. There have been alterations to the standard Ewald techniques, via corrections and reformulations, to compensate for these systems; but the pairwise techniques discussed here require no modifications, making them natural tools to tackle these problems. Additionally, this transferability gives them benefits over other pairwise methods, like reaction field, because estimations of physical properties (e.g. the dielectric constant) are unnecessary. |
| 798 |
> |
Aside from the computational cost benefit, these techniques have |
| 799 |
> |
applicability in situations where the use of the Ewald sum can prove |
| 800 |
> |
problematic. Primary among them is their use in interfacial systems, |
| 801 |
> |
where the unmodified lattice sum techniques artificially accentuate |
| 802 |
> |
the periodicity of the system in an undesirable manner. There have |
| 803 |
> |
been alterations to the standard Ewald techniques, via corrections and |
| 804 |
> |
reformulations, to compensate for these systems; but the pairwise |
| 805 |
> |
techniques discussed here require no modifications, making them |
| 806 |
> |
natural tools to tackle these problems. Additionally, this |
| 807 |
> |
transferability gives them benefits over other pairwise methods, like |
| 808 |
> |
reaction field, because estimations of physical properties (e.g. the |
| 809 |
> |
dielectric constant) are unnecessary. |
| 810 |
|
|
| 811 |
< |
We are not suggesting any flaw with the Ewald sum; in fact, it is the standard by which these simple pairwise sums are judged. However, these results do suggest that in the typical simulations performed today, the Ewald summation may no longer be required to obtain the level of accuracy most researcher have come to expect |
| 811 |
> |
We are not suggesting any flaw with the Ewald sum; in fact, it is the |
| 812 |
> |
standard by which these simple pairwise sums are judged. However, |
| 813 |
> |
these results do suggest that in the typical simulations performed |
| 814 |
> |
today, the Ewald summation may no longer be required to obtain the |
| 815 |
> |
level of accuracy most researcher have come to expect |
| 816 |
|
|
| 817 |
|
\section{Acknowledgments} |
| 818 |
|
\newpage |
