--- trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/03/15 15:51:44 2623 +++ trunk/electrostaticMethodsPaper/electrostaticMethods.tex 2006/03/15 17:09:09 2624 @@ -98,7 +98,7 @@ for an accurate accumulation of electrostatic interact \subsection{The Wolf and Zahn Methods} In a recent paper by Wolf \textit{et al.}, a procedure was outlined -for an accurate accumulation of electrostatic interactions in an +for the accurate accumulation of electrostatic interactions in an efficient pairwise fashion.\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 @@ -111,14 +111,14 @@ V^{\textrm{Wolf}}(r_{ij})= \frac{q_iq_j \textrm{erfc}( 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_iq_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\}. +V_{\textrm{Wolf}}(r_{ij})= \frac{q_iq_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 enables Wolf {\it et al.} to obtain excellent estimates of +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 @@ -126,7 +126,7 @@ F^{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{-\left[\frac derivative of this potential prior to evaluation of the limit. This procedure gives an expression for the forces, \begin{equation} -F^{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{-\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right]+\left[\frac{\textrm{erfc}\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]\right\}, +F_{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{-\left[\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right]+\left[\frac{\textrm{erfc}\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]\right\}, \label{eq:WolfForces} \end{equation} that incorporates both image charges and damping of the electrostatic @@ -134,12 +134,13 @@ In their work, they pointed out that the method that t 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 method 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 is mathematically invalid. +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.} proposed a modified form of this ``Wolf summation method'' as a way to use this technique in Molecular Dynamics @@ -147,7 +148,7 @@ V^{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\frac{\mathrm \ref{eq:WolfForces}, they proposed a new damped Coulomb potential, \begin{equation} -V^{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\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)\right\}. +V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\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)\right\}. \label{eq:ZahnPot} \end{equation} They showed that this potential does fairly well at capturing the @@ -158,10 +159,10 @@ tricks: \begin{itemize} The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et al.} are constructed using two different (and separable) computational -tricks: \begin{itemize} +tricks: \begin{enumerate} \item shifting through the use of image charges, and \item damping the electrostatic interaction. -\end{itemize} Wolf \textit{et al.} treated the +\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 @@ -181,7 +182,7 @@ v^\textrm{SP}(r) = \begin{cases} \textit{et al.} and Zahn \textit{et al.} by considering the standard shifted potential, \begin{equation} -v^\textrm{SP}(r) = \begin{cases} +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}, @@ -189,8 +190,8 @@ v^\textrm{SF}(r) = \begin{cases} \end{equation} and shifted force, \begin{equation} -v^\textrm{SF}(r) = \begin{cases} -v(r)-v_\textrm{c}-\left(\frac{\textrm{d}v(r)}{\textrm{d}r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) +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} @@ -202,82 +203,107 @@ potential is smooth at the cutoff radius potential is smooth at the cutoff radius ($R_\textrm{c}$).\cite{Allen87} - - - -If the derivative term is taken to be zero, we are left with the shifted Coulomb potential devised by Wolf \textit{et al.},\cite{Wolf99} +The forces associated with the shifted potential are simply the forces +of the unshifted potential itself (when inside the cutoff sphere), \begin{equation} -V^\textrm{SP}(r_{ij}) = q_iq_j\left(\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}. \label{eq:WolfSP} +F_{\textrm{SP}} = \left( \frac{d v(r)}{dr} \right), \end{equation} -The forces associated with this potential are obtained by taking the derivative, resulting in the following, +and are zero outside. Inside the cutoff sphere, the forces associated +with the shifted force form can be written, \begin{equation} -F^\textrm{SP}(r_{ij}) = q_iq_j\left(-\frac{1}{r_{ij}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}. +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:WolfSP} +\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:FWolfSP} \end{equation} -These forces are identical to the forces of the standard electrostatic interaction, and this was addressed by Wolf \textit{et al.} as undesirable. They pointed out that the effect of the image charges is neglected in the forces when they would expect there to be some pressure exerted due to their presence.\cite{Wolf99} As a consequence the forces, though mathematically valid, may not be physically correct due to this missing component. Additionally, there is a discontinuity in the forces. This can be remedied with the use of a switching function to zero the potential and forces smoothly as particles near $R_\textrm{c}$. +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. -If the derivative term in equation \ref{eq:shiftingForm} is evaluated, we obtain an hitherto undiscussed shifted force Coulomb potential, +The shifted force ({\sc sf}) form using the normal Coulomb potential +will give, \begin{equation} -V^\textrm{SF}(r_{ij}) = q_iq_j\left[\frac{1}{r_{ij}}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r_{ij}-R_\textrm{c})\right] \quad r_{ij}\leqslant R_\textrm{c}. +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} -Taking the derivative of this shifted force potential gives the -following forces, +with associated forces, \begin{equation} -F^\textrm{SF}(r_{ij} = q_iq_j\left(-\frac{1}{r_{ij}^2}+\frac{1}{R_\textrm{c}^2}\right) \quad r_{ij}\leqslant R_\textrm{c}. +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} -Using this formulation rather than the simple shifted potential -(Eq. \ref{eq:WolfSP}) means that there are no discontinuities in the -forces in addition to the potential. This form also has the benefit -that the image charges have a force presence, addressing the concerns -about a missing physical component. One side effect of this treatment -is a slight alteration in the shape of the potential that comes about -from the derivative term. Thus, a degree of clarity about the -original formulation of the potential is lost in order to gain -functionality in dynamics simulations. +This formulation has the benefits that there are no discontinuities at +the cutoff distance, 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:WolfSP}), and they found that -it was still insufficient for accurate determination of the energy. -The energy would fluctuate around the expected value with increasing -cutoff radius, but the oscillations appeared to be converging toward -the correct value.\cite{Wolf99} A damping function was incorporated to -accelerate convergence; and though alternative functional forms could -be used,\cite{Jones56,Heyes81} the complimentary error function was -chosen to draw parallels to the Ewald summation. Incorporating -damping into the simple Coulomb potential, +it was still insufficient for accurate determination of the energy +with reasonable cutoff distances. The calculated Madelung energies +fluctuate around the expected value with increasing cutoff radius, but +the oscillations converge toward the correct value.\cite{Wolf99} A +damping function was incorporated to accelerate the convergence; and +though alternative functional forms 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_{ij}) = \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}, +v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, \label{eq:dampCoulomb} \end{equation} -the shifted potential (Eq. \ref{eq:WolfSP}) can be rederived +the shifted potential (Eq. \ref{eq:WolfSP}) can be recovered \textit{via} equation \ref{eq:shiftingForm}, \begin{equation} -V^{\textrm{DSP}}(r_{ij}) = q_iq_j\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r_{ij}\leqslant R_\textrm{c}. +v_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}(\alpha r)}{r}-\frac{\textrm{erfc}(\alpha R_\textrm{c})}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c}. \label{eq:DSPPot} -\end{equation} -The derivative of this Shifted-Potential can be taken to obtain forces -for use in MD, +\end{equation}, +with associated forces, \begin{equation} -F^{\textrm{DSP}}(r_{ij}) = q_iq_j\left(-\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}-\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right) \quad r_{ij}\leqslant R_\textrm{c}. +f_{\textrm{DSP}}(r) = q_iq_j\left(-\frac{\textrm{erfc}(\alpha r)}{r^2}-\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r^2)}}{r}\right) \quad r\leqslant R_\textrm{c}. \label{eq:DSPForces} \end{equation} -Again, this Shifted-Potential suffers from a discontinuity in the -forces, and a lack of an image-charge component in the forces. To -remedy these concerns, a Shifted-Force variant is obtained by -inclusion of the derivative term in equation \ref{eq:shiftingForm} to -give, +Again, this damped shifted potential suffers from a discontinuity and +a lack of the image charges in the forces. To remedy these concerns, +one may derive a Shifted-Force variant by including the derivative +term in equation \ref{eq:shiftingForm}, \begin{equation} \begin{split} -V^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ &\left.+\left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right)\left(r_{ij}-R_\mathrm{c}\right)\ \right] \quad r_{ij}\leqslant R_\textrm{c}. +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.+\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)\ \right] \quad r\leqslant R_\textrm{c}. \label{eq:DSFPot} \end{split} \end{equation} The derivative of the above potential gives the following forces, \begin{equation} \begin{split} -F^\mathrm{DSF}(r_{ij}) = q_iq_j\Biggr{[}&-\left(\frac{\textrm{erfc}(\alpha r_{ij})}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r_{ij}^2)}}{r_{ij}}\right) \\ &\left.+\left(\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right] \quad r_{ij}\leqslant R_\textrm{c}. +f_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&-\left(\frac{\textrm{erfc}(\alpha r)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{(-\alpha^2r^2)}}{r}\right) \\ &\left.+\left(\frac{\textrm{erfc}(\alpha R_{\textrm{c}})}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right] \quad r\leqslant R_\textrm{c}. \label{eq:DSFForces} \end{split} \end{equation} @@ -285,23 +311,36 @@ from equation \ref{eq:shiftingForm} is equal to equati This new Shifted-Force potential is similar to equation \ref{eq:ZahnPot} derived by Zahn \textit{et al.}; however, there are two important differences.\cite{Zahn02} First, the $v_\textrm{c}$ term -from equation \ref{eq:shiftingForm} is equal to equation -\ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This -term is not present in the Zahn potential, resulting in a -discontinuity as particles cross $R_\textrm{c}$. Second, the sign of -the derivative portion is different. The constant $v_\textrm{c}$ term -is not going to have a presence in the forces after performing the -derivative, but the negative sign does effect the derivative. In -fact, it introduces a discontinuity in the forces at the cutoff, -because the force function is shifted in the wrong direction and -doesn't cross zero at $R_\textrm{c}$. Thus, these alterations make -for an electrostatic summation method that is continuous in both the -potential and forces and incorporates the pairwise sum considerations -stressed by Wolf \textit{et al.}\cite{Wolf99} +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 would be 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 that is continuous in both the +potential and forces 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. + +\subsection{Other alternatives} + +Reaction Field blah + +Group-based cutoff blah + + \section{Methods} -\subsection{What Qualities are Important?}\label{sec:Qualities} 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 @@ -310,39 +349,39 @@ the origins of this method, the Canonical ensemble acc In MC, the potential energy difference between two subsequent configurations dictates the progression of MC sampling. Going back to -the origins of this method, the Canonical ensemble acceptance criteria -laid out by Metropolis \textit{et al.} states that a subsequent -configuration is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta -E/kT)$, where $\xi$ is a random number between 0 and -1.\cite{Metropolis53} Maintaining a consistent $\Delta E$ when using -an alternate method for handling the long-range electrostatics ensures -proper sampling within the ensemble. +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 directs how the system will +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 it develops as a whole. If the -magnitude and direction of these vectors are similar when using -alternate electrostatic summation techniques, the dynamics in the near -term will be indistinguishable. Because error in MD calculations is -cumulative, one should expect greater deviation in the long term -trajectories with greater differences in these vectors between -configurations using different long-range electrostatics. +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} -Evaluation of the pairwise summation techniques (outlined in section -\ref{sec:ESMethods}) for use in MC simulations was performed through -study of the energy differences between conformations. Considering -the SPME results to be the correct or desired behavior, ideal -performance of a tested method was taken to be agreement between the -energy differences calculated. Linear least squares regression of the -$\Delta E$ values between configurations using SPME against $\Delta E$ -values using tested methods provides a quantitative comparison of this -agreement. Unitary results for both the correlation and correlation -coefficient for these regressions indicate equivalent energetic -results between the methods. The correlation is the slope of the -plotted data while the correlation coefficient ($R^2$) is a measure of -the of the data scatter around the fitted line and tells about the -quality of the fit (Fig. \ref{fig:linearFit}). +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 +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. Since none of the methods +provide exact energy differences, we used linear least squares +regressions of the $\Delta E$ values between configurations using SPME +against $\Delta E$ values using tested methods provides a quantitative +comparison of this agreement. Unitary results for both the +correlation and correlation coefficient for these regressions indicate +equivalent energetic results between the method under consideration +and electrostatics handled using SPME. Sample correlation plots for +two alternate methods are shown in Fig. \ref{fig:linearFit}. \begin{figure} \centering @@ -351,38 +390,47 @@ Each system type (detailed in section \ref{sec:RepSims \label{fig:linearFit} \end{figure} -Each system type (detailed in section \ref{sec:RepSims}) studied -consisted of 500 independent configurations, each equilibrated from -higher temperature trajectories. Thus, 124,750 $\Delta E$ data points -are used in a regression of a single system type. Results and -discussion for the individual analysis of each of the system types -appear in the supporting information, while the cumulative results -over all the investigated systems appears below in section -\ref{sec:EnergyResults}. +Each system type (detailed in section \ref{sec:RepSims}) was +represented using 500 independent configurations. Additionally, we +used seven different system types, so each of the alternate +(non-Ewald) electrostatic summation methods was evaluated using +873,250 configurational energy differences. -\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} -Evaluation of the pairwise methods (outlined in section -\ref{sec:ESMethods}) for use in MD simulations was performed through -comparison of the force and torque vectors obtained with those from -SPME. Both the magnitude and the direction of these vectors on each -of the bodies in the system were analyzed. For the magnitude of these -vectors, linear least squares regression analysis can be performed as -described previously for comparing $\Delta E$ values. Instead of a -single value between two system configurations, there is a value for -each particle in each configuration. For a system of 1000 water -molecules and 40 ions, there are 1040 force vectors and 1000 torque -vectors. With 500 configurations, this results in 520,000 force and -500,000 torque vector comparisons samples for each system type. +Results and discussion for the individual analysis of each of the +system types appear in the supporting information, while the +cumulative results over all the investigated systems appears below in +section \ref{sec:EnergyResults}. -The force and torque vector directions were investigated through -measurement of the angle ($\theta$) formed between those from the -particular method and those from SPME +\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 (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 SPME, \begin{equation} -\theta_F = \frac{\vec{F}_\textrm{SPME}}{|\vec{F}_\textrm{SPME}|}\cdot\frac{\vec{F}_\textrm{Method}}{|\vec{F}_\textrm{Method}|}. +\theta_f = \cos^{-1} \hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}, \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, weighted by the area on the unit sphere. Non-linear fits -were used to measure the shape of the resulting distributions. +function, weighted by the area on the unit sphere. Non-linear +Gaussian fits were used to measure the width of the resulting +distributions. \begin{figure} \centering @@ -395,14 +443,15 @@ between two different electrostatic summation methods, non-linear fits. The solid line is a Gaussian profile, while the dotted line is a Voigt profile, a convolution of a Gaussian and a Lorentzian. Since this distribution is a measure of angular error -between two different electrostatic summation methods, there is -particular reason for the profile to adhere to a specific shape. -Because of this and the Gaussian profile's more statistically -meaningful properties, Gaussian fits was used to compare all the -tested methods. The variance ($\sigma^2$) was extracted from each of -these fits and was used to compare distribution widths. Values of -$\sigma^2$ near zero indicate vector directions indistinguishable from -those calculated when using SPME. +between two different electrostatic summation methods, there is no +{\it a priori} reason for the profile to adhere to any specific shape. +Gaussian fits was used to compare all the tested methods. The +variance ($\sigma^2$) was extracted from each of these fits and was +used to compare distribution widths. Values of $\sigma^2$ near zero +indicate vector directions indistinguishable from those calculated +when using the reference method (SPME). + +\subsection{Short-time Dynamics} \subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} Evaluation of the long-time dynamics of charged systems was performed @@ -752,7 +801,7 @@ V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\ clearly, we show how damping strength affects a simple real-space electrostatic potential, \begin{equation} -V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r_{ij}})}{r_{ij}}\right]S(r), +V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r), \end{equation} where $S(r)$ is a switching function that smoothly zeroes the potential at the cutoff radius. Figure \ref{fig:dampInc} shows how