| 11 |
|
\usepackage{tabularx} |
| 12 |
|
\usepackage{graphicx} |
| 13 |
|
\usepackage{booktabs} |
| 14 |
+ |
\usepackage{bibentry} |
| 15 |
+ |
\usepackage{mathrsfs} |
| 16 |
|
%\usepackage{berkeley} |
| 17 |
|
\usepackage[ref]{overcite} |
| 18 |
|
\pagestyle{plain} |
| 37 |
|
|
| 38 |
|
\maketitle |
| 39 |
|
%\doublespacing |
| 40 |
< |
|
| 40 |
> |
\nobibliography{} |
| 41 |
|
\begin{abstract} |
| 42 |
+ |
A new method for accumulating electrostatic interactions was derived from the previous efforts described in \bibentry{Wolf99} and \bibentry{Zahn02} as a possible replacement for lattice sum methods in molecular simulations. Comparisons were performed with this and other pairwise electrostatic summation techniques against the smooth particle mesh Ewald (SPME) summation to see how well they reproduce the energetics and dynamics of a variety of simulation types. The newly derived Shifted-Force technique shows a remarkable ability to reproduce the behavior exhibited in simulations using SPME with an $\mathscr{O}(N)$ computational cost, equivalent to merely the real-space portion of the lattice summation. |
| 43 |
|
\end{abstract} |
| 44 |
|
|
| 45 |
|
%\narrowtext |
| 50 |
|
|
| 51 |
|
\section{Introduction} |
| 52 |
|
|
| 53 |
+ |
In molecular simulations, proper accumulation of the electrostatic interactions is considered one of the most essential and computationally demanding tasks. |
| 54 |
+ |
|
| 55 |
+ |
blah blah blah Ewald Sum Important blah blah blah |
| 56 |
+ |
|
| 57 |
|
In a recent paper by Wolf \textit{et al.}, a procedure was outlined for accumulation of electrostatic interactions in a simple pairwise fashion.\cite{Wolf99} They took the observation that the electrostatic interaction is short-ranged in systems of charges and that charge neutralization within the cutoff spheres is crucial for potential stability, and they devised a pairwise summation method that ensures charge neutrality and gives results similar to those obtained using the Ewald summation. The resulting shifted Coulomb potential (Eq. \ref{eq:WolfPot}) 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 hasten energetic convergence |
| 58 |
|
\begin{equation} |
| 59 |
|
V(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\}. |
| 86 |
|
V^\mathrm{SF}\left(r_{ij}\right)=q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\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_{ij}-R_\mathrm{c}\right)\right\}. |
| 87 |
|
\label{eq:SFPot} |
| 88 |
|
\end{equation} |
| 89 |
< |
Equation \ref{eq:SFPot} 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 is simply equation \ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This term is not present in equation \ref{eq:ZahnPot}, resulting in a discontinuity in the potential 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 results in a discontinuity in the forces at the cutoff, since the function is pushed in the opposite direction and doesn't cross to zero at $R_\textrm{c}$. |
| 89 |
> |
Equation \ref{eq:SFPot} 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 is simply equation \ref{eq:dampCoulomb} with $R_\textrm{c}$ supplied for $r_{ij}$. This term is not present in equation \ref{eq:ZahnPot}, resulting in a discontinuity in the potential 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} |
| 90 |
|
|
| 91 |
< |
In this paper, a variety of simulation situations were analyzed to determine the relative effectiveness of the adapted Wolf spherical truncation schemes at reproducing the results obtained using a smooth particle mesh Ewald (SPME) summation technique.\cite{Essmann95} In addition to the Shifted-Potential and Shifted-Force adapted Wolf methods, both reaction field and uncorrected cutoff methods were included for comparison purposes. The general usability of these methods in both Monte Carlo and Molecular Dynamics calculations was assessed through statistical analysis over the combined results from all of the following studied systems: |
| 91 |
> |
It is important to note that shifted force techniques have a drawback in that they alter the shape of the original potential. We thereby lose a degree of clarity about the original formulation of the potential in order to gain functionality in dynamics simulations. An alternative direction would be use the derivatives of the original potential for the forces. This was addressed by Wolf \textit{et al.} as undesirable, because 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. In Monte Carlo simulations, this argument is mute, because forces are not evaluated. We decided to consider both the Shifted-Force technique described above and this Shifted-Potential technique to determine their usability in the evaluation of both energetic and dynamic results in simulations with electrostatics. |
| 92 |
> |
|
| 93 |
> |
A variety of simulation situations were assembled and analyzed to determine the relative effectiveness of the adapted Wolf spherical truncation schemes at reproducing the results obtained using a smooth particle mesh Ewald (SPME) summation technique.\cite{Essmann95} In addition to the Shifted-Potential and Shifted-Force adapted Wolf methods, both reaction field and uncorrected cutoff methods were included for comparison purposes. The general usability of these methods in both Monte Carlo and Molecular Dynamics calculations was assessed through statistical analysis over the combined results from all of the following studied systems: |
| 94 |
|
\begin{enumerate} |
| 95 |
|
\item Liquid Water |
| 96 |
|
\item Crystalline Water (Ice I$_\textrm{c}$) |
| 100 |
|
\item High Ionic Strength Solution of NaCl in Water |
| 101 |
|
\item 6 \AA\ Radius Sphere of Argon in Water |
| 102 |
|
\end{enumerate} |
| 103 |
< |
Additional discussion on the results from the individual systems was also performed to identify limitations of the considered methods in specific systems. |
| 103 |
> |
By studying these methods in systems composed entirely of neutral groups, charged particles, and mixtures of the two, we can either comment on possible system dependence or universal applicability of the techniques. |
| 104 |
|
|
| 105 |
|
\section{Methods} |
| 106 |
|
|
| 107 |
< |
In each of the simulated systems, 500 distinct configurations were generated, and the electrostatic summation methods were compared via sequential application on each of these fixed configurations. The methods compared include SPME, the aforementioned Shifted Potential and Shifted Force methods - both with damping parameters ($\alpha$) of 0, 0.1, 0.2, and 0.3 \AA$^{-1}$, reaction field with an infinite dielectric constant, and an unmodified cutoff. Group-based cutoffs with a fifth-order polynomial switching function were necessary for the reaction field simulations and were utilized in the SP, SF, and pure cutoff methods for comparison to the standard lack of group-based cutoffs with a hard truncation. The SPME calculations were performed using the TINKER implementation of SPME, while all all other method calculations were performed using the OOPSE molecular mechanics package.\cite{Ponder87,Meineke05} |
| 107 |
> |
In each of the simulated systems, 500 distinct configurations were generated, and the electrostatic summation methods were compared via sequential application on each of these fixed configurations. The methods compared include SPME, the aforementioned Shifted Potential and Shifted Force methods - both with damping parameters ($\alpha$) of 0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, light, moderate, and heavy damping respectively), reaction field with an infinite dielectric constant, and an unmodified cutoff. Group-based cutoffs with a fifth-order polynomial switching function were necessary for the reaction field simulations and were utilized in the SP, SF, and pure cutoff methods for comparison to the standard lack of group-based cutoffs with a hard truncation. The SPME calculations were performed using the TINKER implementation of SPME, while all all other method calculations were performed using the OOPSE molecular mechanics package.\cite{Ponder87,Meineke05} |
| 108 |
|
|
| 109 |
|
Generation of the system configurations was dependent on the system type. For the solid and liquid water configurations, configuration snapshots were taken at regular intervals from higher temperature 1000 SPC/E water molecule trajectories and each equilibrated individually. The solid and liquid NaCl systems consisted of 500 Na+ and 500 Cl- ions and were selected and equilibrated in the same fashion as the water systems. For the low and high ionic strength NaCl solutions, 4 and 40 ions were first solvated in a 1000 water molecule boxes respectively. 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 (Fig. \ref{fig:argonSlice}). |
| 110 |
|
|
| 115 |
|
\label{fig:argonSlice} |
| 116 |
|
\end{figure} |
| 117 |
|
|
| 118 |
< |
All of these comparisons were performed with three different cutoff radii (9, 12, and 15 \AA) to investigate the cutoff radius dependence of the various techniques. It should be noted that the damping parameter chosen in SPME, or so called ``Ewald Coefficient", has a significant effect on the energies and forces calculated. Typical molecular mechanics packages default 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 increased accuracy.\cite{Essmann95} The default TINKER tolerance of $1 \times 10^{-8}$ was used in all 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. |
| 118 |
> |
All of these comparisons were performed with three different cutoff radii (9, 12, and 15 \AA) to investigate the cutoff radius dependence of the various techniques. It should be noted that the damping parameter chosen in SPME, or so called ``Ewald Coefficient", has a significant effect on the energies and forces calculated. Typical molecular mechanics packages default 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 increased accuracy.\cite{Essmann95} The default TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used in all 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. |
| 119 |
|
|
| 120 |
|
\section{Results and Discussion} |
| 121 |
|
|
| 151 |
|
\label{fig:frcMag} |
| 152 |
|
\end{figure} |
| 153 |
|
|
| 154 |
< |
The results in figure \ref{fig:frcMag} for the most part parallel those seen in the previous look at the $\Delta E$ results. The unmodified cutoff results are poor, but using group based cutoffs and a switching function provides a improvement much more significant than what was seen with $\Delta E$. Looking at the Shifted Potential sets, the slope and R$^2$ improve with the use of damping to an optimal result of 0.2 \AA $^{-1}$ for the 12 and 15 \AA\ cutoffs. Further increases in damping, while beneficial for simulations with a cutoff radius of 9 \AA\ , is detrimental to simulations with larger cutoff radii. The undamped Shifted Force method gives forces in line with those obtained using SPME, and use of a damping function results in minor improvement. 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 it correlates well in general. |
| 154 |
> |
The results in figure \ref{fig:frcMag} for the most part parallel those seen in the previous look at the $\Delta E$ results. The unmodified cutoff results are poor, but using group based cutoffs and a switching function provides a improvement much more significant than what was seen with $\Delta E$. Looking at the Shifted Potential sets, the slope and $R^2$ improve with the use of damping to an optimal result of 0.2 \AA $^{-1}$ for the 12 and 15 \AA\ cutoffs. Further increases in damping, while beneficial for simulations with a cutoff radius of 9 \AA\ , is detrimental to simulations with larger cutoff radii. The undamped Shifted Force method gives forces in line with those obtained using SPME, and use of a damping function results in minor improvement. 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 it correlates well in general. |
| 155 |
|
|
| 156 |
|
\subsection{Torque Magnitude Comparison} |
| 157 |
|
|
| 219 |
|
\label{tab:groupAngle} |
| 220 |
|
\end{table} |
| 221 |
|
|
| 222 |
< |
Although not discussed previously, group based cutoffs can be applied to both the Shifted Potential and Force methods. Use off a switching function corrects for the discontinuities that arise when atoms of a group exit the cutoff before the group's center of mass. Though there are no significant benefit or drawbacks observed in $\Delta E$ and vector magnitude results when doing this, there is a measurable improvement in the vector angle results. Table \ref{tab:groupAngle} shows the angular variance values obtained using group based cutoffs and a switching function alongside the standard results seen in figure \ref{frcTrqAng} for comparison purposes. The Shifted Potential shows much narrower angular distributions for both the force and torque vectors when using an $\alpha$ of 0.2 \AA$^{-1}$ or less, while Shifted Force shows improvements in the undamped and lightly damped cases. Thus, by calculating the electrostatic interactions in terms of molecular pairs rather than atomic pairs, the direction of the force and torque vectors are determined more accurately. |
| 222 |
> |
Although not discussed previously, group based cutoffs can be applied to both the Shifted Potential and Force methods. Use off a switching function corrects for the discontinuities that arise when atoms of a group exit the cutoff before the group's center of mass. Though there are no significant benefit or drawbacks observed in $\Delta E$ and vector magnitude results when doing this, there is a measurable improvement in the vector angle results. Table \ref{tab:groupAngle} shows the angular variance values obtained using group based cutoffs and a switching function alongside the standard results seen in figure \ref{fig:frcTrqAng} for comparison purposes. The Shifted Potential shows much narrower angular distributions for both the force and torque vectors when using an $\alpha$ of 0.2 \AA$^{-1}$ or less, while Shifted Force shows improvements in the undamped and lightly damped cases. Thus, by calculating the electrostatic interactions in terms of molecular pairs rather than atomic pairs, the direction of the force and torque vectors are determined more accurately. |
| 223 |
|
|
| 224 |
|
One additional trend to recognize in table \ref{tab:groupAngle} is that the $\sigma^2$ values for both Shifted Potential and Shifted Force converge as $\alpha$ increases, something that is easier to see when using group based cutoffs. Looking back on figures \ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this behavior clearly at large $\alpha$ and cutoff values. The reason for this is that the complimentary error function inserted into the potential weakens the electrostatic interaction as $\alpha$ increases. Thus, at larger values of $\alpha$, both the summation method types progress toward non-interacting functions, so care is required in choosing large damping functions lest one generate an undesirable loss in the pair interaction. 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 the above findings, empirical damping up to 0.2 \AA$^{-1}$ proves to be beneficial, but is arguably unnecessary when using the Shifted-Force method. |
| 225 |
|
|
| 232 |
|
C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle. |
| 233 |
|
\label{eq:vCorr} |
| 234 |
|
\end{equation} |
| 235 |
< |
Velocity autocorrelation functions require detailed short time data and long trajectories for good statistics, thus velocity information was saved every 5 fs over 100 ps trajectories. The power spectrum ($F_t$) is obtained via discrete Fourier transform of the autocorrelation function |
| 235 |
> |
Velocity autocorrelation functions require detailed short time data and long trajectories for good statistics, thus velocity information was saved every 5 fs over 100 ps trajectories. The power spectrum ($I(\omega)$) is obtained via discrete Fourier transform of the autocorrelation function |
| 236 |
|
\begin{equation} |
| 237 |
< |
F_t = \sum^{N-1}_{\omega=0}C_v(t)e^{-i\omega t/N}, |
| 237 |
> |
I(\omega) = \sum^{N-1}_{\omega=0}C_v(t)e^{-i\omega t/N}, |
| 238 |
|
\label{eq:powerSpec} |
| 239 |
|
\end{equation} |
| 240 |
|
where $N$ is the number of time samples in $C_v(t)$ and the frequency, $\omega=0,\ 1,\ ...,\ N-1$. The resulting spectra (Fig. \ref{fig:normalModes}) show the normal mode frequencies for the crystal under the simulated conditions. |
| 254 |
|
\begin{figure} |
| 255 |
|
\centering |
| 256 |
|
\includegraphics[width = 3.25in]{./alphaCompare.pdf} |
| 257 |
< |
\caption{} |
| 257 |
> |
\caption{Normal modes for an NaCl crystal at 1000 K when using SPME and a simple damped Coulombic sum with damping coefficients ($\alpha$)ranging from 0.15 to 0.3 \AA$^{-1}$. As $\alpha$ increases, the normal modes are red-shifted towards and eventually beyond the values given by SPME. The larger $\alpha$ values weaken the real-space electrostatics, explaining this shift towards less strongly correlated motions in the crystal.} |
| 258 |
|
\label{fig:dampInc} |
| 259 |
|
\end{figure} |
| 260 |
|
|
| 261 |
|
\section{Conclusions} |
| 262 |
|
|
| 263 |
+ |
This investigation of pairwise electrostatic summation techniques shows that there are viable and more computationally efficient electrostatic summation techniques than the Ewald summation, chiefly methods derived from the damped Coulombic sum originally proposed by Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the Shifted-Force method, reformulated above, 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 Shifted-Force and moderately damped Shifted-Potential methods produced results nearly identical to SPME. Similarly for the dynamic features, the un- to moderately damped Shifted-Force and moderately damped Shifted-Potential 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 SPME. |
| 264 |
+ |
|
| 265 |
+ |
Aside from the computational cost benefit, these techniques have applicability in situations where the use of the Ewald sum can prove problematic. Primary among them is their 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 these pairwise techniques require no modifications, making them natural tools to tackle these problems. Additionally, this transferability gives it benefits over other pairwise methods, like reaction field, because estimations of physical properties, like the dielectric constant, are unnecessary. |
| 266 |
+ |
|
| 267 |
+ |
These results don't deprecate the use of the Ewald summation; in fact, it is the standard to which these simple pairwise sums are judged. However, these results do speak to the necessity of the Ewald summation in all molecular simulations. That a simple pairwise technique can be substituted to gain nearly all the physical effects provided by the full lattice sum makes us question whether the minimal perturbations bestowed through added complexity and increased cost are worth it. |
| 268 |
+ |
|
| 269 |
|
\section{Acknowledgments} |
| 270 |
|
|
| 271 |
|
\newpage |
