| 25 |
|
|
| 26 |
|
\begin{document} |
| 27 |
|
|
| 28 |
< |
\title{Is the Ewald Summation necessary? Pairwise alternatives to the accepted standard for long-range electrostatics} |
| 28 |
> |
\title{Is the Ewald summation still necessary? \\ |
| 29 |
> |
Pairwise alternatives to the accepted standard for \\ |
| 30 |
> |
long-range electrostatics in molecular simulations} |
| 31 |
|
|
| 32 |
|
\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
| 33 |
|
gezelter@nd.edu} \\ |
| 40 |
|
\maketitle |
| 41 |
|
\doublespacing |
| 42 |
|
|
| 41 |
– |
\nobibliography{} |
| 43 |
|
\begin{abstract} |
| 44 |
< |
A new method for accumulating electrostatic interactions was derived |
| 45 |
< |
from the previous efforts described in \bibentry{Wolf99} and |
| 46 |
< |
\bibentry{Zahn02} as a possible replacement for lattice sum methods in |
| 47 |
< |
molecular simulations. Comparisons were performed with this and other |
| 48 |
< |
pairwise electrostatic summation techniques against the smooth |
| 49 |
< |
particle mesh Ewald (SPME) summation to see how well they reproduce |
| 50 |
< |
the energetics and dynamics of a variety of simulation types. The |
| 51 |
< |
newly derived Shifted-Force technique shows a remarkable ability to |
| 52 |
< |
reproduce the behavior exhibited in simulations using SPME with an |
| 53 |
< |
$\mathscr{O}(N)$ computational cost, equivalent to merely the |
| 54 |
< |
real-space portion of the lattice summation. |
| 55 |
< |
|
| 44 |
> |
We investigate pairwise electrostatic interaction methods and show |
| 45 |
> |
that there are viable and computationally efficient $(\mathscr{O}(N))$ |
| 46 |
> |
alternatives to the Ewald summation for typical modern molecular |
| 47 |
> |
simulations. These methods are extended from the damped and |
| 48 |
> |
cutoff-neutralized Coulombic sum originally proposed by |
| 49 |
> |
[D. Wolf, P. Keblinski, S.~R. Phillpot, and J. Eggebrecht, {\it J. Chem. Phys.} {\bf 110}, 8255 (1999)] One of these, the damped shifted force method, shows |
| 50 |
> |
a remarkable ability to reproduce the energetic and dynamic |
| 51 |
> |
characteristics exhibited by simulations employing lattice summation |
| 52 |
> |
techniques. Comparisons were performed with this and other pairwise |
| 53 |
> |
methods against the smooth particle mesh Ewald ({\sc spme}) summation |
| 54 |
> |
to see how well they reproduce the energetics and dynamics of a |
| 55 |
> |
variety of simulation types. |
| 56 |
|
\end{abstract} |
| 57 |
|
|
| 58 |
|
\newpage |
| 95 |
|
regarding possible artifacts caused by the inherent periodicity of the |
| 96 |
|
explicit Ewald summation.\cite{Tobias01} |
| 97 |
|
|
| 98 |
< |
In this paper, we focus on a new set of shifted methods devised by |
| 98 |
> |
In this paper, we focus on a new set of pairwise methods devised by |
| 99 |
|
Wolf {\it et al.},\cite{Wolf99} which we further extend. These |
| 100 |
|
methods along with a few other mixed methods (i.e. reaction field) are |
| 101 |
|
compared with the smooth particle mesh Ewald |
| 106 |
|
to the direct pairwise sum. They also lack the added periodicity of |
| 107 |
|
the Ewald sum, so they can be used for systems which are non-periodic |
| 108 |
|
or which have one- or two-dimensional periodicity. Below, these |
| 109 |
< |
methods are evaluated using a variety of model systems to establish |
| 110 |
< |
their usability in molecular simulations. |
| 109 |
> |
methods are evaluated using a variety of model systems to |
| 110 |
> |
establish their usability in molecular simulations. |
| 111 |
|
|
| 112 |
|
\subsection{The Ewald Sum} |
| 113 |
< |
The complete accumulation electrostatic interactions in a system with |
| 113 |
> |
The complete accumulation of the electrostatic interactions in a system with |
| 114 |
|
periodic boundary conditions (PBC) requires the consideration of the |
| 115 |
|
effect of all charges within a (cubic) simulation box as well as those |
| 116 |
|
in the periodic replicas, |
| 155 |
|
system is said to be using conducting (or ``tin-foil'') boundary |
| 156 |
|
conditions, $\epsilon_{\rm S} = \infty$. Figure |
| 157 |
|
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
| 158 |
< |
time. Initially, due to the small sizes of the systems that could be |
| 159 |
< |
feasibly simulated, the entire simulation box was replicated to |
| 160 |
< |
convergence. In more modern simulations, the simulation boxes have |
| 161 |
< |
grown large enough that a real-space cutoff could potentially give |
| 162 |
< |
convergent behavior. Indeed, it has often been observed that the |
| 163 |
< |
reciprocal-space portion of the Ewald sum can be small and rapidly |
| 164 |
< |
convergent compared to the real-space portion with the choice of small |
| 165 |
< |
$\alpha$.\cite{Karasawa89,Kolafa92} |
| 158 |
> |
time. Initially, due to the small system sizes that could be |
| 159 |
> |
simulated feasibly, the entire simulation box was replicated to |
| 160 |
> |
convergence. In more modern simulations, the systems have grown large |
| 161 |
> |
enough that a real-space cutoff could potentially give convergent |
| 162 |
> |
behavior. Indeed, it has been observed that with the choice of a |
| 163 |
> |
small $\alpha$, the reciprocal-space portion of the Ewald sum can be |
| 164 |
> |
rapidly convergent and small relative to the real-space |
| 165 |
> |
portion.\cite{Karasawa89,Kolafa92} |
| 166 |
|
|
| 167 |
|
\begin{figure} |
| 168 |
|
\centering |
| 169 |
|
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 170 |
< |
\caption{How the application of the Ewald summation has changed with |
| 171 |
< |
the increase in computer power. Initially, only small numbers of |
| 172 |
< |
particles could be studied, and the Ewald sum acted to replicate the |
| 173 |
< |
unit cell charge distribution out to convergence. Now, much larger |
| 174 |
< |
systems of charges are investigated with fixed distance cutoffs. The |
| 174 |
< |
calculated structure factor is used to sum out to great distance, and |
| 175 |
< |
a surrounding dielectric term is included.} |
| 170 |
> |
\caption{The change in the need for the Ewald sum with |
| 171 |
> |
increasing computational power. A:~Initially, only small systems |
| 172 |
> |
could be studied, and the Ewald sum replicated the simulation box to |
| 173 |
> |
convergence. B:~Now, radial cutoff methods should be able to reach |
| 174 |
> |
convergence for the larger systems of charges that are common today.} |
| 175 |
|
\label{fig:ewaldTime} |
| 176 |
|
\end{figure} |
| 177 |
|
|
| 201 |
|
interfaces and membranes, the intrinsic three-dimensional periodicity |
| 202 |
|
can prove problematic. The Ewald sum has been reformulated to handle |
| 203 |
|
2D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the |
| 204 |
< |
new methods are computationally expensive.\cite{Spohr97,Yeh99} |
| 205 |
< |
Inclusion of a correction term in the Ewald summation is a possible |
| 206 |
< |
direction for handling 2D systems while still enabling the use of the |
| 207 |
< |
modern optimizations.\cite{Yeh99} |
| 204 |
> |
new methods are computationally expensive.\cite{Spohr97,Yeh99} More |
| 205 |
> |
recently, there have been several successful efforts toward reducing |
| 206 |
> |
the computational cost of 2D lattice summations, often enabling the |
| 207 |
> |
use of the mentioned |
| 208 |
> |
optimizations.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} |
| 209 |
|
|
| 210 |
|
Several studies have recognized that the inherent periodicity in the |
| 211 |
|
Ewald sum can also have an effect on three-dimensional |
| 228 |
|
charge contained within the cutoff radius is crucial for potential |
| 229 |
|
stability. They devised a pairwise summation method that ensures |
| 230 |
|
charge neutrality and gives results similar to those obtained with the |
| 231 |
< |
Ewald summation. The resulting shifted Coulomb potential |
| 232 |
< |
(Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through |
| 233 |
< |
placement on the cutoff sphere and a distance-dependent damping |
| 234 |
< |
function (identical to that seen in the real-space portion of the |
| 235 |
< |
Ewald sum) to aid convergence |
| 231 |
> |
Ewald summation. The resulting shifted Coulomb potential includes |
| 232 |
> |
image-charges subtracted out through placement on the cutoff sphere |
| 233 |
> |
and a distance-dependent damping function (identical to that seen in |
| 234 |
> |
the real-space portion of the Ewald sum) to aid convergence |
| 235 |
|
\begin{equation} |
| 236 |
|
V_{\textrm{Wolf}}(r_{ij})= \frac{q_i q_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}-\lim_{r_{ij}\rightarrow R_\textrm{c}}\left\{\frac{q_iq_j \textrm{erfc}(\alpha r_{ij})}{r_{ij}}\right\}. |
| 237 |
|
\label{eq:WolfPot} |
| 517 |
|
The pairwise summation techniques (outlined in section |
| 518 |
|
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 519 |
|
studying the energy differences between conformations. We took the |
| 520 |
< |
SPME-computed energy difference between two conformations to be the |
| 520 |
> |
{\sc spme}-computed energy difference between two conformations to be the |
| 521 |
|
correct behavior. An ideal performance by an alternative method would |
| 522 |
|
reproduce these energy differences exactly (even if the absolute |
| 523 |
|
energies calculated by the methods are different). Since none of the |
| 525 |
|
regressions of energy gap data to evaluate how closely the methods |
| 526 |
|
mimicked the Ewald energy gaps. Unitary results for both the |
| 527 |
|
correlation (slope) and correlation coefficient for these regressions |
| 528 |
< |
indicate perfect agreement between the alternative method and SPME. |
| 528 |
> |
indicate perfect agreement between the alternative method and {\sc spme}. |
| 529 |
|
Sample correlation plots for two alternate methods are shown in |
| 530 |
|
Fig. \ref{fig:linearFit}. |
| 531 |
|
|
| 539 |
|
\label{fig:linearFit} |
| 540 |
|
\end{figure} |
| 541 |
|
|
| 542 |
< |
Each system type (detailed in section \ref{sec:RepSims}) was |
| 543 |
< |
represented using 500 independent configurations. Additionally, we |
| 544 |
< |
used seven different system types, so each of the alternative |
| 545 |
< |
(non-Ewald) electrostatic summation methods was evaluated using |
| 546 |
< |
873,250 configurational energy differences. |
| 542 |
> |
Each of the seven system types (detailed in section \ref{sec:RepSims}) |
| 543 |
> |
were represented using 500 independent configurations. Thus, each of |
| 544 |
> |
the alternative (non-Ewald) electrostatic summation methods was |
| 545 |
> |
evaluated using an accumulated 873,250 configurational energy |
| 546 |
> |
differences. |
| 547 |
|
|
| 548 |
|
Results and discussion for the individual analysis of each of the |
| 549 |
|
system types appear in the supporting information, while the |
| 554 |
|
We evaluated the pairwise methods (outlined in section |
| 555 |
|
\ref{sec:ESMethods}) for use in MD simulations by |
| 556 |
|
comparing the force and torque vectors with those obtained using the |
| 557 |
< |
reference Ewald summation (SPME). Both the magnitude and the |
| 557 |
> |
reference Ewald summation ({\sc spme}). Both the magnitude and the |
| 558 |
|
direction of these vectors on each of the bodies in the system were |
| 559 |
|
analyzed. For the magnitude of these vectors, linear least squares |
| 560 |
|
regression analyses were performed as described previously for |
| 569 |
|
|
| 570 |
|
The {\it directionality} of the force and torque vectors was |
| 571 |
|
investigated through measurement of the angle ($\theta$) formed |
| 572 |
< |
between those computed from the particular method and those from SPME, |
| 572 |
> |
between those computed from the particular method and those from {\sc spme}, |
| 573 |
|
\begin{equation} |
| 574 |
|
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right), |
| 575 |
|
\end{equation} |
| 576 |
< |
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force |
| 576 |
> |
where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force |
| 577 |
|
vector computed using method M. Each of these $\theta$ values was |
| 578 |
|
accumulated in a distribution function and weighted by the area on the |
| 579 |
|
unit sphere. Since this distribution is a measure of angular error |
| 580 |
|
between two different electrostatic summation methods, there is no |
| 581 |
|
{\it a priori} reason for the profile to adhere to any specific |
| 582 |
|
shape. Thus, gaussian fits were used to measure the width of the |
| 583 |
< |
resulting distributions. |
| 584 |
< |
% |
| 585 |
< |
%\begin{figure} |
| 586 |
< |
%\centering |
| 587 |
< |
%\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 589 |
< |
%\caption{Sample fit of the angular distribution of the force vectors |
| 590 |
< |
%accumulated using all of the studied systems. Gaussian fits were used |
| 591 |
< |
%to obtain values for the variance in force and torque vectors.} |
| 592 |
< |
%\label{fig:gaussian} |
| 593 |
< |
%\end{figure} |
| 594 |
< |
% |
| 595 |
< |
%Figure \ref{fig:gaussian} shows an example distribution with applied |
| 596 |
< |
%non-linear fits. The solid line is a Gaussian profile, while the |
| 597 |
< |
%dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 598 |
< |
%Lorentzian. |
| 599 |
< |
%Since this distribution is a measure of angular error between two |
| 600 |
< |
%different electrostatic summation methods, there is no {\it a priori} |
| 601 |
< |
%reason for the profile to adhere to any specific shape. |
| 602 |
< |
%Gaussian fits was used to compare all the tested methods. |
| 603 |
< |
The variance ($\sigma^2$) was extracted from each of these fits and |
| 604 |
< |
was used to compare distribution widths. Values of $\sigma^2$ near |
| 605 |
< |
zero indicate vector directions indistinguishable from those |
| 606 |
< |
calculated when using the reference method (SPME). |
| 583 |
> |
resulting distributions. The variance ($\sigma^2$) was extracted from |
| 584 |
> |
each of these fits and was used to compare distribution widths. |
| 585 |
> |
Values of $\sigma^2$ near zero indicate vector directions |
| 586 |
> |
indistinguishable from those calculated when using the reference |
| 587 |
> |
method ({\sc spme}). |
| 588 |
|
|
| 589 |
|
\subsection{Short-time Dynamics} |
| 590 |
|
|
| 599 |
|
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
| 600 |
|
of the trajectories, |
| 601 |
|
\begin{equation} |
| 602 |
< |
C_v(t) = \frac{\langle v_i(0)\cdot v_i(t)\rangle}{\langle v_i(0)\cdot v_i(0)\rangle}. |
| 602 |
> |
C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}. |
| 603 |
|
\label{eq:vCorr} |
| 604 |
|
\end{equation} |
| 605 |
|
Velocity autocorrelation functions require detailed short time data, |
| 612 |
|
|
| 613 |
|
The effects of the same subset of alternative electrostatic methods on |
| 614 |
|
the {\it long-time} dynamics of charged systems were evaluated using |
| 615 |
< |
the same model system (NaCl crystals at 1000K). The power spectrum |
| 615 |
> |
the same model system (NaCl crystals at 1000~K). The power spectrum |
| 616 |
|
($I(\omega)$) was obtained via Fourier transform of the velocity |
| 617 |
|
autocorrelation function, \begin{equation} I(\omega) = |
| 618 |
|
\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 622 |
|
NaCl crystal is composed of two different atom types, the average of |
| 623 |
|
the two resulting power spectra was used for comparisons. Simulations |
| 624 |
|
were performed under the microcanonical ensemble, and velocity |
| 625 |
< |
information was saved every 5 fs over 100 ps trajectories. |
| 625 |
> |
information was saved every 5~fs over 100~ps trajectories. |
| 626 |
|
|
| 627 |
|
\subsection{Representative Simulations}\label{sec:RepSims} |
| 628 |
< |
A variety of representative simulations were analyzed to determine the |
| 629 |
< |
relative effectiveness of the pairwise summation techniques in |
| 630 |
< |
reproducing the energetics and dynamics exhibited by SPME. We wanted |
| 631 |
< |
to span the space of modern simulations (i.e. from liquids of neutral |
| 632 |
< |
molecules to ionic crystals), so the systems studied were: |
| 628 |
> |
A variety of representative molecular simulations were analyzed to |
| 629 |
> |
determine the relative effectiveness of the pairwise summation |
| 630 |
> |
techniques in reproducing the energetics and dynamics exhibited by |
| 631 |
> |
{\sc spme}. We wanted to span the space of typical molecular |
| 632 |
> |
simulations (i.e. from liquids of neutral molecules to ionic |
| 633 |
> |
crystals), so the systems studied were: |
| 634 |
|
\begin{enumerate} |
| 635 |
|
\item liquid water (SPC/E),\cite{Berendsen87} |
| 636 |
|
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), |
| 653 |
|
the crystal). The solid and liquid NaCl systems consisted of 500 |
| 654 |
|
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
| 655 |
|
these systems were selected and equilibrated in the same manner as the |
| 656 |
< |
water systems. The equilibrated temperatures were 1000~K for the NaCl |
| 657 |
< |
crystal and 7000~K for the liquid. The ionic solutions were made by |
| 658 |
< |
solvating 4 (or 40) ions in a periodic box containing 1000 SPC/E water |
| 659 |
< |
molecules. Ion and water positions were then randomly swapped, and |
| 660 |
< |
the resulting configurations were again equilibrated individually. |
| 661 |
< |
Finally, for the Argon / Water ``charge void'' systems, the identities |
| 662 |
< |
of all the SPC/E waters within 6 \AA\ of the center of the |
| 663 |
< |
equilibrated water configurations were converted to argon. |
| 664 |
< |
%(Fig. \ref{fig:argonSlice}). |
| 656 |
> |
water systems. In order to introduce measurable fluctuations in the |
| 657 |
> |
configuration energy differences, the crystalline simulations were |
| 658 |
> |
equilibrated at 1000~K, near the $T_\textrm{m}$ for NaCl. The liquid |
| 659 |
> |
NaCl configurations needed to represent a fully disordered array of |
| 660 |
> |
point charges, so the high temperature of 7000~K was selected for |
| 661 |
> |
equilibration. The ionic solutions were made by solvating 4 (or 40) |
| 662 |
> |
ions in a periodic box containing 1000 SPC/E water molecules. Ion and |
| 663 |
> |
water positions were then randomly swapped, and the resulting |
| 664 |
> |
configurations were again equilibrated individually. Finally, for the |
| 665 |
> |
Argon / Water ``charge void'' systems, the identities of all the SPC/E |
| 666 |
> |
waters within 6 \AA\ of the center of the equilibrated water |
| 667 |
> |
configurations were converted to argon. |
| 668 |
|
|
| 669 |
|
These procedures guaranteed us a set of representative configurations |
| 670 |
< |
from chemically-relevant systems sampled from an appropriate |
| 671 |
< |
ensemble. Force field parameters for the ions and Argon were taken |
| 670 |
> |
from chemically-relevant systems sampled from appropriate |
| 671 |
> |
ensembles. Force field parameters for the ions and Argon were taken |
| 672 |
|
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 688 |
– |
|
| 689 |
– |
%\begin{figure} |
| 690 |
– |
%\centering |
| 691 |
– |
%\includegraphics[width = \linewidth]{./slice.pdf} |
| 692 |
– |
%\caption{A slice from the center of a water box used in a charge void |
| 693 |
– |
%simulation. The darkened region represents the boundary sphere within |
| 694 |
– |
%which the water molecules were converted to argon atoms.} |
| 695 |
– |
%\label{fig:argonSlice} |
| 696 |
– |
%\end{figure} |
| 673 |
|
|
| 674 |
|
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 675 |
|
We compared the following alternative summation methods with results |
| 676 |
< |
from the reference method (SPME): |
| 676 |
> |
from the reference method ({\sc spme}): |
| 677 |
|
\begin{itemize} |
| 678 |
|
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 679 |
|
and 0.3 \AA$^{-1}$, |
| 684 |
|
\end{itemize} |
| 685 |
|
Group-based cutoffs with a fifth-order polynomial switching function |
| 686 |
|
were utilized for the reaction field simulations. Additionally, we |
| 687 |
< |
investigated the use of these cutoffs with the SP, SF, and pure |
| 688 |
< |
cutoff. The SPME electrostatics were performed using the TINKER |
| 689 |
< |
implementation of SPME,\cite{Ponder87} while all other method |
| 690 |
< |
calculations were performed using the OOPSE molecular mechanics |
| 687 |
> |
investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure |
| 688 |
> |
cutoff. The {\sc spme} electrostatics were performed using the {\sc tinker} |
| 689 |
> |
implementation of {\sc spme},\cite{Ponder87} while all other calculations |
| 690 |
> |
were performed using the {\sc oopse} molecular mechanics |
| 691 |
|
package.\cite{Meineke05} All other portions of the energy calculation |
| 692 |
|
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 693 |
|
manner across all systems and configurations. |
| 694 |
|
|
| 695 |
< |
The althernative methods were also evaluated with three different |
| 695 |
> |
The alternative methods were also evaluated with three different |
| 696 |
|
cutoff radii (9, 12, and 15 \AA). As noted previously, the |
| 697 |
|
convergence parameter ($\alpha$) plays a role in the balance of the |
| 698 |
|
real-space and reciprocal-space portions of the Ewald calculation. |
| 699 |
|
Typical molecular mechanics packages set this to a value dependent on |
| 700 |
|
the cutoff radius and a tolerance (typically less than $1 \times |
| 701 |
|
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 702 |
< |
increased accuracy at the expense of increased time spent calculating |
| 703 |
< |
the reciprocal-space portion of the |
| 704 |
< |
summation.\cite{Perram88,Essmann95} The default TINKER tolerance of $1 |
| 705 |
< |
\times 10^{-8}$ kcal/mol was used in all SPME calculations, resulting |
| 706 |
< |
in Ewald Coefficients of 0.4200, 0.3119, and 0.2476 \AA$^{-1}$ for |
| 707 |
< |
cutoff radii of 9, 12, and 15 \AA\ respectively. |
| 702 |
> |
increasing accuracy at the expense of computational time spent on the |
| 703 |
> |
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 704 |
> |
The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 705 |
> |
in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200, |
| 706 |
> |
0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ |
| 707 |
> |
respectively. |
| 708 |
|
|
| 709 |
|
\section{Results and Discussion} |
| 710 |
|
|
| 712 |
|
In order to evaluate the performance of the pairwise electrostatic |
| 713 |
|
summation methods for Monte Carlo simulations, the energy differences |
| 714 |
|
between configurations were compared to the values obtained when using |
| 715 |
< |
SPME. The results for the subsequent regression analysis are shown in |
| 715 |
> |
{\sc spme}. The results for the subsequent regression analysis are shown in |
| 716 |
|
figure \ref{fig:delE}. |
| 717 |
|
|
| 718 |
|
\begin{figure} |
| 722 |
|
differences for a given electrostatic method compared with the |
| 723 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 724 |
|
indicate $\Delta E$ values indistinguishable from those obtained using |
| 725 |
< |
SPME. Different values of the cutoff radius are indicated with |
| 725 |
> |
{\sc spme}. Different values of the cutoff radius are indicated with |
| 726 |
|
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 727 |
|
inverted triangles).} |
| 728 |
|
\label{fig:delE} |
| 746 |
|
readers can consult the accompanying supporting information for a |
| 747 |
|
comparison where all groups are neutral. |
| 748 |
|
|
| 749 |
< |
For the {\sc sp} method, inclusion of potential damping improves the |
| 750 |
< |
agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ shows |
| 751 |
< |
an excellent correlation and quality of fit with the SPME results, |
| 752 |
< |
particularly with a cutoff radius greater than 12 |
| 749 |
> |
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 750 |
> |
the agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ |
| 751 |
> |
shows an excellent correlation and quality of fit with the {\sc spme} |
| 752 |
> |
results, particularly with a cutoff radius greater than 12 |
| 753 |
|
\AA . Use of a larger damping parameter is more helpful for the |
| 754 |
|
shortest cutoff shown, but it has a detrimental effect on simulations |
| 755 |
|
with larger cutoffs. |
| 756 |
|
|
| 757 |
< |
In the {\sc sf} sets, increasing damping results in progressively |
| 758 |
< |
worse correlation with Ewald. Overall, the undamped case is the best |
| 757 |
> |
In the {\sc sf} sets, increasing damping results in progressively {\it |
| 758 |
> |
worse} correlation with Ewald. Overall, the undamped case is the best |
| 759 |
|
performing set, as the correlation and quality of fits are |
| 760 |
|
consistently superior regardless of the cutoff distance. The undamped |
| 761 |
|
case is also less computationally demanding (because no evaluation of |
| 770 |
|
|
| 771 |
|
Evaluation of pairwise methods for use in Molecular Dynamics |
| 772 |
|
simulations requires consideration of effects on the forces and |
| 773 |
< |
torques. Investigation of the force and torque vector magnitudes |
| 774 |
< |
provides a measure of the strength of these values relative to SPME. |
| 775 |
< |
Figures \ref{fig:frcMag} and \ref{fig:trqMag} respectively show the |
| 776 |
< |
force and torque vector magnitude regression results for the |
| 801 |
< |
accumulated analysis over all the system types. |
| 773 |
> |
torques. Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the |
| 774 |
> |
regression results for the force and torque vector magnitudes, |
| 775 |
> |
respectively. The data in these figures was generated from an |
| 776 |
> |
accumulation of the statistics from all of the system types. |
| 777 |
|
|
| 778 |
|
\begin{figure} |
| 779 |
|
\centering |
| 782 |
|
magnitudes for a given electrostatic method compared with the |
| 783 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 784 |
|
indicate force magnitude values indistinguishable from those obtained |
| 785 |
< |
using SPME. Different values of the cutoff radius are indicated with |
| 785 |
> |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 786 |
|
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 787 |
|
inverted triangles).} |
| 788 |
|
\label{fig:frcMag} |
| 789 |
|
\end{figure} |
| 790 |
|
|
| 791 |
+ |
Again, it is striking how well the Shifted Potential and Shifted Force |
| 792 |
+ |
methods are doing at reproducing the {\sc spme} forces. The undamped and |
| 793 |
+ |
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
| 794 |
+ |
This is perhaps expected because this method explicitly incorporates a |
| 795 |
+ |
smooth transition in the forces at the cutoff radius as well as the |
| 796 |
+ |
neutralizing image charges. |
| 797 |
+ |
|
| 798 |
|
Figure \ref{fig:frcMag}, for the most part, parallels the results seen |
| 799 |
|
in the previous $\Delta E$ section. The unmodified cutoff results are |
| 800 |
|
poor, but using group based cutoffs and a switching function provides |
| 801 |
< |
a improvement much more significant than what was seen with $\Delta |
| 802 |
< |
E$. Looking at the {\sc sp} sets, the slope and $R^2$ |
| 803 |
< |
improve with the use of damping to an optimal result of 0.2 \AA |
| 804 |
< |
$^{-1}$ for the 12 and 15 \AA\ cutoffs. Further increases in damping, |
| 801 |
> |
an improvement much more significant than what was seen with $\Delta |
| 802 |
> |
E$. |
| 803 |
> |
|
| 804 |
> |
With moderate damping and a large enough cutoff radius, the {\sc sp} |
| 805 |
> |
method is generating usable forces. Further increases in damping, |
| 806 |
|
while beneficial for simulations with a cutoff radius of 9 \AA\ , is |
| 807 |
< |
detrimental to simulations with larger cutoff radii. The undamped |
| 808 |
< |
{\sc sf} method gives forces in line with those obtained using |
| 809 |
< |
SPME, and use of a damping function results in minor improvement. The |
| 827 |
< |
reaction field results are surprisingly good, considering the poor |
| 807 |
> |
detrimental to simulations with larger cutoff radii. |
| 808 |
> |
|
| 809 |
> |
The reaction field results are surprisingly good, considering the poor |
| 810 |
|
quality of the fits for the $\Delta E$ results. There is still a |
| 811 |
< |
considerable degree of scatter in the data, but it correlates well in |
| 812 |
< |
general. To be fair, we again note that the reaction field |
| 813 |
< |
calculations do not encompass NaCl crystal and melt systems, so these |
| 811 |
> |
considerable degree of scatter in the data, but the forces correlate |
| 812 |
> |
well with the Ewald forces in general. We note that the reaction |
| 813 |
> |
field calculations do not include the pure NaCl systems, so these |
| 814 |
|
results are partly biased towards conditions in which the method |
| 815 |
|
performs more favorably. |
| 816 |
|
|
| 821 |
|
magnitudes for a given electrostatic method compared with the |
| 822 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 823 |
|
indicate torque magnitude values indistinguishable from those obtained |
| 824 |
< |
using SPME. Different values of the cutoff radius are indicated with |
| 824 |
> |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 825 |
|
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 826 |
|
inverted triangles).} |
| 827 |
|
\label{fig:trqMag} |
| 828 |
|
\end{figure} |
| 829 |
|
|
| 830 |
< |
To evaluate the torque vector magnitudes, the data set from which |
| 831 |
< |
values are drawn is limited to rigid molecules in the systems |
| 832 |
< |
(i.e. water molecules). In spite of this smaller sampling pool, the |
| 851 |
< |
torque vector magnitude results in figure \ref{fig:trqMag} are still |
| 852 |
< |
similar to those seen for the forces; however, they more clearly show |
| 853 |
< |
the improved behavior that comes with increasing the cutoff radius. |
| 854 |
< |
Moderate damping is beneficial to the {\sc sp} and helpful |
| 855 |
< |
yet possibly unnecessary with the {\sc sf} method, and they also |
| 856 |
< |
show that over-damping adversely effects all cutoff radii rather than |
| 857 |
< |
showing an improvement for systems with short cutoffs. The reaction |
| 858 |
< |
field method performs well when calculating the torques, better than |
| 859 |
< |
the Shifted Force method over this limited data set. |
| 830 |
> |
Molecular torques were only available from the systems which contained |
| 831 |
> |
rigid molecules (i.e. the systems containing water). The data in |
| 832 |
> |
fig. \ref{fig:trqMag} is taken from this smaller sampling pool. |
| 833 |
|
|
| 834 |
+ |
Torques appear to be much more sensitive to charges at a longer |
| 835 |
+ |
distance. The striking feature in comparing the new electrostatic |
| 836 |
+ |
methods with {\sc spme} is how much the agreement improves with increasing |
| 837 |
+ |
cutoff radius. Again, the weakly damped and undamped {\sc sf} method |
| 838 |
+ |
appears to be reproducing the {\sc spme} torques most accurately. |
| 839 |
+ |
|
| 840 |
+ |
Water molecules are dipolar, and the reaction field method reproduces |
| 841 |
+ |
the effect of the surrounding polarized medium on each of the |
| 842 |
+ |
molecular bodies. Therefore it is not surprising that reaction field |
| 843 |
+ |
performs best of all of the methods on molecular torques. |
| 844 |
+ |
|
| 845 |
|
\subsection{Directionality of the Force and Torque Vectors} |
| 846 |
|
|
| 847 |
< |
Having force and torque vectors with magnitudes that are well |
| 848 |
< |
correlated to SPME is good, but if they are not pointing in the proper |
| 849 |
< |
direction the results will be incorrect. These vector directions were |
| 850 |
< |
investigated through measurement of the angle formed between them and |
| 851 |
< |
those from SPME. The results (Fig. \ref{fig:frcTrqAng}) are compared |
| 852 |
< |
through the variance ($\sigma^2$) of the Gaussian fits of the angle |
| 853 |
< |
error distributions of the combined set over all system types. |
| 847 |
> |
It is clearly important that a new electrostatic method can reproduce |
| 848 |
> |
the magnitudes of the force and torque vectors obtained via the Ewald |
| 849 |
> |
sum. However, the {\it directionality} of these vectors will also be |
| 850 |
> |
vital in calculating dynamical quantities accurately. Force and |
| 851 |
> |
torque directionalities were investigated by measuring the angles |
| 852 |
> |
formed between these vectors and the same vectors calculated using |
| 853 |
> |
{\sc spme}. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 854 |
> |
variance ($\sigma^2$) of the Gaussian fits of the angle error |
| 855 |
> |
distributions of the combined set over all system types. |
| 856 |
|
|
| 857 |
|
\begin{figure} |
| 858 |
|
\centering |
| 859 |
|
\includegraphics[width=5.5in]{./frcTrqAngplot.pdf} |
| 860 |
< |
\caption{Statistical analysis of the quality of the Gaussian fit of |
| 861 |
< |
the force and torque vector angular distributions for a given |
| 862 |
< |
electrostatic method compared with the reference Ewald sum. Results |
| 863 |
< |
with a variance ($\sigma^2$) equal to zero (dashed line) indicate |
| 864 |
< |
force and torque directions indistinguishable from those obtained |
| 865 |
< |
using SPME. Different values of the cutoff radius are indicated with |
| 866 |
< |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 867 |
< |
inverted triangles).} |
| 860 |
> |
\caption{Statistical analysis of the width of the angular distribution |
| 861 |
> |
that the force and torque vectors from a given electrostatic method |
| 862 |
> |
make with their counterparts obtained using the reference Ewald sum. |
| 863 |
> |
Results with a variance ($\sigma^2$) equal to zero (dashed line) |
| 864 |
> |
indicate force and torque directions indistinguishable from those |
| 865 |
> |
obtained using {\sc spme}. Different values of the cutoff radius are |
| 866 |
> |
indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, |
| 867 |
> |
and 15\AA\ = inverted triangles).} |
| 868 |
|
\label{fig:frcTrqAng} |
| 869 |
|
\end{figure} |
| 870 |
|
|
| 871 |
|
Both the force and torque $\sigma^2$ results from the analysis of the |
| 872 |
|
total accumulated system data are tabulated in figure |
| 873 |
< |
\ref{fig:frcTrqAng}. All of the sets, aside from the over-damped case |
| 874 |
< |
show the improvement afforded by choosing a longer simulation cutoff. |
| 875 |
< |
Increasing the cutoff from 9 to 12 \AA\ typically results in a halving |
| 876 |
< |
of the distribution widths, with a similar improvement going from 12 |
| 877 |
< |
to 15 \AA . The undamped {\sc sf}, Group Based Cutoff, and |
| 878 |
< |
Reaction Field methods all do equivalently well at capturing the |
| 893 |
< |
direction of both the force and torque vectors. Using damping |
| 894 |
< |
improves the angular behavior significantly for the {\sc sp} |
| 895 |
< |
and moderately for the {\sc sf} methods. Increasing the damping |
| 896 |
< |
too far is destructive for both methods, particularly to the torque |
| 897 |
< |
vectors. Again it is important to recognize that the force vectors |
| 898 |
< |
cover all particles in the systems, while torque vectors are only |
| 899 |
< |
available for neutral molecular groups. Damping appears to have a |
| 900 |
< |
more beneficial effect on non-neutral bodies, and this observation is |
| 901 |
< |
investigated further in the accompanying supporting information. |
| 873 |
> |
\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc |
| 874 |
> |
sp}) method would be essentially unusable for molecular dynamics |
| 875 |
> |
unless the damping function is added. The Shifted Force ({\sc sf}) |
| 876 |
> |
method, however, is generating force and torque vectors which are |
| 877 |
> |
within a few degrees of the Ewald results even with weak (or no) |
| 878 |
> |
damping. |
| 879 |
|
|
| 880 |
+ |
All of the sets (aside from the over-damped case) show the improvement |
| 881 |
+ |
afforded by choosing a larger cutoff radius. Increasing the cutoff |
| 882 |
+ |
from 9 to 12 \AA\ typically results in a halving of the width of the |
| 883 |
+ |
distribution, with a similar improvement when going from 12 to 15 |
| 884 |
+ |
\AA . |
| 885 |
+ |
|
| 886 |
+ |
The undamped {\sc sf}, group-based cutoff, and reaction field methods |
| 887 |
+ |
all do equivalently well at capturing the direction of both the force |
| 888 |
+ |
and torque vectors. Using the electrostatic damping improves the |
| 889 |
+ |
angular behavior significantly for the {\sc sp} and moderately for the |
| 890 |
+ |
{\sc sf} methods. Overdamping is detrimental to both methods. Again |
| 891 |
+ |
it is important to recognize that the force vectors cover all |
| 892 |
+ |
particles in all seven systems, while torque vectors are only |
| 893 |
+ |
available for neutral molecular groups. Damping is more beneficial to |
| 894 |
+ |
charged bodies, and this observation is investigated further in the |
| 895 |
+ |
accompanying supporting information. |
| 896 |
+ |
|
| 897 |
+ |
Although not discussed previously, group based cutoffs can be applied |
| 898 |
+ |
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 899 |
+ |
will reintroduce small discontinuities at the cutoff radius, but the |
| 900 |
+ |
effects of these can be minimized by utilizing a switching function. |
| 901 |
+ |
Though there are no significant benefits or drawbacks observed in |
| 902 |
+ |
$\Delta E$ and the force and torque magnitudes when doing this, there |
| 903 |
+ |
is a measurable improvement in the directionality of the forces and |
| 904 |
+ |
torques. Table \ref{tab:groupAngle} shows the angular variances |
| 905 |
+ |
obtained using group based cutoffs along with the results seen in |
| 906 |
+ |
figure \ref{fig:frcTrqAng}. The {\sc sp} (with an $\alpha$ of 0.2 |
| 907 |
+ |
\AA$^{-1}$ or smaller) shows much narrower angular distributions when |
| 908 |
+ |
using group-based cutoffs. The {\sc sf} method likewise shows |
| 909 |
+ |
improvement in the undamped and lightly damped cases. |
| 910 |
+ |
|
| 911 |
|
\begin{table}[htbp] |
| 912 |
< |
\centering |
| 913 |
< |
\caption{Variance ($\sigma^2$) of the force (top set) and torque |
| 914 |
< |
(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$.} |
| 912 |
> |
\centering |
| 913 |
> |
\caption{Statistical analysis of the angular |
| 914 |
> |
distributions that the force (upper) and torque (lower) vectors |
| 915 |
> |
from a given electrostatic method make with their counterparts |
| 916 |
> |
obtained using the reference Ewald sum. Calculations were |
| 917 |
> |
performed both with (Y) and without (N) group based cutoffs and a |
| 918 |
> |
switching function. The $\alpha$ values have units of \AA$^{-1}$ |
| 919 |
> |
and the variance values have units of degrees$^2$.} |
| 920 |
> |
|
| 921 |
|
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 922 |
|
\\ |
| 923 |
|
\toprule |
| 948 |
|
\label{tab:groupAngle} |
| 949 |
|
\end{table} |
| 950 |
|
|
| 951 |
< |
Although not discussed previously, group based cutoffs can be applied |
| 952 |
< |
to both the {\sc sp} and {\sc sf} methods. Use off a |
| 953 |
< |
switching function corrects for the discontinuities that arise when |
| 954 |
< |
atoms of a group exit the cutoff before the group's center of mass. |
| 955 |
< |
Though there are no significant benefit or drawbacks observed in |
| 956 |
< |
$\Delta E$ and vector magnitude results when doing this, there is a |
| 957 |
< |
measurable improvement in the vector angle results. Table |
| 958 |
< |
\ref{tab:groupAngle} shows the angular variance values obtained using |
| 959 |
< |
group based cutoffs and a switching function alongside the standard |
| 960 |
< |
results seen in figure \ref{fig:frcTrqAng} for comparison purposes. |
| 961 |
< |
The {\sc sp} shows much narrower angular distributions for |
| 962 |
< |
both the force and torque vectors when using an $\alpha$ of 0.2 |
| 963 |
< |
\AA$^{-1}$ or less, while {\sc sf} shows improvements in the |
| 964 |
< |
undamped and lightly damped cases. Thus, by calculating the |
| 965 |
< |
electrostatic interactions in terms of molecular pairs rather than |
| 966 |
< |
atomic pairs, the direction of the force and torque vectors are |
| 967 |
< |
determined more accurately. |
| 951 |
> |
One additional trend in table \ref{tab:groupAngle} is that the |
| 952 |
> |
$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$ |
| 953 |
> |
increases, something that is more obvious with group-based cutoffs. |
| 954 |
> |
The complimentary error function inserted into the potential weakens |
| 955 |
> |
the electrostatic interaction as the value of $\alpha$ is increased. |
| 956 |
> |
However, at larger values of $\alpha$, it is possible to overdamp the |
| 957 |
> |
electrostatic interaction and to remove it completely. Kast |
| 958 |
> |
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 959 |
> |
values for these types of electrostatic summation methods by fitting |
| 960 |
> |
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
| 961 |
> |
0.25, and 0.16 \AA$^{-1}$ for cutoff values of 9, 12, and 15 \AA\ |
| 962 |
> |
respectively.\cite{Kast03} These appear to be reasonable choices to |
| 963 |
> |
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
| 964 |
> |
these findings, choices this high would introduce error in the |
| 965 |
> |
molecular torques, particularly for the shorter cutoffs. Based on our |
| 966 |
> |
observations, empirical damping up to 0.2 \AA$^{-1}$ is beneficial, |
| 967 |
> |
but damping may be unnecessary when using the {\sc sf} method. |
| 968 |
|
|
| 955 |
– |
One additional trend to recognize in table \ref{tab:groupAngle} is |
| 956 |
– |
that the $\sigma^2$ values for both {\sc sp} and |
| 957 |
– |
{\sc sf} converge as $\alpha$ increases, something that is easier |
| 958 |
– |
to see when using group based cutoffs. Looking back on figures |
| 959 |
– |
\ref{fig:delE}, \ref{fig:frcMag}, and \ref{fig:trqMag}, show this |
| 960 |
– |
behavior clearly at large $\alpha$ and cutoff values. The reason for |
| 961 |
– |
this is that the complimentary error function inserted into the |
| 962 |
– |
potential weakens the electrostatic interaction as $\alpha$ increases. |
| 963 |
– |
Thus, at larger values of $\alpha$, both the summation method types |
| 964 |
– |
progress toward non-interacting functions, so care is required in |
| 965 |
– |
choosing large damping functions lest one generate an undesirable loss |
| 966 |
– |
in the pair interaction. Kast \textit{et al.} developed a method for |
| 967 |
– |
choosing appropriate $\alpha$ values for these types of electrostatic |
| 968 |
– |
summation methods by fitting to $g(r)$ data, and their methods |
| 969 |
– |
indicate optimal values of 0.34, 0.25, and 0.16 \AA$^{-1}$ for cutoff |
| 970 |
– |
values of 9, 12, and 15 \AA\ respectively.\cite{Kast03} These appear |
| 971 |
– |
to be reasonable choices to obtain proper MC behavior |
| 972 |
– |
(Fig. \ref{fig:delE}); however, based on these findings, choices this |
| 973 |
– |
high would introduce error in the molecular torques, particularly for |
| 974 |
– |
the shorter cutoffs. Based on the above findings, empirical damping |
| 975 |
– |
up to 0.2 \AA$^{-1}$ proves to be beneficial, but damping is arguably |
| 976 |
– |
unnecessary when using the {\sc sf} method. |
| 977 |
– |
|
| 969 |
|
\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals} |
| 970 |
|
|
| 971 |
< |
In the previous studies using a {\sc sf} variant of the damped |
| 972 |
< |
Wolf coulomb potential, the structure and dynamics of water were |
| 973 |
< |
investigated rather extensively.\cite{Zahn02,Kast03} Their results |
| 974 |
< |
indicated that the damped {\sc sf} method results in properties |
| 975 |
< |
very similar to those obtained when using the Ewald summation. |
| 976 |
< |
Considering the statistical results shown above, the good performance |
| 977 |
< |
of this method is not that surprising. Rather than consider the same |
| 978 |
< |
systems and simply recapitulate their results, we decided to look at |
| 979 |
< |
the solid state dynamical behavior obtained using the best performing |
| 980 |
< |
summation methods from the above results. |
| 971 |
> |
Zahn {\it et al.} investigated the structure and dynamics of water |
| 972 |
> |
using eqs. (\ref{eq:ZahnPot}) and |
| 973 |
> |
(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated |
| 974 |
> |
that a method similar (but not identical with) the damped {\sc sf} |
| 975 |
> |
method resulted in properties very similar to those obtained when |
| 976 |
> |
using the Ewald summation. The properties they studied (pair |
| 977 |
> |
distribution functions, diffusion constants, and velocity and |
| 978 |
> |
orientational correlation functions) may not be particularly sensitive |
| 979 |
> |
to the long-range and collective behavior that governs the |
| 980 |
> |
low-frequency behavior in crystalline systems. Additionally, the |
| 981 |
> |
ionic crystals are the worst case scenario for the pairwise methods |
| 982 |
> |
because they lack the reciprocal space contribution contained in the |
| 983 |
> |
Ewald summation. |
| 984 |
|
|
| 985 |
+ |
We are using two separate measures to probe the effects of these |
| 986 |
+ |
alternative electrostatic methods on the dynamics in crystalline |
| 987 |
+ |
materials. For short- and intermediate-time dynamics, we are |
| 988 |
+ |
computing the velocity autocorrelation function, and for long-time |
| 989 |
+ |
and large length-scale collective motions, we are looking at the |
| 990 |
+ |
low-frequency portion of the power spectrum. |
| 991 |
+ |
|
| 992 |
|
\begin{figure} |
| 993 |
|
\centering |
| 994 |
|
\includegraphics[width = \linewidth]{./vCorrPlot.pdf} |
| 995 |
< |
\caption{Velocity auto-correlation functions of NaCl crystals at |
| 996 |
< |
1000 K while using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and |
| 997 |
< |
{\sc sp} ($\alpha$ = 0.2). The inset is a magnification of the first |
| 998 |
< |
trough. The times to first collision are nearly identical, but the |
| 999 |
< |
differences can be seen in the peaks and troughs, where the undamped |
| 1000 |
< |
to weakly damped methods are stiffer than the moderately damped and |
| 1001 |
< |
SPME methods.} |
| 995 |
> |
\caption{Velocity autocorrelation functions of NaCl crystals at |
| 996 |
> |
1000 K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc |
| 997 |
> |
sp} ($\alpha$ = 0.2). The inset is a magnification of the area around |
| 998 |
> |
the first minimum. The times to first collision are nearly identical, |
| 999 |
> |
but differences can be seen in the peaks and troughs, where the |
| 1000 |
> |
undamped and weakly damped methods are stiffer than the moderately |
| 1001 |
> |
damped and {\sc spme} methods.} |
| 1002 |
|
\label{fig:vCorrPlot} |
| 1003 |
|
\end{figure} |
| 1004 |
|
|
| 1005 |
< |
The short-time decays through the first collision are nearly identical |
| 1006 |
< |
in figure \ref{fig:vCorrPlot}, but the peaks and troughs of the |
| 1007 |
< |
functions show how the methods differ. The undamped {\sc sf} method |
| 1008 |
< |
has deeper troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher |
| 1009 |
< |
peaks than any of the other methods. As the damping function is |
| 1010 |
< |
increased, these peaks are smoothed out, and approach the SPME |
| 1011 |
< |
curve. The damping acts as a distance dependent Gaussian screening of |
| 1012 |
< |
the point charges for the pairwise summation methods; thus, the |
| 1013 |
< |
collisions are more elastic in the undamped {\sc sf} potential, and the |
| 1014 |
< |
stiffness of the potential is diminished as the electrostatic |
| 1015 |
< |
interactions are softened by the damping function. With $\alpha$ |
| 1016 |
< |
values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are |
| 1017 |
< |
nearly identical and track the SPME features quite well. This is not |
| 1018 |
< |
too surprising in that the differences between the {\sc sf} and {\sc |
| 1019 |
< |
sp} potentials are mitigated with increased damping. However, this |
| 1019 |
< |
appears to indicate that once damping is utilized, the form of the |
| 1020 |
< |
potential seems to play a lesser role in the crystal dynamics. |
| 1005 |
> |
The short-time decay of the velocity autocorrelation function through |
| 1006 |
> |
the first collision are nearly identical in figure |
| 1007 |
> |
\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show |
| 1008 |
> |
how the methods differ. The undamped {\sc sf} method has deeper |
| 1009 |
> |
troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than |
| 1010 |
> |
any of the other methods. As the damping parameter ($\alpha$) is |
| 1011 |
> |
increased, these peaks are smoothed out, and the {\sc sf} method |
| 1012 |
> |
approaches the {\sc spme} results. With $\alpha$ values of 0.2 \AA$^{-1}$, |
| 1013 |
> |
the {\sc sf} and {\sc sp} functions are nearly identical and track the |
| 1014 |
> |
{\sc spme} features quite well. This is not surprising because the {\sc sf} |
| 1015 |
> |
and {\sc sp} potentials become nearly identical with increased |
| 1016 |
> |
damping. However, this appears to indicate that once damping is |
| 1017 |
> |
utilized, the details of the form of the potential (and forces) |
| 1018 |
> |
constructed out of the damped electrostatic interaction are less |
| 1019 |
> |
important. |
| 1020 |
|
|
| 1021 |
|
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 1022 |
|
|
| 1023 |
< |
The short time dynamics were extended to evaluate how the differences |
| 1024 |
< |
between the methods affect the collective long-time motion. The same |
| 1025 |
< |
electrostatic summation methods were used as in the short time |
| 1026 |
< |
velocity autocorrelation function evaluation, but the trajectories |
| 1027 |
< |
were sampled over a much longer time. The power spectra of the |
| 1028 |
< |
resulting velocity autocorrelation functions were calculated and are |
| 1029 |
< |
displayed in figure \ref{fig:methodPS}. |
| 1023 |
> |
To evaluate how the differences between the methods affect the |
| 1024 |
> |
collective long-time motion, we computed power spectra from long-time |
| 1025 |
> |
traces of the velocity autocorrelation function. The power spectra for |
| 1026 |
> |
the best-performing alternative methods are shown in |
| 1027 |
> |
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
| 1028 |
> |
a cubic switching function between 40 and 50 ps was used to reduce the |
| 1029 |
> |
ringing resulting from data truncation. This procedure had no |
| 1030 |
> |
noticeable effect on peak location or magnitude. |
| 1031 |
|
|
| 1032 |
|
\begin{figure} |
| 1033 |
|
\centering |
| 1034 |
|
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 1035 |
|
\caption{Power spectra obtained from the velocity auto-correlation |
| 1036 |
< |
functions of NaCl crystals at 1000 K while using SPME, {\sc sf} |
| 1037 |
< |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). |
| 1038 |
< |
Apodization of the correlation functions via a cubic switching |
| 1039 |
< |
function between 40 and 50 ps was used to clear up the spectral noise |
| 1040 |
< |
resulting from data truncation, and had no noticeable effect on peak |
| 1041 |
< |
location or magnitude. The inset shows the frequency region below 100 |
| 1042 |
< |
cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 1036 |
> |
functions of NaCl crystals at 1000 K while using {\sc spme}, {\sc sf} |
| 1037 |
> |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset |
| 1038 |
> |
shows the frequency region below 100 cm$^{-1}$ to highlight where the |
| 1039 |
> |
spectra differ.} |
| 1040 |
|
\label{fig:methodPS} |
| 1041 |
|
\end{figure} |
| 1042 |
|
|
| 1043 |
< |
While high frequency peaks of the spectra in this figure overlap, |
| 1044 |
< |
showing the same general features, the low frequency region shows how |
| 1045 |
< |
the summation methods differ. Considering the low-frequency inset |
| 1046 |
< |
(expanded in the upper frame of figure \ref{fig:dampInc}), at |
| 1047 |
< |
frequencies below 100 cm$^{-1}$, the correlated motions are |
| 1048 |
< |
blue-shifted when using undamped or weakly damped {\sc sf}. When |
| 1049 |
< |
using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the {\sc sf} |
| 1050 |
< |
and {\sc sp} methods give near identical correlated motion behavior as |
| 1051 |
< |
the Ewald method (which has a damping value of 0.3119). This |
| 1052 |
< |
weakening of the electrostatic interaction with increased damping |
| 1053 |
< |
explains why the long-ranged correlated motions are at lower |
| 1054 |
< |
frequencies for the moderately damped methods than for undamped or |
| 1055 |
< |
weakly damped methods. To see this effect more clearly, we show how |
| 1056 |
< |
damping strength alone affects a simple real-space electrostatic |
| 1057 |
< |
potential, |
| 1058 |
< |
\begin{equation} |
| 1059 |
< |
V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r), |
| 1060 |
< |
\end{equation} |
| 1061 |
< |
where $S(r)$ is a switching function that smoothly zeroes the |
| 1062 |
< |
potential at the cutoff radius. Figure \ref{fig:dampInc} shows how |
| 1063 |
< |
the low frequency motions are dependent on the damping used in the |
| 1064 |
< |
direct electrostatic sum. As the damping increases, the peaks drop to |
| 1065 |
< |
lower frequencies. Incidentally, use of an $\alpha$ of 0.25 |
| 1066 |
< |
\AA$^{-1}$ on a simple electrostatic summation results in low |
| 1067 |
< |
frequency correlated dynamics equivalent to a simulation using SPME. |
| 1068 |
< |
When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks |
| 1069 |
< |
shift to higher frequency in exponential fashion. Though not shown, |
| 1070 |
< |
the spectrum for the simple undamped electrostatic potential is |
| 1074 |
< |
blue-shifted such that the lowest frequency peak resides near 325 |
| 1075 |
< |
cm$^{-1}$. In light of these results, the undamped {\sc sf} method |
| 1076 |
< |
producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite |
| 1077 |
< |
respectable and shows that the shifted force procedure accounts for |
| 1078 |
< |
most of the effect afforded through use of the Ewald summation. |
| 1079 |
< |
However, it appears as though moderate damping is required for |
| 1080 |
< |
accurate reproduction of crystal dynamics. |
| 1043 |
> |
While the high frequency regions of the power spectra for the |
| 1044 |
> |
alternative methods are quantitatively identical with Ewald spectrum, |
| 1045 |
> |
the low frequency region shows how the summation methods differ. |
| 1046 |
> |
Considering the low-frequency inset (expanded in the upper frame of |
| 1047 |
> |
figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the |
| 1048 |
> |
correlated motions are blue-shifted when using undamped or weakly |
| 1049 |
> |
damped {\sc sf}. When using moderate damping ($\alpha = 0.2$ |
| 1050 |
> |
\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical |
| 1051 |
> |
correlated motion to the Ewald method (which has a convergence |
| 1052 |
> |
parameter of 0.3119 \AA$^{-1}$). This weakening of the electrostatic |
| 1053 |
> |
interaction with increased damping explains why the long-ranged |
| 1054 |
> |
correlated motions are at lower frequencies for the moderately damped |
| 1055 |
> |
methods than for undamped or weakly damped methods. |
| 1056 |
> |
|
| 1057 |
> |
To isolate the role of the damping constant, we have computed the |
| 1058 |
> |
spectra for a single method ({\sc sf}) with a range of damping |
| 1059 |
> |
constants and compared this with the {\sc spme} spectrum. |
| 1060 |
> |
Fig. \ref{fig:dampInc} shows more clearly that increasing the |
| 1061 |
> |
electrostatic damping red-shifts the lowest frequency phonon modes. |
| 1062 |
> |
However, even without any electrostatic damping, the {\sc sf} method |
| 1063 |
> |
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
| 1064 |
> |
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
| 1065 |
> |
would predict the lowest frequency peak near 325 cm$^{-1}$. {\it |
| 1066 |
> |
Most} of the collective behavior in the crystal is accurately captured |
| 1067 |
> |
using the {\sc sf} method. Quantitative agreement with Ewald can be |
| 1068 |
> |
obtained using moderate damping in addition to the shifting at the |
| 1069 |
> |
cutoff distance. |
| 1070 |
> |
|
| 1071 |
|
\begin{figure} |
| 1072 |
|
\centering |
| 1073 |
< |
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 1074 |
< |
\caption{Regions of spectra showing the low-frequency correlated |
| 1075 |
< |
motions for NaCl crystals at 1000 K using various electrostatic |
| 1076 |
< |
summation methods. The upper plot is a zoomed inset from figure |
| 1077 |
< |
\ref{fig:methodPS}. As the damping value for the {\sc sf} potential |
| 1078 |
< |
increases, the low-frequency peaks red-shift. The lower plot is of |
| 1079 |
< |
spectra when using SPME and a simple damped Coulombic sum with damping |
| 1080 |
< |
coefficients ($\alpha$) ranging from 0.15 to 0.3 \AA$^{-1}$. As |
| 1091 |
< |
$\alpha$ increases, the peaks are red-shifted toward and eventually |
| 1092 |
< |
beyond the values given by SPME. The larger $\alpha$ values weaken |
| 1093 |
< |
the real-space electrostatics, explaining this shift towards less |
| 1094 |
< |
strongly correlated motions in the crystal.} |
| 1073 |
> |
\includegraphics[width = \linewidth]{./increasedDamping.pdf} |
| 1074 |
> |
\caption{Effect of damping on the two lowest-frequency phonon modes in |
| 1075 |
> |
the NaCl crystal at 1000~K. The undamped shifted force ({\sc sf}) |
| 1076 |
> |
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1077 |
> |
electrostatic damping to 0.25 \AA$^{-1}$ gives quantitative agreement |
| 1078 |
> |
with the power spectrum obtained using the Ewald sum. Overdamping can |
| 1079 |
> |
result in underestimates of frequencies of the long-wavelength |
| 1080 |
> |
motions.} |
| 1081 |
|
\label{fig:dampInc} |
| 1082 |
|
\end{figure} |
| 1083 |
|
|
| 1084 |
|
\section{Conclusions} |
| 1085 |
|
|
| 1086 |
|
This investigation of pairwise electrostatic summation techniques |
| 1087 |
< |
shows that there are viable and more computationally efficient |
| 1088 |
< |
electrostatic summation techniques than the Ewald summation, chiefly |
| 1089 |
< |
methods derived from the damped Coulombic sum originally proposed by |
| 1090 |
< |
Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the |
| 1091 |
< |
{\sc sf} method, reformulated above as eq. (\ref{eq:DSFPot}), |
| 1092 |
< |
shows a remarkable ability to reproduce the energetic and dynamic |
| 1093 |
< |
characteristics exhibited by simulations employing lattice summation |
| 1094 |
< |
techniques. The cumulative energy difference results showed the |
| 1095 |
< |
undamped {\sc sf} and moderately damped {\sc sp} methods |
| 1096 |
< |
produced results nearly identical to SPME. Similarly for the dynamic |
| 1097 |
< |
features, the undamped or moderately damped {\sc sf} and |
| 1098 |
< |
moderately damped {\sc sp} methods produce force and torque |
| 1099 |
< |
vector magnitude and directions very similar to the expected values. |
| 1100 |
< |
These results translate into long-time dynamic behavior equivalent to |
| 1101 |
< |
that produced in simulations using SPME. |
| 1087 |
> |
shows that there are viable and computationally efficient alternatives |
| 1088 |
> |
to the Ewald summation. These methods are derived from the damped and |
| 1089 |
> |
cutoff-neutralized Coulombic sum originally proposed by Wolf |
| 1090 |
> |
\textit{et al.}\cite{Wolf99} In particular, the {\sc sf} |
| 1091 |
> |
method, reformulated above as eqs. (\ref{eq:DSFPot}) and |
| 1092 |
> |
(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the |
| 1093 |
> |
energetic and dynamic characteristics exhibited by simulations |
| 1094 |
> |
employing lattice summation techniques. The cumulative energy |
| 1095 |
> |
difference results showed the undamped {\sc sf} and moderately damped |
| 1096 |
> |
{\sc sp} methods produced results nearly identical to {\sc spme}. Similarly |
| 1097 |
> |
for the dynamic features, the undamped or moderately damped {\sc sf} |
| 1098 |
> |
and moderately damped {\sc sp} methods produce force and torque vector |
| 1099 |
> |
magnitude and directions very similar to the expected values. These |
| 1100 |
> |
results translate into long-time dynamic behavior equivalent to that |
| 1101 |
> |
produced in simulations using {\sc spme}. |
| 1102 |
|
|
| 1103 |
+ |
As in all purely-pairwise cutoff methods, these methods are expected |
| 1104 |
+ |
to scale approximately {\it linearly} with system size, and they are |
| 1105 |
+ |
easily parallelizable. This should result in substantial reductions |
| 1106 |
+ |
in the computational cost of performing large simulations. |
| 1107 |
+ |
|
| 1108 |
|
Aside from the computational cost benefit, these techniques have |
| 1109 |
|
applicability in situations where the use of the Ewald sum can prove |
| 1110 |
< |
problematic. Primary among them is their use in interfacial systems, |
| 1111 |
< |
where the unmodified lattice sum techniques artificially accentuate |
| 1112 |
< |
the periodicity of the system in an undesirable manner. There have |
| 1113 |
< |
been alterations to the standard Ewald techniques, via corrections and |
| 1114 |
< |
reformulations, to compensate for these systems; but the pairwise |
| 1115 |
< |
techniques discussed here require no modifications, making them |
| 1116 |
< |
natural tools to tackle these problems. Additionally, this |
| 1117 |
< |
transferability gives them benefits over other pairwise methods, like |
| 1118 |
< |
reaction field, because estimations of physical properties (e.g. the |
| 1119 |
< |
dielectric constant) are unnecessary. |
| 1110 |
> |
problematic. Of greatest interest is their potential use in |
| 1111 |
> |
interfacial systems, where the unmodified lattice sum techniques |
| 1112 |
> |
artificially accentuate the periodicity of the system in an |
| 1113 |
> |
undesirable manner. There have been alterations to the standard Ewald |
| 1114 |
> |
techniques, via corrections and reformulations, to compensate for |
| 1115 |
> |
these systems; but the pairwise techniques discussed here require no |
| 1116 |
> |
modifications, making them natural tools to tackle these problems. |
| 1117 |
> |
Additionally, this transferability gives them benefits over other |
| 1118 |
> |
pairwise methods, like reaction field, because estimations of physical |
| 1119 |
> |
properties (e.g. the dielectric constant) are unnecessary. |
| 1120 |
|
|
| 1121 |
< |
We are not suggesting any flaw with the Ewald sum; in fact, it is the |
| 1122 |
< |
standard by which these simple pairwise sums are judged. However, |
| 1123 |
< |
these results do suggest that in the typical simulations performed |
| 1124 |
< |
today, the Ewald summation may no longer be required to obtain the |
| 1125 |
< |
level of accuracy most researchers have come to expect |
| 1121 |
> |
If a researcher is using Monte Carlo simulations of large chemical |
| 1122 |
> |
systems containing point charges, most structural features will be |
| 1123 |
> |
accurately captured using the undamped {\sc sf} method or the {\sc sp} |
| 1124 |
> |
method with an electrostatic damping of 0.2 \AA$^{-1}$. These methods |
| 1125 |
> |
would also be appropriate for molecular dynamics simulations where the |
| 1126 |
> |
data of interest is either structural or short-time dynamical |
| 1127 |
> |
quantities. For long-time dynamics and collective motions, the safest |
| 1128 |
> |
pairwise method we have evaluated is the {\sc sf} method with an |
| 1129 |
> |
electrostatic damping between 0.2 and 0.25 |
| 1130 |
> |
\AA$^{-1}$. |
| 1131 |
|
|
| 1132 |
+ |
We are not suggesting that there is any flaw with the Ewald sum; in |
| 1133 |
+ |
fact, it is the standard by which these simple pairwise sums have been |
| 1134 |
+ |
judged. However, these results do suggest that in the typical |
| 1135 |
+ |
simulations performed today, the Ewald summation may no longer be |
| 1136 |
+ |
required to obtain the level of accuracy most researchers have come to |
| 1137 |
+ |
expect. |
| 1138 |
+ |
|
| 1139 |
|
\section{Acknowledgments} |
| 1140 |
+ |
Support for this project was provided by the National Science |
| 1141 |
+ |
Foundation under grant CHE-0134881. The authors would like to thank |
| 1142 |
+ |
Steve Corcelli and Ed Maginn for helpful discussions and comments. |
| 1143 |
+ |
|
| 1144 |
|
\newpage |
| 1145 |
|
|
| 1146 |
|
\bibliographystyle{jcp2} |
