| 96 |
|
|
| 97 |
|
\section{Methods} |
| 98 |
|
|
| 99 |
< |
\subsection{What Qualities are Important?} |
| 99 |
> |
\subsection{What Qualities are Important?}\label{sec:Qualities} |
| 100 |
> |
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. |
| 101 |
|
|
| 102 |
< |
\subsection{Monte Carlo and the Energy Gap} |
| 102 |
> |
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. |
| 103 |
|
|
| 104 |
< |
\subsection{Molecular Dynamics and the Force and Torque Vectors} |
| 104 |
> |
In MD, the derivative of the potential directs 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. |
| 105 |
|
|
| 106 |
< |
\subsection{Long-Time and Collective Motion} |
| 106 |
> |
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 107 |
> |
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}). |
| 108 |
|
|
| 109 |
< |
\subsection{Representative Simulations} |
| 110 |
< |
A variety of common and representative simulations 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: |
| 109 |
> |
\begin{figure} |
| 110 |
> |
\centering |
| 111 |
> |
\includegraphics[width=3.25in]{./linearFit.pdf} |
| 112 |
> |
\caption{Example least squares regression of the $\Delta E$ between configurations for the SF method against SPME in the pure water system. } |
| 113 |
> |
\label{fig:linearFit} |
| 114 |
> |
\end{figure} |
| 115 |
> |
|
| 116 |
> |
Each system type (detailed in section \ref{sec:Simulations}) 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}. |
| 117 |
> |
|
| 118 |
> |
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
| 119 |
> |
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. |
| 120 |
> |
|
| 121 |
> |
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. 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. |
| 122 |
> |
|
| 123 |
> |
\begin{figure} |
| 124 |
> |
\centering |
| 125 |
> |
\includegraphics[width=3.25in]{./gaussFit.pdf} |
| 126 |
> |
\caption{Sample fit of the angular distribution of the force vectors over all of the studied systems. Gaussian fits were used to obtain values for the variance in force and torque vectors used in the following figure.} |
| 127 |
> |
\label{fig:gaussian} |
| 128 |
> |
\end{figure} |
| 129 |
> |
|
| 130 |
> |
Figure \ref{fig:gaussian} shows an example distribution with applied 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. |
| 131 |
> |
|
| 132 |
> |
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 133 |
> |
Evaluation of the long-time dynamics of charged systems was performed by considering the NaCl crystal system while using a subset of the best performing pairwise methods. The NaCl crystal was chosen to avoid possible complications involving the propagation techniques of orientational motion in molecular systems. To enhance the atomic motion, these crystals were equilibrated at 1000 K, near the experimental $T_m$ for NaCl. Simulations were performed under the microcanonical ensemble, and velocity autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each of the trajectories, |
| 134 |
> |
\begin{equation} |
| 135 |
> |
C_v(t) = \langle v_i(0)\cdot v_i(t)\rangle. |
| 136 |
> |
\label{eq:vCorr} |
| 137 |
> |
\end{equation} |
| 138 |
> |
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 Fourier transform of the autocorrelation function |
| 139 |
> |
\begin{equation} |
| 140 |
> |
I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 141 |
> |
\label{eq:powerSpec} |
| 142 |
> |
\end{equation} |
| 143 |
> |
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. |
| 144 |
> |
|
| 145 |
> |
\subsection{Representative Simulations}\label{sec:RepSims} |
| 146 |
> |
A variety of common and representative simulations were analyzed to determine the relative effectiveness of the pairwise summation techniques in reproducing the energetics and dynamics exhibited by SPME. The studied systems were as follows: |
| 147 |
|
\begin{enumerate} |
| 148 |
|
\item Liquid Water |
| 149 |
|
\item Crystalline Water (Ice I$_\textrm{c}$) |
| 153 |
|
\item High Ionic Strength Solution of NaCl in Water |
| 154 |
|
\item 6 \AA\ Radius Sphere of Argon in Water |
| 155 |
|
\end{enumerate} |
| 156 |
< |
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. |
| 156 |
> |
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 can comment on possible system dependence and/or universal applicability of the techniques. |
| 157 |
|
|
| 120 |
– |
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,\cite{Ponder87} while all all other method calculations were performed using the OOPSE molecular mechanics package.\cite{Meineke05} |
| 121 |
– |
|
| 158 |
|
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}). |
| 159 |
|
|
| 160 |
|
\begin{figure} |
| 164 |
|
\label{fig:argonSlice} |
| 165 |
|
\end{figure} |
| 166 |
|
|
| 167 |
< |
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. |
| 167 |
> |
\subsection{Electrostatic Summation Methods}\label{sec:ESMethods} |
| 168 |
> |
Electrostatic summation method comparisons were performed using SPME, the Shifted-Potential and Shifted-Force methods - both with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak, moderate, and strong 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,\cite{Ponder87} while all other method calculations were performed using the OOPSE molecular mechanics package.\cite{Meineke05} |
| 169 |
|
|
| 170 |
+ |
These methods were additionally evaluated with three different cutoff radii (9, 12, and 15 \AA) to investigate possible cutoff radius dependence. 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 in the real-space portion of the summation.\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. |
| 171 |
+ |
|
| 172 |
|
\section{Results and Discussion} |
| 173 |
|
|
| 174 |
< |
\subsection{Configuration Energy Differences} |
| 175 |
< |
In order to evaluate the performance of the adapted Wolf Shifted Potential and Shifted Force electrostatic summation methods for Monte Carlo simulations, the energy differences between configurations need to be compared to the results using SPME. Considering the SPME results to be the correct or desired behavior, ideal performance of a tested method is 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 gives an idea of the quality of the fit (Fig. \ref{fig:linearFit}). |
| 174 |
> |
\subsection{Configuration Energy Differences}\label{sec:EnergyResults} |
| 175 |
> |
In order to evaluate the performance of the pairwise electrostatic summation methods for Monte Carlo simulations, the energy differences between configurations were compared to the values obtained when using SPME. The results for the subsequent regression analysis are shown in figure \ref{fig:delE}. |
| 176 |
|
|
| 177 |
|
\begin{figure} |
| 178 |
|
\centering |
| 140 |
– |
\includegraphics[width=3.25in]{./linearFit.pdf} |
| 141 |
– |
\caption{Example least squares regression of the $\Delta E$ between configurations for the SF method against SPME in the pure water system. } |
| 142 |
– |
\label{fig:linearFit} |
| 143 |
– |
\end{figure} |
| 144 |
– |
|
| 145 |
– |
With 500 independent configurations, 124,750 $\Delta E$ data points are used in a regression of a single system. Results and discussion for the individual analysis of each of the system types appear in the supporting information. To probe the applicability of each method in the general case, all the different system types were included in a single regression. The results for this regression are shown in figure \ref{fig:delE}. |
| 146 |
– |
|
| 147 |
– |
\begin{figure} |
| 148 |
– |
\centering |
| 179 |
|
\includegraphics[width=3.25in]{./delEplot.pdf} |
| 180 |
|
\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 SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 181 |
|
\label{fig:delE} |
| 182 |
|
\end{figure} |
| 183 |
|
|
| 184 |
< |
In figure \ref{fig:delE}, it is apparent that it is unreasonable to expect realistic results using an unmodified cutoff. This is not all that surprising since this results in large energy fluctuations as atoms move in and out of the cutoff radius. These fluctuations can be alleviated to some degree by using group based cutoffs with a switching function. The Group Switch Cutoff row doesn't show a significant improvement in this plot because the salt and salt solution systems contain non-neutral groups, see the accompanying supporting information for a comparison where all groups are neutral. Correcting the resulting charged cutoff sphere is one of the purposes of the shifted potential proposed by Wolf \textit{et al.}, and this correction indeed improves the results as seen in the Shifted Potental rows. While the undamped case of this method is a significant improvement over the pure cutoff, it still doesn't correlate that well with SPME. Inclusion of potential damping improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows an excellent correlation and quality of fit with the 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 Shifted Force sets, increasing damping results in progressively poorer correlation. Overall, the undamped case is the best performing set, as the correlation and quality of fits are consistently superior regardless of the cutoff distance. This result is beneficial in that the undamped case is less computationally prohibitive do to the lack of complimentary error function calculation when performing the electrostatic pair interaction. 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. |
| 184 |
> |
In this figure, it is apparent that it is unreasonable to expect realistic results using an unmodified cutoff. This is not all that surprising since this results in large energy fluctuations as atoms move in and out of the cutoff radius. These fluctuations can be alleviated to some degree by using group based cutoffs with a switching function. The Group Switch Cutoff row doesn't show a significant improvement in this plot because the salt and salt solution systems contain non-neutral groups, see the accompanying supporting information for a comparison where all groups are neutral. |
| 185 |
|
|
| 186 |
+ |
Correcting the resulting charged cutoff sphere is one of the purposes of the shifted potential proposed by Wolf \textit{et al.}, and this correction indeed improves the results as seen in the Shifted Potental rows. While the undamped case of this method is a significant improvement over the pure cutoff, it still doesn't correlate that well with SPME. Inclusion of potential damping improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows an excellent correlation and quality of fit with the 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 Shifted Force sets, increasing damping results in progressively poorer correlation. Overall, the undamped case is the best performing set, as the correlation and quality of fits are consistently superior regardless of the cutoff distance. This result is beneficial in that the undamped case is less computationally prohibitive do to the lack of complimentary error function calculation when performing the electrostatic pair interaction. 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. |
| 187 |
+ |
|
| 188 |
|
\subsection{Magnitudes of the Force and Torque Vectors} |
| 189 |
|
|
| 190 |
|
While studying the energy differences provides insight into how comparable these methods are energetically, if we want to use these methods in Molecular Dynamics simulations, we also need to consider their effect on forces and torques. Both the magnitude and the direction of the force and torque vectors of each of the bodies in the system can be compared to those observed while using SPME. Analysis of the magnitude of these vectors can be performed in the manner described previously for comparing $\Delta E$ values, only 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 excess of 500,000 data samples for each system type. Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the force and torque vector magnitude results for the accumulated analysis over all the system types. |
