--- trunk/fennellDissertation/dissertation.tex 2006/07/08 13:13:27 2927 +++ trunk/fennellDissertation/dissertation.tex 2006/08/18 15:49:29 2973 @@ -1,12 +1,15 @@ -\documentclass[12pt]{ndthesis} +\documentclass[11pt]{ndthesis} % some packages for things like equations and graphics +\usepackage[tbtags]{amsmath} \usepackage{amsmath,bm} \usepackage{amssymb} \usepackage{mathrsfs} \usepackage{tabularx} \usepackage{graphicx} \usepackage{booktabs} +\usepackage{cite} +\usepackage{enumitem} \begin{document} @@ -39,1983 +42,9 @@ STUDY OF WATER} \mainmatter \chapter{\label{chap:intro}INTRODUCTION AND BACKGROUND} - -\chapter{\label{chap:electrostatics}ELECTROSTATIC INTERACTION CORRECTION -TECHNIQUES} - -In molecular simulations, proper accumulation of the electrostatic -interactions is essential and is one of the most -computationally-demanding tasks. The common molecular mechanics force -fields represent atomic sites with full or partial charges protected -by Lennard-Jones (short range) interactions. This means that nearly -every pair interaction involves a calculation of charge-charge forces. -Coupled with relatively long-ranged $r^{-1}$ decay, the monopole -interactions quickly become the most expensive part of molecular -simulations. Historically, the electrostatic pair interaction would -not have decayed appreciably within the typical box lengths that could -be feasibly simulated. In the larger systems that are more typical of -modern simulations, large cutoffs should be used to incorporate -electrostatics correctly. - -There have been many efforts to address the proper and practical -handling of electrostatic interactions, and these have resulted in a -variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are -typically classified as implicit methods (i.e., continuum dielectrics, -static dipolar fields),\cite{Born20,Grossfield00} explicit methods -(i.e., Ewald summations, interaction shifting or -truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., -reaction field type methods, fast multipole -methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are -often preferred because they physically incorporate solvent molecules -in the system of interest, but these methods are sometimes difficult -to utilize because of their high computational cost.\cite{Roux99} In -addition to the computational cost, there have been some questions -regarding possible artifacts caused by the inherent periodicity of the -explicit Ewald summation.\cite{Tobias01} - -In this chapter, we focus on a new set of pairwise methods devised by -Wolf {\it et al.},\cite{Wolf99} which we further extend. These -methods along with a few other mixed methods (i.e. reaction field) are -compared with the smooth particle mesh Ewald -sum,\cite{Onsager36,Essmann99} which is our reference method for -handling long-range electrostatic interactions. The new methods for -handling electrostatics have the potential to scale linearly with -increasing system size since they involve only a simple modification -to the direct pairwise sum. They also lack the added periodicity of -the Ewald sum, so they can be used for systems which are non-periodic -or which have one- or two-dimensional periodicity. Below, these -methods are evaluated using a variety of model systems to -establish their usability in molecular simulations. - -\section{The Ewald Sum} - -The complete accumulation of the electrostatic interactions in a system with -periodic boundary conditions (PBC) requires the consideration of the -effect of all charges within a (cubic) simulation box as well as those -in the periodic replicas, -\begin{equation} -V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime -\left[ \sum_{i=1}^N\sum_{j=1}^N \phi -\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right) -\right], -\label{eq:PBCSum} -\end{equation} -where the sum over $\mathbf{n}$ is a sum over all periodic box -replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the -prime indicates $i = j$ are neglected for $\mathbf{n} = -0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic -particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is -the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and -$j$, and $\phi$ is the solution to Poisson's equation -($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for -charge-charge interactions). In the case of monopole electrostatics, -eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for -non-neutral systems. - -The electrostatic summation problem was originally studied by Ewald -for the case of an infinite crystal.\cite{Ewald21}. The approach he -took was to convert this conditionally convergent sum into two -absolutely convergent summations: a short-ranged real-space summation -and a long-ranged reciprocal-space summation, -\begin{equation} -\begin{split} -V_\textrm{elec} = \frac{1}{2}& -\sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime -\frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)} -{|\mathbf{r}_{ij}+\mathbf{n}|} \\ -&+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} -\exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) -\cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ -&- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 -+ \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3} -\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2, -\end{split} -\label{eq:EwaldSum} -\end{equation} -where $\alpha$ is the damping or convergence parameter with units of -\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to -$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric -constant of the surrounding medium. The final two terms of -eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term -for interacting with a surrounding dielectric.\cite{Allen87} This -dipolar term was neglected in early applications in molecular -simulations,\cite{Brush66,Woodcock71} until it was introduced by de -Leeuw {\it et al.} to address situations where the unit cell has a -dipole moment which is magnified through replication of the periodic -images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the -system is said to be using conducting (or ``tin-foil'') boundary -conditions, $\epsilon_{\rm S} = \infty$. Figure -\ref{fig:ewaldTime} shows how the Ewald sum has been applied over -time. Initially, due to the small system sizes that could be -simulated feasibly, the entire simulation box was replicated to -convergence. In more modern simulations, the systems have grown large -enough that a real-space cutoff could potentially give convergent -behavior. Indeed, it has been observed that with the choice of a -small $\alpha$, the reciprocal-space portion of the Ewald sum can be -rapidly convergent and small relative to the real-space -portion.\cite{Karasawa89,Kolafa92} - -\begin{figure} -\includegraphics[width=\linewidth]{./figures/ewaldProgression.pdf} -\caption{The change in the need for the Ewald sum with -increasing computational power. A:~Initially, only small systems -could be studied, and the Ewald sum replicated the simulation box to -convergence. B:~Now, radial cutoff methods should be able to reach -convergence for the larger systems of charges that are common today.} -\label{fig:ewaldTime} -\end{figure} - -The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The -convergence parameter $(\alpha)$ plays an important role in balancing -the computational cost between the direct and reciprocal-space -portions of the summation. The choice of this value allows one to -select whether the real-space or reciprocal space portion of the -summation is an $\mathscr{O}(N^2)$ calculation (with the other being -$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of -$\alpha$ and thoughtful algorithm development, this cost can be -reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route -taken to reduce the cost of the Ewald summation even further is to set -$\alpha$ such that the real-space interactions decay rapidly, allowing -for a short spherical cutoff. Then the reciprocal space summation is -optimized. These optimizations usually involve utilization of the -fast Fourier transform (FFT),\cite{Hockney81} leading to the -particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) -methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these -methods, the cost of the reciprocal-space portion of the Ewald -summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N -\log N)$. - -These developments and optimizations have made the use of the Ewald -summation routine in simulations with periodic boundary -conditions. However, in certain systems, such as vapor-liquid -interfaces and membranes, the intrinsic three-dimensional periodicity -can prove problematic. The Ewald sum has been reformulated to handle -2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} but these -methods are computationally expensive.\cite{Spohr97,Yeh99} More -recently, there have been several successful efforts toward reducing -the computational cost of 2-D lattice -summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} -bringing them more in line with the cost of the full 3-D summation. - -Several studies have recognized that the inherent periodicity in the -Ewald sum can also have an effect on three-dimensional -systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} -Solvated proteins are essentially kept at high concentration due to -the periodicity of the electrostatic summation method. In these -systems, the more compact folded states of a protein can be -artificially stabilized by the periodic replicas introduced by the -Ewald summation.\cite{Weber00} Thus, care must be taken when -considering the use of the Ewald summation where the assumed -periodicity would introduce spurious effects in the system dynamics. - - -\section{The Wolf and Zahn Methods} - -In a recent paper by Wolf \textit{et al.}, a procedure was outlined -for the accurate accumulation of electrostatic interactions in an -efficient pairwise fashion. This procedure lacks the inherent -periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.} -observed that the electrostatic interaction is effectively -short-ranged in condensed phase systems and that neutralization of the -charge contained within the cutoff radius is crucial for potential -stability. They devised a pairwise summation method that ensures -charge neutrality and gives results similar to those obtained with the -Ewald summation. The resulting shifted Coulomb potential includes -image-charges subtracted out through placement on the cutoff sphere -and a distance-dependent damping function (identical to that seen in -the real-space portion of the Ewald sum) to aid convergence -\begin{equation} -V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}} -- \lim_{r_{ij}\rightarrow R_\textrm{c}} -\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. -\label{eq:WolfPot} -\end{equation} -Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted -potential. However, neutralizing the charge contained within each -cutoff sphere requires the placement of a self-image charge on the -surface of the cutoff sphere. This additional self-term in the total -potential enabled Wolf {\it et al.} to obtain excellent estimates of -Madelung energies for many crystals. - -In order to use their charge-neutralized potential in molecular -dynamics simulations, Wolf \textit{et al.} suggested taking the -derivative of this potential prior to evaluation of the limit. This -procedure gives an expression for the forces, -\begin{equation} -\begin{split} -F_{\textrm{Wolf}}(r_{ij}) = q_i q_j \Biggr\{& -\Biggr[\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r^2_{ij}} -+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{r_{ij}} -\Biggr]\\ -&-\Biggr[ -\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2} -+ \frac{2\alpha}{\pi^{1/2}} -\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}} -\Biggr]\Biggr\}, -\end{split} -\label{eq:WolfForces} -\end{equation} -that incorporates both image charges and damping of the electrostatic -interaction. - -More recently, Zahn \textit{et al.} investigated these potential and -force expressions for use in simulations involving water.\cite{Zahn02} -In their work, they pointed out that the forces and derivative of -the potential are not commensurate. Attempts to use both -eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead -to poor energy conservation. They correctly observed that taking the -limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the -derivatives gives forces for a different potential energy function -than the one shown in eq. (\ref{eq:WolfPot}). - -Zahn \textit{et al.} introduced a modified form of this summation -method as a way to use the technique in Molecular Dynamics -simulations. They proposed a new damped Coulomb potential, -\begin{equation} -\begin{split} -V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\Biggr\{& -\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} \\ -&- \left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} -+ \frac{2\alpha}{\pi^{1/2}} -\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}} -\right]\left(r_{ij}-R_\mathrm{c}\right)\Biggr\}, -\end{split} -\label{eq:ZahnPot} -\end{equation} -and showed that this potential does fairly well at capturing the -structural and dynamic properties of water compared the same -properties obtained using the Ewald sum. - -\section{Simple Forms for Pairwise Electrostatics} - -The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et -al.} are constructed using two different (and separable) computational -tricks: - -\begin{enumerate} -\item shifting through the use of image charges, and -\item damping the electrostatic interaction. -\end{enumerate} -Wolf \textit{et al.} treated the -development of their summation method as a progressive application of -these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded -their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the -post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using -both techniques. It is possible, however, to separate these -tricks and study their effects independently. - -Starting with the original observation that the effective range of the -electrostatic interaction in condensed phases is considerably less -than $r^{-1}$, either the cutoff sphere neutralization or the -distance-dependent damping technique could be used as a foundation for -a new pairwise summation method. Wolf \textit{et al.} made the -observation that charge neutralization within the cutoff sphere plays -a significant role in energy convergence; therefore we will begin our -analysis with the various shifted forms that maintain this charge -neutralization. We can evaluate the methods of Wolf -\textit{et al.} and Zahn \textit{et al.} by considering the standard -shifted potential, -\begin{equation} -V_\textrm{SP}(r) = \begin{cases} -v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > -R_\textrm{c} -\end{cases}, -\label{eq:shiftingPotForm} -\end{equation} -and shifted force, -\begin{equation} -V_\textrm{SF}(r) = \begin{cases} -v(r) - v_\textrm{c} -- \left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) -&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} -\end{cases}, -\label{eq:shiftingForm} -\end{equation} -functions where $v(r)$ is the unshifted form of the potential, and -$v_c$ is $v(R_\textrm{c})$. The Shifted Force ({\sc sf}) form ensures -that both the potential and the forces goes to zero at the cutoff -radius, while the Shifted Potential ({\sc sp}) form only ensures the -potential is smooth at the cutoff radius -($R_\textrm{c}$).\cite{Allen87} - -The forces associated with the shifted potential are simply the forces -of the unshifted potential itself (when inside the cutoff sphere), -\begin{equation} -F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), -\end{equation} -and are zero outside. Inside the cutoff sphere, the forces associated -with the shifted force form can be written, -\begin{equation} -F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d -v(r)}{dr} \right)_{r=R_\textrm{c}}. -\end{equation} - -If the potential, $v(r)$, is taken to be the normal Coulomb potential, -\begin{equation} -v(r) = \frac{q_i q_j}{r}, -\label{eq:Coulomb} -\end{equation} -then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et -al.}'s undamped prescription: -\begin{equation} -V_\textrm{SP}(r) = -q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad -r\leqslant R_\textrm{c}, -\label{eq:SPPot} -\end{equation} -with associated forces, -\begin{equation} -F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) -\quad r\leqslant R_\textrm{c}. -\label{eq:SPForces} -\end{equation} -These forces are identical to the forces of the standard Coulomb -interaction, and cutting these off at $R_c$ was addressed by Wolf -\textit{et al.} as undesirable. They pointed out that the effect of -the image charges is neglected in the forces when this form is -used,\cite{Wolf99} thereby eliminating any benefit from the method in -molecular dynamics. Additionally, there is a discontinuity in the -forces at the cutoff radius which results in energy drift during MD -simulations. - -The shifted force ({\sc sf}) form using the normal Coulomb potential -will give, -\begin{equation} -V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}} -+ \left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] -\quad r\leqslant R_\textrm{c}. -\label{eq:SFPot} -\end{equation} -with associated forces, -\begin{equation} -F_\textrm{SF}(r) = q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) -\quad r\leqslant R_\textrm{c}. -\label{eq:SFForces} -\end{equation} -This formulation has the benefits that there are no discontinuities at -the cutoff radius, while the neutralizing image charges are present in -both the energy and force expressions. It would be simple to add the -self-neutralizing term back when computing the total energy of the -system, thereby maintaining the agreement with the Madelung energies. -A side effect of this treatment is the alteration in the shape of the -potential that comes from the derivative term. Thus, a degree of -clarity about agreement with the empirical potential is lost in order -to gain functionality in dynamics simulations. - -Wolf \textit{et al.} originally discussed the energetics of the -shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was -insufficient for accurate determination of the energy with reasonable -cutoff distances. The calculated Madelung energies fluctuated around -the expected value as the cutoff radius was increased, but the -oscillations converged toward the correct value.\cite{Wolf99} A -damping function was incorporated to accelerate the convergence; and -though alternative forms for the damping function could be -used,\cite{Jones56,Heyes81} the complimentary error function was -chosen to mirror the effective screening used in the Ewald summation. -Incorporating this error function damping into the simple Coulomb -potential, -\begin{equation} -v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, -\label{eq:dampCoulomb} -\end{equation} -the shifted potential (eq. (\ref{eq:SPPot})) becomes -\begin{equation} -V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r} -- \frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) -\quad r\leqslant R_\textrm{c}, -\label{eq:DSPPot} -\end{equation} -with associated forces, -\begin{equation} -F_{\textrm{DSP}}(r) = q_iq_j -\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2} -+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) -\quad r\leqslant R_\textrm{c}. -\label{eq:DSPForces} -\end{equation} -Again, this damped shifted potential suffers from a -force-discontinuity at the cutoff radius, and the image charges play -no role in the forces. To remedy these concerns, one may derive a -{\sc sf} variant by including the derivative term in -eq. (\ref{eq:shiftingForm}), -\begin{equation} -\begin{split} -V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}& -\frac{\mathrm{erfc}\left(\alpha r\right)}{r} -- \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ -&+ \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} -+ \frac{2\alpha}{\pi^{1/2}} -\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}} -\right)\left(r-R_\mathrm{c}\right)\ \Biggr{]} -\quad r\leqslant R_\textrm{c}. -\label{eq:DSFPot} -\end{split} -\end{equation} -The derivative of the above potential will lead to the following forces, -\begin{equation} -\begin{split} -F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}& -\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2} -+ \frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ -&- \left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)} -{R_{\textrm{c}}^2} -+ \frac{2\alpha}{\pi^{1/2}} -\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}} -\right)\Biggr{]} -\quad r\leqslant R_\textrm{c}. -\label{eq:DSFForces} -\end{split} -\end{equation} -If the damping parameter $(\alpha)$ is set to zero, the undamped case, -eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly -recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}). - -This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} -derived by Zahn \textit{et al.}; however, there are two important -differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from -eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb}) -with $r$ replaced by $R_\textrm{c}$. This term is {\it not} present -in the Zahn potential, resulting in a potential discontinuity as -particles cross $R_\textrm{c}$. Second, the sign of the derivative -portion is different. The missing $v_\textrm{c}$ term would not -affect molecular dynamics simulations (although the computed energy -would be expected to have sudden jumps as particle distances crossed -$R_c$). The sign problem is a potential source of errors, however. -In fact, it introduces a discontinuity in the forces at the cutoff, -because the force function is shifted in the wrong direction and -doesn't cross zero at $R_\textrm{c}$. - -Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an -electrostatic summation method in which the potential and forces are -continuous at the cutoff radius and which incorporates the damping -function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of -this paper, we will evaluate exactly how good these methods ({\sc sp}, -{\sc sf}, damping) are at reproducing the correct electrostatic -summation performed by the Ewald sum. - - -\section{Evaluating Pairwise Summation Techniques} - -In classical molecular mechanics simulations, there are two primary -techniques utilized to obtain information about the system of -interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these -techniques utilize pairwise summations of interactions between -particle sites, but they use these summations in different ways. - -In MC, the potential energy difference between configurations dictates -the progression of MC sampling. Going back to the origins of this -method, the acceptance criterion for the canonical ensemble laid out -by Metropolis \textit{et al.} states that a subsequent configuration -is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where -$\xi$ is a random number between 0 and 1.\cite{Metropolis53} -Maintaining the correct $\Delta E$ when using an alternate method for -handling the long-range electrostatics will ensure proper sampling -from the ensemble. - -In MD, the derivative of the potential governs how the system will -progress in time. Consequently, the force and torque vectors on each -body in the system dictate how the system evolves. If the magnitude -and direction of these vectors are similar when using alternate -electrostatic summation techniques, the dynamics in the short term -will be indistinguishable. Because error in MD calculations is -cumulative, one should expect greater deviation at longer times, -although methods which have large differences in the force and torque -vectors will diverge from each other more rapidly. - -\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} - -The pairwise summation techniques (outlined in section -\ref{sec:ESMethods}) were evaluated for use in MC simulations by -studying the energy differences between conformations. We took the -{\sc spme}-computed energy difference between two conformations to be the -correct behavior. An ideal performance by an alternative method would -reproduce these energy differences exactly (even if the absolute -energies calculated by the methods are different). Since none of the -methods provide exact energy differences, we used linear least squares -regressions of energy gap data to evaluate how closely the methods -mimicked the Ewald energy gaps. Unitary results for both the -correlation (slope) and correlation coefficient for these regressions -indicate perfect agreement between the alternative method and {\sc spme}. -Sample correlation plots for two alternate methods are shown in -Fig. \ref{fig:linearFit}. -\begin{figure} -\centering -\includegraphics[width = \linewidth]{./figures/dualLinear.pdf} -\caption{Example least squares regressions of the configuration energy -differences for SPC/E water systems. The upper plot shows a data set -with a poor correlation coefficient ($R^2$), while the lower plot -shows a data set with a good correlation coefficient.} -\label{fig:linearFit} -\end{figure} - -Each of the seven system types (detailed in section \ref{sec:RepSims}) -were represented using 500 independent configurations. Thus, each of -the alternative (non-Ewald) electrostatic summation methods was -evaluated using an accumulated 873,250 configurational energy -differences. - -Results and discussion for the individual analysis of each of the -system types appear in sections \ref{sec:IndividualResults}, while the -cumulative results over all the investigated systems appear below in -sections \ref{sec:EnergyResults}. - -\subsection{Molecular Dynamics and the Force and Torque -Vectors}\label{sec:MDMethods} We evaluated the pairwise methods -(outlined in section \ref{sec:ESMethods}) for use in MD simulations by -comparing the force and torque vectors with those obtained using the -reference Ewald summation ({\sc spme}). Both the magnitude and the -direction of these vectors on each of the bodies in the system were -analyzed. For the magnitude of these vectors, linear least squares -regression analyses were performed as described previously for -comparing $\Delta E$ values. Instead of a single energy difference -between two system configurations, we compared the magnitudes of the -forces (and torques) on each molecule in each configuration. For a -system of 1000 water molecules and 40 ions, there are 1040 force -vectors and 1000 torque vectors. With 500 configurations, this -results in 520,000 force and 500,000 torque vector comparisons. -Additionally, data from seven different system types was aggregated -before the comparison was made. - -The {\it directionality} of the force and torque vectors was -investigated through measurement of the angle ($\theta$) formed -between those computed from the particular method and those from {\sc spme}, -\begin{equation} -\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} -\cdot \hat{F}_\textrm{M}\right), -\end{equation} -where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force -vector computed using method M. Each of these $\theta$ values was -accumulated in a distribution function and weighted by the area on the -unit sphere. Since this distribution is a measure of angular error -between two different electrostatic summation methods, there is no -{\it a priori} reason for the profile to adhere to any specific -shape. Thus, gaussian fits were used to measure the width of the -resulting distributions. The variance ($\sigma^2$) was extracted from -each of these fits and was used to compare distribution widths. -Values of $\sigma^2$ near zero indicate vector directions -indistinguishable from those calculated when using the reference -method ({\sc spme}). - -\subsection{Short-time Dynamics} - -The effects of the alternative electrostatic summation methods on the -short-time dynamics of charged systems were evaluated by considering a -NaCl crystal at a temperature of 1000 K. A subset of the best -performing pairwise methods was used in this comparison. The NaCl -crystal was chosen to avoid possible complications from the treatment -of orientational motion in molecular systems. All systems were -started with the same initial positions and velocities. Simulations -were performed under the microcanonical ensemble, and velocity -autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each -of the trajectories, -\begin{equation} -C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}. -\label{eq:vCorr} -\end{equation} -Velocity autocorrelation functions require detailed short time data, -thus velocity information was saved every 2 fs over 10 ps -trajectories. Because the NaCl crystal is composed of two different -atom types, the average of the two resulting velocity autocorrelation -functions was used for comparisons. - -\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} - -The effects of the same subset of alternative electrostatic methods on -the {\it long-time} dynamics of charged systems were evaluated using -the same model system (NaCl crystals at 1000K). The power spectrum -($I(\omega)$) was obtained via Fourier transform of the velocity -autocorrelation function, \begin{equation} I(\omega) = -\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, -\label{eq:powerSpec} -\end{equation} -where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the -NaCl crystal is composed of two different atom types, the average of -the two resulting power spectra was used for comparisons. Simulations -were performed under the microcanonical ensemble, and velocity -information was saved every 5~fs over 100~ps trajectories. - -\subsection{Representative Simulations}\label{sec:RepSims} -A variety of representative molecular simulations were analyzed to -determine the relative effectiveness of the pairwise summation -techniques in reproducing the energetics and dynamics exhibited by -{\sc spme}. We wanted to span the space of typical molecular -simulations (i.e. from liquids of neutral molecules to ionic -crystals), so the systems studied were: - -\begin{enumerate} -\item liquid water (SPC/E),\cite{Berendsen87} -\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), -\item NaCl crystals, -\item NaCl melts, -\item a low ionic strength solution of NaCl in water (0.11 M), -\item a high ionic strength solution of NaCl in water (1.1 M), and -\item a 6\AA\ radius sphere of Argon in water. -\end{enumerate} - -By utilizing the pairwise techniques (outlined in section -\ref{sec:ESMethods}) in systems composed entirely of neutral groups, -charged particles, and mixtures of the two, we hope to discern under -which conditions it will be possible to use one of the alternative -summation methodologies instead of the Ewald sum. - -For the solid and liquid water configurations, configurations were -taken at regular intervals from high temperature trajectories of 1000 -SPC/E water molecules. Each configuration was equilibrated -independently at a lower temperature (300K for the liquid, 200K for -the crystal). The solid and liquid NaCl systems consisted of 500 -$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for -these systems were selected and equilibrated in the same manner as the -water systems. In order to introduce measurable fluctuations in the -configuration energy differences, the crystalline simulations were -equilibrated at 1000K, near the $T_\textrm{m}$ for NaCl. The liquid -NaCl configurations needed to represent a fully disordered array of -point charges, so the high temperature of 7000K was selected for -equilibration. The ionic solutions were made by solvating 4 (or 40) -ions in a periodic box containing 1000 SPC/E water molecules. Ion and -water positions were then randomly swapped, and the resulting -configurations were again equilibrated individually. Finally, for the -Argon / Water ``charge void'' systems, the identities of all the SPC/E -waters within 6\AA\ of the center of the equilibrated water -configurations were converted to argon. - -These procedures guaranteed us a set of representative configurations -from chemically-relevant systems sampled from appropriate -ensembles. Force field parameters for the ions and Argon were taken -from the force field utilized by {\sc oopse}.\cite{Meineke05} - -\subsection{Comparison of Summation Methods}\label{sec:ESMethods} -We compared the following alternative summation methods with results -from the reference method ({\sc spme}): - -\begin{enumerate} -\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, -and 0.3\AA$^{-1}$, -\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, -and 0.3\AA$^{-1}$, -\item reaction field with an infinite dielectric constant, and -\item an unmodified cutoff. -\end{enumerate} - -Group-based cutoffs with a fifth-order polynomial switching function -were utilized for the reaction field simulations. Additionally, we -investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure -cutoff. The {\sc spme} electrostatics were performed using the {\sc tinker} -implementation of {\sc spme},\cite{Ponder87} while all other calculations -were performed using the {\sc oopse} molecular mechanics -package.\cite{Meineke05} All other portions of the energy calculation -(i.e. Lennard-Jones interactions) were handled in exactly the same -manner across all systems and configurations. - -The alternative methods were also evaluated with three different -cutoff radii (9, 12, and 15\AA). As noted previously, the -convergence parameter ($\alpha$) plays a role in the balance of the -real-space and reciprocal-space portions of the Ewald calculation. -Typical molecular mechanics packages set this to a value dependent on -the cutoff radius and a tolerance (typically less than $1 \times -10^{-4}$ kcal/mol). Smaller tolerances are typically associated with -increasing accuracy at the expense of computational time spent on the -reciprocal-space portion of the summation.\cite{Perram88,Essmann95} -The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used -in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200, -0.3119, and 0.2476\AA$^{-1}$ for cutoff radii of 9, 12, and 15\AA\ -respectively. - -\section{Configuration Energy Difference Results}\label{sec:EnergyResults} -In order to evaluate the performance of the pairwise electrostatic -summation methods for Monte Carlo (MC) simulations, the energy -differences between configurations were compared to the values -obtained when using {\sc spme}. The results for the combined -regression analysis of all of the systems are shown in figure -\ref{fig:delE}. - -\begin{figure} -\centering -\includegraphics[width=4.75in]{./figures/delEplot.pdf} -\caption{Statistical analysis of the quality of configurational energy -differences for a given electrostatic method compared with the -reference Ewald sum. Results with a value equal to 1 (dashed line) -indicate $\Delta E$ values indistinguishable from those obtained using -{\sc spme}. Different values of the cutoff radius are indicated with -different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = -inverted triangles).} -\label{fig:delE} -\end{figure} - -The most striking feature of this plot is how well the Shifted Force -({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy -differences. For the undamped {\sc sf} method, and the -moderately-damped {\sc sp} methods, the results are nearly -indistinguishable from the Ewald results. The other common methods do -significantly less well. - -The unmodified cutoff method is essentially unusable. This is not -surprising since hard cutoffs give large energy fluctuations as atoms -or molecules move in and out of the cutoff -radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to -some degree by using group based cutoffs with a switching -function.\cite{Adams79,Steinbach94,Leach01} However, we do not see -significant improvement using the group-switched cutoff because the -salt and salt solution systems contain non-neutral groups. Section -\ref{sec:IndividualResults} includes results for systems comprised entirely -of neutral groups. - -For the {\sc sp} method, inclusion of electrostatic damping improves -the agreement with Ewald, and using an $\alpha$ of 0.2\AA $^{-1}$ -shows an excellent correlation and quality of fit with the {\sc spme} -results, particularly with a cutoff radius greater than 12 -\AA . Use of a larger damping parameter is more helpful for the -shortest cutoff shown, but it has a detrimental effect on simulations -with larger cutoffs. - -In the {\sc sf} sets, increasing damping results in progressively {\it -worse} correlation with Ewald. Overall, the undamped case is the best -performing set, as the correlation and quality of fits are -consistently superior regardless of the cutoff distance. The undamped -case is also less computationally demanding (because no evaluation of -the complementary error function is required). - -The reaction field results illustrates some of that method's -limitations, primarily that it was developed for use in homogenous -systems; although it does provide results that are an improvement over -those from an unmodified cutoff. - -\section{Magnitude of the Force and Torque Vector Results}\label{sec:FTMagResults} - -Evaluation of pairwise methods for use in Molecular Dynamics -simulations requires consideration of effects on the forces and -torques. Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the -regression results for the force and torque vector magnitudes, -respectively. The data in these figures was generated from an -accumulation of the statistics from all of the system types. - -\begin{figure} -\centering -\includegraphics[width=4.75in]{./figures/frcMagplot.pdf} -\caption{Statistical analysis of the quality of the force vector -magnitudes for a given electrostatic method compared with the -reference Ewald sum. Results with a value equal to 1 (dashed line) -indicate force magnitude values indistinguishable from those obtained -using {\sc spme}. Different values of the cutoff radius are indicated with -different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = -inverted triangles).} -\label{fig:frcMag} -\end{figure} - -Again, it is striking how well the Shifted Potential and Shifted Force -methods are doing at reproducing the {\sc spme} forces. The undamped and -weakly-damped {\sc sf} method gives the best agreement with Ewald. -This is perhaps expected because this method explicitly incorporates a -smooth transition in the forces at the cutoff radius as well as the -neutralizing image charges. - -Figure \ref{fig:frcMag}, for the most part, parallels the results seen -in the previous $\Delta E$ section. The unmodified cutoff results are -poor, but using group based cutoffs and a switching function provides -an improvement much more significant than what was seen with $\Delta -E$. - -With moderate damping and a large enough cutoff radius, the {\sc sp} -method is generating usable forces. Further increases in damping, -while beneficial for simulations with a cutoff radius of 9\AA\ , is -detrimental to simulations with larger cutoff radii. - -The reaction field results are surprisingly good, considering the poor -quality of the fits for the $\Delta E$ results. There is still a -considerable degree of scatter in the data, but the forces correlate -well with the Ewald forces in general. We note that the reaction -field calculations do not include the pure NaCl systems, so these -results are partly biased towards conditions in which the method -performs more favorably. - -\begin{figure} -\centering -\includegraphics[width=4.75in]{./figures/trqMagplot.pdf} -\caption{Statistical analysis of the quality of the torque vector -magnitudes for a given electrostatic method compared with the -reference Ewald sum. Results with a value equal to 1 (dashed line) -indicate torque magnitude values indistinguishable from those obtained -using {\sc spme}. Different values of the cutoff radius are indicated with -different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = -inverted triangles).} -\label{fig:trqMag} -\end{figure} - -Molecular torques were only available from the systems which contained -rigid molecules (i.e. the systems containing water). The data in -fig. \ref{fig:trqMag} is taken from this smaller sampling pool. - -Torques appear to be much more sensitive to charges at a longer -distance. The striking feature in comparing the new electrostatic -methods with {\sc spme} is how much the agreement improves with increasing -cutoff radius. Again, the weakly damped and undamped {\sc sf} method -appears to be reproducing the {\sc spme} torques most accurately. - -Water molecules are dipolar, and the reaction field method reproduces -the effect of the surrounding polarized medium on each of the -molecular bodies. Therefore it is not surprising that reaction field -performs best of all of the methods on molecular torques. - -\section{Directionality of the Force and Torque Vector Results}\label{sec:FTDirResults} - -It is clearly important that a new electrostatic method can reproduce -the magnitudes of the force and torque vectors obtained via the Ewald -sum. However, the {\it directionality} of these vectors will also be -vital in calculating dynamical quantities accurately. Force and -torque directionalities were investigated by measuring the angles -formed between these vectors and the same vectors calculated using -{\sc spme}. The results (Fig. \ref{fig:frcTrqAng}) are compared through the -variance ($\sigma^2$) of the Gaussian fits of the angle error -distributions of the combined set over all system types. - -\begin{figure} -\centering -\includegraphics[width=4.75in]{./figures/frcTrqAngplot.pdf} -\caption{Statistical analysis of the width of the angular distribution -that the force and torque vectors from a given electrostatic method -make with their counterparts obtained using the reference Ewald sum. -Results with a variance ($\sigma^2$) equal to zero (dashed line) -indicate force and torque directions indistinguishable from those -obtained using {\sc spme}. Different values of the cutoff radius are -indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, -and 15\AA\ = inverted triangles).} -\label{fig:frcTrqAng} -\end{figure} - -Both the force and torque $\sigma^2$ results from the analysis of the -total accumulated system data are tabulated in figure -\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc -sp}) method would be essentially unusable for molecular dynamics -unless the damping function is added. The Shifted Force ({\sc sf}) -method, however, is generating force and torque vectors which are -within a few degrees of the Ewald results even with weak (or no) -damping. - -All of the sets (aside from the over-damped case) show the improvement -afforded by choosing a larger cutoff radius. Increasing the cutoff -from 9 to 12\AA\ typically results in a halving of the width of the -distribution, with a similar improvement when going from 12 to 15 -\AA . - -The undamped {\sc sf}, group-based cutoff, and reaction field methods -all do equivalently well at capturing the direction of both the force -and torque vectors. Using the electrostatic damping improves the -angular behavior significantly for the {\sc sp} and moderately for the -{\sc sf} methods. Overdamping is detrimental to both methods. Again -it is important to recognize that the force vectors cover all -particles in all seven systems, while torque vectors are only -available for neutral molecular groups. Damping is more beneficial to -charged bodies, and this observation is investigated further in -section \ref{IndividualResults}. - -Although not discussed previously, group based cutoffs can be applied -to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs -will reintroduce small discontinuities at the cutoff radius, but the -effects of these can be minimized by utilizing a switching function. -Though there are no significant benefits or drawbacks observed in -$\Delta E$ and the force and torque magnitudes when doing this, there -is a measurable improvement in the directionality of the forces and -torques. Table \ref{tab:groupAngle} shows the angular variances -obtained both without (N) and with (Y) group based cutoffs and a -switching function. Note that the $\alpha$ values have units of -\AA$^{-1}$ and the variance values have units of degrees$^2$. The -{\sc sp} (with an $\alpha$ of 0.2\AA$^{-1}$ or smaller) shows much -narrower angular distributions when using group-based cutoffs. The -{\sc sf} method likewise shows improvement in the undamped and lightly -damped cases. - -\begin{table} -\caption{STATISTICAL ANALYSIS OF THE ANGULAR DISTRIBUTIONS ($\sigma^2$) -THAT THE FORCE ({\it upper}) AND TORQUE ({\it lower}) VECTORS FROM A -GIVEN ELECTROSTATIC METHOD MAKE WITH THEIR COUNTERPARTS OBTAINED USING -THE REFERENCE EWALD SUMMATION} - -\footnotesize -\begin{center} -\begin{tabular}{@{} ccrrrrrrrr @{}} \\ -\toprule -\toprule - -& &\multicolumn{4}{c}{Shifted Potential} & \multicolumn{4}{c}{Shifted -Force} \\ -\cmidrule(lr){3-6} \cmidrule(l){7-10} $R_\textrm{c}$ & Groups & -$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$ & -$\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$\\ - -\midrule - -9\AA & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\ - & \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} \\ -12\AA & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\ - & \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} \\ -15\AA & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\ - & \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} \\ - -\midrule - -9\AA & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\ - & \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} \\ -12\AA & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\ - & \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} \\ -15\AA & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\ - & \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} \\ - -\bottomrule -\end{tabular} -\end{center} -\label{tab:groupAngle} -\end{table} - -One additional trend in table \ref{tab:groupAngle} is that the -$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$ -increases, something that is more obvious with group-based cutoffs. -The complimentary error function inserted into the potential weakens -the electrostatic interaction as the value of $\alpha$ is increased. -However, at larger values of $\alpha$, it is possible to overdamp the -electrostatic interaction and to remove it completely. Kast -\textit{et al.} developed a method for choosing appropriate $\alpha$ -values for these types of electrostatic summation methods by fitting -to $g(r)$ data, and their methods indicate optimal values of 0.34, -0.25, and 0.16\AA$^{-1}$ for cutoff values of 9, 12, and 15\AA\ -respectively.\cite{Kast03} These appear to be reasonable choices to -obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on -these findings, choices this high would introduce error in the -molecular torques, particularly for the shorter cutoffs. Based on our -observations, empirical damping up to 0.2\AA$^{-1}$ is beneficial, -but damping may be unnecessary when using the {\sc sf} method. - -\section{Individual System Analysis Results}\label{sec:IndividualResults} - -The combined results of the previous sections show how the pairwise -methods compare to the Ewald summation in the general sense over all -of the system types. It is also useful to consider each of the -studied systems in an individual fashion, so that we can identify -conditions that are particularly difficult for a selected pairwise -method to address. This allows us to further establish the limitations -of these pairwise techniques. Below, the energy difference, force -vector, and torque vector analyses are presented on an individual -system basis. - -\subsection{SPC/E Water Results}\label{sec:WaterResults} - -The first system considered was liquid water at 300K using the SPC/E -model of water.\cite{Berendsen87} The results for the energy gap -comparisons and the force and torque vector magnitude comparisons are -shown in table \ref{tab:spce}. The force and torque vector -directionality results are displayed separately in table -\ref{tab:spceAng}, where the effect of group-based cutoffs and -switching functions on the {\sc sp} and {\sc sf} potentials are also -investigated. In all of the individual results table, the method -abbreviations are as follows: - -\begin{itemize} -\item PC = Pure Cutoff, -\item SP = Shifted Potential, -\item SF = Shifted Force, -\item GSC = Group Switched Cutoff, -\item RF = Reaction Field (where $\varepsilon \approx\infty$), -\item GSSP = Group Switched Shifted Potential, and -\item GSSF = Group Switched Shifted Force. -\end{itemize} - -\begin{table}[htbp] -\centering -\caption{REGRESSION RESULTS OF THE LIQUID WATER SYSTEM FOR THE -$\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it middle}) -AND TORQUE VECTOR MAGNITUDES ({\it lower})} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ -\cmidrule(lr){3-4} -\cmidrule(lr){5-6} -\cmidrule(l){7-8} -Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ -\midrule -PC & & 3.046 & 0.002 & -3.018 & 0.002 & 4.719 & 0.005 \\ -SP & 0.0 & 1.035 & 0.218 & 0.908 & 0.313 & 1.037 & 0.470 \\ - & 0.1 & 1.021 & 0.387 & 0.965 & 0.752 & 1.006 & 0.947 \\ - & 0.2 & 0.997 & 0.962 & 1.001 & 0.994 & 0.994 & 0.996 \\ - & 0.3 & 0.984 & 0.980 & 0.997 & 0.985 & 0.982 & 0.987 \\ -SF & 0.0 & 0.977 & 0.974 & 0.996 & 0.992 & 0.991 & 0.997 \\ - & 0.1 & 0.983 & 0.974 & 1.001 & 0.994 & 0.996 & 0.998 \\ - & 0.2 & 0.992 & 0.989 & 1.001 & 0.995 & 0.994 & 0.996 \\ - & 0.3 & 0.984 & 0.980 & 0.996 & 0.985 & 0.982 & 0.987 \\ -GSC & & 0.918 & 0.862 & 0.852 & 0.756 & 0.801 & 0.700 \\ -RF & & 0.971 & 0.958 & 0.975 & 0.987 & 0.959 & 0.983 \\ -\midrule -PC & & -1.647 & 0.000 & -0.127 & 0.000 & -0.979 & 0.000 \\ -SP & 0.0 & 0.735 & 0.368 & 0.813 & 0.537 & 0.865 & 0.659 \\ - & 0.1 & 0.850 & 0.612 & 0.956 & 0.887 & 0.992 & 0.979 \\ - & 0.2 & 0.996 & 0.989 & 1.000 & 1.000 & 1.000 & 1.000 \\ - & 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\ -SF & 0.0 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 0.999 \\ - & 0.1 & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\ - & 0.2 & 0.999 & 0.998 & 1.000 & 1.000 & 1.000 & 1.000 \\ - & 0.3 & 0.996 & 0.998 & 0.997 & 0.998 & 0.996 & 0.998 \\ -GSC & & 0.998 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\ -RF & & 0.999 & 0.995 & 1.000 & 0.999 & 1.000 & 1.000 \\ -\midrule -PC & & 2.387 & 0.000 & 0.183 & 0.000 & 1.282 & 0.000 \\ -SP & 0.0 & 0.847 & 0.543 & 0.904 & 0.694 & 0.935 & 0.786 \\ - & 0.1 & 0.922 & 0.749 & 0.980 & 0.934 & 0.996 & 0.988 \\ - & 0.2 & 0.987 & 0.985 & 0.989 & 0.992 & 0.990 & 0.993 \\ - & 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\ -SF & 0.0 & 0.978 & 0.990 & 0.988 & 0.997 & 0.993 & 0.999 \\ - & 0.1 & 0.983 & 0.991 & 0.993 & 0.997 & 0.997 & 0.999 \\ - & 0.2 & 0.986 & 0.989 & 0.989 & 0.992 & 0.990 & 0.993 \\ - & 0.3 & 0.965 & 0.973 & 0.967 & 0.975 & 0.967 & 0.976 \\ -GSC & & 0.995 & 0.981 & 0.999 & 0.991 & 1.001 & 0.994 \\ -RF & & 0.993 & 0.989 & 0.998 & 0.996 & 1.000 & 0.999 \\ -\bottomrule -\end{tabular} -\label{tab:spce} -\end{table} - -\begin{table}[htbp] -\centering -\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR -DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE LIQUID WATER -SYSTEM} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ -\cmidrule(lr){3-5} -\cmidrule(l){6-8} -Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ -\midrule -PC & & 783.759 & 481.353 & 332.677 & 248.674 & 144.382 & 98.535 \\ -SP & 0.0 & 659.440 & 380.699 & 250.002 & 235.151 & 134.661 & 88.135 \\ - & 0.1 & 293.849 & 67.772 & 11.609 & 105.090 & 23.813 & 4.369 \\ - & 0.2 & 5.975 & 0.136 & 0.094 & 5.553 & 1.784 & 1.536 \\ - & 0.3 & 0.725 & 0.707 & 0.693 & 7.293 & 6.933 & 6.748 \\ -SF & 0.0 & 2.238 & 0.713 & 0.292 & 3.290 & 1.090 & 0.416 \\ - & 0.1 & 2.238 & 0.524 & 0.115 & 3.184 & 0.945 & 0.326 \\ - & 0.2 & 0.374 & 0.102 & 0.094 & 2.598 & 1.755 & 1.537 \\ - & 0.3 & 0.721 & 0.707 & 0.693 & 7.322 & 6.933 & 6.748 \\ -GSC & & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\ -RF & & 2.091 & 0.403 & 0.113 & 3.583 & 1.071 & 0.399 \\ -\midrule -GSSP & 0.0 & 2.431 & 0.614 & 0.274 & 5.135 & 2.133 & 1.339 \\ - & 0.1 & 1.879 & 0.291 & 0.057 & 3.983 & 1.117 & 0.370 \\ - & 0.2 & 0.443 & 0.103 & 0.093 & 2.821 & 1.794 & 1.532 \\ - & 0.3 & 0.728 & 0.694 & 0.692 & 7.387 & 6.942 & 6.748 \\ -GSSF & 0.0 & 1.298 & 0.270 & 0.083 & 3.098 & 0.992 & 0.375 \\ - & 0.1 & 1.296 & 0.210 & 0.044 & 3.055 & 0.922 & 0.330 \\ - & 0.2 & 0.433 & 0.104 & 0.093 & 2.895 & 1.797 & 1.532 \\ - & 0.3 & 0.728 & 0.694 & 0.692 & 7.410 & 6.942 & 6.748 \\ -\bottomrule -\end{tabular} -\label{tab:spceAng} -\end{table} - -The water results parallel the combined results seen in sections -\ref{sec:EnergyResults} through \ref{sec:FTDirResults}. There is good -agreement with {\sc spme} in both energetic and dynamic behavior when -using the {\sc sf} method with and without damping. The {\sc sp} -method does well with an $\alpha$ around 0.2\AA$^{-1}$, particularly -with cutoff radii greater than 12\AA. Overdamping the electrostatics -reduces the agreement between both these methods and {\sc spme}. - -The pure cutoff ({\sc pc}) method performs poorly, again mirroring the -observations from the combined results. In contrast to these results, however, the use of a switching function and group -based cutoffs greatly improves the results for these neutral water -molecules. The group switched cutoff ({\sc gsc}) does not mimic the -energetics of {\sc spme} as well as the {\sc sp} (with moderate -damping) and {\sc sf} methods, but the dynamics are quite good. The -switching functions correct discontinuities in the potential and -forces, leading to these improved results. Such improvements with the -use of a switching function have been recognized in previous -studies,\cite{Andrea83,Steinbach94} and this proves to be a useful -tactic for stably incorporating local area electrostatic effects. - -The reaction field ({\sc rf}) method simply extends upon the results -observed in the {\sc gsc} case. Both methods are similar in form -(i.e. neutral groups, switching function), but {\sc rf} incorporates -an added effect from the external dielectric. This similarity -translates into the same good dynamic results and improved energetic -agreement with {\sc spme}. Though this agreement is not to the level -of the moderately damped {\sc sp} and {\sc sf} methods, these results -show how incorporating some implicit properties of the surroundings -(i.e. $\epsilon_\textrm{S}$) can improve the solvent depiction. - -As a final note for the liquid water system, use of group cutoffs and a -switching function leads to noticeable improvements in the {\sc sp} -and {\sc sf} methods, primarily in directionality of the force and -torque vectors (table \ref{tab:spceAng}). The {\sc sp} method shows -significant narrowing of the angle distribution when using little to -no damping and only modest improvement for the recommended conditions -($\alpha = 0.2$\AA$^{-1}$ and $R_\textrm{c}~\geqslant~12$\AA). The -{\sc sf} method shows modest narrowing across all damping and cutoff -ranges of interest. When overdamping these methods, group cutoffs and -the switching function do not improve the force and torque -directionalities. - -\subsection{SPC/E Ice I$_\textrm{c}$ Results}\label{sec:IceResults} - -In addition to the disordered molecular system above, the ordered -molecular system of ice I$_\textrm{c}$ was also considered. Ice -polymorph could have been used to fit this role; however, ice -I$_\textrm{c}$ was chosen because it can form an ideal periodic -lattice with the same number of water molecules used in the disordered -liquid state case. The results for the energy gap comparisons and the -force and torque vector magnitude comparisons are shown in table -\ref{tab:ice}. The force and torque vector directionality results are -displayed separately in table \ref{tab:iceAng}, where the effect of -group-based cutoffs and switching functions on the {\sc sp} and {\sc -sf} potentials are also displayed. - -\begin{table}[htbp] -\centering -\caption{REGRESSION RESULTS OF THE ICE I$_\textrm{c}$ SYSTEM FOR -$\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES ({\it -middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ -\cmidrule(lr){3-4} -\cmidrule(lr){5-6} -\cmidrule(l){7-8} -Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ -\midrule -PC & & 19.897 & 0.047 & -29.214 & 0.048 & -3.771 & 0.001 \\ -SP & 0.0 & -0.014 & 0.000 & 2.135 & 0.347 & 0.457 & 0.045 \\ - & 0.1 & 0.321 & 0.017 & 1.490 & 0.584 & 0.886 & 0.796 \\ - & 0.2 & 0.896 & 0.872 & 1.011 & 0.998 & 0.997 & 0.999 \\ - & 0.3 & 0.983 & 0.997 & 0.992 & 0.997 & 0.991 & 0.997 \\ -SF & 0.0 & 0.943 & 0.979 & 1.048 & 0.978 & 0.995 & 0.999 \\ - & 0.1 & 0.948 & 0.979 & 1.044 & 0.983 & 1.000 & 0.999 \\ - & 0.2 & 0.982 & 0.997 & 0.969 & 0.960 & 0.997 & 0.999 \\ - & 0.3 & 0.985 & 0.997 & 0.961 & 0.961 & 0.991 & 0.997 \\ -GSC & & 0.983 & 0.985 & 0.966 & 0.994 & 1.003 & 0.999 \\ -RF & & 0.924 & 0.944 & 0.990 & 0.996 & 0.991 & 0.998 \\ -\midrule -PC & & -4.375 & 0.000 & 6.781 & 0.000 & -3.369 & 0.000 \\ -SP & 0.0 & 0.515 & 0.164 & 0.856 & 0.426 & 0.743 & 0.478 \\ - & 0.1 & 0.696 & 0.405 & 0.977 & 0.817 & 0.974 & 0.964 \\ - & 0.2 & 0.981 & 0.980 & 1.001 & 1.000 & 1.000 & 1.000 \\ - & 0.3 & 0.996 & 0.998 & 0.997 & 0.999 & 0.997 & 0.999 \\ -SF & 0.0 & 0.991 & 0.995 & 1.003 & 0.998 & 0.999 & 1.000 \\ - & 0.1 & 0.992 & 0.995 & 1.003 & 0.998 & 1.000 & 1.000 \\ - & 0.2 & 0.998 & 0.998 & 0.981 & 0.962 & 1.000 & 1.000 \\ - & 0.3 & 0.996 & 0.998 & 0.976 & 0.957 & 0.997 & 0.999 \\ -GSC & & 0.997 & 0.996 & 0.998 & 0.999 & 1.000 & 1.000 \\ -RF & & 0.988 & 0.989 & 1.000 & 0.999 & 1.000 & 1.000 \\ -\midrule -PC & & -6.367 & 0.000 & -3.552 & 0.000 & -3.447 & 0.000 \\ -SP & 0.0 & 0.643 & 0.409 & 0.833 & 0.607 & 0.961 & 0.805 \\ - & 0.1 & 0.791 & 0.683 & 0.957 & 0.914 & 1.000 & 0.989 \\ - & 0.2 & 0.974 & 0.991 & 0.993 & 0.998 & 0.993 & 0.998 \\ - & 0.3 & 0.976 & 0.992 & 0.977 & 0.992 & 0.977 & 0.992 \\ -SF & 0.0 & 0.979 & 0.997 & 0.992 & 0.999 & 0.994 & 1.000 \\ - & 0.1 & 0.984 & 0.997 & 0.996 & 0.999 & 0.998 & 1.000 \\ - & 0.2 & 0.991 & 0.997 & 0.974 & 0.958 & 0.993 & 0.998 \\ - & 0.3 & 0.977 & 0.992 & 0.956 & 0.948 & 0.977 & 0.992 \\ -GSC & & 0.999 & 0.997 & 0.996 & 0.999 & 1.002 & 1.000 \\ -RF & & 0.994 & 0.997 & 0.997 & 0.999 & 1.000 & 1.000 \\ -\bottomrule -\end{tabular} -\label{tab:ice} -\end{table} - -\begin{table}[htbp] -\centering -\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS -OF THE FORCE AND TORQUE VECTORS IN THE ICE I$_\textrm{c}$ SYSTEM} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque -$\sigma^2$} \\ -\cmidrule(lr){3-5} -\cmidrule(l){6-8} -Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ -\midrule -PC & & 2128.921 & 603.197 & 715.579 & 329.056 & 221.397 & 81.042 \\ -SP & 0.0 & 1429.341 & 470.320 & 447.557 & 301.678 & 197.437 & 73.840 \\ - & 0.1 & 590.008 & 107.510 & 18.883 & 118.201 & 32.472 & 3.599 \\ - & 0.2 & 10.057 & 0.105 & 0.038 & 2.875 & 0.572 & 0.518 \\ - & 0.3 & 0.245 & 0.260 & 0.262 & 2.365 & 2.396 & 2.327 \\ -SF & 0.0 & 1.745 & 1.161 & 0.212 & 1.135 & 0.426 & 0.155 \\ - & 0.1 & 1.721 & 0.868 & 0.082 & 1.118 & 0.358 & 0.118 \\ - & 0.2 & 0.201 & 0.040 & 0.038 & 0.786 & 0.555 & 0.518 \\ - & 0.3 & 0.241 & 0.260 & 0.262 & 2.368 & 2.400 & 2.327 \\ -GSC & & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\ -RF & & 2.887 & 0.217 & 0.107 & 1.006 & 0.281 & 0.085 \\ -\midrule -GSSP & 0.0 & 1.483 & 0.261 & 0.099 & 0.926 & 0.295 & 0.095 \\ - & 0.1 & 1.341 & 0.123 & 0.037 & 0.835 & 0.234 & 0.085 \\ - & 0.2 & 0.558 & 0.040 & 0.037 & 0.823 & 0.557 & 0.519 \\ - & 0.3 & 0.250 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\ -GSSF & 0.0 & 2.124 & 0.132 & 0.069 & 0.919 & 0.263 & 0.099 \\ - & 0.1 & 2.165 & 0.101 & 0.035 & 0.895 & 0.244 & 0.096 \\ - & 0.2 & 0.706 & 0.040 & 0.037 & 0.870 & 0.559 & 0.519 \\ - & 0.3 & 0.251 & 0.251 & 0.259 & 2.387 & 2.395 & 2.328 \\ -\bottomrule -\end{tabular} -\label{tab:iceAng} -\end{table} - -Highly ordered systems are a difficult test for the pairwise methods -in that they lack the implicit periodicity of the Ewald summation. As -expected, the energy gap agreement with {\sc spme} is reduced for the -{\sc sp} and {\sc sf} methods with parameters that were ideal for the -disordered liquid system. Moving to higher $R_\textrm{c}$ helps -improve the agreement, though at an increase in computational cost. -The dynamics of this crystalline system (both in magnitude and -direction) are little affected. Both methods still reproduce the Ewald -behavior with the same parameter recommendations from the previous -section. - -It is also worth noting that {\sc rf} exhibits improved energy gap -results over the liquid water system. One possible explanation is -that the ice I$_\textrm{c}$ crystal is ordered such that the net -dipole moment of the crystal is zero. With $\epsilon_\textrm{S} = -\infty$, the reaction field incorporates this structural organization -by actively enforcing a zeroed dipole moment within each cutoff -sphere. - -\subsection{NaCl Melt Results}\label{sec:SaltMeltResults} - -A high temperature NaCl melt was tested to gauge the accuracy of the -pairwise summation methods in a disordered system of charges. The -results for the energy gap comparisons and the force vector magnitude -comparisons are shown in table \ref{tab:melt}. The force vector -directionality results are displayed separately in table -\ref{tab:meltAng}. - -\begin{table}[htbp] -\centering -\caption{REGRESSION RESULTS OF THE MOLTEN SODIUM CHLORIDE SYSTEM FOR -$\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES ({\it -lower})} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ -\cmidrule(lr){3-4} -\cmidrule(lr){5-6} -\cmidrule(l){7-8} -Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ -\midrule -PC & & -0.008 & 0.000 & -0.049 & 0.005 & -0.136 & 0.020 \\ -SP & 0.0 & 0.928 & 0.996 & 0.931 & 0.998 & 0.950 & 0.999 \\ - & 0.1 & 0.977 & 0.998 & 0.998 & 1.000 & 0.997 & 1.000 \\ - & 0.2 & 0.960 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\ - & 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\ -SF & 0.0 & 0.996 & 1.000 & 0.995 & 1.000 & 0.997 & 1.000 \\ - & 0.1 & 1.021 & 1.000 & 1.024 & 1.000 & 1.007 & 1.000 \\ - & 0.2 & 0.966 & 1.000 & 0.813 & 0.996 & 0.811 & 0.954 \\ - & 0.3 & 0.671 & 0.994 & 0.439 & 0.929 & 0.535 & 0.831 \\ - \midrule -PC & & 1.103 & 0.000 & 0.989 & 0.000 & 0.802 & 0.000 \\ -SP & 0.0 & 0.973 & 0.981 & 0.975 & 0.988 & 0.979 & 0.992 \\ - & 0.1 & 0.987 & 0.992 & 0.993 & 0.998 & 0.997 & 0.999 \\ - & 0.2 & 0.993 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\ - & 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\ -SF & 0.0 & 0.996 & 0.997 & 0.997 & 0.999 & 0.998 & 1.000 \\ - & 0.1 & 1.000 & 0.997 & 1.001 & 0.999 & 1.000 & 1.000 \\ - & 0.2 & 0.994 & 0.996 & 0.985 & 0.988 & 0.986 & 0.981 \\ - & 0.3 & 0.956 & 0.956 & 0.940 & 0.912 & 0.948 & 0.929 \\ -\bottomrule -\end{tabular} -\label{tab:melt} -\end{table} - -\begin{table}[htbp] -\centering -\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS -OF THE FORCE VECTORS IN THE MOLTEN SODIUM CHLORIDE SYSTEM} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{3}{c}{Force $\sigma^2$} \\ -\cmidrule(lr){3-5} -\cmidrule(l){6-8} -Method & $\alpha$ & 9\AA & 12\AA & 15\AA \\ -\midrule -PC & & 13.294 & 8.035 & 5.366 \\ -SP & 0.0 & 13.316 & 8.037 & 5.385 \\ - & 0.1 & 5.705 & 1.391 & 0.360 \\ - & 0.2 & 2.415 & 7.534 & 13.927 \\ - & 0.3 & 23.769 & 67.306 & 57.252 \\ -SF & 0.0 & 1.693 & 0.603 & 0.256 \\ - & 0.1 & 1.687 & 0.653 & 0.272 \\ - & 0.2 & 2.598 & 7.523 & 13.930 \\ - & 0.3 & 23.734 & 67.305 & 57.252 \\ -\bottomrule -\end{tabular} -\label{tab:meltAng} -\end{table} - -The molten NaCl system shows more sensitivity to the electrostatic -damping than the water systems. The most noticeable point is that the -undamped {\sc sf} method does very well at replicating the {\sc spme} -configurational energy differences and forces. Light damping appears -to minimally improve the dynamics, but this comes with a deterioration -of the energy gap results. In contrast, this light damping improves -the {\sc sp} energy gaps and forces. Moderate and heavy electrostatic -damping reduce the agreement with {\sc spme} for both methods. From -these observations, the undamped {\sc sf} method is the best choice -for disordered systems of charges. - -\subsection{NaCl Crystal Results}\label{sec:SaltCrystalResults} - -Similar to the use of ice I$_\textrm{c}$ to investigate the role of -order in molecular systems on the effectiveness of the pairwise -methods, the 1000K NaCl crystal system was used to investigate the -accuracy of the pairwise summation methods in an ordered system of -charged particles. The results for the energy gap comparisons and the -force vector magnitude comparisons are shown in table \ref{tab:salt}. -The force vector directionality results are displayed separately in -table \ref{tab:saltAng}. - -\begin{table}[htbp] -\centering -\caption{REGRESSION RESULTS OF THE CRYSTALLINE SODIUM CHLORIDE -SYSTEM FOR $\Delta E$ VALUES ({\it upper}) AND FORCE VECTOR MAGNITUDES -({\it lower})} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ -\cmidrule(lr){3-4} -\cmidrule(lr){5-6} -\cmidrule(l){7-8} -Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ -\midrule -PC & & -20.241 & 0.228 & -20.248 & 0.229 & -20.239 & 0.228 \\ -SP & 0.0 & 1.039 & 0.733 & 2.037 & 0.565 & 1.225 & 0.743 \\ - & 0.1 & 1.049 & 0.865 & 1.424 & 0.784 & 1.029 & 0.980 \\ - & 0.2 & 0.982 & 0.976 & 0.969 & 0.980 & 0.960 & 0.980 \\ - & 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.945 \\ -SF & 0.0 & 1.041 & 0.967 & 0.994 & 0.989 & 0.957 & 0.993 \\ - & 0.1 & 1.050 & 0.968 & 0.996 & 0.991 & 0.972 & 0.995 \\ - & 0.2 & 0.982 & 0.975 & 0.959 & 0.980 & 0.960 & 0.980 \\ - & 0.3 & 0.873 & 0.944 & 0.872 & 0.945 & 0.872 & 0.944 \\ -\midrule -PC & & 0.795 & 0.000 & 0.792 & 0.000 & 0.793 & 0.000 \\ -SP & 0.0 & 0.916 & 0.829 & 1.086 & 0.791 & 1.010 & 0.936 \\ - & 0.1 & 0.958 & 0.917 & 1.049 & 0.943 & 1.001 & 0.995 \\ - & 0.2 & 0.981 & 0.981 & 0.982 & 0.984 & 0.981 & 0.984 \\ - & 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\ -SF & 0.0 & 1.002 & 0.983 & 0.997 & 0.994 & 0.991 & 0.997 \\ - & 0.1 & 1.003 & 0.984 & 0.996 & 0.995 & 0.993 & 0.997 \\ - & 0.2 & 0.983 & 0.980 & 0.981 & 0.984 & 0.981 & 0.984 \\ - & 0.3 & 0.950 & 0.952 & 0.950 & 0.953 & 0.950 & 0.953 \\ -\bottomrule -\end{tabular} -\label{tab:salt} -\end{table} - -\begin{table}[htbp] -\centering -\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR -DISTRIBUTIONS OF THE FORCE VECTORS IN THE CRYSTALLINE SODIUM CHLORIDE -SYSTEM} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{3}{c}{Force $\sigma^2$} \\ -\cmidrule(lr){3-5} -\cmidrule(l){6-8} -Method & $\alpha$ & 9\AA & 12\AA & 15\AA \\ -\midrule -PC & & 111.945 & 111.824 & 111.866 \\ -SP & 0.0 & 112.414 & 152.215 & 38.087 \\ - & 0.1 & 52.361 & 42.574 & 2.819 \\ - & 0.2 & 10.847 & 9.709 & 9.686 \\ - & 0.3 & 31.128 & 31.104 & 31.029 \\ -SF & 0.0 & 10.025 & 3.555 & 1.648 \\ - & 0.1 & 9.462 & 3.303 & 1.721 \\ - & 0.2 & 11.454 & 9.813 & 9.701 \\ - & 0.3 & 31.120 & 31.105 & 31.029 \\ -\bottomrule -\end{tabular} -\label{tab:saltAng} -\end{table} - -The crystalline NaCl system is the most challenging test case for the -pairwise summation methods, as evidenced by the results in tables -\ref{tab:salt} and \ref{tab:saltAng}. The undamped and weakly damped -{\sc sf} methods seem to be the best choices. These methods match well -with {\sc spme} across the energy gap, force magnitude, and force -directionality tests. The {\sc sp} method struggles in all cases, -with the exception of good dynamics reproduction when using weak -electrostatic damping with a large cutoff radius. - -The moderate electrostatic damping case is not as good as we would -expect given the long-time dynamics results observed for this system -(see section \ref{sec:LongTimeDynamics}). Since the data tabulated in -tables \ref{tab:salt} and \ref{tab:saltAng} are a test of -instantaneous dynamics, this indicates that good long-time dynamics -comes in part at the expense of short-time dynamics. - -\subsection{0.11M NaCl Solution Results} - -In an effort to bridge the charged atomic and neutral molecular -systems, Na$^+$ and Cl$^-$ ion charge defects were incorporated into -the liquid water system. This low ionic strength system consists of 4 -ions in the 1000 SPC/E water solvent ($\approx$0.11 M). The results -for the energy gap comparisons and the force and torque vector -magnitude comparisons are shown in table \ref{tab:solnWeak}. The -force and torque vector directionality results are displayed -separately in table \ref{tab:solnWeakAng}, where the effect of -group-based cutoffs and switching functions on the {\sc sp} and {\sc -sf} potentials are investigated. - -\begin{table}[htbp] -\centering -\caption{REGRESSION RESULTS OF THE WEAK SODIUM CHLORIDE SOLUTION -SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES -({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ -\cmidrule(lr){3-4} -\cmidrule(lr){5-6} -\cmidrule(l){7-8} -Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ -\midrule -PC & & 0.247 & 0.000 & -1.103 & 0.001 & 5.480 & 0.015 \\ -SP & 0.0 & 0.935 & 0.388 & 0.984 & 0.541 & 1.010 & 0.685 \\ - & 0.1 & 0.951 & 0.603 & 0.993 & 0.875 & 1.001 & 0.979 \\ - & 0.2 & 0.969 & 0.968 & 0.996 & 0.997 & 0.994 & 0.997 \\ - & 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\ -SF & 0.0 & 0.963 & 0.971 & 0.989 & 0.996 & 0.991 & 0.998 \\ - & 0.1 & 0.970 & 0.971 & 0.995 & 0.997 & 0.997 & 0.999 \\ - & 0.2 & 0.972 & 0.975 & 0.996 & 0.997 & 0.994 & 0.997 \\ - & 0.3 & 0.955 & 0.966 & 0.984 & 0.992 & 0.978 & 0.991 \\ -GSC & & 0.964 & 0.731 & 0.984 & 0.704 & 1.005 & 0.770 \\ -RF & & 0.968 & 0.605 & 0.974 & 0.541 & 1.014 & 0.614 \\ -\midrule -PC & & 1.354 & 0.000 & -1.190 & 0.000 & -0.314 & 0.000 \\ -SP & 0.0 & 0.720 & 0.338 & 0.808 & 0.523 & 0.860 & 0.643 \\ - & 0.1 & 0.839 & 0.583 & 0.955 & 0.882 & 0.992 & 0.978 \\ - & 0.2 & 0.995 & 0.987 & 0.999 & 1.000 & 0.999 & 1.000 \\ - & 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\ -SF & 0.0 & 0.998 & 0.994 & 1.000 & 0.998 & 1.000 & 0.999 \\ - & 0.1 & 0.997 & 0.994 & 1.000 & 0.999 & 1.000 & 1.000 \\ - & 0.2 & 0.999 & 0.998 & 0.999 & 1.000 & 0.999 & 1.000 \\ - & 0.3 & 0.995 & 0.996 & 0.996 & 0.998 & 0.996 & 0.998 \\ -GSC & & 0.995 & 0.990 & 0.998 & 0.997 & 0.998 & 0.996 \\ -RF & & 0.998 & 0.993 & 0.999 & 0.998 & 0.999 & 0.996 \\ -\midrule -PC & & 2.437 & 0.000 & -1.872 & 0.000 & 2.138 & 0.000 \\ -SP & 0.0 & 0.838 & 0.525 & 0.901 & 0.686 & 0.932 & 0.779 \\ - & 0.1 & 0.914 & 0.733 & 0.979 & 0.932 & 0.995 & 0.987 \\ - & 0.2 & 0.977 & 0.969 & 0.988 & 0.990 & 0.989 & 0.990 \\ - & 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\ -SF & 0.0 & 0.969 & 0.977 & 0.987 & 0.996 & 0.993 & 0.998 \\ - & 0.1 & 0.975 & 0.978 & 0.993 & 0.996 & 0.997 & 0.998 \\ - & 0.2 & 0.976 & 0.973 & 0.988 & 0.990 & 0.989 & 0.990 \\ - & 0.3 & 0.952 & 0.950 & 0.964 & 0.971 & 0.965 & 0.970 \\ -GSC & & 0.980 & 0.959 & 0.990 & 0.983 & 0.992 & 0.989 \\ -RF & & 0.984 & 0.975 & 0.996 & 0.995 & 0.998 & 0.998 \\ -\bottomrule -\end{tabular} -\label{tab:solnWeak} -\end{table} - -\begin{table}[htbp] -\centering -\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR -DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE WEAK SODIUM -CHLORIDE SOLUTION SYSTEM} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ -\cmidrule(lr){3-5} -\cmidrule(l){6-8} -Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ -\midrule -PC & & 882.863 & 510.435 & 344.201 & 277.691 & 154.231 & 100.131 \\ -SP & 0.0 & 732.569 & 405.704 & 257.756 & 261.445 & 142.245 & 91.497 \\ - & 0.1 & 329.031 & 70.746 & 12.014 & 118.496 & 25.218 & 4.711 \\ - & 0.2 & 6.772 & 0.153 & 0.118 & 9.780 & 2.101 & 2.102 \\ - & 0.3 & 0.951 & 0.774 & 0.784 & 12.108 & 7.673 & 7.851 \\ -SF & 0.0 & 2.555 & 0.762 & 0.313 & 6.590 & 1.328 & 0.558 \\ - & 0.1 & 2.561 & 0.560 & 0.123 & 6.464 & 1.162 & 0.457 \\ - & 0.2 & 0.501 & 0.118 & 0.118 & 5.698 & 2.074 & 2.099 \\ - & 0.3 & 0.943 & 0.774 & 0.784 & 12.118 & 7.674 & 7.851 \\ -GSC & & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\ -RF & & 2.415 & 0.452 & 0.130 & 6.915 & 1.423 & 0.507 \\ -\midrule -GSSP & 0.0 & 2.915 & 0.643 & 0.261 & 9.576 & 3.133 & 1.812 \\ - & 0.1 & 2.251 & 0.324 & 0.064 & 7.628 & 1.639 & 0.497 \\ - & 0.2 & 0.590 & 0.118 & 0.116 & 6.080 & 2.096 & 2.103 \\ - & 0.3 & 0.953 & 0.759 & 0.780 & 12.347 & 7.683 & 7.849 \\ -GSSF & 0.0 & 1.541 & 0.301 & 0.096 & 6.407 & 1.316 & 0.496 \\ - & 0.1 & 1.541 & 0.237 & 0.050 & 6.356 & 1.202 & 0.457 \\ - & 0.2 & 0.568 & 0.118 & 0.116 & 6.166 & 2.105 & 2.105 \\ - & 0.3 & 0.954 & 0.759 & 0.780 & 12.337 & 7.684 & 7.849 \\ -\bottomrule -\end{tabular} -\label{tab:solnWeakAng} -\end{table} - -Because this system is a perturbation of the pure liquid water system, -comparisons are best drawn between these two sets. The {\sc sp} and -{\sc sf} methods are not significantly affected by the inclusion of a -few ions. The aspect of cutoff sphere neutralization aids in the -smooth incorporation of these ions; thus, all of the observations -regarding these methods carry over from section -\ref{sec:WaterResults}. The differences between these systems are more -visible for the {\sc rf} method. Though good force agreement is still -maintained, the energy gaps show a significant increase in the scatter -of the data. - -\subsection{1.1M NaCl Solution Results} - -The bridging of the charged atomic and neutral molecular systems was -further developed by considering a high ionic strength system -consisting of 40 ions in the 1000 SPC/E water solvent ($\approx$1.1 -M). The results for the energy gap comparisons and the force and -torque vector magnitude comparisons are shown in table -\ref{tab:solnStr}. The force and torque vector directionality -results are displayed separately in table \ref{tab:solnStrAng}, where -the effect of group-based cutoffs and switching functions on the {\sc -sp} and {\sc sf} potentials are investigated. - -\begin{table}[htbp] -\centering -\caption{REGRESSION RESULTS OF THE STRONG SODIUM CHLORIDE SOLUTION -SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR MAGNITUDES -({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ -\cmidrule(lr){3-4} -\cmidrule(lr){5-6} -\cmidrule(l){7-8} -Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ -\midrule -PC & & -0.081 & 0.000 & 0.945 & 0.001 & 0.073 & 0.000 \\ -SP & 0.0 & 0.978 & 0.469 & 0.996 & 0.672 & 0.975 & 0.668 \\ - & 0.1 & 0.944 & 0.645 & 0.997 & 0.886 & 0.991 & 0.978 \\ - & 0.2 & 0.873 & 0.896 & 0.985 & 0.993 & 0.980 & 0.993 \\ - & 0.3 & 0.831 & 0.860 & 0.960 & 0.979 & 0.955 & 0.977 \\ -SF & 0.0 & 0.858 & 0.905 & 0.985 & 0.970 & 0.990 & 0.998 \\ - & 0.1 & 0.865 & 0.907 & 0.992 & 0.974 & 0.994 & 0.999 \\ - & 0.2 & 0.862 & 0.894 & 0.985 & 0.993 & 0.980 & 0.993 \\ - & 0.3 & 0.831 & 0.859 & 0.960 & 0.979 & 0.955 & 0.977 \\ -GSC & & 1.985 & 0.152 & 0.760 & 0.031 & 1.106 & 0.062 \\ -RF & & 2.414 & 0.116 & 0.813 & 0.017 & 1.434 & 0.047 \\ -\midrule -PC & & -7.028 & 0.000 & -9.364 & 0.000 & 0.925 & 0.865 \\ -SP & 0.0 & 0.701 & 0.319 & 0.909 & 0.773 & 0.861 & 0.665 \\ - & 0.1 & 0.824 & 0.565 & 0.970 & 0.930 & 0.990 & 0.979 \\ - & 0.2 & 0.988 & 0.981 & 0.995 & 0.998 & 0.991 & 0.998 \\ - & 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\ -SF & 0.0 & 0.993 & 0.988 & 0.992 & 0.984 & 0.998 & 0.999 \\ - & 0.1 & 0.993 & 0.989 & 0.993 & 0.986 & 0.998 & 1.000 \\ - & 0.2 & 0.993 & 0.992 & 0.995 & 0.998 & 0.991 & 0.998 \\ - & 0.3 & 0.983 & 0.985 & 0.985 & 0.991 & 0.978 & 0.990 \\ -GSC & & 0.964 & 0.897 & 0.970 & 0.917 & 0.925 & 0.865 \\ -RF & & 0.994 & 0.864 & 0.988 & 0.865 & 0.980 & 0.784 \\ -\midrule -PC & & -2.212 & 0.000 & -0.588 & 0.000 & 0.953 & 0.925 \\ -SP & 0.0 & 0.800 & 0.479 & 0.930 & 0.804 & 0.924 & 0.759 \\ - & 0.1 & 0.883 & 0.694 & 0.976 & 0.942 & 0.993 & 0.986 \\ - & 0.2 & 0.952 & 0.943 & 0.980 & 0.984 & 0.980 & 0.983 \\ - & 0.3 & 0.914 & 0.909 & 0.943 & 0.948 & 0.944 & 0.946 \\ -SF & 0.0 & 0.945 & 0.953 & 0.980 & 0.984 & 0.991 & 0.998 \\ - & 0.1 & 0.951 & 0.954 & 0.987 & 0.986 & 0.995 & 0.998 \\ - & 0.2 & 0.951 & 0.946 & 0.980 & 0.984 & 0.980 & 0.983 \\ - & 0.3 & 0.914 & 0.908 & 0.943 & 0.948 & 0.944 & 0.946 \\ -GSC & & 0.882 & 0.818 & 0.939 & 0.902 & 0.953 & 0.925 \\ -RF & & 0.949 & 0.939 & 0.988 & 0.988 & 0.992 & 0.993 \\ -\bottomrule -\end{tabular} -\label{tab:solnStr} -\end{table} - -\begin{table}[htbp] -\centering -\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR DISTRIBUTIONS -OF THE FORCE AND TORQUE VECTORS IN THE STRONG SODIUM CHLORIDE SOLUTION -SYSTEM} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ -\cmidrule(lr){3-5} -\cmidrule(l){6-8} -Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ -\midrule -PC & & 957.784 & 513.373 & 2.260 & 340.043 & 179.443 & 13.079 \\ -SP & 0.0 & 786.244 & 139.985 & 259.289 & 311.519 & 90.280 & 105.187 \\ - & 0.1 & 354.697 & 38.614 & 12.274 & 144.531 & 23.787 & 5.401 \\ - & 0.2 & 7.674 & 0.363 & 0.215 & 16.655 & 3.601 & 3.634 \\ - & 0.3 & 1.745 & 1.456 & 1.449 & 23.669 & 14.376 & 14.240 \\ -SF & 0.0 & 3.282 & 8.567 & 0.369 & 11.904 & 6.589 & 0.717 \\ - & 0.1 & 3.263 & 7.479 & 0.142 & 11.634 & 5.750 & 0.591 \\ - & 0.2 & 0.686 & 0.324 & 0.215 & 10.809 & 3.580 & 3.635 \\ - & 0.3 & 1.749 & 1.456 & 1.449 & 23.635 & 14.375 & 14.240 \\ -GSC & & 6.181 & 2.904 & 2.263 & 44.349 & 19.442 & 12.873 \\ -RF & & 3.891 & 0.847 & 0.323 & 18.628 & 3.995 & 2.072 \\ -\midrule -GSSP & 0.0 & 6.197 & 2.929 & 2.290 & 44.441 & 19.442 & 12.873 \\ - & 0.1 & 4.688 & 1.064 & 0.260 & 31.208 & 6.967 & 2.303 \\ - & 0.2 & 1.021 & 0.218 & 0.213 & 14.425 & 3.629 & 3.649 \\ - & 0.3 & 1.752 & 1.454 & 1.451 & 23.540 & 14.390 & 14.245 \\ -GSSF & 0.0 & 2.494 & 0.546 & 0.217 & 16.391 & 3.230 & 1.613 \\ - & 0.1 & 2.448 & 0.429 & 0.106 & 16.390 & 2.827 & 1.159 \\ - & 0.2 & 0.899 & 0.214 & 0.213 & 13.542 & 3.583 & 3.645 \\ - & 0.3 & 1.752 & 1.454 & 1.451 & 23.587 & 14.390 & 14.245 \\ -\bottomrule -\end{tabular} -\label{tab:solnStrAng} -\end{table} - -The {\sc rf} method struggles with the jump in ionic strength. The -configuration energy differences degrade to unusable levels while the -forces and torques show a more modest reduction in the agreement with -{\sc spme}. The {\sc rf} method was designed for homogeneous systems, -and this attribute is apparent in these results. - -The {\sc sp} and {\sc sf} methods require larger cutoffs to maintain -their agreement with {\sc spme}. With these results, we still -recommend undamped to moderate damping for the {\sc sf} method and -moderate damping for the {\sc sp} method, both with cutoffs greater -than 12\AA. - -\subsection{6\AA\ Argon Sphere in SPC/E Water Results} - -The final model system studied was a 6\AA\ sphere of Argon solvated -by SPC/E water. This serves as a test case of a specifically sized -electrostatic defect in a disordered molecular system. The results for -the energy gap comparisons and the force and torque vector magnitude -comparisons are shown in table \ref{tab:argon}. The force and torque -vector directionality results are displayed separately in table -\ref{tab:argonAng}, where the effect of group-based cutoffs and -switching functions on the {\sc sp} and {\sc sf} potentials are -investigated. - -\begin{table}[htbp] -\centering -\caption{REGRESSION RESULTS OF THE 6\AA\ ARGON SPHERE IN LIQUID -WATER SYSTEM FOR $\Delta E$ VALUES ({\it upper}), FORCE VECTOR -MAGNITUDES ({\it middle}) AND TORQUE VECTOR MAGNITUDES ({\it lower})} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{2}{c}{9\AA} & \multicolumn{2}{c}{12\AA} & \multicolumn{2}{c}{15\AA}\\ -\cmidrule(lr){3-4} -\cmidrule(lr){5-6} -\cmidrule(l){7-8} -Method & $\alpha$ & slope & $R^2$ & slope & $R^2$ & slope & $R^2$ \\ -\midrule -PC & & 2.320 & 0.008 & -0.650 & 0.001 & 3.848 & 0.029 \\ -SP & 0.0 & 1.053 & 0.711 & 0.977 & 0.820 & 0.974 & 0.882 \\ - & 0.1 & 1.032 & 0.846 & 0.989 & 0.965 & 0.992 & 0.994 \\ - & 0.2 & 0.993 & 0.995 & 0.982 & 0.998 & 0.986 & 0.998 \\ - & 0.3 & 0.968 & 0.995 & 0.954 & 0.992 & 0.961 & 0.994 \\ -SF & 0.0 & 0.982 & 0.996 & 0.992 & 0.999 & 0.993 & 1.000 \\ - & 0.1 & 0.987 & 0.996 & 0.996 & 0.999 & 0.997 & 1.000 \\ - & 0.2 & 0.989 & 0.998 & 0.984 & 0.998 & 0.989 & 0.998 \\ - & 0.3 & 0.971 & 0.995 & 0.957 & 0.992 & 0.965 & 0.994 \\ -GSC & & 1.002 & 0.983 & 0.992 & 0.973 & 0.996 & 0.971 \\ -RF & & 0.998 & 0.995 & 0.999 & 0.998 & 0.998 & 0.998 \\ -\midrule -PC & & -36.559 & 0.002 & -44.917 & 0.004 & -52.945 & 0.006 \\ -SP & 0.0 & 0.890 & 0.786 & 0.927 & 0.867 & 0.949 & 0.909 \\ - & 0.1 & 0.942 & 0.895 & 0.984 & 0.974 & 0.997 & 0.995 \\ - & 0.2 & 0.999 & 0.997 & 1.000 & 1.000 & 1.000 & 1.000 \\ - & 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\ -SF & 0.0 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ - & 0.1 & 1.000 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ - & 0.2 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 & 1.000 \\ - & 0.3 & 1.001 & 0.999 & 1.001 & 1.000 & 1.001 & 1.000 \\ -GSC & & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ -RF & & 0.999 & 0.999 & 1.000 & 1.000 & 1.000 & 1.000 \\ -\midrule -PC & & 1.984 & 0.000 & 0.012 & 0.000 & 1.357 & 0.000 \\ -SP & 0.0 & 0.850 & 0.552 & 0.907 & 0.703 & 0.938 & 0.793 \\ - & 0.1 & 0.924 & 0.755 & 0.980 & 0.936 & 0.995 & 0.988 \\ - & 0.2 & 0.985 & 0.983 & 0.986 & 0.988 & 0.987 & 0.988 \\ - & 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\ -SF & 0.0 & 0.977 & 0.989 & 0.987 & 0.995 & 0.992 & 0.998 \\ - & 0.1 & 0.982 & 0.989 & 0.992 & 0.996 & 0.997 & 0.998 \\ - & 0.2 & 0.984 & 0.987 & 0.986 & 0.987 & 0.987 & 0.988 \\ - & 0.3 & 0.961 & 0.966 & 0.959 & 0.964 & 0.960 & 0.966 \\ -GSC & & 0.995 & 0.981 & 0.999 & 0.990 & 1.000 & 0.993 \\ -RF & & 0.993 & 0.988 & 0.997 & 0.995 & 0.999 & 0.998 \\ -\bottomrule -\end{tabular} -\label{tab:argon} -\end{table} +\input{electrostaticsChapter} -\begin{table}[htbp] -\centering -\caption{VARIANCE RESULTS FROM GAUSSIAN FITS TO ANGULAR -DISTRIBUTIONS OF THE FORCE AND TORQUE VECTORS IN THE 6\AA\ SPHERE OF -ARGON IN LIQUID WATER SYSTEM} - -\footnotesize -\begin{tabular}{@{} ccrrrrrr @{}} -\\ -\toprule -\toprule -& & \multicolumn{3}{c}{Force $\sigma^2$} & \multicolumn{3}{c}{Torque $\sigma^2$} \\ -\cmidrule(lr){3-5} -\cmidrule(l){6-8} -Method & $\alpha$ & 9\AA & 12\AA & 15\AA & 9\AA & 12\AA & 15\AA \\ -\midrule -PC & & 568.025 & 265.993 & 195.099 & 246.626 & 138.600 & 91.654 \\ -SP & 0.0 & 504.578 & 251.694 & 179.932 & 231.568 & 131.444 & 85.119 \\ - & 0.1 & 224.886 & 49.746 & 9.346 & 104.482 & 23.683 & 4.480 \\ - & 0.2 & 4.889 & 0.197 & 0.155 & 6.029 & 2.507 & 2.269 \\ - & 0.3 & 0.817 & 0.833 & 0.812 & 8.286 & 8.436 & 8.135 \\ -SF & 0.0 & 1.924 & 0.675 & 0.304 & 3.658 & 1.448 & 0.600 \\ - & 0.1 & 1.937 & 0.515 & 0.143 & 3.565 & 1.308 & 0.546 \\ - & 0.2 & 0.407 & 0.166 & 0.156 & 3.086 & 2.501 & 2.274 \\ - & 0.3 & 0.815 & 0.833 & 0.812 & 8.330 & 8.437 & 8.135 \\ -GSC & & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\ -RF & & 1.822 & 0.408 & 0.142 & 3.799 & 1.362 & 0.550 \\ -\midrule -GSSP & 0.0 & 2.098 & 0.584 & 0.284 & 5.391 & 2.414 & 1.501 \\ - & 0.1 & 1.652 & 0.309 & 0.087 & 4.197 & 1.401 & 0.590 \\ - & 0.2 & 0.465 & 0.165 & 0.153 & 3.323 & 2.529 & 2.273 \\ - & 0.3 & 0.813 & 0.825 & 0.816 & 8.316 & 8.447 & 8.132 \\ -GSSF & 0.0 & 1.173 & 0.292 & 0.113 & 3.452 & 1.347 & 0.583 \\ - & 0.1 & 1.166 & 0.240 & 0.076 & 3.381 & 1.281 & 0.575 \\ - & 0.2 & 0.459 & 0.165 & 0.153 & 3.430 & 2.542 & 2.273 \\ - & 0.3 & 0.814 & 0.825 & 0.816 & 8.325 & 8.447 & 8.132 \\ -\bottomrule -\end{tabular} -\label{tab:argonAng} -\end{table} - -This system does not appear to show any significant deviations from -the previously observed results. The {\sc sp} and {\sc sf} methods -have aggrements similar to those observed in section -\ref{sec:WaterResults}. The only significant difference is the -improvement in the configuration energy differences for the {\sc rf} -method. This is surprising in that we are introducing an inhomogeneity -to the system; however, this inhomogeneity is charge-neutral and does -not result in charged cutoff spheres. The charge-neutrality of the -cutoff spheres, which the {\sc sp} and {\sc sf} methods explicitly -enforce, seems to play a greater role in the stability of the {\sc rf} -method than the required homogeneity of the environment. - - -\section{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals}\label{sec:ShortTimeDynamics} - -Zahn {\it et al.} investigated the structure and dynamics of water -using eqs. (\ref{eq:ZahnPot}) and -(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated -that a method similar (but not identical with) the damped {\sc sf} -method resulted in properties very similar to those obtained when -using the Ewald summation. The properties they studied (pair -distribution functions, diffusion constants, and velocity and -orientational correlation functions) may not be particularly sensitive -to the long-range and collective behavior that governs the -low-frequency behavior in crystalline systems. Additionally, the -ionic crystals are the worst case scenario for the pairwise methods -because they lack the reciprocal space contribution contained in the -Ewald summation. - -We are using two separate measures to probe the effects of these -alternative electrostatic methods on the dynamics in crystalline -materials. For short- and intermediate-time dynamics, we are -computing the velocity autocorrelation function, and for long-time -and large length-scale collective motions, we are looking at the -low-frequency portion of the power spectrum. - -\begin{figure} -\centering -\includegraphics[width = \linewidth]{./figures/vCorrPlot.pdf} -\caption{Velocity autocorrelation functions of NaCl crystals at -1000K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& -0.2\AA$^{-1}$), and {\sc sp} ($\alpha$ = 0.2\AA$^{-1}$). The inset is -a magnification of the area around the first minimum. The times to -first collision are nearly identical, but differences can be seen in -the peaks and troughs, where the undamped and weakly damped methods -are stiffer than the moderately damped and {\sc spme} methods.} -\label{fig:vCorrPlot} -\end{figure} - -The short-time decay of the velocity autocorrelation function through -the first collision are nearly identical in figure -\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show -how the methods differ. The undamped {\sc sf} method has deeper -troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than -any of the other methods. As the damping parameter ($\alpha$) is -increased, these peaks are smoothed out, and the {\sc sf} method -approaches the {\sc spme} results. With $\alpha$ values of 0.2\AA$^{-1}$, -the {\sc sf} and {\sc sp} functions are nearly identical and track the -{\sc spme} features quite well. This is not surprising because the {\sc sf} -and {\sc sp} potentials become nearly identical with increased -damping. However, this appears to indicate that once damping is -utilized, the details of the form of the potential (and forces) -constructed out of the damped electrostatic interaction are less -important. - -\section{Collective Motion: Power Spectra of NaCl Crystals}\label{sec:LongTimeDynamics} - -To evaluate how the differences between the methods affect the -collective long-time motion, we computed power spectra from long-time -traces of the velocity autocorrelation function. The power spectra for -the best-performing alternative methods are shown in -fig. \ref{fig:methodPS}. Apodization of the correlation functions via -a cubic switching function between 40 and 50 ps was used to reduce the -ringing resulting from data truncation. This procedure had no -noticeable effect on peak location or magnitude. - -\begin{figure} -\centering -\includegraphics[width = \linewidth]{./figures/spectraSquare.pdf} -\caption{Power spectra obtained from the velocity auto-correlation -functions of NaCl crystals at 1000K while using {\sc spme}, {\sc sf} -($\alpha$ = 0, 0.1, \& 0.2\AA$^{-1}$), and {\sc sp} ($\alpha$ = -0.2\AA$^{-1}$). The inset shows the frequency region below 100 -cm$^{-1}$ to highlight where the spectra differ.} -\label{fig:methodPS} -\end{figure} - -While the high frequency regions of the power spectra for the -alternative methods are quantitatively identical with Ewald spectrum, -the low frequency region shows how the summation methods differ. -Considering the low-frequency inset (expanded in the upper frame of -figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the -correlated motions are blue-shifted when using undamped or weakly -damped {\sc sf}. When using moderate damping ($\alpha = 0.2$ -\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical -correlated motion to the Ewald method (which has a convergence -parameter of 0.3119\AA$^{-1}$). This weakening of the electrostatic -interaction with increased damping explains why the long-ranged -correlated motions are at lower frequencies for the moderately damped -methods than for undamped or weakly damped methods. - -To isolate the role of the damping constant, we have computed the -spectra for a single method ({\sc sf}) with a range of damping -constants and compared this with the {\sc spme} spectrum. -Fig. \ref{fig:dampInc} shows more clearly that increasing the -electrostatic damping red-shifts the lowest frequency phonon modes. -However, even without any electrostatic damping, the {\sc sf} method -has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. -Without the {\sc sf} modifications, an undamped (pure cutoff) method -would predict the lowest frequency peak near 325 cm$^{-1}$. {\it -Most} of the collective behavior in the crystal is accurately captured -using the {\sc sf} method. Quantitative agreement with Ewald can be -obtained using moderate damping in addition to the shifting at the -cutoff distance. - -\begin{figure} -\centering -\includegraphics[width = \linewidth]{./figures/increasedDamping.pdf} -\caption{Effect of damping on the two lowest-frequency phonon modes in -the NaCl crystal at 1000K. The undamped shifted force ({\sc sf}) -method is off by less than 10 cm$^{-1}$, and increasing the -electrostatic damping to 0.25\AA$^{-1}$ gives quantitative agreement -with the power spectrum obtained using the Ewald sum. Overdamping can -result in underestimates of frequencies of the long-wavelength -motions.} -\label{fig:dampInc} -\end{figure} - -\section{Synopsis of the Pairwise Method Evaluation}\label{sec:PairwiseSynopsis} - -The above investigation of pairwise electrostatic summation techniques -shows that there are viable and computationally efficient alternatives -to the Ewald summation. These methods are derived from the damped and -cutoff-neutralized Coulombic sum originally proposed by Wolf -\textit{et al.}\cite{Wolf99} In particular, the {\sc sf} -method, reformulated above as eqs. (\ref{eq:DSFPot}) and -(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the -energetic and dynamic characteristics exhibited by simulations -employing lattice summation techniques. The cumulative energy -difference results showed the undamped {\sc sf} and moderately damped -{\sc sp} methods produced results nearly identical to {\sc spme}. -Similarly for the dynamic features, the undamped or moderately damped -{\sc sf} and moderately damped {\sc sp} methods produce force and -torque vector magnitude and directions very similar to the expected -values. These results translate into long-time dynamic behavior -equivalent to that produced in simulations using {\sc spme}. - -As in all purely-pairwise cutoff methods, these methods are expected -to scale approximately {\it linearly} with system size, and they are -easily parallelizable. This should result in substantial reductions -in the computational cost of performing large simulations. - -Aside from the computational cost benefit, these techniques have -applicability in situations where the use of the Ewald sum can prove -problematic. Of greatest interest is their potential use in -interfacial systems, where the unmodified lattice sum techniques -artificially accentuate the periodicity of the system in an -undesirable manner. There have been alterations to the standard Ewald -techniques, via corrections and reformulations, to compensate for -these systems; but the pairwise techniques discussed here require no -modifications, making them natural tools to tackle these problems. -Additionally, this transferability gives them benefits over other -pairwise methods, like reaction field, because estimations of physical -properties (e.g. the dielectric constant) are unnecessary. - -If a researcher is using Monte Carlo simulations of large chemical -systems containing point charges, most structural features will be -accurately captured using the undamped {\sc sf} method or the {\sc sp} -method with an electrostatic damping of 0.2\AA$^{-1}$. These methods -would also be appropriate for molecular dynamics simulations where the -data of interest is either structural or short-time dynamical -quantities. For long-time dynamics and collective motions, the safest -pairwise method we have evaluated is the {\sc sf} method with an -electrostatic damping between 0.2 and 0.25\AA$^{-1}$. - -We are not suggesting that there is any flaw with the Ewald sum; in -fact, it is the standard by which these simple pairwise sums have been -judged. However, these results do suggest that in the typical -simulations performed today, the Ewald summation may no longer be -required to obtain the level of accuracy most researchers have come to -expect. - -\section{An Application: TIP5P-E Water} - - \chapter{\label{chap:water}SIMPLE MODELS FOR WATER} \chapter{\label{chap:ice}PHASE BEHAVIOR OF WATER IN COMPUTER SIMULATIONS}