| 558 |
|
|
| 559 |
|
The effects of the alternative electrostatic summation methods on the |
| 560 |
|
short-time dynamics of charged systems were evaluated by considering a |
| 561 |
< |
NaCl crystal at a temperature of 1000 K. A subset of the best |
| 561 |
> |
NaCl crystal at a temperature of 1000~K. A subset of the best |
| 562 |
|
performing pairwise methods was used in this comparison. The NaCl |
| 563 |
|
crystal was chosen to avoid possible complications from the treatment |
| 564 |
|
of orientational motion in molecular systems. All systems were |
| 571 |
|
\label{eq:vCorr} |
| 572 |
|
\end{equation} |
| 573 |
|
Velocity autocorrelation functions require detailed short time data, |
| 574 |
< |
thus velocity information was saved every 2fs over 10ps |
| 574 |
> |
thus velocity information was saved every 2~fs over 10~ps |
| 575 |
|
trajectories. Because the NaCl crystal is composed of two different |
| 576 |
|
atom types, the average of the two resulting velocity autocorrelation |
| 577 |
|
functions was used for comparisons. |
| 591 |
|
NaCl crystal is composed of two different atom types, the average of |
| 592 |
|
the two resulting power spectra was used for comparisons. Simulations |
| 593 |
|
were performed under the microcanonical ensemble, and velocity |
| 594 |
< |
information was saved every 5fs over 100ps trajectories. |
| 594 |
> |
information was saved every 5~fs over 100~ps trajectories. |
| 595 |
|
|
| 596 |
|
\subsection{Representative Simulations}\label{sec:RepSims} |
| 597 |
|
A variety of representative molecular simulations were analyzed to |
| 608 |
|
\item NaCl melts, |
| 609 |
|
\item a low ionic strength solution of NaCl in water (0.11 M), |
| 610 |
|
\item a high ionic strength solution of NaCl in water (1.1 M), and |
| 611 |
< |
\item a 6\AA\ radius sphere of Argon in water. |
| 611 |
> |
\item a 6~\AA\ radius sphere of Argon in water. |
| 612 |
|
\end{enumerate} |
| 613 |
|
|
| 614 |
|
By utilizing the pairwise techniques (outlined in section |
| 620 |
|
For the solid and liquid water configurations, configurations were |
| 621 |
|
taken at regular intervals from high temperature trajectories of 1000 |
| 622 |
|
SPC/E water molecules. Each configuration was equilibrated |
| 623 |
< |
independently at a lower temperature (300K for the liquid, 200K for |
| 623 |
> |
independently at a lower temperature (300~K for the liquid, 200~K for |
| 624 |
|
the crystal). The solid and liquid NaCl systems consisted of 500 |
| 625 |
|
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
| 626 |
|
these systems were selected and equilibrated in the same manner as the |
| 627 |
|
water systems. In order to introduce measurable fluctuations in the |
| 628 |
|
configuration energy differences, the crystalline simulations were |
| 629 |
< |
equilibrated at 1000K, near the $T_\textrm{m}$ for NaCl. The liquid |
| 629 |
> |
equilibrated at 1000~K, near the $T_\textrm{m}$ for NaCl. The liquid |
| 630 |
|
NaCl configurations needed to represent a fully disordered array of |
| 631 |
< |
point charges, so the high temperature of 7000K was selected for |
| 631 |
> |
point charges, so the high temperature of 7000~K was selected for |
| 632 |
|
equilibration. The ionic solutions were made by solvating 4 (or 40) |
| 633 |
|
ions in a periodic box containing 1000 SPC/E water molecules. Ion and |
| 634 |
|
water positions were then randomly swapped, and the resulting |
| 635 |
|
configurations were again equilibrated individually. Finally, for the |
| 636 |
|
Argon / Water ``charge void'' systems, the identities of all the SPC/E |
| 637 |
< |
waters within 6\AA\ of the center of the equilibrated water |
| 637 |
> |
waters within 6~\AA\ of the center of the equilibrated water |
| 638 |
|
configurations were converted to argon. |
| 639 |
|
|
| 640 |
|
These procedures guaranteed us a set of representative configurations |
| 648 |
|
|
| 649 |
|
\begin{enumerate}[itemsep=0pt] |
| 650 |
|
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 651 |
< |
and 0.3\AA$^{-1}$, |
| 651 |
> |
and 0.3~\AA$^{-1}$, |
| 652 |
|
\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 653 |
< |
and 0.3\AA$^{-1}$, |
| 653 |
> |
and 0.3~\AA$^{-1}$, |
| 654 |
|
\item reaction field with an infinite dielectric constant, and |
| 655 |
|
\item an unmodified cutoff. |
| 656 |
|
\end{enumerate} |
| 666 |
|
manner across all systems and configurations. |
| 667 |
|
|
| 668 |
|
The alternative methods were also evaluated with three different |
| 669 |
< |
cutoff radii (9, 12, and 15\AA). As noted previously, the |
| 669 |
> |
cutoff radii (9, 12, and 15~\AA). As noted previously, the |
| 670 |
|
convergence parameter ($\alpha$) plays a role in the balance of the |
| 671 |
|
real-space and reciprocal-space portions of the Ewald calculation. |
| 672 |
|
Typical molecular mechanics packages set this to a value dependent on |
| 673 |
|
the cutoff radius and a tolerance (typically less than $1 \times |
| 674 |
< |
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 674 |
> |
10^{-4}$~kcal/mol). Smaller tolerances are typically associated with |
| 675 |
|
increasing accuracy at the expense of computational time spent on the |
| 676 |
|
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 677 |
< |
The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 677 |
> |
The default {\sc tinker} tolerance of $1 \times 10^{-8}$~kcal/mol was used |
| 678 |
|
in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200, |
| 679 |
< |
0.3119, and 0.2476\AA$^{-1}$ for cutoff radii of 9, 12, and 15\AA\ |
| 679 |
> |
0.3119, and 0.2476~\AA$^{-1}$ for cutoff radii of 9, 12, and 15~\AA\ |
| 680 |
|
respectively. |
| 681 |
|
|
| 682 |
|
\section{Configuration Energy Difference Results}\label{sec:EnergyResults} |
| 695 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 696 |
|
indicate $\Delta E$ values indistinguishable from those obtained using |
| 697 |
|
{\sc spme}. Different values of the cutoff radius are indicated with |
| 698 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 698 |
> |
different symbols (9~\AA\ = circles, 12~\AA\ = squares, and 15~\AA\ = |
| 699 |
|
inverted triangles).} |
| 700 |
|
\label{fig:delE} |
| 701 |
|
\end{figure} |
| 719 |
|
of neutral groups. |
| 720 |
|
|
| 721 |
|
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 722 |
< |
the agreement with Ewald, and using an $\alpha$ of 0.2\AA $^{-1}$ |
| 722 |
> |
the agreement with Ewald, and using an $\alpha$ of 0.2~\AA $^{-1}$ |
| 723 |
|
shows an excellent correlation and quality of fit with the {\sc spme} |
| 724 |
< |
results, particularly with a cutoff radius greater than 12 |
| 725 |
< |
\AA . Use of a larger damping parameter is more helpful for the |
| 726 |
< |
shortest cutoff shown, but it has a detrimental effect on simulations |
| 727 |
< |
with larger cutoffs. |
| 724 |
> |
results, particularly with a cutoff radius greater than 12~\AA\. Use |
| 725 |
> |
of a larger damping parameter is more helpful for the shortest cutoff |
| 726 |
> |
shown, but it has a detrimental effect on simulations with larger |
| 727 |
> |
cutoffs. |
| 728 |
|
|
| 729 |
|
In the {\sc sf} sets, increasing damping results in progressively {\it |
| 730 |
|
worse} correlation with Ewald. Overall, the undamped case is the best |
| 755 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 756 |
|
indicate force magnitude values indistinguishable from those obtained |
| 757 |
|
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 758 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 758 |
> |
different symbols (9~\AA\ = circles, 12~\AA\ = squares, and 15~\AA\ = |
| 759 |
|
inverted triangles).} |
| 760 |
|
\label{fig:frcMag} |
| 761 |
|
\end{figure} |
| 775 |
|
|
| 776 |
|
With moderate damping and a large enough cutoff radius, the {\sc sp} |
| 777 |
|
method is generating usable forces. Further increases in damping, |
| 778 |
< |
while beneficial for simulations with a cutoff radius of 9\AA\ , is |
| 778 |
> |
while beneficial for simulations with a cutoff radius of 9~\AA\ , is |
| 779 |
|
detrimental to simulations with larger cutoff radii. |
| 780 |
|
|
| 781 |
|
The reaction field results are surprisingly good, considering the poor |
| 794 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 795 |
|
indicate torque magnitude values indistinguishable from those obtained |
| 796 |
|
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 797 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 797 |
> |
different symbols (9~\AA\ = circles, 12~\AA\ = squares, and 15~\AA\ = |
| 798 |
|
inverted triangles).} |
| 799 |
|
\label{fig:trqMag} |
| 800 |
|
\end{figure} |
| 835 |
|
Results with a variance ($\sigma^2$) equal to zero (dashed line) |
| 836 |
|
indicate force and torque directions indistinguishable from those |
| 837 |
|
obtained using {\sc spme}. Different values of the cutoff radius are |
| 838 |
< |
indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, |
| 839 |
< |
and 15\AA\ = inverted triangles).} |
| 838 |
> |
indicated with different symbols (9~\AA\ = circles, 12~\AA\ = squares, |
| 839 |
> |
and 15~\AA\ = inverted triangles).} |
| 840 |
|
\label{fig:frcTrqAng} |
| 841 |
|
\end{figure} |
| 842 |
|
|
| 851 |
|
|
| 852 |
|
All of the sets (aside from the over-damped case) show the improvement |
| 853 |
|
afforded by choosing a larger cutoff radius. Increasing the cutoff |
| 854 |
< |
from 9 to 12\AA\ typically results in a halving of the width of the |
| 855 |
< |
distribution, with a similar improvement when going from 12 to 15 |
| 856 |
< |
\AA . |
| 854 |
> |
from 9 to 12~\AA\ typically results in a halving of the width of the |
| 855 |
> |
distribution, with a similar improvement when going from 12 to |
| 856 |
> |
15~\AA . |
| 857 |
|
|
| 858 |
|
The undamped {\sc sf}, group-based cutoff, and reaction field methods |
| 859 |
|
all do equivalently well at capturing the direction of both the force |
| 877 |
|
obtained both without (N) and with (Y) group based cutoffs and a |
| 878 |
|
switching function. Note that the $\alpha$ values have units of |
| 879 |
|
\AA$^{-1}$ and the variance values have units of degrees$^2$. The |
| 880 |
< |
{\sc sp} (with an $\alpha$ of 0.2\AA$^{-1}$ or smaller) shows much |
| 880 |
> |
{\sc sp} (with an $\alpha$ of 0.2~\AA$^{-1}$ or smaller) shows much |
| 881 |
|
narrower angular distributions when using group-based cutoffs. The |
| 882 |
|
{\sc sf} method likewise shows improvement in the undamped and lightly |
| 883 |
|
damped cases. |
| 901 |
|
|
| 902 |
|
\midrule |
| 903 |
|
|
| 904 |
< |
9\AA & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\ |
| 904 |
> |
9~\AA & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\ |
| 905 |
|
& \textbf{Y} & \textbf{2.486} & \textbf{2.160} & \textbf{0.667} & \textbf{0.608} & \textbf{1.768} & \textbf{1.766} & \textbf{0.676} & \textbf{0.609} \\ |
| 906 |
< |
12\AA & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\ |
| 906 |
> |
12~\AA & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\ |
| 907 |
|
& \textbf{Y} & \textbf{0.515} & \textbf{0.288} & \textbf{0.127} & \textbf{0.586} & \textbf{0.308} & \textbf{0.249} & \textbf{0.127} & \textbf{0.586} \\ |
| 908 |
< |
15\AA & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\ |
| 908 |
> |
15~\AA & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\ |
| 909 |
|
& \textbf{Y} & \textbf{0.228} & \textbf{0.099} & \textbf{0.121} & \textbf{0.598} & \textbf{0.144} & \textbf{0.090} & \textbf{0.121} & \textbf{0.598} \\ |
| 910 |
|
|
| 911 |
|
\midrule |
| 912 |
|
|
| 913 |
< |
9\AA & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\ |
| 913 |
> |
9~\AA & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\ |
| 914 |
|
& \textbf{Y} & \textbf{2.115} & \textbf{1.914} & \textbf{1.878} & \textbf{5.142} & \textbf{2.076} & \textbf{2.039} & \textbf{1.972} & \textbf{5.146} \\ |
| 915 |
< |
12\AA & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\ |
| 915 |
> |
12~\AA & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\ |
| 916 |
|
& \textbf{Y} & \textbf{0.810} & \textbf{0.685} & \textbf{1.352} & \textbf{5.082} & \textbf{0.765} & \textbf{0.714} & \textbf{1.360} & \textbf{5.082} \\ |
| 917 |
< |
15\AA & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\ |
| 917 |
> |
15~\AA & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\ |
| 918 |
|
& \textbf{Y} & \textbf{0.282} & \textbf{0.294} & \textbf{1.272} & \textbf{4.999} & \textbf{0.324} & \textbf{0.318} & \textbf{1.272} & \textbf{4.999} \\ |
| 919 |
|
|
| 920 |
|
\bottomrule |
| 933 |
|
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 934 |
|
values for these types of electrostatic summation methods by fitting |
| 935 |
|
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
| 936 |
< |
0.25, and 0.16\AA$^{-1}$ for cutoff values of 9, 12, and 15\AA\ |
| 936 |
> |
0.25, and 0.16~\AA$^{-1}$ for cutoff values of 9, 12, and 15~\AA\ |
| 937 |
|
respectively.\cite{Kast03} These appear to be reasonable choices to |
| 938 |
|
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
| 939 |
|
these findings, choices this high would introduce error in the |
| 940 |
|
molecular torques, particularly for the shorter cutoffs. Based on our |
| 941 |
< |
observations, empirical damping up to 0.2\AA$^{-1}$ is beneficial, |
| 941 |
> |
observations, empirical damping up to 0.2~\AA$^{-1}$ is beneficial, |
| 942 |
|
but damping may be unnecessary when using the {\sc sf} method. |
| 943 |
|
|
| 944 |
|
\section{Individual System Analysis Results}\label{sec:IndividualResults} |
| 955 |
|
|
| 956 |
|
\subsection{SPC/E Water Results}\label{sec:WaterResults} |
| 957 |
|
|
| 958 |
< |
The first system considered was liquid water at 300K using the SPC/E |
| 958 |
> |
The first system considered was liquid water at 300~K using the SPC/E |
| 959 |
|
model of water.\cite{Berendsen87} The results for the energy gap |
| 960 |
|
comparisons and the force and torque vector magnitude comparisons are |
| 961 |
|
shown in table \ref{tab:spce}. The force and torque vector |
| 985 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 986 |
|
\toprule |
| 987 |
|
\toprule |
| 988 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 988 |
> |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 989 |
|
\cmidrule(lr){3-4} |
| 990 |
|
\cmidrule(lr){5-6} |
| 991 |
|
\cmidrule(l){7-8} |
| 1044 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1045 |
|
\cmidrule(lr){3-5} |
| 1046 |
|
\cmidrule(l){6-8} |
| 1047 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1047 |
> |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1048 |
|
\midrule |
| 1049 |
|
PC & & 783.759 & 481.353 & 332.677 & 248.674 & 144.382 & 98.535 \\ |
| 1050 |
|
SP & 0.0 & 659.440 & 380.699 & 250.002 & 235.151 & 134.661 & 88.135 \\ |
| 1075 |
|
\ref{sec:EnergyResults} through \ref{sec:FTDirResults}. There is good |
| 1076 |
|
agreement with {\sc spme} in both energetic and dynamic behavior when |
| 1077 |
|
using the {\sc sf} method with and without damping. The {\sc sp} |
| 1078 |
< |
method does well with an $\alpha$ around 0.2\AA$^{-1}$, particularly |
| 1079 |
< |
with cutoff radii greater than 12\AA. Over-damping the electrostatics |
| 1078 |
> |
method does well with an $\alpha$ around 0.2~\AA$^{-1}$, particularly |
| 1079 |
> |
with cutoff radii greater than 12~\AA. Over-damping the electrostatics |
| 1080 |
|
reduces the agreement between both these methods and {\sc spme}. |
| 1081 |
|
|
| 1082 |
|
The pure cutoff ({\sc pc}) method performs poorly, again mirroring the |
| 1107 |
|
torque vectors (table \ref{tab:spceAng}). The {\sc sp} method shows |
| 1108 |
|
significant narrowing of the angle distribution when using little to |
| 1109 |
|
no damping and only modest improvement for the recommended conditions |
| 1110 |
< |
($\alpha = 0.2$\AA$^{-1}$ and $R_\textrm{c}~\geqslant~12$\AA). The |
| 1110 |
> |
($\alpha = 0.2$~\AA$^{-1}$ and $R_\textrm{c}~\geqslant~12$~\AA). The |
| 1111 |
|
{\sc sf} method shows modest narrowing across all damping and cutoff |
| 1112 |
|
ranges of interest. When over-damping these methods, group cutoffs and |
| 1113 |
|
the switching function do not improve the force and torque |
| 1137 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1138 |
|
\toprule |
| 1139 |
|
\toprule |
| 1140 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1140 |
> |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1141 |
|
\cmidrule(lr){3-4} |
| 1142 |
|
\cmidrule(lr){5-6} |
| 1143 |
|
\cmidrule(l){7-8} |
| 1196 |
|
$\sigma^2$} \\ |
| 1197 |
|
\cmidrule(lr){3-5} |
| 1198 |
|
\cmidrule(l){6-8} |
| 1199 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1199 |
> |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1200 |
|
\midrule |
| 1201 |
|
PC & & 2128.921 & 603.197 & 715.579 & 329.056 & 221.397 & 81.042 \\ |
| 1202 |
|
SP & 0.0 & 1429.341 & 470.320 & 447.557 & 301.678 & 197.437 & 73.840 \\ |
| 1261 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1262 |
|
\toprule |
| 1263 |
|
\toprule |
| 1264 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1264 |
> |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1265 |
|
\cmidrule(lr){3-4} |
| 1266 |
|
\cmidrule(lr){5-6} |
| 1267 |
|
\cmidrule(l){7-8} |
| 1303 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1304 |
|
\cmidrule(lr){3-5} |
| 1305 |
|
\cmidrule(l){6-8} |
| 1306 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA \\ |
| 1306 |
> |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA \\ |
| 1307 |
|
\midrule |
| 1308 |
|
PC & & 13.294 & 8.035 & 5.366 \\ |
| 1309 |
|
SP & 0.0 & 13.316 & 8.037 & 5.385 \\ |
| 1334 |
|
|
| 1335 |
|
Similar to the use of ice I$_\textrm{c}$ to investigate the role of |
| 1336 |
|
order in molecular systems on the effectiveness of the pairwise |
| 1337 |
< |
methods, the 1000K NaCl crystal system was used to investigate the |
| 1337 |
> |
methods, the 1000~K NaCl crystal system was used to investigate the |
| 1338 |
|
accuracy of the pairwise summation methods in an ordered system of |
| 1339 |
|
charged particles. The results for the energy gap comparisons and the |
| 1340 |
|
force vector magnitude comparisons are shown in table \ref{tab:salt}. |
| 1351 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1352 |
|
\toprule |
| 1353 |
|
\toprule |
| 1354 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1354 |
> |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1355 |
|
\cmidrule(lr){3-4} |
| 1356 |
|
\cmidrule(lr){5-6} |
| 1357 |
|
\cmidrule(l){7-8} |
| 1394 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} \\ |
| 1395 |
|
\cmidrule(lr){3-5} |
| 1396 |
|
\cmidrule(l){6-8} |
| 1397 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA \\ |
| 1397 |
> |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA \\ |
| 1398 |
|
\midrule |
| 1399 |
|
PC & & 111.945 & 111.824 & 111.866 \\ |
| 1400 |
|
SP & 0.0 & 112.414 & 152.215 & 38.087 \\ |
| 1449 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1450 |
|
\toprule |
| 1451 |
|
\toprule |
| 1452 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1452 |
> |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1453 |
|
\cmidrule(lr){3-4} |
| 1454 |
|
\cmidrule(lr){5-6} |
| 1455 |
|
\cmidrule(l){7-8} |
| 1508 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1509 |
|
\cmidrule(lr){3-5} |
| 1510 |
|
\cmidrule(l){6-8} |
| 1511 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1511 |
> |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1512 |
|
\midrule |
| 1513 |
|
PC & & 882.863 & 510.435 & 344.201 & 277.691 & 154.231 & 100.131 \\ |
| 1514 |
|
SP & 0.0 & 732.569 & 405.704 & 257.756 & 261.445 & 142.245 & 91.497 \\ |
| 1568 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1569 |
|
\toprule |
| 1570 |
|
\toprule |
| 1571 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1571 |
> |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1572 |
|
\cmidrule(lr){3-4} |
| 1573 |
|
\cmidrule(lr){5-6} |
| 1574 |
|
\cmidrule(l){7-8} |
| 1627 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1628 |
|
\cmidrule(lr){3-5} |
| 1629 |
|
\cmidrule(l){6-8} |
| 1630 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1630 |
> |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1631 |
|
\midrule |
| 1632 |
|
PC & & 957.784 & 513.373 & 2.260 & 340.043 & 179.443 & 13.079 \\ |
| 1633 |
|
SP & 0.0 & 786.244 & 139.985 & 259.289 & 311.519 & 90.280 & 105.187 \\ |
| 1664 |
|
their agreement with {\sc spme}. With these results, we still |
| 1665 |
|
recommend undamped to moderate damping for the {\sc sf} method and |
| 1666 |
|
moderate damping for the {\sc sp} method, both with cutoffs greater |
| 1667 |
< |
than 12\AA. |
| 1667 |
> |
than 12~\AA. |
| 1668 |
|
|
| 1669 |
< |
\subsection{6\AA\ Argon Sphere in SPC/E Water Results} |
| 1669 |
> |
\subsection{6~\AA\ Argon Sphere in SPC/E Water Results} |
| 1670 |
|
|
| 1671 |
< |
The final model system studied was a 6\AA\ sphere of Argon solvated |
| 1671 |
> |
The final model system studied was a 6~\AA\ sphere of Argon solvated |
| 1672 |
|
by SPC/E water. This serves as a test case of a specifically sized |
| 1673 |
|
electrostatic defect in a disordered molecular system. The results for |
| 1674 |
|
the energy gap comparisons and the force and torque vector magnitude |
| 1680 |
|
|
| 1681 |
|
\begin{table}[htbp] |
| 1682 |
|
\centering |
| 1683 |
< |
\caption{REGRESSION RESULTS OF THE 6\AA\ ARGON SPHERE IN LIQUID |
| 1683 |
> |
\caption{REGRESSION RESULTS OF THE 6~\AA\ ARGON SPHERE IN LIQUID |
| 1684 |
|
WATER SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR |
| 1685 |
|
MAGNITUDES ({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} |
| 1686 |
|
|
| 1688 |
|
\begin{tabular}{@{} ccrrrrrr @{}} |
| 1689 |
|
\toprule |
| 1690 |
|
\toprule |
| 1691 |
< |
& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ |
| 1691 |
> |
& & \multicolumn{2}{c}{9~\AA} & \multicolumn{2}{c}{12~\AA} & \multicolumn{2}{c}{15~\AA}\\ |
| 1692 |
|
\cmidrule(lr){3-4} |
| 1693 |
|
\cmidrule(lr){5-6} |
| 1694 |
|
\cmidrule(l){7-8} |
| 1737 |
|
\begin{table}[htbp] |
| 1738 |
|
\centering |
| 1739 |
|
\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR |
| 1740 |
< |
DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE 6\AA\ SPHERE OF |
| 1740 |
> |
DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE 6~\AA\ SPHERE OF |
| 1741 |
|
ARGON IN LIQUID WATER SYSTEM} |
| 1742 |
|
|
| 1743 |
|
\footnotesize |
| 1747 |
|
& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ |
| 1748 |
|
\cmidrule(lr){3-5} |
| 1749 |
|
\cmidrule(l){6-8} |
| 1750 |
< |
Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ |
| 1750 |
> |
Method & $\alpha$ & 9~\AA & 12~\AA & 15~\AA & 9~\AA & 12~\AA & 15~\AA \\ |
| 1751 |
|
\midrule |
| 1752 |
|
PC & & 568.025 & 265.993 & 195.099 & 246.626 & 138.600 & 91.654 \\ |
| 1753 |
|
SP & 0.0 & 504.578 & 251.694 & 179.932 & 231.568 & 131.444 & 85.119 \\ |
| 1814 |
|
\centering |
| 1815 |
|
\includegraphics[width = \linewidth]{./figures/vCorrPlot.pdf} |
| 1816 |
|
\caption{Velocity autocorrelation functions of NaCl crystals at |
| 1817 |
< |
1000K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& |
| 1818 |
< |
0.2\AA$^{-1}$), and {\sc sp} ($\alpha$ = 0.2\AA$^{-1}$). The inset is |
| 1817 |
> |
1000~K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& |
| 1818 |
> |
0.2~\AA$^{-1}$), and {\sc sp} ($\alpha$ = 0.2~\AA$^{-1}$). The inset is |
| 1819 |
|
a magnification of the area around the first minimum. The times to |
| 1820 |
|
first collision are nearly identical, but differences can be seen in |
| 1821 |
|
the peaks and troughs, where the undamped and weakly damped methods |
| 1830 |
|
troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than |
| 1831 |
|
any of the other methods. As the damping parameter ($\alpha$) is |
| 1832 |
|
increased, these peaks are smoothed out, and the {\sc sf} method |
| 1833 |
< |
approaches the {\sc spme} results. With $\alpha$ values of 0.2\AA$^{-1}$, |
| 1833 |
> |
approaches the {\sc spme} results. With $\alpha$ values of 0.2~\AA$^{-1}$, |
| 1834 |
|
the {\sc sf} and {\sc sp} functions are nearly identical and track the |
| 1835 |
|
{\sc spme} features quite well. This is not surprising because the {\sc sf} |
| 1836 |
|
and {\sc sp} potentials become nearly identical with increased |
| 1846 |
|
traces of the velocity autocorrelation function. The power spectra for |
| 1847 |
|
the best-performing alternative methods are shown in |
| 1848 |
|
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
| 1849 |
< |
a cubic switching function between 40 and 50ps was used to reduce the |
| 1849 |
> |
a cubic switching function between 40 and 50~ps was used to reduce the |
| 1850 |
|
ringing resulting from data truncation. This procedure had no |
| 1851 |
|
noticeable effect on peak location or magnitude. |
| 1852 |
|
|
| 1854 |
|
\centering |
| 1855 |
|
\includegraphics[width = \linewidth]{./figures/spectraSquare.pdf} |
| 1856 |
|
\caption{Power spectra obtained from the velocity auto-correlation |
| 1857 |
< |
functions of NaCl crystals at 1000K while using {\sc spme}, {\sc sf} |
| 1858 |
< |
($\alpha$ = 0, 0.1, \& 0.2\AA$^{-1}$), and {\sc sp} ($\alpha$ = |
| 1859 |
< |
0.2\AA$^{-1}$). The inset shows the frequency region below 100 |
| 1860 |
< |
cm$^{-1}$ to highlight where the spectra differ.} |
| 1857 |
> |
functions of NaCl crystals at 1000~K while using {\sc spme}, {\sc sf} |
| 1858 |
> |
($\alpha$ = 0, 0.1, \& 0.2~\AA$^{-1}$), and {\sc sp} ($\alpha$ = |
| 1859 |
> |
0.2~\AA$^{-1}$). The inset shows the frequency region below |
| 1860 |
> |
100~cm$^{-1}$ to highlight where the spectra differ.} |
| 1861 |
|
\label{fig:methodPS} |
| 1862 |
|
\end{figure} |
| 1863 |
|
|
| 1865 |
|
alternative methods are quantitatively identical with Ewald spectrum, |
| 1866 |
|
the low frequency region shows how the summation methods differ. |
| 1867 |
|
Considering the low-frequency inset (expanded in the upper frame of |
| 1868 |
< |
figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the |
| 1868 |
> |
figure \ref{fig:dampInc}), at frequencies below 100~cm$^{-1}$, the |
| 1869 |
|
correlated motions are blue-shifted when using undamped or weakly |
| 1870 |
< |
damped {\sc sf}. When using moderate damping ($\alpha = 0.2$ |
| 1871 |
< |
\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical |
| 1872 |
< |
correlated motion to the Ewald method (which has a convergence |
| 1873 |
< |
parameter of 0.3119\AA$^{-1}$). This weakening of the electrostatic |
| 1874 |
< |
interaction with increased damping explains why the long-ranged |
| 1875 |
< |
correlated motions are at lower frequencies for the moderately damped |
| 1876 |
< |
methods than for undamped or weakly damped methods. |
| 1870 |
> |
damped {\sc sf}. When using moderate damping ($\alpha = |
| 1871 |
> |
0.2$~\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly |
| 1872 |
> |
identical correlated motion to the Ewald method (which has a |
| 1873 |
> |
convergence parameter of 0.3119~\AA$^{-1}$). This weakening of the |
| 1874 |
> |
electrostatic interaction with increased damping explains why the |
| 1875 |
> |
long-ranged correlated motions are at lower frequencies for the |
| 1876 |
> |
moderately damped methods than for undamped or weakly damped methods. |
| 1877 |
|
|
| 1878 |
|
To isolate the role of the damping constant, we have computed the |
| 1879 |
|
spectra for a single method ({\sc sf}) with a range of damping |
| 1883 |
|
However, even without any electrostatic damping, the {\sc sf} method |
| 1884 |
|
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
| 1885 |
|
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
| 1886 |
< |
would predict the lowest frequency peak near 325 cm$^{-1}$. {\it |
| 1886 |
> |
would predict the lowest frequency peak near 325~cm$^{-1}$. {\it |
| 1887 |
|
Most} of the collective behavior in the crystal is accurately captured |
| 1888 |
|
using the {\sc sf} method. Quantitative agreement with Ewald can be |
| 1889 |
|
obtained using moderate damping in addition to the shifting at the |
| 1893 |
|
\centering |
| 1894 |
|
\includegraphics[width = \linewidth]{./figures/increasedDamping.pdf} |
| 1895 |
|
\caption{Effect of damping on the two lowest-frequency phonon modes in |
| 1896 |
< |
the NaCl crystal at 1000K. The undamped shifted force ({\sc sf}) |
| 1897 |
< |
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1898 |
< |
electrostatic damping to 0.25\AA$^{-1}$ gives quantitative agreement |
| 1896 |
> |
the NaCl crystal at 1000~K. The undamped shifted force ({\sc sf}) |
| 1897 |
> |
method is off by less than 10~cm$^{-1}$, and increasing the |
| 1898 |
> |
electrostatic damping to 0.25~\AA$^{-1}$ gives quantitative agreement |
| 1899 |
|
with the power spectrum obtained using the Ewald sum. Over-damping can |
| 1900 |
|
result in underestimates of frequencies of the long-wavelength |
| 1901 |
|
motions.} |
| 1934 |
|
the Lennard-Jones parameters to correct the density at |
| 1935 |
|
298K.\cite{Rick04} With the density corrected, he compared common |
| 1936 |
|
water properties for TIP5P-E using the Ewald sum with TIP5P using a |
| 1937 |
< |
9\AA\ cutoff. |
| 1937 |
> |
9~\AA\ cutoff. |
| 1938 |
|
|
| 1939 |
|
In the following sections, we compared these same water properties |
| 1940 |
|
calculated from TIP5P-E using the Ewald sum with TIP5P-E using the |
| 1941 |
|
{\sc sf} technique. In the above evaluation of the pairwise |
| 1942 |
|
techniques, we observed some flexibility in the choice of parameters. |
| 1943 |
|
Because of this, the following comparisons include the {\sc sf} |
| 1944 |
< |
technique with a 12\AA\ cutoff and an $\alpha$ = 0.0, 0.1, and |
| 1945 |
< |
0.2\AA$^{-1}$, as well as a 9\AA\ cutoff with an $\alpha$ = |
| 1946 |
< |
0.2\AA$^{-1}$. |
| 1944 |
> |
technique with a 12~\AA\ cutoff and an $\alpha$ = 0.0, 0.1, and |
| 1945 |
> |
0.2~\AA$^{-1}$, as well as a 9~\AA\ cutoff with an $\alpha$ = |
| 1946 |
> |
0.2~\AA$^{-1}$. |
| 1947 |
|
|
| 1948 |
|
\subsection{Density}\label{sec:t5peDensity} |
| 1949 |
|
|
| 1975 |
|
To compare densities, $NPT$ simulations were performed with the same |
| 1976 |
|
temperatures as those selected by Rick in his Ewald summation |
| 1977 |
|
simulations.\cite{Rick04} In order to improve statistics around the |
| 1978 |
< |
density maximum, 3ns trajectories were accumulated at 0, 12.5, and |
| 1979 |
< |
25$^\circ$C, while 2ns trajectories were obtained at all other |
| 1978 |
> |
density maximum, 3~ns trajectories were accumulated at 0, 12.5, and |
| 1979 |
> |
25$^\circ$C, while 2~ns trajectories were obtained at all other |
| 1980 |
|
temperatures. The average densities were calculated from the later |
| 1981 |
|
three-fourths of each trajectory. Similar to Mahoney and Jorgensen's |
| 1982 |
|
method for accumulating statistics, these sequences were spliced into |
| 1990 |
|
with various parameters. The pressure term from the image-charge shell |
| 1991 |
|
is larger than that provided by the reciprocal-space portion of the |
| 1992 |
|
Ewald summation, leading to slightly lower densities. This effect is |
| 1993 |
< |
more visible with the 9\AA\ cutoff, where the image charges exert a |
| 1993 |
> |
more visible with the 9~\AA\ cutoff, where the image charges exert a |
| 1994 |
|
greater force on the central particle. The error bars for the {\sc sf} |
| 1995 |
|
methods show plus or minus the standard deviation of the density |
| 1996 |
|
measurement at each temperature.} |
| 2005 |
|
indicate that the pressure component from the image charges at |
| 2006 |
|
R$_\textrm{c}$ is larger than that exerted by the reciprocal-space |
| 2007 |
|
portion of the Ewald summation. Bringing the image charges closer to |
| 2008 |
< |
the central particle by choosing a 9\AA\ R$_\textrm{c}$ (rather than |
| 2009 |
< |
the preferred 12\AA\ R$_\textrm{c}$) increases the strength of their |
| 2008 |
> |
the central particle by choosing a 9~\AA\ R$_\textrm{c}$ (rather than |
| 2009 |
> |
the preferred 12~\AA\ R$_\textrm{c}$) increases the strength of their |
| 2010 |
|
interactions, resulting in a further reduction of the densities. |
| 2011 |
|
|
| 2012 |
|
Because the strength of the image charge interactions has a noticeable |
| 2025 |
|
dependence of the computed densities when using the Ewald summation, |
| 2026 |
|
most likely due to his tying of the convergence parameter to the box |
| 2027 |
|
dimensions.\cite{Rick04} For systems of 256 water molecules, the |
| 2028 |
< |
calculated densities were on average 0.002 to 0.003 g/cm$^3$ lower. A |
| 2028 |
> |
calculated densities were on average 0.002 to 0.003~g/cm$^3$ lower. A |
| 2029 |
|
system size of 256 molecules would force the use of a shorter |
| 2030 |
|
R$_\textrm{c}$ when using the {\sc sf} technique, and this would also |
| 2031 |
|
lower the densities. Moving to larger systems, as long as the |
| 2069 |
|
|
| 2070 |
|
\begin{figure} |
| 2071 |
|
\includegraphics[width=\linewidth]{./figures/tip5peGofR.pdf} |
| 2072 |
< |
\caption{The $g_\textrm{OO}(r)$s calculated for TIP5P-E at 298K and |
| 2072 |
> |
\caption{The $g_\textrm{OO}(r)$s calculated for TIP5P-E at 298~K and |
| 2073 |
|
1atm while using the Ewald summation (Ref. \cite{Rick04}) and the {\sc |
| 2074 |
|
sf} technique with varying parameters. Even with the reduced densities |
| 2075 |
|
using the {\sc sf} technique, the $g_\textrm{OO}(r)$s are essentially |
| 2191 |
|
All of the above properties were calculated from the same trajectories |
| 2192 |
|
used to determine the densities in section \ref{sec:t5peDensity} |
| 2193 |
|
except for the static dielectric constants. The $\epsilon$ values were |
| 2194 |
< |
accumulated from 2ns $NVE$ ensemble trajectories with system densities |
| 2194 |
> |
accumulated from 2~ns $NVE$ ensemble trajectories with system densities |
| 2195 |
|
fixed at the average values from the $NPT$ simulations at each of the |
| 2196 |
|
temperatures. The resulting values are displayed in figure |
| 2197 |
|
\ref{fig:t5peThermo}. |
| 2243 |
|
indicate a more pronounced transition in the supercooled regime, |
| 2244 |
|
particularly in the case of {\sc sf} without damping. This points to |
| 2245 |
|
the onset of a more frustrated or glassy behavior for TIP5P-E at |
| 2246 |
< |
temperatures below 250K in these simulations. Because the systems are |
| 2246 |
> |
temperatures below 250~K in these simulations. Because the systems are |
| 2247 |
|
locked in different regions of phase-space, comparisons between |
| 2248 |
|
properties at these temperatures are not exactly fair. This |
| 2249 |
|
observation is explored in more detail in section |
| 2258 |
|
converged $\epsilon$ values accumulated for the {\sc sf} |
| 2259 |
|
simulations. Lack of a damping function results in dielectric |
| 2260 |
|
constants significantly smaller than that obtained using the Ewald |
| 2261 |
< |
sum. Increasing the damping coefficient to 0.2\AA$^{-1}$ improves the |
| 2261 |
> |
sum. Increasing the damping coefficient to 0.2~\AA$^{-1}$ improves the |
| 2262 |
|
agreement considerably. It should be noted that the choice of the |
| 2263 |
|
``Ewald coefficient'' value also has a significant effect on the |
| 2264 |
|
calculated value when using the Ewald summation. In the simulations of |
| 2275 |
|
\subsection{Dynamic Properties}\label{sec:t5peDynamics} |
| 2276 |
|
|
| 2277 |
|
To look at the dynamic properties of TIP5P-E when using the {\sc sf} |
| 2278 |
< |
method, 200ps $NVE$ simulations were performed for each temperature at |
| 2278 |
> |
method, 200~ps $NVE$ simulations were performed for each temperature at |
| 2279 |
|
the average density reported by the $NPT$ simulations. The |
| 2280 |
|
self-diffusion constants ($D$) were calculated with the Einstein |
| 2281 |
|
relation using the mean square displacement (MSD), |
| 2314 |
|
\end{figure} |
| 2315 |
|
In addition to translational diffusion, reorientational time constants |
| 2316 |
|
were calculated for comparisons with the Ewald simulations and with |
| 2317 |
< |
experiments. These values were determined from 25ps $NVE$ trajectories |
| 2317 |
> |
experiments. These values were determined from 25~ps $NVE$ trajectories |
| 2318 |
|
through calculation of the orientational time correlation function, |
| 2319 |
|
\begin{equation} |
| 2320 |
|
C_l^\alpha(t) = \left\langle P_l\left[\hat{\mathbf{u}}_i^\alpha(t) |
| 2340 |
|
constants for rotational relaxation. Figure \ref{fig:OrientCorr} shows |
| 2341 |
|
some example plots of orientational autocorrelation functions for the |
| 2342 |
|
first and second Legendre polynomials. The relatively short time |
| 2343 |
< |
portions (between 1 and 3ps for water) of these curves can be fit to |
| 2343 |
> |
portions (between 1 and 3~ps for water) of these curves can be fit to |
| 2344 |
|
an exponential decay to obtain these constants, and they are directly |
| 2345 |
|
comparable to water orientational relaxation times from nuclear |
| 2346 |
|
magnetic resonance (NMR). The relaxation constant obtained from |
| 2384 |
|
of the choice of damping constant. Their are several possible reasons |
| 2385 |
|
for this deviation between techniques. The Ewald results were taken |
| 2386 |
|
shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick |
| 2387 |
< |
calculation from a 10ps trajectory with {\sc sf} with an $\alpha$ of |
| 2388 |
< |
0.2\AA$^-1$ at 25$^\circ$C showed a 0.4ps drop in $\tau_2^y$, placing |
| 2387 |
> |
calculation from a 10~ps trajectory with {\sc sf} with an $\alpha$ of |
| 2388 |
> |
0.2~\AA$^-1$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$, placing |
| 2389 |
|
the result more in line with that obtained using the Ewald sum. These |
| 2390 |
|
results support this explanation; however, recomputing the results to |
| 2391 |
|
meet a poorer statistical standard is counter-productive. Assuming the |
| 2586 |
|
In order to find these optimal values, we mapped out the static |
| 2587 |
|
dielectric constant as a function of both the damping parameter and |
| 2588 |
|
cutoff radius for several different water models. To calculate the |
| 2589 |
< |
static dielectric constant, we performed 5ns $NPT$ calculations on |
| 2589 |
> |
static dielectric constant, we performed 5~ns $NPT$ calculations on |
| 2590 |
|
systems of 512 water molecules, using the TIP5P-E, TIP4P-Ew, SPC/E, |
| 2591 |
|
and SSD/RF water models. TIP4P-Ew is a reparametrized version of the |
| 2592 |
|
four-point transferable intermolecular potential (TIP4P) for water |
| 2598 |
|
charge-charge interactions like the other three models. Damping of the |
| 2599 |
|
dipole-dipole interaction was handled as described in section |
| 2600 |
|
\ref{sec:dampingMultipoles}. Each of these systems were studied with |
| 2601 |
< |
cutoff radii of 9, 10, 11, and 12\AA\ and with damping parameter values |
| 2602 |
< |
ranging from 0 to 0.35\AA$^{-1}$. |
| 2601 |
> |
cutoff radii of 9, 10, 11, and 12~\AA\ and with damping parameter values |
| 2602 |
> |
ranging from 0 to 0.35~\AA$^{-1}$. |
| 2603 |
|
\begin{figure} |
| 2604 |
|
\centering |
| 2605 |
|
\includegraphics[width=\linewidth]{./figures/dielectricMap.pdf} |
| 2614 |
|
interesting aspect of all four contour plots is that the dielectric |
| 2615 |
|
constant is effectively linear with respect to $\alpha$ and |
| 2616 |
|
$R_\textrm{c}$ in the low to moderate damping regions, and the slope |
| 2617 |
< |
is 0.025\AA$^{-1}\cdot R_\textrm{c}^{-1}$. Another point to note is |
| 2617 |
> |
is 0.025~\AA$^{-1}\cdot R_\textrm{c}^{-1}$. Another point to note is |
| 2618 |
|
that choosing $\alpha$ and $R_\textrm{c}$ identical to those used in |
| 2619 |
|
studies with the Ewald summation results in the same calculated |
| 2620 |
|
dielectric constant. As an example, in the paper outlining the |
| 2621 |
|
development of TIP5P-E, the real-space cutoff and Ewald coefficient |
| 2622 |
|
were tethered to the system size, and for a 512 molecule system are |
| 2623 |
< |
approximately 12\AA\ and 0.25\AA$^{-1}$ respectively.\cite{Rick04} |
| 2623 |
> |
approximately 12~\AA\ and 0.25~\AA$^{-1}$ respectively.\cite{Rick04} |
| 2624 |
|
These parameters resulted in a dielectric constant of 92$\pm$14, while |
| 2625 |
|
with {\sc sf} these parameters give a dielectric constant of |
| 2626 |
|
90.8$\pm$0.9. Another example comes from the TIP4P-Ew paper where |
| 2627 |
< |
$\alpha$ and $R_\textrm{c}$ were chosen to be 9.5\AA\ and |
| 2628 |
< |
0.35\AA$^{-1}$, and these parameters resulted in a $\epsilon_0$ equal |
| 2627 |
> |
$\alpha$ and $R_\textrm{c}$ were chosen to be 9.5~\AA\ and |
| 2628 |
> |
0.35~\AA$^{-1}$, and these parameters resulted in a $\epsilon_0$ equal |
| 2629 |
|
to 63$\pm$1.\cite{Horn04} We did not perform calculations with these |
| 2630 |
|
exact parameters, but interpolating between surrounding values gives a |
| 2631 |
|
$\epsilon_0$ of 61$\pm$1. Seeing a dependence of the dielectric |
| 2637 |
|
these Ewald examples, the results discussed in sections |
| 2638 |
|
\ref{sec:EnergyResults} through \ref{sec:IndividualResults} indicate |
| 2639 |
|
that values this high are destructive to both the energetics and |
| 2640 |
< |
dynamics. Ideally, $\alpha$ should not exceed 0.3\AA$^{-1}$ for any of |
| 2640 |
> |
dynamics. Ideally, $\alpha$ should not exceed 0.3~\AA$^{-1}$ for any of |
| 2641 |
|
the cutoff values in this range. If the optimal damping parameter is |
| 2642 |
< |
chosen to be midway between 0.275 and 0.3\AA$^{-1}$ (0.2875\AA$^{-1}$) |
| 2643 |
< |
at the 9\AA\ cutoff, then the linear trend with $R_\textrm{c}$ will |
| 2644 |
< |
always keep $\alpha$ below 0.3\AA$^{-1}$. This linear progression |
| 2645 |
< |
would give values of 0.2875, 0.2625, 0.2375, and 0.2125\AA$^{-1}$ for |
| 2646 |
< |
cutoff radii of 9, 10, 11, and 12\AA. Setting this to be the default |
| 2642 |
> |
chosen to be midway between 0.275 and 0.3~\AA$^{-1}$ (0.2875~\AA$^{-1}$) |
| 2643 |
> |
at the 9~\AA\ cutoff, then the linear trend with $R_\textrm{c}$ will |
| 2644 |
> |
always keep $\alpha$ below 0.3~\AA$^{-1}$. This linear progression |
| 2645 |
> |
would give values of 0.2875, 0.2625, 0.2375, and 0.2125~\AA$^{-1}$ for |
| 2646 |
> |
cutoff radii of 9, 10, 11, and 12~\AA. Setting this to be the default |
| 2647 |
|
behavior for the damped {\sc sf} technique will result in consistent |
| 2648 |
|
dielectric behavior for these and other condensed molecular systems, |
| 2649 |
|
regardless of the chosen cutoff radius. The static dielectric |
| 2696 |
|
If a researcher is using Monte Carlo simulations of large chemical |
| 2697 |
|
systems containing point charges, most structural features will be |
| 2698 |
|
accurately captured using the undamped {\sc sf} method or the {\sc sp} |
| 2699 |
< |
method with an electrostatic damping of 0.2\AA$^{-1}$. These methods |
| 2699 |
> |
method with an electrostatic damping of 0.2~\AA$^{-1}$. These methods |
| 2700 |
|
would also be appropriate for molecular dynamics simulations where the |
| 2701 |
|
data of interest is either structural or short-time dynamical |
| 2702 |
|
quantities. For long-time dynamics and collective motions, the safest |
| 2703 |
|
pairwise method we have evaluated is the {\sc sf} method with an |
| 2704 |
< |
electrostatic damping between 0.2 and 0.25\AA$^{-1}$. It is also |
| 2704 |
> |
electrostatic damping between 0.2 and 0.25~\AA$^{-1}$. It is also |
| 2705 |
|
important to note that the static dielectric constant in water |
| 2706 |
|
simulations is highly dependent on both $\alpha$ and |
| 2707 |
|
$R_\textrm{c}$. For consistent dielectric behavior, the damped {\sc |
| 2708 |
< |
sf} method should use an $\alpha$ of 0.2175\AA$^{-1}$ for an |
| 2709 |
< |
$R_\textrm{c}$ of 12\AA, and $\alpha$ should decrease by |
| 2710 |
< |
0.025\AA$^{-1}$ for every 1\AA\ increase in cutoff radius. |
| 2708 |
> |
sf} method should use an $\alpha$ of 0.2175~\AA$^{-1}$ for an |
| 2709 |
> |
$R_\textrm{c}$ of 12~\AA, and $\alpha$ should decrease by |
| 2710 |
> |
0.025~\AA$^{-1}$ for every 1~\AA\ increase in cutoff radius. |
| 2711 |
|
|
| 2712 |
|
We are not suggesting that there is any flaw with the Ewald sum; in |
| 2713 |
|
fact, it is the standard by which these simple pairwise sums have been |
