| 140 |
|
\end{split} |
| 141 |
|
\label{eq:EwaldSum} |
| 142 |
|
\end{equation} |
| 143 |
< |
where $\alpha$ is a damping parameter, or separation constant, with |
| 144 |
< |
units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are |
| 145 |
< |
equal to $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the |
| 146 |
< |
dielectric constant of the surrounding medium. The final two terms of |
| 143 |
> |
where $\alpha$ is the damping or convergence parameter with units of |
| 144 |
> |
\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to |
| 145 |
> |
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 146 |
> |
constant of the surrounding medium. The final two terms of |
| 147 |
|
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 148 |
|
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 149 |
|
dipolar term was neglected in early applications in molecular |
| 159 |
|
convergence. In more modern simulations, the simulation boxes have |
| 160 |
|
grown large enough that a real-space cutoff could potentially give |
| 161 |
|
convergent behavior. Indeed, it has often been observed that the |
| 162 |
< |
reciprocal-space portion of the Ewald sum can be vanishingly |
| 163 |
< |
small compared to the real-space portion.\cite{XXX} |
| 162 |
> |
reciprocal-space portion of the Ewald sum can be small and rapidly |
| 163 |
> |
convergent compared to the real-space portion with the choice of small |
| 164 |
> |
$\alpha$.\cite{Karasawa89,Kolafa92} |
| 165 |
|
|
| 166 |
|
\begin{figure} |
| 167 |
|
\centering |
| 177 |
|
\end{figure} |
| 178 |
|
|
| 179 |
|
The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The |
| 180 |
< |
separation constant $(\alpha)$ plays an important role in balancing |
| 180 |
> |
convergence parameter $(\alpha)$ plays an important role in balancing |
| 181 |
|
the computational cost between the direct and reciprocal-space |
| 182 |
|
portions of the summation. The choice of this value allows one to |
| 183 |
|
select whether the real-space or reciprocal space portion of the |
| 575 |
|
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right), |
| 576 |
|
\end{equation} |
| 577 |
|
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force |
| 578 |
< |
vector computed using method M. |
| 579 |
< |
|
| 580 |
< |
Each of these $\theta$ values was accumulated in a distribution |
| 581 |
< |
function and weighted by the area on the unit sphere. Non-linear |
| 582 |
< |
Gaussian fits were used to measure the width of the resulting |
| 583 |
< |
distributions. |
| 584 |
< |
|
| 585 |
< |
\begin{figure} |
| 586 |
< |
\centering |
| 587 |
< |
\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 588 |
< |
\caption{Sample fit of the angular distribution of the force vectors |
| 589 |
< |
accumulated using all of the studied systems. Gaussian fits were used |
| 590 |
< |
to obtain values for the variance in force and torque vectors.} |
| 591 |
< |
\label{fig:gaussian} |
| 592 |
< |
\end{figure} |
| 593 |
< |
|
| 594 |
< |
Figure \ref{fig:gaussian} shows an example distribution with applied |
| 595 |
< |
non-linear fits. The solid line is a Gaussian profile, while the |
| 596 |
< |
dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 597 |
< |
Lorentzian. Since this distribution is a measure of angular error |
| 598 |
< |
between two different electrostatic summation methods, there is no |
| 599 |
< |
{\it a priori} reason for the profile to adhere to any specific shape. |
| 600 |
< |
Gaussian fits was used to compare all the tested methods. The |
| 601 |
< |
variance ($\sigma^2$) was extracted from each of these fits and was |
| 602 |
< |
used to compare distribution widths. Values of $\sigma^2$ near zero |
| 603 |
< |
indicate vector directions indistinguishable from those calculated |
| 603 |
< |
when using the reference method (SPME). |
| 578 |
> |
vector computed using method M. Each of these $\theta$ values was |
| 579 |
> |
accumulated in a distribution function and weighted by the area on the |
| 580 |
> |
unit sphere. Non-linear Gaussian fits were used to measure the width |
| 581 |
> |
of the resulting distributions. |
| 582 |
> |
% |
| 583 |
> |
%\begin{figure} |
| 584 |
> |
%\centering |
| 585 |
> |
%\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 586 |
> |
%\caption{Sample fit of the angular distribution of the force vectors |
| 587 |
> |
%accumulated using all of the studied systems. Gaussian fits were used |
| 588 |
> |
%to obtain values for the variance in force and torque vectors.} |
| 589 |
> |
%\label{fig:gaussian} |
| 590 |
> |
%\end{figure} |
| 591 |
> |
% |
| 592 |
> |
%Figure \ref{fig:gaussian} shows an example distribution with applied |
| 593 |
> |
%non-linear fits. The solid line is a Gaussian profile, while the |
| 594 |
> |
%dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 595 |
> |
%Lorentzian. |
| 596 |
> |
%Since this distribution is a measure of angular error between two |
| 597 |
> |
%different electrostatic summation methods, there is no {\it a priori} |
| 598 |
> |
%reason for the profile to adhere to any specific shape. |
| 599 |
> |
%Gaussian fits was used to compare all the tested methods. |
| 600 |
> |
The variance ($\sigma^2$) was extracted from each of these fits and |
| 601 |
> |
was used to compare distribution widths. Values of $\sigma^2$ near |
| 602 |
> |
zero indicate vector directions indistinguishable from those |
| 603 |
> |
calculated when using the reference method (SPME). |
| 604 |
|
|
| 605 |
|
\subsection{Short-time Dynamics} |
| 606 |
|
|
| 675 |
|
the resulting configurations were again equilibrated individually. |
| 676 |
|
Finally, for the Argon / Water ``charge void'' systems, the identities |
| 677 |
|
of all the SPC/E waters within 6 \AA\ of the center of the |
| 678 |
< |
equilibrated water configurations were converted to argon |
| 679 |
< |
(Fig. \ref{fig:argonSlice}). |
| 678 |
> |
equilibrated water configurations were converted to argon. |
| 679 |
> |
%(Fig. \ref{fig:argonSlice}). |
| 680 |
|
|
| 681 |
|
These procedures guaranteed us a set of representative configurations |
| 682 |
|
from chemically-relevant systems sampled from an appropriate |
| 683 |
|
ensemble. Force field parameters for the ions and Argon were taken |
| 684 |
|
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 685 |
|
|
| 686 |
< |
\begin{figure} |
| 687 |
< |
\centering |
| 688 |
< |
\includegraphics[width = \linewidth]{./slice.pdf} |
| 689 |
< |
\caption{A slice from the center of a water box used in a charge void |
| 690 |
< |
simulation. The darkened region represents the boundary sphere within |
| 691 |
< |
which the water molecules were converted to argon atoms.} |
| 692 |
< |
\label{fig:argonSlice} |
| 693 |
< |
\end{figure} |
| 686 |
> |
%\begin{figure} |
| 687 |
> |
%\centering |
| 688 |
> |
%\includegraphics[width = \linewidth]{./slice.pdf} |
| 689 |
> |
%\caption{A slice from the center of a water box used in a charge void |
| 690 |
> |
%simulation. The darkened region represents the boundary sphere within |
| 691 |
> |
%which the water molecules were converted to argon atoms.} |
| 692 |
> |
%\label{fig:argonSlice} |
| 693 |
> |
%\end{figure} |
| 694 |
|
|
| 695 |
|
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 696 |
|
We compared the following alternative summation methods with results |
| 714 |
|
manner across all systems and configurations. |
| 715 |
|
|
| 716 |
|
The althernative methods were also evaluated with three different |
| 717 |
< |
cutoff radii (9, 12, and 15 \AA). It should be noted that the damping |
| 718 |
< |
parameter chosen in SPME, or so called ``Ewald Coefficient'', has a |
| 719 |
< |
significant effect on the energies and forces calculated. Typical |
| 720 |
< |
molecular mechanics packages set this to a value dependent on the |
| 721 |
< |
cutoff radius and a tolerance (typically less than $1 \times 10^{-4}$ |
| 722 |
< |
kcal/mol). Smaller tolerances are typically associated with increased |
| 723 |
< |
accuracy at the expense of increased time spent calculating the |
| 724 |
< |
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 725 |
< |
The default TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 726 |
< |
in all SPME calculations, resulting in Ewald Coefficients of 0.4200, |
| 727 |
< |
0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ |
| 728 |
< |
respectively. |
| 717 |
> |
cutoff radii (9, 12, and 15 \AA). As noted previously, the |
| 718 |
> |
convergence parameter ($\alpha$) plays a role in the balance of the |
| 719 |
> |
real-space and reciprocal-space portions of the Ewald calculation. |
| 720 |
> |
Typical molecular mechanics packages set this to a value dependent on |
| 721 |
> |
the cutoff radius and a tolerance (typically less than $1 \times |
| 722 |
> |
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 723 |
> |
increased accuracy at the expense of increased time spent calculating |
| 724 |
> |
the reciprocal-space portion of the |
| 725 |
> |
summation.\cite{Perram88,Essmann95} The default TINKER tolerance of $1 |
| 726 |
> |
\times 10^{-8}$ kcal/mol was used in all SPME calculations, resulting |
| 727 |
> |
in Ewald Coefficients of 0.4200, 0.3119, and 0.2476 \AA$^{-1}$ for |
| 728 |
> |
cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 729 |
|
|
| 730 |
|
\section{Results and Discussion} |
| 731 |
|
|
| 800 |
|
\begin{figure} |
| 801 |
|
\centering |
| 802 |
|
\includegraphics[width=5.5in]{./frcMagplot.pdf} |
| 803 |
< |
\caption{Statistical analysis of the quality of the force vector magnitudes for a given electrostatic method compared with the reference Ewald sum. Results with a value equal to 1 (dashed line) indicate force magnitude values indistinguishable from those obtained using SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 803 |
> |
\caption{Statistical analysis of the quality of the force vector |
| 804 |
> |
magnitudes for a given electrostatic method compared with the |
| 805 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 806 |
> |
indicate force magnitude values indistinguishable from those obtained |
| 807 |
> |
using SPME. Different values of the cutoff radius are indicated with |
| 808 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 809 |
> |
inverted triangles).} |
| 810 |
|
\label{fig:frcMag} |
| 811 |
|
\end{figure} |
| 812 |
|
|
| 832 |
|
\begin{figure} |
| 833 |
|
\centering |
| 834 |
|
\includegraphics[width=5.5in]{./trqMagplot.pdf} |
| 835 |
< |
\caption{Statistical analysis of the quality of the torque vector magnitudes for a given electrostatic method compared with the reference Ewald sum. Results with a value equal to 1 (dashed line) indicate torque magnitude values indistinguishable from those obtained using SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 835 |
> |
\caption{Statistical analysis of the quality of the torque vector |
| 836 |
> |
magnitudes for a given electrostatic method compared with the |
| 837 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 838 |
> |
indicate torque magnitude values indistinguishable from those obtained |
| 839 |
> |
using SPME. Different values of the cutoff radius are indicated with |
| 840 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 841 |
> |
inverted triangles).} |
| 842 |
|
\label{fig:trqMag} |
| 843 |
|
\end{figure} |
| 844 |
|
|
| 868 |
|
\begin{figure} |
| 869 |
|
\centering |
| 870 |
|
\includegraphics[width=5.5in]{./frcTrqAngplot.pdf} |
| 871 |
< |
\caption{Statistical analysis of the quality of the Gaussian fit of the force and torque vector angular distributions for a given electrostatic method compared with the reference Ewald sum. Results with a variance ($\sigma^2$) equal to zero (dashed line) indicate force and torque directions indistinguishable from those obtained using SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 871 |
> |
\caption{Statistical analysis of the quality of the Gaussian fit of |
| 872 |
> |
the force and torque vector angular distributions for a given |
| 873 |
> |
electrostatic method compared with the reference Ewald sum. Results |
| 874 |
> |
with a variance ($\sigma^2$) equal to zero (dashed line) indicate |
| 875 |
> |
force and torque directions indistinguishable from those obtained |
| 876 |
> |
using SPME. Different values of the cutoff radius are indicated with |
| 877 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 878 |
> |
inverted triangles).} |
| 879 |
|
\label{fig:frcTrqAng} |
| 880 |
|
\end{figure} |
| 881 |
|
|
| 899 |
|
|
| 900 |
|
\begin{table}[htbp] |
| 901 |
|
\centering |
| 902 |
< |
\caption{Variance ($\sigma^2$) of the force (top set) and torque (bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods. Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function. The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.} |
| 902 |
> |
\caption{Variance ($\sigma^2$) of the force (top set) and torque |
| 903 |
> |
(bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods. Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function. The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.} |
| 904 |
|
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 905 |
|
\\ |
| 906 |
|
\toprule |
| 988 |
|
\begin{figure} |
| 989 |
|
\centering |
| 990 |
|
\includegraphics[width = \linewidth]{./vCorrPlot.pdf} |
| 991 |
< |
\caption{Velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset is a magnification of the first trough. The times to first collision are nearly identical, but the differences can be seen in the peaks and troughs, where the undamped to weakly damped methods are stiffer than the moderately damped and SPME methods.} |
| 991 |
> |
\caption{Velocity auto-correlation functions of NaCl crystals at |
| 992 |
> |
1000 K while using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and |
| 993 |
> |
{\sc sp} ($\alpha$ = 0.2). The inset is a magnification of the first |
| 994 |
> |
trough. The times to first collision are nearly identical, but the |
| 995 |
> |
differences can be seen in the peaks and troughs, where the undamped |
| 996 |
> |
to weakly damped methods are stiffer than the moderately damped and |
| 997 |
> |
SPME methods.} |
| 998 |
|
\label{fig:vCorrPlot} |
| 999 |
|
\end{figure} |
| 1000 |
|
|
| 1029 |
|
\begin{figure} |
| 1030 |
|
\centering |
| 1031 |
|
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 1032 |
< |
\caption{Power spectra obtained from the velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, {\sc sf} ($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). Apodization of the correlation functions via a cubic switching function between 40 and 50 ps was used to clear up the spectral noise resulting from data truncation, and had no noticeable effect on peak location or magnitude. The inset shows the frequency region below 100 cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 1032 |
> |
\caption{Power spectra obtained from the velocity auto-correlation |
| 1033 |
> |
functions of NaCl crystals at 1000 K while using SPME, {\sc sf} |
| 1034 |
> |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). |
| 1035 |
> |
Apodization of the correlation functions via a cubic switching |
| 1036 |
> |
function between 40 and 50 ps was used to clear up the spectral noise |
| 1037 |
> |
resulting from data truncation, and had no noticeable effect on peak |
| 1038 |
> |
location or magnitude. The inset shows the frequency region below 100 |
| 1039 |
> |
cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 1040 |
|
\label{fig:methodPS} |
| 1041 |
|
\end{figure} |
| 1042 |
|
|
| 1078 |
|
\begin{figure} |
| 1079 |
|
\centering |
| 1080 |
|
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 1081 |
< |
\caption{Regions of spectra showing the low-frequency correlated motions for NaCl crystals at 1000 K using various electrostatic summation methods. The upper plot is a zoomed inset from figure \ref{fig:methodPS}. As the damping value for the {\sc sf} potential increases, the low-frequency peaks red-shift. The lower plot is of spectra when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$. As $\alpha$ increases, the peaks are red-shifted toward and eventually beyond the values given by SPME. The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.} |
| 1081 |
> |
\caption{Regions of spectra showing the low-frequency correlated |
| 1082 |
> |
motions for NaCl crystals at 1000 K using various electrostatic |
| 1083 |
> |
summation methods. The upper plot is a zoomed inset from figure |
| 1084 |
> |
\ref{fig:methodPS}. As the damping value for the {\sc sf} potential |
| 1085 |
> |
increases, the low-frequency peaks red-shift. The lower plot is of |
| 1086 |
> |
spectra when using SPME and a simple damped Coulombic sum with damping |
| 1087 |
> |
coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$. As |
| 1088 |
> |
$\alpha$ increases, the peaks are red-shifted toward and eventually |
| 1089 |
> |
beyond the values given by SPME. The larger $\alpha$ values weaken |
| 1090 |
> |
the real-space electrostatics, explaining this shift towards less |
| 1091 |
> |
strongly correlated motions in the crystal.} |
| 1092 |
|
\label{fig:dampInc} |
| 1093 |
|
\end{figure} |
| 1094 |
|
|
