ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/chrisDissertation/Electrostatics.tex
(Generate patch)

Comparing trunk/chrisDissertation/Electrostatics.tex (file contents):
Revision 2987 by chrisfen, Wed Aug 30 22:36:06 2006 UTC vs.
Revision 3023 by chrisfen, Tue Sep 26 03:07:59 2006 UTC

# Line 1 | Line 1
1 < \chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION
2 < TECHNIQUES}
1 > \chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION \\ TECHNIQUES}
2  
3 < In molecular simulations, proper accumulation of the electrostatic
3 > In molecular simulations, proper accumulation of electrostatic
4   interactions is essential and is one of the most
5   computationally-demanding tasks.  The common molecular mechanics force
6   fields represent atomic sites with full or partial charges protected
7 < by repulsive Lennard-Jones interactions.  This means that nearly
8 < every pair interaction involves a calculation of charge-charge forces.
7 > by repulsive Lennard-Jones interactions.  This means that nearly every
8 > pair interaction involves a calculation of charge-charge forces.
9   Coupled with relatively long-ranged $r^{-1}$ decay, the monopole
10   interactions quickly become the most expensive part of molecular
11   simulations.  Historically, the electrostatic pair interaction would
# Line 102 | Line 101 | system is said to be using conducting (or ``tin-foil''
101   dipole moment which is magnified through replication of the periodic
102   images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the
103   system is said to be using conducting (or ``tin-foil'') boundary
104 < conditions, $\epsilon_{\rm S} = \infty$. Figure
106 < \ref{fig:ewaldTime} shows how the Ewald sum has been applied over
107 < time.  Initially, due to the small system sizes that could be
108 < simulated feasibly, the entire simulation box was replicated to
109 < convergence.  In more modern simulations, the systems have grown large
110 < enough that a real-space cutoff could potentially give convergent
111 < behavior.  Indeed, it has been observed that with the choice of a
112 < small $\alpha$, the reciprocal-space portion of the Ewald sum can be
113 < rapidly convergent and small relative to the real-space
114 < portion.\cite{Karasawa89,Kolafa92}
104 > conditions, $\epsilon_{\rm S} = \infty$.
105  
106   \begin{figure}
107   \includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf}
# Line 122 | Line 112 | convergence for the larger systems of charges that are
112   convergence for the larger systems of charges that are common today.}
113   \label{fig:ewaldTime}
114   \end{figure}
115 + Figure \ref{fig:ewaldTime} shows how the Ewald sum has been applied
116 + over time.  Initially, due to the small system sizes that could be
117 + simulated feasibly, the entire simulation box was replicated to
118 + convergence.  In more modern simulations, the systems have grown large
119 + enough that a real-space cutoff could potentially give convergent
120 + behavior.  Indeed, it has been observed that with the choice of a
121 + small $\alpha$, the reciprocal-space portion of the Ewald sum can be
122 + rapidly convergent and small relative to the real-space
123 + portion.\cite{Karasawa89,Kolafa92}
124  
125 < The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm.  The
125 > The original Ewald summation is an $\mathcal{O}(N^2)$ algorithm.  The
126   convergence parameter $(\alpha)$ plays an important role in balancing
127   the computational cost between the direct and reciprocal-space
128   portions of the summation.  The choice of this value allows one to
129   select whether the real-space or reciprocal space portion of the
130 < summation is an $\mathscr{O}(N^2)$ calculation (with the other being
131 < $\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of
130 > summation is an $\mathcal{O}(N^2)$ calculation (with the other being
131 > $\mathcal{O}(N)$).\cite{Sagui99} With the appropriate choice of
132   $\alpha$ and thoughtful algorithm development, this cost can be
133 < reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route
133 > reduced to $\mathcal{O}(N^{3/2})$.\cite{Perram88} The typical route
134   taken to reduce the cost of the Ewald summation even further is to set
135   $\alpha$ such that the real-space interactions decay rapidly, allowing
136   for a short spherical cutoff. Then the reciprocal space summation is
# Line 140 | Line 139 | methods, the cost of the reciprocal-space portion of t
139   particle-particle particle-mesh (P3M) and particle mesh Ewald (PME)
140   methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these
141   methods, the cost of the reciprocal-space portion of the Ewald
142 < summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N
142 > summation is reduced from $\mathcal{O}(N^2)$ down to $\mathcal{O}(N
143   \log N)$.
144  
145   These developments and optimizations have made the use of the Ewald
# Line 164 | Line 163 | considering the use of the Ewald summation where the a
163   artificially stabilized by the periodic replicas introduced by the
164   Ewald summation.\cite{Weber00} Thus, care must be taken when
165   considering the use of the Ewald summation where the assumed
166 < periodicity would introduce spurious effects in the system dynamics.
166 > periodicity would introduce spurious effects.
167  
168  
169   \section{The Wolf and Zahn Methods}
# Line 453 | Line 452 | summation performed by the Ewald sum.
452  
453   \section{Evaluating Pairwise Summation Techniques}
454  
455 < In classical molecular mechanics simulations, there are two primary
456 < techniques utilized to obtain information about the system of
457 < interest: Monte Carlo (MC) and Molecular Dynamics (MD).  Both of these
458 < techniques utilize pairwise summations of interactions between
459 < particle sites, but they use these summations in different ways.
455 > As mentioned in the introduction, there are two primary techniques
456 > utilized to obtain information about the system of interest in
457 > classical molecular mechanics simulations: Monte Carlo (MC) and
458 > Molecular Dynamics (MD).  Both of these techniques utilize pairwise
459 > summations of interactions between particle sites, but they use these
460 > summations in different ways.
461  
462   In MC, the potential energy difference between configurations dictates
463   the progression of MC sampling.  Going back to the origins of this
# Line 481 | Line 481 | vectors will diverge from each other more rapidly.
481  
482   \subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods}
483  
484 + \begin{figure}
485 + \centering
486 + \includegraphics[width = 3.5in]{./figures/dualLinear.pdf}
487 + \caption{Example least squares regressions of the configuration energy
488 + differences for SPC/E water systems. The upper plot shows a data set
489 + with a poor correlation coefficient ($R^2$), while the lower plot
490 + shows a data set with a good correlation coefficient.}
491 + \label{fig:linearFit}
492 + \end{figure}
493   The pairwise summation techniques (outlined in section
494   \ref{sec:ESMethods}) were evaluated for use in MC simulations by
495   studying the energy differences between conformations.  We took the
# Line 496 | Line 505 | Fig. \ref{fig:linearFit}.
505   Sample correlation plots for two alternate methods are shown in
506   Fig. \ref{fig:linearFit}.
507  
499 \begin{figure}
500 \centering
501 \includegraphics[width = 3.5in]{./figures/dualLinear.pdf}
502 \caption{Example least squares regressions of the configuration energy
503 differences for SPC/E water systems. The upper plot shows a data set
504 with a poor correlation coefficient ($R^2$), while the lower plot
505 shows a data set with a good correlation coefficient.}
506 \label{fig:linearFit}
507 \end{figure}
508
508   Each of the seven system types (detailed in section \ref{sec:RepSims})
509   were represented using 500 independent configurations.  Thus, each of
510   the alternative (non-Ewald) electrostatic summation methods was
# Line 513 | Line 512 | Results and discussion for the individual analysis of
512   differences.
513  
514   Results and discussion for the individual analysis of each of the
515 < system types appear in sections \ref{sec:IndividualResults}, while the
515 > system types appear in appendix \ref{app:IndividualResults}, while the
516   cumulative results over all the investigated systems appear below in
517   sections \ref{sec:EnergyResults}.
518  
# Line 699 | Line 698 | inverted triangles).}
698   inverted triangles).}
699   \label{fig:delE}
700   \end{figure}
702
701   The most striking feature of this plot is how well the Shifted Force
702   ({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy
703   differences.  For the undamped {\sc sf} method, and the
# Line 714 | Line 712 | significant improvement using the group-switched cutof
712   some degree by using group based cutoffs with a switching
713   function.\cite{Adams79,Steinbach94,Leach01} However, we do not see
714   significant improvement using the group-switched cutoff because the
715 < salt and salt solution systems contain non-neutral groups.  Section
716 < \ref{sec:IndividualResults} includes results for systems comprised entirely
717 < of neutral groups.
715 > salt and salt solution systems contain non-neutral groups.  Appendix
716 > \ref{app:IndividualResults} includes results for systems comprised
717 > entirely of neutral groups.
718  
719   For the {\sc sp} method, inclusion of electrostatic damping improves
720   the agreement with Ewald, and using an $\alpha$ of 0.2~\AA $^{-1}$
721   shows an excellent correlation and quality of fit with the {\sc spme}
722 < results, particularly with a cutoff radius greater than 12~\AA\.  Use
722 > results, particularly with a cutoff radius greater than 12~\AA . Use
723   of a larger damping parameter is more helpful for the shortest cutoff
724   shown, but it has a detrimental effect on simulations with larger
725   cutoffs.
# Line 759 | Line 757 | inverted triangles).}
757   inverted triangles).}
758   \label{fig:frcMag}
759   \end{figure}
762
760   Again, it is striking how well the Shifted Potential and Shifted Force
761   methods are doing at reproducing the {\sc spme} forces.  The undamped and
762   weakly-damped {\sc sf} method gives the best agreement with Ewald.
# Line 798 | Line 795 | inverted triangles).}
795   inverted triangles).}
796   \label{fig:trqMag}
797   \end{figure}
801
798   Molecular torques were only available from the systems which contained
799   rigid molecules (i.e. the systems containing water).  The data in
800   fig. \ref{fig:trqMag} is taken from this smaller sampling pool.
# Line 807 | Line 803 | cutoff radius.  Again, the weakly damped and undamped
803   distance.   The striking feature in comparing the new electrostatic
804   methods with {\sc spme} is how much the agreement improves with increasing
805   cutoff radius.  Again, the weakly damped and undamped {\sc sf} method
806 < appears to be reproducing the {\sc spme} torques most accurately.  
806 > appears to reproduce the {\sc spme} torques most accurately.  
807  
808   Water molecules are dipolar, and the reaction field method reproduces
809   the effect of the surrounding polarized medium on each of the
# Line 816 | Line 812 | performs best of all of the methods on molecular torqu
812  
813   \section{Directionality of the Force and Torque Vector Results}\label{sec:FTDirResults}
814  
815 < It is clearly important that a new electrostatic method can reproduce
816 < the magnitudes of the force and torque vectors obtained via the Ewald
817 < sum. However, the {\it directionality} of these vectors will also be
818 < vital in calculating dynamical quantities accurately.  Force and
819 < torque directionalities were investigated by measuring the angles
820 < formed between these vectors and the same vectors calculated using
821 < {\sc spme}.  The results (Fig. \ref{fig:frcTrqAng}) are compared through the
822 < variance ($\sigma^2$) of the Gaussian fits of the angle error
823 < distributions of the combined set over all system types.
815 > It is clearly important that a new electrostatic method should be able
816 > to reproduce the magnitudes of the force and torque vectors obtained
817 > via the Ewald sum. However, the {\it directionality} of these vectors
818 > will also be vital in calculating dynamical quantities accurately.
819 > Force and torque directionalities were investigated by measuring the
820 > angles formed between these vectors and the same vectors calculated
821 > using {\sc spme}.  The results (Fig. \ref{fig:frcTrqAng}) are compared
822 > through the variance ($\sigma^2$) of the Gaussian fits of the angle
823 > error distributions of the combined set over all system types.
824  
825   \begin{figure}
826   \centering
# Line 839 | Line 835 | and 15~\AA\ = inverted triangles).}
835   and 15~\AA\ = inverted triangles).}
836   \label{fig:frcTrqAng}
837   \end{figure}
842
838   Both the force and torque $\sigma^2$ results from the analysis of the
839   total accumulated system data are tabulated in figure
840   \ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc
# Line 864 | Line 859 | charged bodies, and this observation is investigated f
859   particles in all seven systems, while torque vectors are only
860   available for neutral molecular groups.  Damping is more beneficial to
861   charged bodies, and this observation is investigated further in
862 < section \ref{sec:IndividualResults}.
862 > appendix \ref{app:IndividualResults}.
863  
864   Although not discussed previously, group based cutoffs can be applied
865   to both the {\sc sp} and {\sc sf} methods.  The group-based cutoffs
# Line 940 | Line 935 | but damping may be unnecessary when using the {\sc sf}
935   molecular torques, particularly for the shorter cutoffs.  Based on our
936   observations, empirical damping up to 0.2~\AA$^{-1}$ is beneficial,
937   but damping may be unnecessary when using the {\sc sf} method.
943
944 \section{Individual System Analysis Results}\label{sec:IndividualResults}
945
946 The combined results of the previous sections show how the pairwise
947 methods compare to the Ewald summation in the general sense over all
948 of the system types.  It is also useful to consider each of the
949 studied systems in an individual fashion, so that we can identify
950 conditions that are particularly difficult for a selected pairwise
951 method to address. This allows us to further establish the limitations
952 of these pairwise techniques.  Below, the energy difference, force
953 vector, and torque vector analyses are presented on an individual
954 system basis.
955
956 \subsection{SPC/E Water Results}\label{sec:WaterResults}
957
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
962 directionality results are displayed separately in table
963 \ref{tab:spceAng}, where the effect of group-based cutoffs and
964 switching functions on the {\sc sp} and {\sc sf} potentials are also
965 investigated.  In all of the individual results table, the method
966 abbreviations are as follows:
967
968 \begin{itemize}[itemsep=0pt]
969 \item PC = Pure Cutoff,
970 \item SP = Shifted Potential,
971 \item SF = Shifted Force,
972 \item GSC = Group Switched Cutoff,
973 \item RF = Reaction Field (where $\varepsilon \approx\infty$),
974 \item GSSP = Group Switched Shifted Potential, and
975 \item GSSF = Group Switched Shifted Force.
976 \end{itemize}
977
978 \begin{table}[htbp]
979 \centering
980 \caption{REGRESSION RESULTS OF THE LIQUID WATER SYSTEM FOR THE
981 $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it middle})
982 AND TORQUE VECTOR MAGNITUDES ({\it lower})}
983
984 \footnotesize
985 \begin{tabular}{@{} ccrrrrrr @{}}
986 \toprule
987 \toprule
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}
992 Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
993 \midrule
994 PC  &     & 3.046 & 0.002 & -3.018 & 0.002 & 4.719 & 0.005 \\
995 SP  & 0.0 & 1.035 & 0.218 & 0.908 & 0.313 & 1.037 & 0.470 \\
996    & 0.1 & 1.021 & 0.387 & 0.965 & 0.752 & 1.006 & 0.947 \\
997    & 0.2 & 0.997 & 0.962 & 1.001 & 0.994 & 0.994 & 0.996 \\
998    & 0.3 & 0.984 & 0.980 & 0.997 & 0.985 & 0.982 & 0.987 \\
999 SF  & 0.0 & 0.977 & 0.974 & 0.996 & 0.992 & 0.991 & 0.997 \\
1000    & 0.1 & 0.983 & 0.974 & 1.001 & 0.994 & 0.996 & 0.998 \\
1001    & 0.2 & 0.992 & 0.989 & 1.001 & 0.995 & 0.994 & 0.996 \\
1002    & 0.3 & 0.984 & 0.980 & 0.996 & 0.985 & 0.982 & 0.987 \\
1003 GSC &     & 0.918 & 0.862 & 0.852 & 0.756 & 0.801 & 0.700 \\
1004 RF  &     & 0.971 & 0.958 & 0.975 & 0.987 & 0.959 & 0.983 \\                
1005 \midrule
1006 PC  &     & -1.647 & 0.000 & -0.127 & 0.000 & -0.979 & 0.000 \\
1007 SP  & 0.0 & 0.735 & 0.368 & 0.813 & 0.537 & 0.865 & 0.659 \\
1008    & 0.1 & 0.850 & 0.612 & 0.956 & 0.887 & 0.992 & 0.979 \\
1009    & 0.2 & 0.996 & 0.989 & 1.000 & 1.000 & 1.000 & 1.000 \\
1010    & 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\
1011 SF  & 0.0 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 0.999 \\
1012    & 0.1 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\
1013    & 0.2 & 0.999 & 0.998 & 1.000 & 1.000 & 1.000 & 1.000 \\
1014    & 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\
1015 GSC &     & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\
1016 RF  &     & 0.999 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\          
1017 \midrule
1018 PC  &     & 2.387 & 0.000 & 0.183 & 0.000 & 1.282 & 0.000 \\
1019 SP  & 0.0 & 0.847 & 0.543 & 0.904 & 0.694 & 0.935 & 0.786 \\
1020    & 0.1 & 0.922 & 0.749 & 0.980 & 0.934 & 0.996 & 0.988 \\
1021    & 0.2 & 0.987 & 0.985 & 0.989 & 0.992 & 0.990 & 0.993 \\
1022    & 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\
1023 SF  & 0.0 & 0.978 & 0.990 & 0.988 & 0.997 & 0.993 & 0.999 \\
1024    & 0.1 & 0.983 & 0.991 & 0.993 & 0.997 & 0.997 & 0.999 \\
1025    & 0.2 & 0.986 & 0.989 & 0.989 & 0.992 & 0.990 & 0.993 \\
1026    & 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\
1027 GSC &     & 0.995 & 0.981 & 0.999 & 0.991 & 1.001 & 0.994 \\
1028 RF  &     & 0.993 & 0.989 & 0.998 & 0.996 & 1.000 & 0.999 \\
1029 \bottomrule
1030 \end{tabular}
1031 \label{tab:spce}
1032 \end{table}
1033
1034 \begin{table}[htbp]
1035 \centering
1036 \caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR
1037 DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE LIQUID WATER
1038 SYSTEM}
1039
1040 \footnotesize
1041 \begin{tabular}{@{} ccrrrrrr @{}}
1042 \toprule
1043 \toprule
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 \\
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 \\
1051    & 0.1 & 293.849 & 67.772 & 11.609 & 105.090 & 23.813 & 4.369 \\
1052    & 0.2 & 5.975 & 0.136 & 0.094 & 5.553 & 1.784 & 1.536 \\
1053    & 0.3 & 0.725 & 0.707 & 0.693 & 7.293 & 6.933 & 6.748 \\
1054 SF  & 0.0 & 2.238 & 0.713 & 0.292 & 3.290 & 1.090 & 0.416 \\
1055    & 0.1 & 2.238 & 0.524 & 0.115 & 3.184 & 0.945 & 0.326 \\
1056    & 0.2 & 0.374 & 0.102 & 0.094 & 2.598 & 1.755 & 1.537 \\
1057    & 0.3 & 0.721 & 0.707 & 0.693 & 7.322 & 6.933 & 6.748 \\
1058 GSC &     & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\
1059 RF  &     & 2.091 & 0.403 & 0.113 & 3.583 & 1.071 & 0.399 \\      
1060 \midrule
1061 GSSP  & 0.0 & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\
1062      & 0.1 & 1.879 & 0.291 & 0.057 & 3.983 & 1.117 & 0.370 \\
1063      & 0.2 & 0.443 & 0.103 & 0.093 & 2.821 & 1.794 & 1.532 \\
1064      & 0.3 & 0.728 & 0.694 & 0.692 & 7.387 & 6.942 & 6.748 \\
1065 GSSF  & 0.0 & 1.298 & 0.270 & 0.083 & 3.098 & 0.992 & 0.375 \\
1066      & 0.1 & 1.296 & 0.210 & 0.044 & 3.055 & 0.922 & 0.330 \\
1067      & 0.2 & 0.433 & 0.104 & 0.093 & 2.895 & 1.797 & 1.532 \\
1068      & 0.3 & 0.728 & 0.694 & 0.692 & 7.410 & 6.942 & 6.748 \\
1069 \bottomrule
1070 \end{tabular}
1071 \label{tab:spceAng}
1072 \end{table}
1073
1074 The water results parallel the combined results seen in sections
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
1080 reduces the agreement between both these methods and {\sc spme}.
1081
1082 The pure cutoff ({\sc pc}) method performs poorly, again mirroring the
1083 observations from the combined results.  In contrast to these results, however, the use of a switching function and group
1084 based cutoffs greatly improves the results for these neutral water
1085 molecules.  The group switched cutoff ({\sc gsc}) does not mimic the
1086 energetics of {\sc spme} as well as the {\sc sp} (with moderate
1087 damping) and {\sc sf} methods, but the dynamics are quite good.  The
1088 switching functions correct discontinuities in the potential and
1089 forces, leading to these improved results.  Such improvements with the
1090 use of a switching function have been recognized in previous
1091 studies,\cite{Andrea83,Steinbach94} and this proves to be a useful
1092 tactic for stably incorporating local area electrostatic effects.
1093
1094 The reaction field ({\sc rf}) method simply extends upon the results
1095 observed in the {\sc gsc} case.  Both methods are similar in form
1096 (i.e. neutral groups, switching function), but {\sc rf} incorporates
1097 an added effect from the external dielectric. This similarity
1098 translates into the same good dynamic results and improved energetic
1099 agreement with {\sc spme}.  Though this agreement is not to the level
1100 of the moderately damped {\sc sp} and {\sc sf} methods, these results
1101 show how incorporating some implicit properties of the surroundings
1102 (i.e. $\epsilon_\textrm{S}$) can improve the solvent depiction.
1103
1104 As a final note for the liquid water system, use of group cutoffs and a
1105 switching function leads to noticeable improvements in the {\sc sp}
1106 and {\sc sf} methods, primarily in directionality of the force and
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
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
1114 directionalities.
1115
1116 \subsection{SPC/E Ice I$_\textrm{c}$ Results}\label{sec:IceResults}
1117
1118 In addition to the disordered molecular system above, the ordered
1119 molecular system of ice I$_\textrm{c}$ was also considered.  Ice
1120 polymorph could have been used to fit this role; however, ice
1121 I$_\textrm{c}$ was chosen because it can form an ideal periodic
1122 lattice with the same number of water molecules used in the disordered
1123 liquid state case.  The results for the energy gap comparisons and the
1124 force and torque vector magnitude comparisons are shown in table
1125 \ref{tab:ice}.  The force and torque vector directionality results are
1126 displayed separately in table \ref{tab:iceAng}, where the effect of
1127 group-based cutoffs and switching functions on the {\sc sp} and {\sc
1128 sf} potentials are also displayed.
938  
1130 \begin{table}[htbp]
1131 \centering
1132 \caption{REGRESSION RESULTS OF THE ICE I$_\textrm{c}$ SYSTEM FOR
1133 $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it
1134 middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})}
939  
1136 \footnotesize
1137 \begin{tabular}{@{} ccrrrrrr @{}}
1138 \toprule
1139 \toprule
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}
1144 Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1145 \midrule
1146 PC  &     & 19.897 & 0.047 & -29.214 & 0.048 & -3.771 & 0.001 \\
1147 SP  & 0.0 & -0.014 & 0.000 & 2.135 & 0.347 & 0.457 & 0.045 \\
1148    & 0.1 & 0.321 & 0.017 & 1.490 & 0.584 & 0.886 & 0.796 \\
1149    & 0.2 & 0.896 & 0.872 & 1.011 & 0.998 & 0.997 & 0.999 \\
1150    & 0.3 & 0.983 & 0.997 & 0.992 & 0.997 & 0.991 & 0.997 \\
1151 SF  & 0.0 & 0.943 & 0.979 & 1.048 & 0.978 & 0.995 & 0.999 \\
1152    & 0.1 & 0.948 & 0.979 & 1.044 & 0.983 & 1.000 & 0.999 \\
1153    & 0.2 & 0.982 & 0.997 & 0.969 & 0.960 & 0.997 & 0.999 \\
1154    & 0.3 & 0.985 & 0.997 & 0.961 & 0.961 & 0.991 & 0.997 \\
1155 GSC &     & 0.983 & 0.985 & 0.966 & 0.994 & 1.003 & 0.999 \\
1156 RF  &     & 0.924 & 0.944 & 0.990 & 0.996 & 0.991 & 0.998 \\
1157 \midrule
1158 PC  &     & -4.375 & 0.000 & 6.781 & 0.000 & -3.369 & 0.000 \\
1159 SP  & 0.0 & 0.515 & 0.164 & 0.856 & 0.426 & 0.743 & 0.478 \\
1160    & 0.1 & 0.696 & 0.405 & 0.977 & 0.817 & 0.974 & 0.964 \\
1161    & 0.2 & 0.981 & 0.980 & 1.001 & 1.000 & 1.000 & 1.000 \\
1162    & 0.3 & 0.996 & 0.998 & 0.997 & 0.999 & 0.997 & 0.999 \\
1163 SF  & 0.0 & 0.991 & 0.995 & 1.003 & 0.998 & 0.999 & 1.000 \\
1164    & 0.1 & 0.992 & 0.995 & 1.003 & 0.998 & 1.000 & 1.000 \\
1165    & 0.2 & 0.998 & 0.998 & 0.981 & 0.962 & 1.000 & 1.000 \\
1166    & 0.3 & 0.996 & 0.998 & 0.976 & 0.957 & 0.997 & 0.999 \\
1167 GSC &     & 0.997 & 0.996 & 0.998 & 0.999 & 1.000 & 1.000 \\
1168 RF  &     & 0.988 & 0.989 & 1.000 & 0.999 & 1.000 & 1.000 \\
1169 \midrule
1170 PC  &     & -6.367 & 0.000 & -3.552 & 0.000 & -3.447 & 0.000 \\
1171 SP  & 0.0 & 0.643 & 0.409 & 0.833 & 0.607 & 0.961 & 0.805 \\
1172    & 0.1 & 0.791 & 0.683 & 0.957 & 0.914 & 1.000 & 0.989 \\
1173    & 0.2 & 0.974 & 0.991 & 0.993 & 0.998 & 0.993 & 0.998 \\
1174    & 0.3 & 0.976 & 0.992 & 0.977 & 0.992 & 0.977 & 0.992 \\
1175 SF  & 0.0 & 0.979 & 0.997 & 0.992 & 0.999 & 0.994 & 1.000 \\
1176    & 0.1 & 0.984 & 0.997 & 0.996 & 0.999 & 0.998 & 1.000 \\
1177    & 0.2 & 0.991 & 0.997 & 0.974 & 0.958 & 0.993 & 0.998 \\
1178    & 0.3 & 0.977 & 0.992 & 0.956 & 0.948 & 0.977 & 0.992 \\
1179 GSC &     & 0.999 & 0.997 & 0.996 & 0.999 & 1.002 & 1.000 \\
1180 RF  &     & 0.994 & 0.997 & 0.997 & 0.999 & 1.000 & 1.000 \\
1181 \bottomrule
1182 \end{tabular}
1183 \label{tab:ice}
1184 \end{table}
1185
1186 \begin{table}[htbp]
1187 \centering
1188 \caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS
1189 OF THE FORCE AND TORQUE VECTORS IN THE ICE I$_\textrm{c}$ SYSTEM}      
1190
1191 \footnotesize
1192 \begin{tabular}{@{} ccrrrrrr @{}}
1193 \toprule
1194 \toprule
1195 & & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque
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 \\
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 \\
1203    & 0.1 & 590.008 & 107.510 & 18.883 & 118.201 & 32.472 & 3.599 \\
1204    & 0.2 & 10.057 & 0.105 & 0.038 & 2.875 & 0.572 & 0.518 \\
1205    & 0.3 & 0.245 & 0.260 & 0.262 & 2.365 & 2.396 & 2.327 \\
1206 SF  & 0.0 & 1.745 & 1.161 & 0.212 & 1.135 & 0.426 & 0.155 \\
1207    & 0.1 & 1.721 & 0.868 & 0.082 & 1.118 & 0.358 & 0.118 \\
1208    & 0.2 & 0.201 & 0.040 & 0.038 & 0.786 & 0.555 & 0.518 \\
1209    & 0.3 & 0.241 & 0.260 & 0.262 & 2.368 & 2.400 & 2.327 \\
1210 GSC &     & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\
1211 RF  &     & 2.887 & 0.217 & 0.107 & 1.006 & 0.281 & 0.085 \\
1212 \midrule
1213 GSSP  & 0.0 & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\
1214      & 0.1 & 1.341 & 0.123 & 0.037 & 0.835 & 0.234 & 0.085 \\
1215      & 0.2 & 0.558 & 0.040 & 0.037 & 0.823 & 0.557 & 0.519 \\
1216      & 0.3 & 0.250 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\
1217 GSSF  & 0.0 & 2.124 & 0.132 & 0.069 & 0.919 & 0.263 & 0.099 \\
1218      & 0.1 & 2.165 & 0.101 & 0.035 & 0.895 & 0.244 & 0.096 \\
1219      & 0.2 & 0.706 & 0.040 & 0.037 & 0.870 & 0.559 & 0.519 \\
1220      & 0.3 & 0.251 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\
1221 \bottomrule
1222 \end{tabular}
1223 \label{tab:iceAng}
1224 \end{table}
1225
1226 Highly ordered systems are a difficult test for the pairwise methods
1227 in that they lack the implicit periodicity of the Ewald summation.  As
1228 expected, the energy gap agreement with {\sc spme} is reduced for the
1229 {\sc sp} and {\sc sf} methods with parameters that were ideal for the
1230 disordered liquid system.  Moving to higher $R_\textrm{c}$ helps
1231 improve the agreement, though at an increase in computational cost.
1232 The dynamics of this crystalline system (both in magnitude and
1233 direction) are little affected. Both methods still reproduce the Ewald
1234 behavior with the same parameter recommendations from the previous
1235 section.
1236
1237 It is also worth noting that {\sc rf} exhibits improved energy gap
1238 results over the liquid water system.  One possible explanation is
1239 that the ice I$_\textrm{c}$ crystal is ordered such that the net
1240 dipole moment of the crystal is zero.  With $\epsilon_\textrm{S} =
1241 \infty$, the reaction field incorporates this structural organization
1242 by actively enforcing a zeroed dipole moment within each cutoff
1243 sphere.
1244
1245 \subsection{NaCl Melt Results}\label{sec:SaltMeltResults}
1246
1247 A high temperature NaCl melt was tested to gauge the accuracy of the
1248 pairwise summation methods in a disordered system of charges. The
1249 results for the energy gap comparisons and the force vector magnitude
1250 comparisons are shown in table \ref{tab:melt}.  The force vector
1251 directionality results are displayed separately in table
1252 \ref{tab:meltAng}.
1253
1254 \begin{table}[htbp]
1255 \centering
1256 \caption{REGRESSION RESULTS OF THE MOLTEN SODIUM CHLORIDE SYSTEM FOR
1257 $\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES ({\it
1258 lower})}
1259
1260 \footnotesize
1261 \begin{tabular}{@{} ccrrrrrr @{}}
1262 \toprule
1263 \toprule
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}
1268 Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1269 \midrule
1270 PC  &     & -0.008 & 0.000 & -0.049 & 0.005 & -0.136 & 0.020 \\
1271 SP  & 0.0 & 0.928 & 0.996 & 0.931 & 0.998 & 0.950 & 0.999 \\
1272    & 0.1 & 0.977 & 0.998 & 0.998 & 1.000 & 0.997 & 1.000 \\
1273    & 0.2 & 0.960 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\
1274    & 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\
1275 SF  & 0.0 & 0.996 & 1.000 & 0.995 & 1.000 & 0.997 & 1.000 \\
1276    & 0.1 & 1.021 & 1.000 & 1.024 & 1.000 & 1.007 & 1.000 \\
1277    & 0.2 & 0.966 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\
1278    & 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\
1279            \midrule
1280 PC  &     & 1.103 & 0.000 & 0.989 & 0.000 & 0.802 & 0.000 \\
1281 SP  & 0.0 & 0.973 & 0.981 & 0.975 & 0.988 & 0.979 & 0.992 \\
1282    & 0.1 & 0.987 & 0.992 & 0.993 & 0.998 & 0.997 & 0.999 \\
1283    & 0.2 & 0.993 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\
1284    & 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\
1285 SF  & 0.0 & 0.996 & 0.997 & 0.997 & 0.999 & 0.998 & 1.000 \\
1286    & 0.1 & 1.000 & 0.997 & 1.001 & 0.999 & 1.000 & 1.000 \\
1287    & 0.2 & 0.994 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\
1288    & 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\
1289 \bottomrule
1290 \end{tabular}
1291 \label{tab:melt}
1292 \end{table}
1293
1294 \begin{table}[htbp]
1295 \centering
1296 \caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS
1297 OF THE FORCE VECTORS IN THE MOLTEN SODIUM CHLORIDE SYSTEM}      
1298
1299 \footnotesize
1300 \begin{tabular}{@{} ccrrrrrr @{}}
1301 \toprule
1302 \toprule
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 \\
1307 \midrule
1308 PC  &     & 13.294 & 8.035 & 5.366 \\
1309 SP  & 0.0 & 13.316 & 8.037 & 5.385 \\
1310    & 0.1 & 5.705 & 1.391 & 0.360 \\
1311    & 0.2 & 2.415 & 7.534 & 13.927 \\
1312    & 0.3 & 23.769 & 67.306 & 57.252 \\
1313 SF  & 0.0 & 1.693 & 0.603 & 0.256 \\
1314    & 0.1 & 1.687 & 0.653 & 0.272 \\
1315    & 0.2 & 2.598 & 7.523 & 13.930 \\
1316    & 0.3 & 23.734 & 67.305 & 57.252 \\
1317 \bottomrule
1318 \end{tabular}
1319 \label{tab:meltAng}
1320 \end{table}
1321
1322 The molten NaCl system shows more sensitivity to the electrostatic
1323 damping than the water systems. The most noticeable point is that the
1324 undamped {\sc sf} method does very well at replicating the {\sc spme}
1325 configurational energy differences and forces. Light damping appears
1326 to minimally improve the dynamics, but this comes with a deterioration
1327 of the energy gap results. In contrast, this light damping improves
1328 the {\sc sp} energy gaps and forces. Moderate and heavy electrostatic
1329 damping reduce the agreement with {\sc spme} for both methods. From
1330 these observations, the undamped {\sc sf} method is the best choice
1331 for disordered systems of charges.
1332
1333 \subsection{NaCl Crystal Results}\label{sec:SaltCrystalResults}
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 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}.
1341 The force vector directionality results are displayed separately in
1342 table \ref{tab:saltAng}.
1343
1344 \begin{table}[htbp]
1345 \centering
1346 \caption{REGRESSION RESULTS OF THE CRYSTALLINE SODIUM CHLORIDE
1347 SYSTEM FOR $\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES
1348 ({\it lower})}
1349
1350 \footnotesize
1351 \begin{tabular}{@{} ccrrrrrr @{}}
1352 \toprule
1353 \toprule
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}
1358 Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1359 \midrule
1360 PC  &     & -20.241 & 0.228 & -20.248 & 0.229 & -20.239 & 0.228 \\
1361 SP  & 0.0 & 1.039 & 0.733 & 2.037 & 0.565 & 1.225 & 0.743 \\
1362    & 0.1 & 1.049 & 0.865 & 1.424 & 0.784 & 1.029 & 0.980 \\
1363    & 0.2 & 0.982 & 0.976 & 0.969 & 0.980 & 0.960 & 0.980 \\
1364    & 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.945 \\
1365 SF  & 0.0 & 1.041 & 0.967 & 0.994 & 0.989 & 0.957 & 0.993 \\
1366    & 0.1 & 1.050 & 0.968 & 0.996 & 0.991 & 0.972 & 0.995 \\
1367    & 0.2 & 0.982 & 0.975 & 0.959 & 0.980 & 0.960 & 0.980 \\
1368    & 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.944 \\
1369 \midrule
1370 PC  &     & 0.795 & 0.000 & 0.792 & 0.000 & 0.793 & 0.000 \\
1371 SP  & 0.0 & 0.916 & 0.829 & 1.086 & 0.791 & 1.010 & 0.936 \\
1372    & 0.1 & 0.958 & 0.917 & 1.049 & 0.943 & 1.001 & 0.995 \\
1373    & 0.2 & 0.981 & 0.981 & 0.982 & 0.984 & 0.981 & 0.984 \\
1374    & 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\
1375 SF  & 0.0 & 1.002 & 0.983 & 0.997 & 0.994 & 0.991 & 0.997 \\
1376    & 0.1 & 1.003 & 0.984 & 0.996 & 0.995 & 0.993 & 0.997 \\
1377    & 0.2 & 0.983 & 0.980 & 0.981 & 0.984 & 0.981 & 0.984 \\
1378    & 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\
1379 \bottomrule
1380 \end{tabular}
1381 \label{tab:salt}
1382 \end{table}
1383
1384 \begin{table}[htbp]
1385 \centering
1386 \caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR
1387 DISTRIBUTIONS OF THE FORCE VECTORS IN THE CRYSTALLINE SODIUM CHLORIDE
1388 SYSTEM}
1389
1390 \footnotesize
1391 \begin{tabular}{@{} ccrrrrrr @{}}
1392 \toprule
1393 \toprule
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 \\
1398 \midrule
1399 PC  &     & 111.945 & 111.824 & 111.866 \\
1400 SP  & 0.0 & 112.414 & 152.215 & 38.087 \\
1401    & 0.1 & 52.361 & 42.574 & 2.819 \\
1402    & 0.2 & 10.847 & 9.709 & 9.686 \\
1403    & 0.3 & 31.128 & 31.104 & 31.029 \\
1404 SF  & 0.0 & 10.025 & 3.555 & 1.648 \\
1405    & 0.1 & 9.462 & 3.303 & 1.721 \\
1406    & 0.2 & 11.454 & 9.813 & 9.701 \\
1407    & 0.3 & 31.120 & 31.105 & 31.029 \\
1408 \bottomrule
1409 \end{tabular}
1410 \label{tab:saltAng}
1411 \end{table}
1412
1413 The crystalline NaCl system is the most challenging test case for the
1414 pairwise summation methods, as evidenced by the results in tables
1415 \ref{tab:salt} and \ref{tab:saltAng}. The undamped and weakly damped
1416 {\sc sf} methods seem to be the best choices. These methods match well
1417 with {\sc spme} across the energy gap, force magnitude, and force
1418 directionality tests.  The {\sc sp} method struggles in all cases,
1419 with the exception of good dynamics reproduction when using weak
1420 electrostatic damping with a large cutoff radius.
1421
1422 The moderate electrostatic damping case is not as good as we would
1423 expect given the long-time dynamics results observed for this system
1424 (see section \ref{sec:LongTimeDynamics}). Since the data tabulated in
1425 tables \ref{tab:salt} and \ref{tab:saltAng} are a test of
1426 instantaneous dynamics, this indicates that good long-time dynamics
1427 comes in part at the expense of short-time dynamics.
1428
1429 \subsection{0.11M NaCl Solution Results}
1430
1431 In an effort to bridge the charged atomic and neutral molecular
1432 systems, Na$^+$ and Cl$^-$ ion charge defects were incorporated into
1433 the liquid water system. This low ionic strength system consists of 4
1434 ions in the 1000 SPC/E water solvent ($\approx$0.11 M). The results
1435 for the energy gap comparisons and the force and torque vector
1436 magnitude comparisons are shown in table \ref{tab:solnWeak}.  The
1437 force and torque vector directionality results are displayed
1438 separately in table \ref{tab:solnWeakAng}, where the effect of
1439 group-based cutoffs and switching functions on the {\sc sp} and {\sc
1440 sf} potentials are investigated.
1441
1442 \begin{table}[htbp]
1443 \centering
1444 \caption{REGRESSION RESULTS OF THE WEAK SODIUM CHLORIDE SOLUTION
1445 SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES
1446 ({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})}
1447
1448 \footnotesize
1449 \begin{tabular}{@{} ccrrrrrr @{}}
1450 \toprule
1451 \toprule
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}
1456 Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1457 \midrule
1458 PC  &     & 0.247 & 0.000 & -1.103 & 0.001 & 5.480 & 0.015 \\
1459 SP  & 0.0 & 0.935 & 0.388 & 0.984 & 0.541 & 1.010 & 0.685 \\
1460    & 0.1 & 0.951 & 0.603 & 0.993 & 0.875 & 1.001 & 0.979 \\
1461    & 0.2 & 0.969 & 0.968 & 0.996 & 0.997 & 0.994 & 0.997 \\
1462    & 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\
1463 SF  & 0.0 & 0.963 & 0.971 & 0.989 & 0.996 & 0.991 & 0.998 \\
1464    & 0.1 & 0.970 & 0.971 & 0.995 & 0.997 & 0.997 & 0.999 \\
1465    & 0.2 & 0.972 & 0.975 & 0.996 & 0.997 & 0.994 & 0.997 \\
1466    & 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\
1467 GSC &     & 0.964 & 0.731 & 0.984 & 0.704 & 1.005 & 0.770 \\
1468 RF  &     & 0.968 & 0.605 & 0.974 & 0.541 & 1.014 & 0.614 \\
1469 \midrule
1470 PC  &     & 1.354 & 0.000 & -1.190 & 0.000 & -0.314 & 0.000 \\
1471 SP  & 0.0 & 0.720 & 0.338 & 0.808 & 0.523 & 0.860 & 0.643 \\
1472    & 0.1 & 0.839 & 0.583 & 0.955 & 0.882 & 0.992 & 0.978 \\
1473    & 0.2 & 0.995 & 0.987 & 0.999 & 1.000 & 0.999 & 1.000 \\
1474    & 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\
1475 SF  & 0.0 & 0.998 & 0.994 & 1.000 & 0.998 & 1.000 & 0.999 \\
1476    & 0.1 & 0.997 & 0.994 & 1.000 & 0.999 & 1.000 & 1.000 \\
1477    & 0.2 & 0.999 & 0.998 & 0.999 & 1.000 & 0.999 & 1.000 \\
1478    & 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\
1479 GSC &     & 0.995 & 0.990 & 0.998 & 0.997 & 0.998 & 0.996 \\
1480 RF  &     & 0.998 & 0.993 & 0.999 & 0.998 & 0.999 & 0.996 \\
1481 \midrule
1482 PC  &     & 2.437 & 0.000 & -1.872 & 0.000 & 2.138 & 0.000 \\
1483 SP  & 0.0 & 0.838 & 0.525 & 0.901 & 0.686 & 0.932 & 0.779 \\
1484    & 0.1 & 0.914 & 0.733 & 0.979 & 0.932 & 0.995 & 0.987 \\
1485    & 0.2 & 0.977 & 0.969 & 0.988 & 0.990 & 0.989 & 0.990 \\
1486    & 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\
1487 SF  & 0.0 & 0.969 & 0.977 & 0.987 & 0.996 & 0.993 & 0.998 \\
1488    & 0.1 & 0.975 & 0.978 & 0.993 & 0.996 & 0.997 & 0.998 \\
1489    & 0.2 & 0.976 & 0.973 & 0.988 & 0.990 & 0.989 & 0.990 \\
1490    & 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\
1491 GSC &     & 0.980 & 0.959 & 0.990 & 0.983 & 0.992 & 0.989 \\
1492 RF  &     & 0.984 & 0.975 & 0.996 & 0.995 & 0.998 & 0.998 \\
1493 \bottomrule
1494 \end{tabular}
1495 \label{tab:solnWeak}
1496 \end{table}
1497
1498 \begin{table}[htbp]
1499 \centering
1500 \caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR
1501 DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE WEAK SODIUM
1502 CHLORIDE SOLUTION SYSTEM}
1503
1504 \footnotesize
1505 \begin{tabular}{@{} ccrrrrrr @{}}
1506 \toprule
1507 \toprule
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 \\
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 \\
1515    & 0.1 & 329.031 & 70.746 & 12.014 & 118.496 & 25.218 & 4.711 \\
1516    & 0.2 & 6.772 & 0.153 & 0.118 & 9.780 & 2.101 & 2.102 \\
1517    & 0.3 & 0.951 & 0.774 & 0.784 & 12.108 & 7.673 & 7.851 \\
1518 SF  & 0.0 & 2.555 & 0.762 & 0.313 & 6.590 & 1.328 & 0.558 \\
1519    & 0.1 & 2.561 & 0.560 & 0.123 & 6.464 & 1.162 & 0.457 \\
1520    & 0.2 & 0.501 & 0.118 & 0.118 & 5.698 & 2.074 & 2.099 \\
1521    & 0.3 & 0.943 & 0.774 & 0.784 & 12.118 & 7.674 & 7.851 \\
1522 GSC &     & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\
1523 RF  &     & 2.415 & 0.452 & 0.130 & 6.915 & 1.423 & 0.507 \\
1524 \midrule
1525 GSSP  & 0.0 & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\
1526      & 0.1 & 2.251 & 0.324 & 0.064 & 7.628 & 1.639 & 0.497 \\
1527      & 0.2 & 0.590 & 0.118 & 0.116 & 6.080 & 2.096 & 2.103 \\
1528      & 0.3 & 0.953 & 0.759 & 0.780 & 12.347 & 7.683 & 7.849 \\
1529 GSSF  & 0.0 & 1.541 & 0.301 & 0.096 & 6.407 & 1.316 & 0.496 \\
1530      & 0.1 & 1.541 & 0.237 & 0.050 & 6.356 & 1.202 & 0.457 \\
1531      & 0.2 & 0.568 & 0.118 & 0.116 & 6.166 & 2.105 & 2.105 \\
1532      & 0.3 & 0.954 & 0.759 & 0.780 & 12.337 & 7.684 & 7.849 \\
1533 \bottomrule
1534 \end{tabular}
1535 \label{tab:solnWeakAng}
1536 \end{table}
1537
1538 Because this system is a perturbation of the pure liquid water system,
1539 comparisons are best drawn between these two sets. The {\sc sp} and
1540 {\sc sf} methods are not significantly affected by the inclusion of a
1541 few ions. The aspect of cutoff sphere neutralization aids in the
1542 smooth incorporation of these ions; thus, all of the observations
1543 regarding these methods carry over from section
1544 \ref{sec:WaterResults}. The differences between these systems are more
1545 visible for the {\sc rf} method. Though good force agreement is still
1546 maintained, the energy gaps show a significant increase in the scatter
1547 of the data.
1548
1549 \subsection{1.1M NaCl Solution Results}
1550
1551 The bridging of the charged atomic and neutral molecular systems was
1552 further developed by considering a high ionic strength system
1553 consisting of 40 ions in the 1000 SPC/E water solvent ($\approx$1.1
1554 M). The results for the energy gap comparisons and the force and
1555 torque vector magnitude comparisons are shown in table
1556 \ref{tab:solnStr}.  The force and torque vector directionality
1557 results are displayed separately in table \ref{tab:solnStrAng}, where
1558 the effect of group-based cutoffs and switching functions on the {\sc
1559 sp} and {\sc sf} potentials are investigated.
1560
1561 \begin{table}[htbp]
1562 \centering
1563 \caption{REGRESSION RESULTS OF THE STRONG SODIUM CHLORIDE SOLUTION
1564 SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES
1565 ({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})}
1566
1567 \footnotesize
1568 \begin{tabular}{@{} ccrrrrrr @{}}
1569 \toprule
1570 \toprule
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}
1575 Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1576 \midrule
1577 PC  &     & -0.081 & 0.000 & 0.945 & 0.001 & 0.073 & 0.000 \\
1578 SP  & 0.0 & 0.978 & 0.469 & 0.996 & 0.672 & 0.975 & 0.668 \\
1579    & 0.1 & 0.944 & 0.645 & 0.997 & 0.886 & 0.991 & 0.978 \\
1580    & 0.2 & 0.873 & 0.896 & 0.985 & 0.993 & 0.980 & 0.993 \\
1581    & 0.3 & 0.831 & 0.860 & 0.960 & 0.979 & 0.955 & 0.977 \\
1582 SF  & 0.0 & 0.858 & 0.905 & 0.985 & 0.970 & 0.990 & 0.998 \\
1583    & 0.1 & 0.865 & 0.907 & 0.992 & 0.974 & 0.994 & 0.999 \\
1584    & 0.2 & 0.862 & 0.894 & 0.985 & 0.993 & 0.980 & 0.993 \\
1585    & 0.3 & 0.831 & 0.859 & 0.960 & 0.979 & 0.955 & 0.977 \\
1586 GSC &     & 1.985 & 0.152 & 0.760 & 0.031 & 1.106 & 0.062 \\
1587 RF  &     & 2.414 & 0.116 & 0.813 & 0.017 & 1.434 & 0.047 \\
1588 \midrule
1589 PC  &     & -7.028 & 0.000 & -9.364 & 0.000 & 0.925 & 0.865 \\
1590 SP  & 0.0 & 0.701 & 0.319 & 0.909 & 0.773 & 0.861 & 0.665 \\
1591    & 0.1 & 0.824 & 0.565 & 0.970 & 0.930 & 0.990 & 0.979 \\
1592    & 0.2 & 0.988 & 0.981 & 0.995 & 0.998 & 0.991 & 0.998 \\
1593    & 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\
1594 SF  & 0.0 & 0.993 & 0.988 & 0.992 & 0.984 & 0.998 & 0.999 \\
1595    & 0.1 & 0.993 & 0.989 & 0.993 & 0.986 & 0.998 & 1.000 \\
1596    & 0.2 & 0.993 & 0.992 & 0.995 & 0.998 & 0.991 & 0.998 \\
1597    & 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\
1598 GSC &     & 0.964 & 0.897 & 0.970 & 0.917 & 0.925 & 0.865 \\
1599 RF  &     & 0.994 & 0.864 & 0.988 & 0.865 & 0.980 & 0.784 \\
1600 \midrule
1601 PC  &     & -2.212 & 0.000 & -0.588 & 0.000 & 0.953 & 0.925 \\
1602 SP  & 0.0 & 0.800 & 0.479 & 0.930 & 0.804 & 0.924 & 0.759 \\
1603    & 0.1 & 0.883 & 0.694 & 0.976 & 0.942 & 0.993 & 0.986 \\
1604    & 0.2 & 0.952 & 0.943 & 0.980 & 0.984 & 0.980 & 0.983 \\
1605    & 0.3 & 0.914 & 0.909 & 0.943 & 0.948 & 0.944 & 0.946 \\
1606 SF  & 0.0 & 0.945 & 0.953 & 0.980 & 0.984 & 0.991 & 0.998 \\
1607    & 0.1 & 0.951 & 0.954 & 0.987 & 0.986 & 0.995 & 0.998 \\
1608    & 0.2 & 0.951 & 0.946 & 0.980 & 0.984 & 0.980 & 0.983 \\
1609    & 0.3 & 0.914 & 0.908 & 0.943 & 0.948 & 0.944 & 0.946 \\
1610 GSC &     & 0.882 & 0.818 & 0.939 & 0.902 & 0.953 & 0.925 \\
1611 RF  &     & 0.949 & 0.939 & 0.988 & 0.988 & 0.992 & 0.993 \\
1612 \bottomrule
1613 \end{tabular}
1614 \label{tab:solnStr}
1615 \end{table}
1616
1617 \begin{table}[htbp]
1618 \centering
1619 \caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS
1620 OF THE FORCE AND TORQUE VECTORS IN THE STRONG SODIUM CHLORIDE SOLUTION
1621 SYSTEM}
1622
1623 \footnotesize
1624 \begin{tabular}{@{} ccrrrrrr @{}}
1625 \toprule
1626 \toprule
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 \\
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 \\
1634    & 0.1 & 354.697 & 38.614 & 12.274 & 144.531 & 23.787 & 5.401 \\
1635    & 0.2 & 7.674 & 0.363 & 0.215 & 16.655 & 3.601 & 3.634 \\
1636    & 0.3 & 1.745 & 1.456 & 1.449 & 23.669 & 14.376 & 14.240 \\
1637 SF  & 0.0 & 3.282 & 8.567 & 0.369 & 11.904 & 6.589 & 0.717 \\
1638    & 0.1 & 3.263 & 7.479 & 0.142 & 11.634 & 5.750 & 0.591 \\
1639    & 0.2 & 0.686 & 0.324 & 0.215 & 10.809 & 3.580 & 3.635 \\
1640    & 0.3 & 1.749 & 1.456 & 1.449 & 23.635 & 14.375 & 14.240 \\
1641 GSC &     & 6.181 & 2.904 & 2.263 & 44.349 & 19.442 & 12.873 \\
1642 RF  &     & 3.891 & 0.847 & 0.323 & 18.628 & 3.995 & 2.072 \\
1643 \midrule
1644 GSSP  & 0.0 & 6.197 & 2.929 & 2.290 & 44.441 & 19.442 & 12.873 \\
1645      & 0.1 & 4.688 & 1.064 & 0.260 & 31.208 & 6.967 & 2.303 \\
1646      & 0.2 & 1.021 & 0.218 & 0.213 & 14.425 & 3.629 & 3.649 \\
1647      & 0.3 & 1.752 & 1.454 & 1.451 & 23.540 & 14.390 & 14.245 \\
1648 GSSF  & 0.0 & 2.494 & 0.546 & 0.217 & 16.391 & 3.230 & 1.613 \\
1649      & 0.1 & 2.448 & 0.429 & 0.106 & 16.390 & 2.827 & 1.159 \\
1650      & 0.2 & 0.899 & 0.214 & 0.213 & 13.542 & 3.583 & 3.645 \\
1651      & 0.3 & 1.752 & 1.454 & 1.451 & 23.587 & 14.390 & 14.245 \\
1652 \bottomrule
1653 \end{tabular}
1654 \label{tab:solnStrAng}
1655 \end{table}
1656
1657 The {\sc rf} method struggles with the jump in ionic strength. The
1658 configuration energy differences degrade to unusable levels while the
1659 forces and torques show a more modest reduction in the agreement with
1660 {\sc spme}. The {\sc rf} method was designed for homogeneous systems,
1661 and this attribute is apparent in these results.
1662
1663 The {\sc sp} and {\sc sf} methods require larger cutoffs to maintain
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.
1668
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
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
1675 comparisons are shown in table \ref{tab:argon}.  The force and torque
1676 vector directionality results are displayed separately in table
1677 \ref{tab:argonAng}, where the effect of group-based cutoffs and
1678 switching functions on the {\sc sp} and {\sc sf} potentials are
1679 investigated.
1680
1681 \begin{table}[htbp]
1682 \centering
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
1687 \footnotesize
1688 \begin{tabular}{@{} ccrrrrrr @{}}
1689 \toprule
1690 \toprule
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}
1695 Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\
1696 \midrule
1697 PC  &     & 2.320 & 0.008 & -0.650 & 0.001 & 3.848 & 0.029 \\
1698 SP  & 0.0 & 1.053 & 0.711 & 0.977 & 0.820 & 0.974 & 0.882 \\
1699    & 0.1 & 1.032 & 0.846 & 0.989 & 0.965 & 0.992 & 0.994 \\
1700    & 0.2 & 0.993 & 0.995 & 0.982 & 0.998 & 0.986 & 0.998 \\
1701    & 0.3 & 0.968 & 0.995 & 0.954 & 0.992 & 0.961 & 0.994 \\
1702 SF  & 0.0 & 0.982 & 0.996 & 0.992 & 0.999 & 0.993 & 1.000 \\
1703    & 0.1 & 0.987 & 0.996 & 0.996 & 0.999 & 0.997 & 1.000 \\
1704    & 0.2 & 0.989 & 0.998 & 0.984 & 0.998 & 0.989 & 0.998 \\
1705    & 0.3 & 0.971 & 0.995 & 0.957 & 0.992 & 0.965 & 0.994 \\
1706 GSC &     & 1.002 & 0.983 & 0.992 & 0.973 & 0.996 & 0.971 \\
1707 RF  &     & 0.998 & 0.995 & 0.999 & 0.998 & 0.998 & 0.998 \\
1708 \midrule
1709 PC  &     & -36.559 & 0.002 & -44.917 & 0.004 & -52.945 & 0.006 \\
1710 SP  & 0.0 & 0.890 & 0.786 & 0.927 & 0.867 & 0.949 & 0.909 \\
1711    & 0.1 & 0.942 & 0.895 & 0.984 & 0.974 & 0.997 & 0.995 \\
1712    & 0.2 & 0.999 & 0.997 & 1.000 & 1.000 & 1.000 & 1.000 \\
1713    & 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\
1714 SF  & 0.0 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\
1715    & 0.1 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\
1716    & 0.2 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\
1717    & 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\
1718 GSC &     & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\
1719 RF  &     & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\
1720 \midrule
1721 PC  &     & 1.984 & 0.000 & 0.012 & 0.000 & 1.357 & 0.000 \\
1722 SP  & 0.0 & 0.850 & 0.552 & 0.907 & 0.703 & 0.938 & 0.793 \\
1723    & 0.1 & 0.924 & 0.755 & 0.980 & 0.936 & 0.995 & 0.988 \\
1724    & 0.2 & 0.985 & 0.983 & 0.986 & 0.988 & 0.987 & 0.988 \\
1725    & 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\
1726 SF  & 0.0 & 0.977 & 0.989 & 0.987 & 0.995 & 0.992 & 0.998 \\
1727    & 0.1 & 0.982 & 0.989 & 0.992 & 0.996 & 0.997 & 0.998 \\
1728    & 0.2 & 0.984 & 0.987 & 0.986 & 0.987 & 0.987 & 0.988 \\
1729    & 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\
1730 GSC &     & 0.995 & 0.981 & 0.999 & 0.990 & 1.000 & 0.993 \\
1731 RF  &     & 0.993 & 0.988 & 0.997 & 0.995 & 0.999 & 0.998 \\
1732 \bottomrule
1733 \end{tabular}
1734 \label{tab:argon}
1735 \end{table}
1736
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
1741 ARGON IN LIQUID WATER SYSTEM}  
1742
1743 \footnotesize
1744 \begin{tabular}{@{} ccrrrrrr @{}}
1745 \toprule
1746 \toprule
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 \\
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 \\
1754    & 0.1 & 224.886 & 49.746 & 9.346 & 104.482 & 23.683 & 4.480 \\
1755    & 0.2 & 4.889 & 0.197 & 0.155 & 6.029 & 2.507 & 2.269 \\
1756    & 0.3 & 0.817 & 0.833 & 0.812 & 8.286 & 8.436 & 8.135 \\
1757 SF  & 0.0 & 1.924 & 0.675 & 0.304 & 3.658 & 1.448 & 0.600 \\
1758    & 0.1 & 1.937 & 0.515 & 0.143 & 3.565 & 1.308 & 0.546 \\
1759    & 0.2 & 0.407 & 0.166 & 0.156 & 3.086 & 2.501 & 2.274 \\
1760    & 0.3 & 0.815 & 0.833 & 0.812 & 8.330 & 8.437 & 8.135 \\
1761 GSC &     & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\
1762 RF  &     & 1.822 & 0.408 & 0.142 & 3.799 & 1.362 & 0.550 \\
1763 \midrule
1764 GSSP  & 0.0 & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\
1765      & 0.1 & 1.652 & 0.309 & 0.087 & 4.197 & 1.401 & 0.590 \\
1766      & 0.2 & 0.465 & 0.165 & 0.153 & 3.323 & 2.529 & 2.273 \\
1767      & 0.3 & 0.813 & 0.825 & 0.816 & 8.316 & 8.447 & 8.132 \\
1768 GSSF  & 0.0 & 1.173 & 0.292 & 0.113 & 3.452 & 1.347 & 0.583 \\
1769      & 0.1 & 1.166 & 0.240 & 0.076 & 3.381 & 1.281 & 0.575 \\
1770      & 0.2 & 0.459 & 0.165 & 0.153 & 3.430 & 2.542 & 2.273 \\
1771      & 0.3 & 0.814 & 0.825 & 0.816 & 8.325 & 8.447 & 8.132 \\
1772 \bottomrule
1773 \end{tabular}
1774 \label{tab:argonAng}
1775 \end{table}
1776
1777 This system does not appear to show any significant deviations from
1778 the previously observed results. The {\sc sp} and {\sc sf} methods
1779 have agreements similar to those observed in section
1780 \ref{sec:WaterResults}. The only significant difference is the
1781 improvement in the configuration energy differences for the {\sc rf}
1782 method. This is surprising in that we are introducing an inhomogeneity
1783 to the system; however, this inhomogeneity is charge-neutral and does
1784 not result in charged cutoff spheres. The charge-neutrality of the
1785 cutoff spheres, which the {\sc sp} and {\sc sf} methods explicitly
1786 enforce, seems to play a greater role in the stability of the {\sc rf}
1787 method than the required homogeneity of the environment.
1788
1789
940   \section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}\label{sec:ShortTimeDynamics}
941  
942   Zahn {\it et al.} investigated the structure and dynamics of water
# Line 1822 | Line 972 | are stiffer than the moderately damped and {\sc spme}
972   are stiffer than the moderately damped and {\sc spme} methods.}
973   \label{fig:vCorrPlot}
974   \end{figure}
1825
975   The short-time decay of the velocity autocorrelation function through
976   the first collision are nearly identical in figure
977   \ref{fig:vCorrPlot}, but the peaks and troughs of the functions show
# Line 1841 | Line 990 | important.
990  
991   \section{Collective Motion: Power Spectra of NaCl Crystals}\label{sec:LongTimeDynamics}
992  
1844 To evaluate how the differences between the methods affect the
1845 collective long-time motion, we computed power spectra from long-time
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 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
993   \begin{figure}
994   \centering
995   \includegraphics[width = \linewidth]{./figures/spectraSquare.pdf}
# Line 1860 | Line 1000 | functions of NaCl crystals at 1000~K while using {\sc
1000   100~cm$^{-1}$ to highlight where the spectra differ.}
1001   \label{fig:methodPS}
1002   \end{figure}
1003 + To evaluate how the differences between the methods affect the
1004 + collective long-time motion, we computed power spectra from long-time
1005 + traces of the velocity autocorrelation function. The power spectra for
1006 + the best-performing alternative methods are shown in
1007 + fig. \ref{fig:methodPS}.  Apodization of the correlation functions via
1008 + a cubic switching function between 40 and 50~ps was used to reduce the
1009 + ringing resulting from data truncation.  This procedure had no
1010 + noticeable effect on peak location or magnitude.
1011  
1012   While the high frequency regions of the power spectra for the
1013   alternative methods are quantitatively identical with Ewald spectrum,
1014   the low frequency region shows how the summation methods differ.
1015   Considering the low-frequency inset (expanded in the upper frame of
1016 < figure \ref{fig:dampInc}), at frequencies below 100~cm$^{-1}$, the
1016 > figure \ref{fig:methodPS}), at frequencies below 100~cm$^{-1}$, the
1017   correlated motions are blue-shifted when using undamped or weakly
1018   damped {\sc sf}.  When using moderate damping ($\alpha =
1019   0.2$~\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly
# Line 1875 | Line 1023 | moderately damped methods than for undamped or weakly
1023   long-ranged correlated motions are at lower frequencies for the
1024   moderately damped methods than for undamped or weakly damped methods.
1025  
1026 + \begin{figure}
1027 + \centering
1028 + \includegraphics[width = \linewidth]{./figures/increasedDamping.pdf}
1029 + \caption{Effect of damping on the two lowest-frequency phonon modes in
1030 + the NaCl crystal at 1000~K.  The undamped shifted force ({\sc sf})
1031 + method is off by less than 10~cm$^{-1}$, and increasing the
1032 + electrostatic damping to 0.25~\AA$^{-1}$ gives quantitative agreement
1033 + with the power spectrum obtained using the Ewald sum.  Over-damping can
1034 + result in underestimates of frequencies of the long-wavelength
1035 + motions.}
1036 + \label{fig:dampInc}
1037 + \end{figure}
1038   To isolate the role of the damping constant, we have computed the
1039   spectra for a single method ({\sc sf}) with a range of damping
1040   constants and compared this with the {\sc spme} spectrum.
# Line 1889 | Line 1049 | cutoff distance.
1049   obtained using moderate damping in addition to the shifting at the
1050   cutoff distance.
1051  
1892 \begin{figure}
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 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.}
1902 \label{fig:dampInc}
1903 \end{figure}
1904
1052   \section{An Application: TIP5P-E Water}\label{sec:t5peApplied}
1053  
1054   The above sections focused on the energetics and dynamics of a variety
# Line 1959 | Line 1106 | p(t) =         \frac{1}{\textrm{d}V}\sum_\mu
1106          + \mathbf{R}_{\mu}\cdot\sum_i\mathbf{F}_{\mu i}\right],
1107   \label{eq:MolecularPressure}
1108   \end{equation}
1109 < where $V$ is the volume, $\mathbf{P}_{\mu}$ is the momentum of
1110 < molecule $\mu$, $\mathbf{R}_\mu$ is the position of the center of mass
1111 < ($M_\mu$) of molecule $\mu$, and $\mathbf{F}_{\mu i}$ is the force on
1112 < atom $i$ of molecule $\mu$.\cite{Melchionna93} The virial term (the
1113 < right term in the brackets of equation \ref{eq:MolecularPressure}) is
1114 < directly dependent on the interatomic forces.  Since the {\sc sp}
1115 < method does not modify the forces (see
1116 < section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp} will
1117 < be identical to that obtained without an electrostatic correction.
1118 < The {\sc sf} method does alter the virial component and, by way of the
1119 < modified pressures, should provide densities more in line with those
1120 < obtained using the Ewald summation.
1109 > where d is the dimensionality of the system, $V$ is the volume,
1110 > $\mathbf{P}_{\mu}$ is the momentum of molecule $\mu$, $\mathbf{R}_\mu$
1111 > is the position of the center of mass ($M_\mu$) of molecule $\mu$, and
1112 > $\mathbf{F}_{\mu i}$ is the force on atom $i$ of molecule
1113 > $\mu$.\cite{Melchionna93} The virial term (the right term in the
1114 > brackets of equation
1115 > \ref{eq:MolecularPressure}) is directly dependent on the interatomic
1116 > forces.  Since the {\sc sp} method does not modify the forces (see
1117 > section. \ref{sec:PairwiseDerivation}), the pressure using {\sc sp}
1118 > will be identical to that obtained without an electrostatic
1119 > correction.  The {\sc sf} method does alter the virial component and,
1120 > by way of the modified pressures, should provide densities more in
1121 > line with those obtained using the Ewald summation.
1122  
1123   To compare densities, $NPT$ simulations were performed with the same
1124   temperatures as those selected by Rick in his Ewald summation
# Line 1992 | Line 1140 | greater force on the central particle. The error bars
1140   Ewald summation, leading to slightly lower densities. This effect is
1141   more visible with the 9~\AA\ cutoff, where the image charges exert a
1142   greater force on the central particle. The error bars for the {\sc sf}
1143 < methods show plus or minus the standard deviation of the density
1144 < measurement at each temperature.}
1143 > methods show the average one-sigma uncertainty of the density
1144 > measurement, and this uncertainty is the same for all the {\sc sf}
1145 > curves.}
1146   \label{fig:t5peDensities}
1147   \end{figure}
1999
1148   Figure \ref{fig:t5peDensities} shows the densities calculated for
1149   TIP5P-E using differing electrostatic corrections overlaid on the
1150   experimental values.\cite{CRC80} The densities when using the {\sc sf}
# Line 2021 | Line 1169 | As a final note, all of the above density calculations
1169   important role in the resulting densities.
1170  
1171   As a final note, all of the above density calculations were performed
1172 < with systems of 512 water molecules. Rick observed a system sized
1172 > with systems of 512 water molecules. Rick observed a system size
1173   dependence of the computed densities when using the Ewald summation,
1174   most likely due to his tying of the convergence parameter to the box
1175   dimensions.\cite{Rick04} For systems of 256 water molecules, the
# Line 2076 | Line 1224 | identical.}
1224   identical.}
1225   \label{fig:t5peGofRs}
1226   \end{figure}
2079
1227   The $g_\textrm{OO}(r)$s calculated for TIP5P-E while using the {\sc
1228   sf} technique with a various parameters are overlaid on the
1229 < $g_\textrm{OO}(r)$ while using the Ewald summation. The differences in
1230 < density do not appear to have any effect on the liquid structure as
1231 < the $g_\textrm{OO}(r)$s are indistinguishable. These results indicate
1232 < that the $g_\textrm{OO}(r)$ is insensitive to the choice of
1233 < electrostatic correction.
1229 > $g_\textrm{OO}(r)$ while using the Ewald summation in figure
1230 > \ref{fig:t5peGofRs}. The differences in density do not appear to have
1231 > any effect on the liquid structure as the $g_\textrm{OO}(r)$s are
1232 > indistinguishable. These results indicate that the $g_\textrm{OO}(r)$
1233 > is insensitive to the choice of electrostatic correction.
1234  
1235   \subsection{Thermodynamic Properties}\label{sec:t5peThermo}
1236  
# Line 2213 | Line 1360 | property trends with temperature seen when using the E
1360  
1361   As observed for the density in section \ref{sec:t5peDensity}, the
1362   property trends with temperature seen when using the Ewald summation
1363 < are reproduced with the {\sc sf} technique. Differences include the
1364 < calculated values of $\Delta H_\textrm{vap}$ under-predicting the Ewald
1365 < values. This is to be expected due to the direct weakening of the
1366 < electrostatic interaction through forced neutralization in {\sc
1367 < sf}. This results in an increase of the intermolecular potential
1368 < producing lower values from equation (\ref{eq:DeltaHVap}). The slopes of
1369 < these values with temperature are similar to that seen using the Ewald
1370 < summation; however, they are both steeper than the experimental trend,
1371 < indirectly resulting in the inflated $C_p$ values at all temperatures.
1363 > are reproduced with the {\sc sf} technique. One noticable difference
1364 > between the properties calculated using the two methods are the lower
1365 > $\Delta H_\textrm{vap}$ values when using {\sc sf}. This is to be
1366 > expected due to the direct weakening of the electrostatic interaction
1367 > through forced neutralization. This results in an increase of the
1368 > intermolecular potential producing lower values from equation
1369 > (\ref{eq:DeltaHVap}). The slopes of these values with temperature are
1370 > similar to that seen using the Ewald summation; however, they are both
1371 > steeper than the experimental trend, indirectly resulting in the
1372 > inflated $C_p$ values at all temperatures.
1373  
1374   Above the supercooled regime, $C_p$, $\kappa_T$, and $\alpha_p$
1375   values all overlap within error. As indicated for the $\Delta
# Line 2243 | Line 1391 | the onset of a more frustrated or glassy behavior for
1391   indicate a more pronounced transition in the supercooled regime,
1392   particularly in the case of {\sc sf} without damping. This points to
1393   the onset of a more frustrated or glassy behavior for TIP5P-E at
1394 < temperatures below 250~K in these simulations. Because the systems are
1395 < locked in different regions of phase-space, comparisons between
1396 < properties at these temperatures are not exactly fair. This
1397 < observation is explored in more detail in section
1398 < \ref{sec:t5peDynamics}.
1394 > temperatures below 250~K in the {\sc sf} simulations, indicating that
1395 > disorder in the reciprical-space term of the Ewald summation might act
1396 > to loosen up the local structure more than the image-charges in {\sc
1397 > sf}. Because the systems are locked in different regions of
1398 > phase-space, comparisons between properties at these temperatures are
1399 > not exactly fair. This observation is explored in more detail in
1400 > section \ref{sec:t5peDynamics}.
1401  
1402   The final thermodynamic property displayed in figure
1403   \ref{fig:t5peThermo}, $\epsilon$, shows the greatest discrepancy
# Line 2280 | Line 1430 | relation using the mean square displacement (MSD),
1430   self-diffusion constants ($D$) were calculated with the Einstein
1431   relation using the mean square displacement (MSD),
1432   \begin{equation}
1433 < D = \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
1433 > D = \lim_{t\rightarrow\infty}
1434 >    \frac{\langle\left|\mathbf{r}_i(t)-\mathbf{r}_i(0)\right|^2\rangle}{6t},
1435   \label{eq:MSD}
1436   \end{equation}
1437   where $t$ is time, and $\mathbf{r}_i$ is the position of particle
# Line 2291 | Line 1442 | regions:
1442   \begin{enumerate}[itemsep=0pt]
1443   \item parabolic short-time ballistic motion,
1444   \item linear diffusive regime, and
1445 < \item poor statistic region at long-time.
1445 > \item a region with poor statistics.
1446   \end{enumerate}
1447   The slope from the linear region (region 2) is used to calculate $D$.
1448   \begin{figure}
# Line 2360 | Line 1511 | easier comparisons in the more relevant temperature re
1511   easier comparisons in the more relevant temperature regime.}
1512   \label{fig:t5peDynamics}
1513   \end{figure}
1514 < Results for the diffusion constants and reorientational time constants
1514 > Results for the diffusion constants and orientational relaxation times
1515   are shown in figure \ref{fig:t5peDynamics}. From this figure, it is
1516   apparent that the trends for both $D$ and $\tau_2^y$ of TIP5P-E using
1517   the Ewald sum are reproduced with the {\sc sf} technique. The enhanced
# Line 2369 | Line 1520 | diffuse a little faster than with the Ewald sum; howev
1520   insight into differences between the electrostatic summation
1521   techniques. With the undamped {\sc sf} technique, TIP5P-E tends to
1522   diffuse a little faster than with the Ewald sum; however, use of light
1523 < to moderate damping results in indistinguishable $D$ values. Though not
1524 < apparent in this figure, {\sc sf} values at the lowest temperature are
1525 < approximately an order of magnitude lower than with Ewald. These
1523 > to moderate damping results in indistinguishable $D$ values. Though
1524 > not apparent in this figure, {\sc sf} values at the lowest temperature
1525 > are approximately an order of magnitude lower than with Ewald. These
1526   values support the observation from section \ref{sec:t5peThermo} that
1527   there appeared to be a change to a more glassy-like phase with the
1528   {\sc sf} technique at these lower temperatures.
# Line 2385 | Line 1536 | calculation from a 10~ps trajectory with {\sc sf} with
1536   for this deviation between techniques. The Ewald results were taken
1537   shorter (10ps) trajectories than the {\sc sf} results (25ps). A quick
1538   calculation from a 10~ps trajectory with {\sc sf} with an $\alpha$ of
1539 < 0.2~\AA$^-1$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$, placing
1540 < the result more in line with that obtained using the Ewald sum. These
1541 < results support this explanation; however, recomputing the results to
1542 < meet a poorer statistical standard is counter-productive. Assuming the
1543 < Ewald results are not the product of poor statistics, differences in
1544 < techniques to integrate the orientational motion could also play a
1545 < role. {\sc shake} is the most commonly used technique for
1546 < approximating rigid-body orientational motion,\cite{Ryckaert77} where
1547 < as in {\sc oopse}, we maintain and integrate the entire rotation
1548 < matrix using the {\sc dlm} method.\cite{Meineke05} Since {\sc shake}
1549 < is an iterative constraint technique, if the convergence tolerances
1550 < are raised for increased performance, error will accumulate in the
1551 < orientational motion. Finally, the Ewald results were calculated using
1552 < the $NVT$ ensemble, while the $NVE$ ensemble was used for {\sc sf}
1539 > 0.2~\AA$^{-1}$ at 25$^\circ$C showed a 0.4~ps drop in $\tau_2^y$,
1540 > placing the result more in line with that obtained using the Ewald
1541 > sum. These results support this explanation; however, recomputing the
1542 > results to meet a poorer statistical standard is
1543 > counter-productive. Assuming the Ewald results are not the product of
1544 > poor statistics, differences in techniques to integrate the
1545 > orientational motion could also play a role. {\sc shake} is the most
1546 > commonly used technique for approximating rigid-body orientational
1547 > motion,\cite{Ryckaert77} where as in {\sc oopse}, we maintain and
1548 > integrate the entire rotation matrix using the {\sc dlm}
1549 > method.\cite{Meineke05} Since {\sc shake} is an iterative constraint
1550 > technique, if the convergence tolerances are raised for increased
1551 > performance, error will accumulate in the orientational
1552 > motion. Finally, the Ewald results were calculated using the $NVT$
1553 > ensemble, while the $NVE$ ensemble was used for {\sc sf}
1554   calculations. The additional mode of motion due to the thermostat will
1555   alter the dynamics, resulting in differences between $NVT$ and $NVE$
1556   results. These differences are increasingly noticeable as the
# Line 2410 | Line 1562 | consider an extension of these techniques to include p
1562   neutralizing the cutoff sphere with charge-charge interaction shifting
1563   and by damping the electrostatic interactions. Now we would like to
1564   consider an extension of these techniques to include point multipole
1565 < interactions. How will the shifting and damping need to develop in
1565 > interactions. How will the shifting and damping need to be modified in
1566   order to accommodate point multipoles?
1567  
1568 < Of the two techniques, the least to vary is shifting. Shifting is
1568 > Of the two techniques, the easiest to adapt is shifting. Shifting is
1569   employed to neutralize the cutoff sphere; however, in a system
1570   composed purely of point multipoles, the cutoff sphere is already
1571   neutralized. This means that shifting is not necessary between point
# Line 2431 | Line 1583 | than considering only the interactions between single
1583   replacing $r^{-1}$ with erfc$(\alpha r)\cdot r^{-1}$ in the multipole
1584   expansion.\cite{Hirschfelder67} In the multipole expansion, rather
1585   than considering only the interactions between single point charges,
1586 < the electrostatic interactions is reformulated such that it describes
1586 > the electrostatic interaction is reformulated such that it describes
1587   the interaction between charge distributions about central sites of
1588   the respective sets of charges. This procedure is what leads to the
1589   familiar charge-dipole,
# Line 2600 | Line 1752 | ranging from 0 to 0.35~\AA$^{-1}$.
1752   \ref{sec:dampingMultipoles}. Each of these systems were studied with
1753   cutoff radii of 9, 10, 11, and 12~\AA\ and with damping parameter values
1754   ranging from 0 to 0.35~\AA$^{-1}$.
1755 +
1756   \begin{figure}
1757   \centering
1758   \includegraphics[width=\linewidth]{./figures/dielectricMap.pdf}
# Line 2608 | Line 1761 | radius ($R_\textrm{c}$) and damping coefficient ($\alp
1761   radius ($R_\textrm{c}$) and damping coefficient ($\alpha$).}
1762   \label{fig:dielectricMap}
1763   \end{figure}
2611
1764   The results of these calculations are displayed in figure
1765   \ref{fig:dielectricMap} in the form of shaded contour plots. An
1766   interesting aspect of all four contour plots is that the dielectric
# Line 2635 | Line 1787 | these Ewald examples, the results discussed in section
1787  
1788   Although it is tempting to choose damping parameters equivalent to
1789   these Ewald examples, the results discussed in sections
1790 < \ref{sec:EnergyResults} through \ref{sec:IndividualResults} indicate
1791 < that values this high are destructive to both the energetics and
1792 < dynamics. Ideally, $\alpha$ should not exceed 0.3~\AA$^{-1}$ for any of
1793 < the cutoff values in this range. If the optimal damping parameter is
1794 < chosen to be midway between 0.275 and 0.3~\AA$^{-1}$ (0.2875~\AA$^{-1}$)
1795 < at the 9~\AA\ cutoff, then the linear trend with $R_\textrm{c}$ will
1796 < always keep $\alpha$ below 0.3~\AA$^{-1}$. This linear progression
1797 < would give values of 0.2875, 0.2625, 0.2375, and 0.2125~\AA$^{-1}$ for
1798 < cutoff radii of 9, 10, 11, and 12~\AA. Setting this to be the default
1799 < behavior for the damped {\sc sf} technique will result in consistent
1800 < dielectric behavior for these and other condensed molecular systems,
1801 < regardless of the chosen cutoff radius. The static dielectric
1802 < constants for TIP5P-E, TIP4P-Ew, SPC/E, and SSD/RF will be
1803 < approximately fixed at 74, 52, 58, and 89 respectively. These values
1804 < are generally lower than the values reported in the literature;
1805 < however, the relative dielectric behavior scales as expected when
1806 < comparing the models to one another.
1790 > \ref{sec:EnergyResults} through \ref{sec:FTDirResults} and appendix
1791 > \ref{app:IndividualResults} indicate that values this high are
1792 > destructive to both the energetics and dynamics. Ideally, $\alpha$
1793 > should not exceed 0.3~\AA$^{-1}$ for any of the cutoff values in this
1794 > range. If the optimal damping parameter is chosen to be midway between
1795 > 0.275 and 0.3~\AA$^{-1}$ (0.2875~\AA$^{-1}$) at the 9~\AA\ cutoff,
1796 > then the linear trend with $R_\textrm{c}$ will always keep $\alpha$
1797 > below 0.3~\AA$^{-1}$. This linear progression would give values of
1798 > 0.2875, 0.2625, 0.2375, and 0.2125~\AA$^{-1}$ for cutoff radii of 9,
1799 > 10, 11, and 12~\AA. Setting this to be the default behavior for the
1800 > damped {\sc sf} technique will result in consistent dielectric
1801 > behavior for these and other condensed molecular systems, regardless
1802 > of the chosen cutoff radius. The static dielectric constants for
1803 > TIP5P-E, TIP4P-Ew, SPC/E, and SSD/RF will be fixed at approximately
1804 > 74, 52, 58, and 89 respectively. These values are generally lower than
1805 > the values reported in the literature; however, the relative
1806 > dielectric behavior scales as expected when comparing the models to
1807 > one another.
1808  
1809   \section{Conclusions}\label{sec:PairwiseConclusions}
1810  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines