| 140 |
|
\end{split} |
| 141 |
|
\label{eq:EwaldSum} |
| 142 |
|
\end{equation} |
| 143 |
< |
where $\alpha$ is a damping parameter, or separation constant, with |
| 144 |
< |
units of \AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are |
| 145 |
< |
equal to $2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the |
| 146 |
< |
dielectric constant of the surrounding medium. The final two terms of |
| 143 |
> |
where $\alpha$ is the damping or convergence parameter with units of |
| 144 |
> |
\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to |
| 145 |
> |
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
| 146 |
> |
constant of the surrounding medium. The final two terms of |
| 147 |
|
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
| 148 |
|
for interacting with a surrounding dielectric.\cite{Allen87} This |
| 149 |
|
dipolar term was neglected in early applications in molecular |
| 159 |
|
convergence. In more modern simulations, the simulation boxes have |
| 160 |
|
grown large enough that a real-space cutoff could potentially give |
| 161 |
|
convergent behavior. Indeed, it has often been observed that the |
| 162 |
< |
reciprocal-space portion of the Ewald sum can be vanishingly |
| 163 |
< |
small compared to the real-space portion.\cite{XXX} |
| 162 |
> |
reciprocal-space portion of the Ewald sum can be small and rapidly |
| 163 |
> |
convergent compared to the real-space portion with the choice of small |
| 164 |
> |
$\alpha$.\cite{Karasawa89,Kolafa92} |
| 165 |
|
|
| 166 |
|
\begin{figure} |
| 167 |
|
\centering |
| 177 |
|
\end{figure} |
| 178 |
|
|
| 179 |
|
The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The |
| 180 |
< |
separation constant $(\alpha)$ plays an important role in balancing |
| 180 |
> |
convergence parameter $(\alpha)$ plays an important role in balancing |
| 181 |
|
the computational cost between the direct and reciprocal-space |
| 182 |
|
portions of the summation. The choice of this value allows one to |
| 183 |
|
select whether the real-space or reciprocal space portion of the |
| 493 |
|
techniques utilize pairwise summations of interactions between |
| 494 |
|
particle sites, but they use these summations in different ways. |
| 495 |
|
|
| 496 |
< |
In MC, the potential energy difference between two subsequent |
| 497 |
< |
configurations dictates the progression of MC sampling. Going back to |
| 498 |
< |
the origins of this method, the acceptance criterion for the canonical |
| 499 |
< |
ensemble laid out by Metropolis \textit{et al.} states that a |
| 500 |
< |
subsequent configuration is accepted if $\Delta E < 0$ or if $\xi < |
| 501 |
< |
\exp(-\Delta E/kT)$, where $\xi$ is a random number between 0 and |
| 502 |
< |
1.\cite{Metropolis53} Maintaining the correct $\Delta E$ when using an |
| 503 |
< |
alternate method for handling the long-range electrostatics will |
| 504 |
< |
ensure proper sampling from the ensemble. |
| 496 |
> |
In MC, the potential energy difference between configurations dictates |
| 497 |
> |
the progression of MC sampling. Going back to the origins of this |
| 498 |
> |
method, the acceptance criterion for the canonical ensemble laid out |
| 499 |
> |
by Metropolis \textit{et al.} states that a subsequent configuration |
| 500 |
> |
is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where |
| 501 |
> |
$\xi$ is a random number between 0 and 1.\cite{Metropolis53} |
| 502 |
> |
Maintaining the correct $\Delta E$ when using an alternate method for |
| 503 |
> |
handling the long-range electrostatics will ensure proper sampling |
| 504 |
> |
from the ensemble. |
| 505 |
|
|
| 506 |
|
In MD, the derivative of the potential governs how the system will |
| 507 |
|
progress in time. Consequently, the force and torque vectors on each |
| 514 |
|
vectors will diverge from each other more rapidly. |
| 515 |
|
|
| 516 |
|
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
| 517 |
+ |
|
| 518 |
|
The pairwise summation techniques (outlined in section |
| 519 |
|
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 520 |
|
studying the energy differences between conformations. We took the |
| 521 |
|
SPME-computed energy difference between two conformations to be the |
| 522 |
|
correct behavior. An ideal performance by an alternative method would |
| 523 |
< |
reproduce these energy differences exactly. Since none of the methods |
| 524 |
< |
provide exact energy differences, we used linear least squares |
| 525 |
< |
regressions of the $\Delta E$ values between configurations using SPME |
| 526 |
< |
against $\Delta E$ values using tested methods provides a quantitative |
| 527 |
< |
comparison of this agreement. Unitary results for both the |
| 528 |
< |
correlation and correlation coefficient for these regressions indicate |
| 529 |
< |
equivalent energetic results between the method under consideration |
| 530 |
< |
and electrostatics handled using SPME. Sample correlation plots for |
| 531 |
< |
two alternate methods are shown in Fig. \ref{fig:linearFit}. |
| 523 |
> |
reproduce these energy differences exactly (even if the absolute |
| 524 |
> |
energies calculated by the methods are different). Since none of the |
| 525 |
> |
methods provide exact energy differences, we used linear least squares |
| 526 |
> |
regressions of energy gap data to evaluate how closely the methods |
| 527 |
> |
mimicked the Ewald energy gaps. Unitary results for both the |
| 528 |
> |
correlation (slope) and correlation coefficient for these regressions |
| 529 |
> |
indicate perfect agreement between the alternative method and SPME. |
| 530 |
> |
Sample correlation plots for two alternate methods are shown in |
| 531 |
> |
Fig. \ref{fig:linearFit}. |
| 532 |
|
|
| 533 |
|
\begin{figure} |
| 534 |
|
\centering |
| 535 |
|
\includegraphics[width = \linewidth]{./dualLinear.pdf} |
| 536 |
< |
\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.} |
| 537 |
< |
\label{fig:linearFit} |
| 536 |
> |
\caption{Example least squares regressions of the configuration energy |
| 537 |
> |
differences for SPC/E water systems. The upper plot shows a data set |
| 538 |
> |
with a poor correlation coefficient ($R^2$), while the lower plot |
| 539 |
> |
shows a data set with a good correlation coefficient.} |
| 540 |
> |
\label{fig:linearFit} |
| 541 |
|
\end{figure} |
| 542 |
|
|
| 543 |
|
Each system type (detailed in section \ref{sec:RepSims}) was |
| 544 |
|
represented using 500 independent configurations. Additionally, we |
| 545 |
< |
used seven different system types, so each of the alternate |
| 545 |
> |
used seven different system types, so each of the alternative |
| 546 |
|
(non-Ewald) electrostatic summation methods was evaluated using |
| 547 |
|
873,250 configurational energy differences. |
| 548 |
|
|
| 572 |
|
investigated through measurement of the angle ($\theta$) formed |
| 573 |
|
between those computed from the particular method and those from SPME, |
| 574 |
|
\begin{equation} |
| 575 |
< |
\theta_f = \cos^{-1} \left(\hat{f}_\textrm{SPME} \cdot \hat{f}_\textrm{Method}\right), |
| 575 |
> |
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right), |
| 576 |
|
\end{equation} |
| 577 |
< |
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the |
| 578 |
< |
force vector computed using method $M$. |
| 579 |
< |
|
| 580 |
< |
Each of these $\theta$ values was accumulated in a distribution |
| 581 |
< |
function, weighted by the area on the unit sphere. Non-linear |
| 582 |
< |
Gaussian fits were used to measure the width of the resulting |
| 583 |
< |
distributions. |
| 584 |
< |
|
| 585 |
< |
\begin{figure} |
| 586 |
< |
\centering |
| 587 |
< |
\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 588 |
< |
\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.} |
| 589 |
< |
\label{fig:gaussian} |
| 590 |
< |
\end{figure} |
| 591 |
< |
|
| 592 |
< |
Figure \ref{fig:gaussian} shows an example distribution with applied |
| 593 |
< |
non-linear fits. The solid line is a Gaussian profile, while the |
| 594 |
< |
dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 595 |
< |
Lorentzian. Since this distribution is a measure of angular error |
| 596 |
< |
between two different electrostatic summation methods, there is no |
| 597 |
< |
{\it a priori} reason for the profile to adhere to any specific shape. |
| 598 |
< |
Gaussian fits was used to compare all the tested methods. The |
| 599 |
< |
variance ($\sigma^2$) was extracted from each of these fits and was |
| 600 |
< |
used to compare distribution widths. Values of $\sigma^2$ near zero |
| 601 |
< |
indicate vector directions indistinguishable from those calculated |
| 602 |
< |
when using the reference method (SPME). |
| 577 |
> |
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force |
| 578 |
> |
vector computed using method M. Each of these $\theta$ values was |
| 579 |
> |
accumulated in a distribution function and weighted by the area on the |
| 580 |
> |
unit sphere. Non-linear Gaussian fits were used to measure the width |
| 581 |
> |
of the resulting distributions. |
| 582 |
> |
% |
| 583 |
> |
%\begin{figure} |
| 584 |
> |
%\centering |
| 585 |
> |
%\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 586 |
> |
%\caption{Sample fit of the angular distribution of the force vectors |
| 587 |
> |
%accumulated using all of the studied systems. Gaussian fits were used |
| 588 |
> |
%to obtain values for the variance in force and torque vectors.} |
| 589 |
> |
%\label{fig:gaussian} |
| 590 |
> |
%\end{figure} |
| 591 |
> |
% |
| 592 |
> |
%Figure \ref{fig:gaussian} shows an example distribution with applied |
| 593 |
> |
%non-linear fits. The solid line is a Gaussian profile, while the |
| 594 |
> |
%dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 595 |
> |
%Lorentzian. |
| 596 |
> |
%Since this distribution is a measure of angular error between two |
| 597 |
> |
%different electrostatic summation methods, there is no {\it a priori} |
| 598 |
> |
%reason for the profile to adhere to any specific shape. |
| 599 |
> |
%Gaussian fits was used to compare all the tested methods. |
| 600 |
> |
The variance ($\sigma^2$) was extracted from each of these fits and |
| 601 |
> |
was used to compare distribution widths. Values of $\sigma^2$ near |
| 602 |
> |
zero indicate vector directions indistinguishable from those |
| 603 |
> |
calculated when using the reference method (SPME). |
| 604 |
|
|
| 605 |
|
\subsection{Short-time Dynamics} |
| 606 |
< |
Evaluation of the short-time dynamics of charged systems was performed |
| 607 |
< |
by considering the 1000 K NaCl crystal system while using a subset of the |
| 608 |
< |
best performing pairwise methods. The NaCl crystal was chosen to |
| 609 |
< |
avoid possible complications involving the propagation techniques of |
| 610 |
< |
orientational motion in molecular systems. All systems were started |
| 611 |
< |
with the same initial positions and velocities. Simulations were |
| 612 |
< |
performed under the microcanonical ensemble, and velocity |
| 606 |
> |
|
| 607 |
> |
The effects of the alternative electrostatic summation methods on the |
| 608 |
> |
short-time dynamics of charged systems were evaluated by considering a |
| 609 |
> |
NaCl crystal at a temperature of 1000 K. A subset of the best |
| 610 |
> |
performing pairwise methods was used in this comparison. The NaCl |
| 611 |
> |
crystal was chosen to avoid possible complications from the treatment |
| 612 |
> |
of orientational motion in molecular systems. All systems were |
| 613 |
> |
started with the same initial positions and velocities. Simulations |
| 614 |
> |
were performed under the microcanonical ensemble, and velocity |
| 615 |
|
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
| 616 |
|
of the trajectories, |
| 617 |
|
\begin{equation} |
| 625 |
|
functions was used for comparisons. |
| 626 |
|
|
| 627 |
|
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
| 628 |
< |
Evaluation of the long-time dynamics of charged systems was performed |
| 629 |
< |
by considering the NaCl crystal system, again while using a subset of |
| 630 |
< |
the best performing pairwise methods. To enhance the atomic motion, |
| 631 |
< |
these crystals were equilibrated at 1000 K, near the experimental |
| 632 |
< |
$T_m$ for NaCl. Simulations were performed under the microcanonical |
| 633 |
< |
ensemble, and velocity information was saved every 5 fs over 100 ps |
| 634 |
< |
trajectories. The power spectrum ($I(\omega)$) was obtained via |
| 627 |
< |
Fourier transform of the velocity autocorrelation function |
| 628 |
< |
\begin{equation} |
| 629 |
< |
I(\omega) = \frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 628 |
> |
|
| 629 |
> |
The effects of the same subset of alternative electrostatic methods on |
| 630 |
> |
the {\it long-time} dynamics of charged systems were evaluated using |
| 631 |
> |
the same model system (NaCl crystals at 1000K). The power spectrum |
| 632 |
> |
($I(\omega)$) was obtained via Fourier transform of the velocity |
| 633 |
> |
autocorrelation function, \begin{equation} I(\omega) = |
| 634 |
> |
\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 635 |
|
\label{eq:powerSpec} |
| 636 |
|
\end{equation} |
| 637 |
|
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the |
| 638 |
|
NaCl crystal is composed of two different atom types, the average of |
| 639 |
< |
the two resulting power spectra was used for comparisons. |
| 639 |
> |
the two resulting power spectra was used for comparisons. Simulations |
| 640 |
> |
were performed under the microcanonical ensemble, and velocity |
| 641 |
> |
information was saved every 5 fs over 100 ps trajectories. |
| 642 |
|
|
| 643 |
|
\subsection{Representative Simulations}\label{sec:RepSims} |
| 644 |
< |
A variety of common and representative simulations were analyzed to |
| 645 |
< |
determine the relative effectiveness of the pairwise summation |
| 646 |
< |
techniques in reproducing the energetics and dynamics exhibited by |
| 647 |
< |
SPME. The studied systems were as follows: |
| 644 |
> |
A variety of representative simulations were analyzed to determine the |
| 645 |
> |
relative effectiveness of the pairwise summation techniques in |
| 646 |
> |
reproducing the energetics and dynamics exhibited by SPME. We wanted |
| 647 |
> |
to span the space of modern simulations (i.e. from liquids of neutral |
| 648 |
> |
molecules to ionic crystals), so the systems studied were: |
| 649 |
|
\begin{enumerate} |
| 650 |
< |
\item Liquid Water |
| 651 |
< |
\item Crystalline Water (Ice I$_\textrm{c}$) |
| 652 |
< |
\item NaCl Crystal |
| 653 |
< |
\item NaCl Melt |
| 654 |
< |
\item Low Ionic Strength Solution of NaCl in Water |
| 655 |
< |
\item High Ionic Strength Solution of NaCl in Water |
| 656 |
< |
\item 6 \AA\ Radius Sphere of Argon in Water |
| 650 |
> |
\item liquid water (SPC/E),\cite{Berendsen87} |
| 651 |
> |
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), |
| 652 |
> |
\item NaCl crystals, |
| 653 |
> |
\item NaCl melts, |
| 654 |
> |
\item a low ionic strength solution of NaCl in water (0.11 M), |
| 655 |
> |
\item a high ionic strength solution of NaCl in water (1.1 M), and |
| 656 |
> |
\item a 6 \AA\ radius sphere of Argon in water. |
| 657 |
|
\end{enumerate} |
| 658 |
|
By utilizing the pairwise techniques (outlined in section |
| 659 |
|
\ref{sec:ESMethods}) in systems composed entirely of neutral groups, |
| 660 |
< |
charged particles, and mixtures of the two, we can comment on possible |
| 661 |
< |
system dependence and/or universal applicability of the techniques. |
| 660 |
> |
charged particles, and mixtures of the two, we hope to discern under |
| 661 |
> |
which conditions it will be possible to use one of the alternative |
| 662 |
> |
summation methodologies instead of the Ewald sum. |
| 663 |
|
|
| 664 |
< |
Generation of the system configurations was dependent on the system |
| 665 |
< |
type. For the solid and liquid water configurations, configuration |
| 666 |
< |
snapshots were taken at regular intervals from higher temperature 1000 |
| 667 |
< |
SPC/E water molecule trajectories and each equilibrated |
| 668 |
< |
individually.\cite{Berendsen87} The solid and liquid NaCl systems |
| 669 |
< |
consisted of 500 Na+ and 500 Cl- ions and were selected and |
| 670 |
< |
equilibrated in the same fashion as the water systems. For the low |
| 671 |
< |
and high ionic strength NaCl solutions, 4 and 40 ions were first |
| 672 |
< |
solvated in a 1000 water molecule boxes respectively. Ion and water |
| 673 |
< |
positions were then randomly swapped, and the resulting configurations |
| 674 |
< |
were again equilibrated individually. Finally, for the Argon/Water |
| 675 |
< |
"charge void" systems, the identities of all the SPC/E waters within 6 |
| 676 |
< |
\AA\ of the center of the equilibrated water configurations were |
| 677 |
< |
converted to argon (Fig. \ref{fig:argonSlice}). |
| 678 |
< |
|
| 679 |
< |
\begin{figure} |
| 671 |
< |
\centering |
| 672 |
< |
\includegraphics[width = \linewidth]{./slice.pdf} |
| 673 |
< |
\caption{A slice from the center of a water box used in a charge void simulation. The darkened region represents the boundary sphere within which the water molecules were converted to argon atoms.} |
| 674 |
< |
\label{fig:argonSlice} |
| 675 |
< |
\end{figure} |
| 664 |
> |
For the solid and liquid water configurations, configurations were |
| 665 |
> |
taken at regular intervals from high temperature trajectories of 1000 |
| 666 |
> |
SPC/E water molecules. Each configuration was equilibrated |
| 667 |
> |
independently at a lower temperature (300~K for the liquid, 200~K for |
| 668 |
> |
the crystal). The solid and liquid NaCl systems consisted of 500 |
| 669 |
> |
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
| 670 |
> |
these systems were selected and equilibrated in the same manner as the |
| 671 |
> |
water systems. The equilibrated temperatures were 1000~K for the NaCl |
| 672 |
> |
crystal and 7000~K for the liquid. The ionic solutions were made by |
| 673 |
> |
solvating 4 (or 40) ions in a periodic box containing 1000 SPC/E water |
| 674 |
> |
molecules. Ion and water positions were then randomly swapped, and |
| 675 |
> |
the resulting configurations were again equilibrated individually. |
| 676 |
> |
Finally, for the Argon / Water ``charge void'' systems, the identities |
| 677 |
> |
of all the SPC/E waters within 6 \AA\ of the center of the |
| 678 |
> |
equilibrated water configurations were converted to argon. |
| 679 |
> |
%(Fig. \ref{fig:argonSlice}). |
| 680 |
|
|
| 681 |
< |
\subsection{Electrostatic Summation Methods}\label{sec:ESMethods} |
| 682 |
< |
Electrostatic summation method comparisons were performed using SPME, |
| 683 |
< |
the {\sc sp} and {\sc sf} methods - both with damping |
| 684 |
< |
parameters ($\alpha$) of 0.0, 0.1, 0.2, and 0.3 \AA$^{-1}$ (no, weak, |
| 681 |
< |
moderate, and strong damping respectively), reaction field with an |
| 682 |
< |
infinite dielectric constant, and an unmodified cutoff. Group-based |
| 683 |
< |
cutoffs with a fifth-order polynomial switching function were |
| 684 |
< |
necessary for the reaction field simulations and were utilized in the |
| 685 |
< |
SP, SF, and pure cutoff methods for comparison to the standard lack of |
| 686 |
< |
group-based cutoffs with a hard truncation. The SPME calculations |
| 687 |
< |
were performed using the TINKER implementation of SPME,\cite{Ponder87} |
| 688 |
< |
while all other method calculations were performed using the OOPSE |
| 689 |
< |
molecular mechanics package.\cite{Meineke05} |
| 681 |
> |
These procedures guaranteed us a set of representative configurations |
| 682 |
> |
from chemically-relevant systems sampled from an appropriate |
| 683 |
> |
ensemble. Force field parameters for the ions and Argon were taken |
| 684 |
> |
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 685 |
|
|
| 686 |
< |
These methods were additionally evaluated with three different cutoff |
| 687 |
< |
radii (9, 12, and 15 \AA) to investigate possible cutoff radius |
| 688 |
< |
dependence. It should be noted that the damping parameter chosen in |
| 689 |
< |
SPME, or so called ``Ewald Coefficient", has a significant effect on |
| 690 |
< |
the energies and forces calculated. Typical molecular mechanics |
| 691 |
< |
packages default this to a value dependent on the cutoff radius and a |
| 692 |
< |
tolerance (typically less than $1 \times 10^{-4}$ kcal/mol). Smaller |
| 693 |
< |
tolerances are typically associated with increased accuracy, but this |
| 699 |
< |
usually means more time spent calculating the reciprocal-space portion |
| 700 |
< |
of the summation.\cite{Perram88,Essmann95} The default TINKER |
| 701 |
< |
tolerance of $1 \times 10^{-8}$ kcal/mol was used in all SPME |
| 702 |
< |
calculations, resulting in Ewald Coefficients of 0.4200, 0.3119, and |
| 703 |
< |
0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 686 |
> |
%\begin{figure} |
| 687 |
> |
%\centering |
| 688 |
> |
%\includegraphics[width = \linewidth]{./slice.pdf} |
| 689 |
> |
%\caption{A slice from the center of a water box used in a charge void |
| 690 |
> |
%simulation. The darkened region represents the boundary sphere within |
| 691 |
> |
%which the water molecules were converted to argon atoms.} |
| 692 |
> |
%\label{fig:argonSlice} |
| 693 |
> |
%\end{figure} |
| 694 |
|
|
| 695 |
+ |
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 696 |
+ |
We compared the following alternative summation methods with results |
| 697 |
+ |
from the reference method (SPME): |
| 698 |
+ |
\begin{itemize} |
| 699 |
+ |
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 700 |
+ |
and 0.3 \AA$^{-1}$, |
| 701 |
+ |
\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 702 |
+ |
and 0.3 \AA$^{-1}$, |
| 703 |
+ |
\item reaction field with an infinite dielectric constant, and |
| 704 |
+ |
\item an unmodified cutoff. |
| 705 |
+ |
\end{itemize} |
| 706 |
+ |
Group-based cutoffs with a fifth-order polynomial switching function |
| 707 |
+ |
were utilized for the reaction field simulations. Additionally, we |
| 708 |
+ |
investigated the use of these cutoffs with the SP, SF, and pure |
| 709 |
+ |
cutoff. The SPME electrostatics were performed using the TINKER |
| 710 |
+ |
implementation of SPME,\cite{Ponder87} while all other method |
| 711 |
+ |
calculations were performed using the OOPSE molecular mechanics |
| 712 |
+ |
package.\cite{Meineke05} All other portions of the energy calculation |
| 713 |
+ |
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 714 |
+ |
manner across all systems and configurations. |
| 715 |
+ |
|
| 716 |
+ |
The althernative methods were also evaluated with three different |
| 717 |
+ |
cutoff radii (9, 12, and 15 \AA). As noted previously, the |
| 718 |
+ |
convergence parameter ($\alpha$) plays a role in the balance of the |
| 719 |
+ |
real-space and reciprocal-space portions of the Ewald calculation. |
| 720 |
+ |
Typical molecular mechanics packages set this to a value dependent on |
| 721 |
+ |
the cutoff radius and a tolerance (typically less than $1 \times |
| 722 |
+ |
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 723 |
+ |
increased accuracy at the expense of increased time spent calculating |
| 724 |
+ |
the reciprocal-space portion of the |
| 725 |
+ |
summation.\cite{Perram88,Essmann95} The default TINKER tolerance of $1 |
| 726 |
+ |
\times 10^{-8}$ kcal/mol was used in all SPME calculations, resulting |
| 727 |
+ |
in Ewald Coefficients of 0.4200, 0.3119, and 0.2476 \AA$^{-1}$ for |
| 728 |
+ |
cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 729 |
+ |
|
| 730 |
|
\section{Results and Discussion} |
| 731 |
|
|
| 732 |
|
\subsection{Configuration Energy Differences}\label{sec:EnergyResults} |
| 739 |
|
\begin{figure} |
| 740 |
|
\centering |
| 741 |
|
\includegraphics[width=5.5in]{./delEplot.pdf} |
| 742 |
< |
\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).} |
| 742 |
> |
\caption{Statistical analysis of the quality of configurational energy |
| 743 |
> |
differences for a given electrostatic method compared with the |
| 744 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 745 |
> |
indicate $\Delta E$ values indistinguishable from those obtained using |
| 746 |
> |
SPME. Different values of the cutoff radius are indicated with |
| 747 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 748 |
> |
inverted triangles).} |
| 749 |
|
\label{fig:delE} |
| 750 |
|
\end{figure} |
| 751 |
|
|
| 752 |
< |
In this figure, it is apparent that it is unreasonable to expect |
| 753 |
< |
realistic results using an unmodified cutoff. This is not all that |
| 754 |
< |
surprising since this results in large energy fluctuations as atoms or |
| 755 |
< |
molecules move in and out of the cutoff radius.\cite{Rahman71,Adams79} |
| 756 |
< |
These fluctuations can be alleviated to some degree by using group |
| 757 |
< |
based cutoffs with a switching |
| 758 |
< |
function.\cite{Adams79,Steinbach94,Leach01} The Group Switch Cutoff |
| 759 |
< |
row doesn't show a significant improvement in this plot because the |
| 760 |
< |
salt and salt solution systems contain non-neutral groups, see the |
| 761 |
< |
accompanying supporting information for a comparison where all groups |
| 762 |
< |
are neutral. |
| 763 |
< |
|
| 764 |
< |
Correcting the resulting charged cutoff sphere is one of the purposes |
| 765 |
< |
of the damped Coulomb summation proposed by Wolf \textit{et |
| 766 |
< |
al.},\cite{Wolf99} and this correction indeed improves the results as |
| 767 |
< |
seen in the {\sc sp} rows. While the undamped case of this |
| 768 |
< |
method is a significant improvement over the pure cutoff, it still |
| 769 |
< |
doesn't correlate that well with SPME. Inclusion of potential damping |
| 770 |
< |
improves the results, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 752 |
> |
The most striking feature of this plot is how well the Shifted Force |
| 753 |
> |
({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy |
| 754 |
> |
differences. For the undamped {\sc sf} method, and the |
| 755 |
> |
moderately-damped {\sc sp} methods, the results are nearly |
| 756 |
> |
indistinguishable from the Ewald results. The other common methods do |
| 757 |
> |
significantly less well. |
| 758 |
> |
|
| 759 |
> |
The unmodified cutoff method is essentially unusable. This is not |
| 760 |
> |
surprising since hard cutoffs give large energy fluctuations as atoms |
| 761 |
> |
or molecules move in and out of the cutoff |
| 762 |
> |
radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to |
| 763 |
> |
some degree by using group based cutoffs with a switching |
| 764 |
> |
function.\cite{Adams79,Steinbach94,Leach01} However, we do not see |
| 765 |
> |
significant improvement using the group-switched cutoff because the |
| 766 |
> |
salt and salt solution systems contain non-neutral groups. Interested |
| 767 |
> |
readers can consult the accompanying supporting information for a |
| 768 |
> |
comparison where all groups are neutral. |
| 769 |
> |
|
| 770 |
> |
For the {\sc sp} method, inclusion of potential damping improves the |
| 771 |
> |
agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 772 |
|
an excellent correlation and quality of fit with the SPME results, |
| 773 |
< |
particularly with a cutoff radius greater than 12 \AA . Use of a |
| 774 |
< |
larger damping parameter is more helpful for the shortest cutoff |
| 775 |
< |
shown, but it has a detrimental effect on simulations with larger |
| 776 |
< |
cutoffs. In the {\sc sf} sets, increasing damping results in |
| 745 |
< |
progressively poorer correlation. Overall, the undamped case is the |
| 746 |
< |
best performing set, as the correlation and quality of fits are |
| 747 |
< |
consistently superior regardless of the cutoff distance. This result |
| 748 |
< |
is beneficial in that the undamped case is less computationally |
| 749 |
< |
prohibitive do to the lack of complimentary error function calculation |
| 750 |
< |
when performing the electrostatic pair interaction. The reaction |
| 751 |
< |
field results illustrates some of that method's limitations, primarily |
| 752 |
< |
that it was developed for use in homogenous systems; although it does |
| 753 |
< |
provide results that are an improvement over those from an unmodified |
| 754 |
< |
cutoff. |
| 773 |
> |
particularly with a cutoff radius greater than 12 |
| 774 |
> |
\AA . Use of a larger damping parameter is more helpful for the |
| 775 |
> |
shortest cutoff shown, but it has a detrimental effect on simulations |
| 776 |
> |
with larger cutoffs. |
| 777 |
|
|
| 778 |
+ |
In the {\sc sf} sets, increasing damping results in progressively |
| 779 |
+ |
worse correlation with Ewald. Overall, the undamped case is the best |
| 780 |
+ |
performing set, as the correlation and quality of fits are |
| 781 |
+ |
consistently superior regardless of the cutoff distance. The undamped |
| 782 |
+ |
case is also less computationally demanding (because no evaluation of |
| 783 |
+ |
the complementary error function is required). |
| 784 |
+ |
|
| 785 |
+ |
The reaction field results illustrates some of that method's |
| 786 |
+ |
limitations, primarily that it was developed for use in homogenous |
| 787 |
+ |
systems; although it does provide results that are an improvement over |
| 788 |
+ |
those from an unmodified cutoff. |
| 789 |
+ |
|
| 790 |
|
\subsection{Magnitudes of the Force and Torque Vectors} |
| 791 |
|
|
| 792 |
|
Evaluation of pairwise methods for use in Molecular Dynamics |
| 800 |
|
\begin{figure} |
| 801 |
|
\centering |
| 802 |
|
\includegraphics[width=5.5in]{./frcMagplot.pdf} |
| 803 |
< |
\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 SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 803 |
> |
\caption{Statistical analysis of the quality of the force vector |
| 804 |
> |
magnitudes for a given electrostatic method compared with the |
| 805 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 806 |
> |
indicate force magnitude values indistinguishable from those obtained |
| 807 |
> |
using SPME. Different values of the cutoff radius are indicated with |
| 808 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 809 |
> |
inverted triangles).} |
| 810 |
|
\label{fig:frcMag} |
| 811 |
|
\end{figure} |
| 812 |
|
|
| 832 |
|
\begin{figure} |
| 833 |
|
\centering |
| 834 |
|
\includegraphics[width=5.5in]{./trqMagplot.pdf} |
| 835 |
< |
\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 SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 835 |
> |
\caption{Statistical analysis of the quality of the torque vector |
| 836 |
> |
magnitudes for a given electrostatic method compared with the |
| 837 |
> |
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 838 |
> |
indicate torque magnitude values indistinguishable from those obtained |
| 839 |
> |
using SPME. Different values of the cutoff radius are indicated with |
| 840 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 841 |
> |
inverted triangles).} |
| 842 |
|
\label{fig:trqMag} |
| 843 |
|
\end{figure} |
| 844 |
|
|
| 868 |
|
\begin{figure} |
| 869 |
|
\centering |
| 870 |
|
\includegraphics[width=5.5in]{./frcTrqAngplot.pdf} |
| 871 |
< |
\caption{Statistical analysis of the quality of the Gaussian fit of the force and torque vector angular distributions for a given electrostatic method compared with 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 SPME. Different values of the cutoff radius are indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = inverted triangles).} |
| 871 |
> |
\caption{Statistical analysis of the quality of the Gaussian fit of |
| 872 |
> |
the force and torque vector angular distributions for a given |
| 873 |
> |
electrostatic method compared with the reference Ewald sum. Results |
| 874 |
> |
with a variance ($\sigma^2$) equal to zero (dashed line) indicate |
| 875 |
> |
force and torque directions indistinguishable from those obtained |
| 876 |
> |
using SPME. Different values of the cutoff radius are indicated with |
| 877 |
> |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 878 |
> |
inverted triangles).} |
| 879 |
|
\label{fig:frcTrqAng} |
| 880 |
|
\end{figure} |
| 881 |
|
|
| 899 |
|
|
| 900 |
|
\begin{table}[htbp] |
| 901 |
|
\centering |
| 902 |
< |
\caption{Variance ($\sigma^2$) of the force (top set) and torque (bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods. Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function. The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.} |
| 902 |
> |
\caption{Variance ($\sigma^2$) of the force (top set) and torque |
| 903 |
> |
(bottom set) vector angle difference distributions for the Shifted Potential and Shifted Force methods. Calculations were performed both with (Y) and without (N) group based cutoffs and a switching function. The $\alpha$ values have units of \AA$^{-1}$ and the variance values have units of degrees$^2$.} |
| 904 |
|
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 905 |
|
\\ |
| 906 |
|
\toprule |
| 988 |
|
\begin{figure} |
| 989 |
|
\centering |
| 990 |
|
\includegraphics[width = \linewidth]{./vCorrPlot.pdf} |
| 991 |
< |
\caption{Velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset is a magnification of the first trough. The times to first collision are nearly identical, but the differences can be seen in the peaks and troughs, where the undamped to weakly damped methods are stiffer than the moderately damped and SPME methods.} |
| 991 |
> |
\caption{Velocity auto-correlation functions of NaCl crystals at |
| 992 |
> |
1000 K while using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and |
| 993 |
> |
{\sc sp} ($\alpha$ = 0.2). The inset is a magnification of the first |
| 994 |
> |
trough. The times to first collision are nearly identical, but the |
| 995 |
> |
differences can be seen in the peaks and troughs, where the undamped |
| 996 |
> |
to weakly damped methods are stiffer than the moderately damped and |
| 997 |
> |
SPME methods.} |
| 998 |
|
\label{fig:vCorrPlot} |
| 999 |
|
\end{figure} |
| 1000 |
|
|
| 1029 |
|
\begin{figure} |
| 1030 |
|
\centering |
| 1031 |
|
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 1032 |
< |
\caption{Power spectra obtained from the velocity auto-correlation functions of NaCl crystals at 1000 K while using SPME, {\sc sf} ($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). Apodization of the correlation functions via a cubic switching function between 40 and 50 ps was used to clear up the spectral noise resulting from data truncation, and had no noticeable effect on peak location or magnitude. The inset shows the frequency region below 100 cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 1032 |
> |
\caption{Power spectra obtained from the velocity auto-correlation |
| 1033 |
> |
functions of NaCl crystals at 1000 K while using SPME, {\sc sf} |
| 1034 |
> |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). |
| 1035 |
> |
Apodization of the correlation functions via a cubic switching |
| 1036 |
> |
function between 40 and 50 ps was used to clear up the spectral noise |
| 1037 |
> |
resulting from data truncation, and had no noticeable effect on peak |
| 1038 |
> |
location or magnitude. The inset shows the frequency region below 100 |
| 1039 |
> |
cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 1040 |
|
\label{fig:methodPS} |
| 1041 |
|
\end{figure} |
| 1042 |
|
|
| 1078 |
|
\begin{figure} |
| 1079 |
|
\centering |
| 1080 |
|
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 1081 |
< |
\caption{Regions of spectra showing the low-frequency correlated motions for NaCl crystals at 1000 K using various electrostatic summation methods. The upper plot is a zoomed inset from figure \ref{fig:methodPS}. As the damping value for the {\sc sf} potential increases, the low-frequency peaks red-shift. The lower plot is of spectra 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 peaks are red-shifted toward 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.} |
| 1081 |
> |
\caption{Regions of spectra showing the low-frequency correlated |
| 1082 |
> |
motions for NaCl crystals at 1000 K using various electrostatic |
| 1083 |
> |
summation methods. The upper plot is a zoomed inset from figure |
| 1084 |
> |
\ref{fig:methodPS}. As the damping value for the {\sc sf} potential |
| 1085 |
> |
increases, the low-frequency peaks red-shift. The lower plot is of |
| 1086 |
> |
spectra when using SPME and a simple damped Coulombic sum with damping |
| 1087 |
> |
coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$. As |
| 1088 |
> |
$\alpha$ increases, the peaks are red-shifted toward and eventually |
| 1089 |
> |
beyond the values given by SPME. The larger $\alpha$ values weaken |
| 1090 |
> |
the real-space electrostatics, explaining this shift towards less |
| 1091 |
> |
strongly correlated motions in the crystal.} |
| 1092 |
|
\label{fig:dampInc} |
| 1093 |
|
\end{figure} |
| 1094 |
|
|
