| 20 |
|
\topmargin -21pt \headsep 10pt |
| 21 |
|
\textheight 9.0in \textwidth 6.5in |
| 22 |
|
\brokenpenalty=10000 |
| 23 |
< |
\renewcommand{\baselinestretch}{1.2} |
| 23 |
> |
%\renewcommand{\baselinestretch}{1.2} |
| 24 |
> |
\renewcommand{\baselinestretch}{2} |
| 25 |
|
\renewcommand\citemid{\ } % no comma in optional reference note |
| 26 |
+ |
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{2}} %doublespace captions |
| 27 |
+ |
\let\Caption\caption |
| 28 |
+ |
\renewcommand\caption[1]{% |
| 29 |
+ |
\Caption[#1]{}% |
| 30 |
+ |
} |
| 31 |
|
|
| 32 |
+ |
|
| 33 |
|
\begin{document} |
| 34 |
|
|
| 35 |
< |
\title{Is the Ewald Summation necessary? Pairwise alternatives to the accepted standard for long-range electrostatics} |
| 35 |
> |
\title{Is the Ewald summation still necessary? \\ |
| 36 |
> |
Pairwise alternatives to the accepted standard for |
| 37 |
> |
long-range electrostatics in molecular simulations} |
| 38 |
|
|
| 39 |
|
\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
| 40 |
|
gezelter@nd.edu} \\ |
| 45 |
|
\date{\today} |
| 46 |
|
|
| 47 |
|
\maketitle |
| 48 |
< |
\doublespacing |
| 48 |
> |
%\doublespacing |
| 49 |
|
|
| 41 |
– |
\nobibliography{} |
| 50 |
|
\begin{abstract} |
| 51 |
< |
A new method for accumulating electrostatic interactions was derived |
| 52 |
< |
from the previous efforts described in \bibentry{Wolf99} and |
| 53 |
< |
\bibentry{Zahn02} as a possible replacement for lattice sum methods in |
| 54 |
< |
molecular simulations. Comparisons were performed with this and other |
| 55 |
< |
pairwise electrostatic summation techniques against the smooth |
| 56 |
< |
particle mesh Ewald (SPME) summation to see how well they reproduce |
| 57 |
< |
the energetics and dynamics of a variety of simulation types. The |
| 58 |
< |
newly derived Shifted-Force technique shows a remarkable ability to |
| 59 |
< |
reproduce the behavior exhibited in simulations using SPME with an |
| 60 |
< |
$\mathscr{O}(N)$ computational cost, equivalent to merely the |
| 61 |
< |
real-space portion of the lattice summation. |
| 62 |
< |
|
| 51 |
> |
We investigate pairwise electrostatic interaction methods and show |
| 52 |
> |
that there are viable and computationally efficient $(\mathscr{O}(N))$ |
| 53 |
> |
alternatives to the Ewald summation for typical modern molecular |
| 54 |
> |
simulations. These methods are extended from the damped and |
| 55 |
> |
cutoff-neutralized Coulombic sum originally proposed by |
| 56 |
> |
[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 |
| 57 |
> |
a remarkable ability to reproduce the energetic and dynamic |
| 58 |
> |
characteristics exhibited by simulations employing lattice summation |
| 59 |
> |
techniques. Comparisons were performed with this and other pairwise |
| 60 |
> |
methods against the smooth particle mesh Ewald ({\sc spme}) summation |
| 61 |
> |
to see how well they reproduce the energetics and dynamics of a |
| 62 |
> |
variety of simulation types. |
| 63 |
|
\end{abstract} |
| 64 |
|
|
| 65 |
|
\newpage |
| 102 |
|
regarding possible artifacts caused by the inherent periodicity of the |
| 103 |
|
explicit Ewald summation.\cite{Tobias01} |
| 104 |
|
|
| 105 |
< |
In this paper, we focus on a new set of shifted methods devised by |
| 105 |
> |
In this paper, we focus on a new set of pairwise methods devised by |
| 106 |
|
Wolf {\it et al.},\cite{Wolf99} which we further extend. These |
| 107 |
|
methods along with a few other mixed methods (i.e. reaction field) are |
| 108 |
|
compared with the smooth particle mesh Ewald |
| 113 |
|
to the direct pairwise sum. They also lack the added periodicity of |
| 114 |
|
the Ewald sum, so they can be used for systems which are non-periodic |
| 115 |
|
or which have one- or two-dimensional periodicity. Below, these |
| 116 |
< |
methods are evaluated using a variety of model systems to establish |
| 117 |
< |
their usability in molecular simulations. |
| 116 |
> |
methods are evaluated using a variety of model systems to |
| 117 |
> |
establish their usability in molecular simulations. |
| 118 |
|
|
| 119 |
|
\subsection{The Ewald Sum} |
| 120 |
< |
The complete accumulation electrostatic interactions in a system with |
| 120 |
> |
The complete accumulation of the electrostatic interactions in a system with |
| 121 |
|
periodic boundary conditions (PBC) requires the consideration of the |
| 122 |
|
effect of all charges within a (cubic) simulation box as well as those |
| 123 |
|
in the periodic replicas, |
| 173 |
|
|
| 174 |
|
\begin{figure} |
| 175 |
|
\centering |
| 176 |
< |
\includegraphics[width = \linewidth]{./ewaldProgression2.pdf} |
| 177 |
< |
\caption{The change in the application of the Ewald sum with |
| 178 |
< |
increasing computational power. Initially, only small systems could |
| 179 |
< |
be studied, and the Ewald sum replicated the simulation box to |
| 180 |
< |
convergence. Now, much larger systems of charges are investigated |
| 181 |
< |
with fixed-distance cutoffs.} |
| 176 |
> |
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 177 |
> |
\caption{The change in the need for the Ewald sum with |
| 178 |
> |
increasing computational power. A:~Initially, only small systems |
| 179 |
> |
could be studied, and the Ewald sum replicated the simulation box to |
| 180 |
> |
convergence. B:~Now, radial cutoff methods should be able to reach |
| 181 |
> |
convergence for the larger systems of charges that are common today.} |
| 182 |
|
\label{fig:ewaldTime} |
| 183 |
|
\end{figure} |
| 184 |
|
|
| 208 |
|
interfaces and membranes, the intrinsic three-dimensional periodicity |
| 209 |
|
can prove problematic. The Ewald sum has been reformulated to handle |
| 210 |
|
2D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89}, but the |
| 211 |
< |
new methods are computationally expensive.\cite{Spohr97,Yeh99} |
| 212 |
< |
Inclusion of a correction term in the Ewald summation is a possible |
| 213 |
< |
direction for handling 2D systems while still enabling the use of the |
| 214 |
< |
modern optimizations.\cite{Yeh99} |
| 211 |
> |
new methods are computationally expensive.\cite{Spohr97,Yeh99} More |
| 212 |
> |
recently, there have been several successful efforts toward reducing |
| 213 |
> |
the computational cost of 2D lattice summations, often enabling the |
| 214 |
> |
use of the mentioned |
| 215 |
> |
optimizations.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} |
| 216 |
|
|
| 217 |
|
Several studies have recognized that the inherent periodicity in the |
| 218 |
|
Ewald sum can also have an effect on three-dimensional |
| 235 |
|
charge contained within the cutoff radius is crucial for potential |
| 236 |
|
stability. They devised a pairwise summation method that ensures |
| 237 |
|
charge neutrality and gives results similar to those obtained with the |
| 238 |
< |
Ewald summation. The resulting shifted Coulomb potential |
| 239 |
< |
(Eq. \ref{eq:WolfPot}) includes image-charges subtracted out through |
| 240 |
< |
placement on the cutoff sphere and a distance-dependent damping |
| 241 |
< |
function (identical to that seen in the real-space portion of the |
| 233 |
< |
Ewald sum) to aid convergence |
| 238 |
> |
Ewald summation. The resulting shifted Coulomb potential includes |
| 239 |
> |
image-charges subtracted out through placement on the cutoff sphere |
| 240 |
> |
and a distance-dependent damping function (identical to that seen in |
| 241 |
> |
the real-space portion of the Ewald sum) to aid convergence |
| 242 |
|
\begin{equation} |
| 243 |
|
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\}. |
| 244 |
|
\label{eq:WolfPot} |
| 524 |
|
The pairwise summation techniques (outlined in section |
| 525 |
|
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 526 |
|
studying the energy differences between conformations. We took the |
| 527 |
< |
SPME-computed energy difference between two conformations to be the |
| 527 |
> |
{\sc spme}-computed energy difference between two conformations to be the |
| 528 |
|
correct behavior. An ideal performance by an alternative method would |
| 529 |
|
reproduce these energy differences exactly (even if the absolute |
| 530 |
|
energies calculated by the methods are different). Since none of the |
| 532 |
|
regressions of energy gap data to evaluate how closely the methods |
| 533 |
|
mimicked the Ewald energy gaps. Unitary results for both the |
| 534 |
|
correlation (slope) and correlation coefficient for these regressions |
| 535 |
< |
indicate perfect agreement between the alternative method and SPME. |
| 535 |
> |
indicate perfect agreement between the alternative method and {\sc spme}. |
| 536 |
|
Sample correlation plots for two alternate methods are shown in |
| 537 |
|
Fig. \ref{fig:linearFit}. |
| 538 |
|
|
| 546 |
|
\label{fig:linearFit} |
| 547 |
|
\end{figure} |
| 548 |
|
|
| 549 |
< |
Each system type (detailed in section \ref{sec:RepSims}) was |
| 550 |
< |
represented using 500 independent configurations. Additionally, we |
| 551 |
< |
used seven different system types, so each of the alternative |
| 552 |
< |
(non-Ewald) electrostatic summation methods was evaluated using |
| 553 |
< |
873,250 configurational energy differences. |
| 549 |
> |
Each of the seven system types (detailed in section \ref{sec:RepSims}) |
| 550 |
> |
were represented using 500 independent configurations. Thus, each of |
| 551 |
> |
the alternative (non-Ewald) electrostatic summation methods was |
| 552 |
> |
evaluated using an accumulated 873,250 configurational energy |
| 553 |
> |
differences. |
| 554 |
|
|
| 555 |
|
Results and discussion for the individual analysis of each of the |
| 556 |
< |
system types appear in the supporting information, while the |
| 557 |
< |
cumulative results over all the investigated systems appears below in |
| 558 |
< |
section \ref{sec:EnergyResults}. |
| 556 |
> |
system types appear in the supporting information,\cite{EPAPSdeposit} |
| 557 |
> |
while the cumulative results over all the investigated systems appears |
| 558 |
> |
below in section \ref{sec:EnergyResults}. |
| 559 |
|
|
| 560 |
|
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
| 561 |
|
We evaluated the pairwise methods (outlined in section |
| 562 |
|
\ref{sec:ESMethods}) for use in MD simulations by |
| 563 |
|
comparing the force and torque vectors with those obtained using the |
| 564 |
< |
reference Ewald summation (SPME). Both the magnitude and the |
| 564 |
> |
reference Ewald summation ({\sc spme}). Both the magnitude and the |
| 565 |
|
direction of these vectors on each of the bodies in the system were |
| 566 |
|
analyzed. For the magnitude of these vectors, linear least squares |
| 567 |
|
regression analyses were performed as described previously for |
| 576 |
|
|
| 577 |
|
The {\it directionality} of the force and torque vectors was |
| 578 |
|
investigated through measurement of the angle ($\theta$) formed |
| 579 |
< |
between those computed from the particular method and those from SPME, |
| 579 |
> |
between those computed from the particular method and those from {\sc spme}, |
| 580 |
|
\begin{equation} |
| 581 |
|
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right), |
| 582 |
|
\end{equation} |
| 583 |
< |
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force |
| 583 |
> |
where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force |
| 584 |
|
vector computed using method M. Each of these $\theta$ values was |
| 585 |
|
accumulated in a distribution function and weighted by the area on the |
| 586 |
|
unit sphere. Since this distribution is a measure of angular error |
| 587 |
|
between two different electrostatic summation methods, there is no |
| 588 |
|
{\it a priori} reason for the profile to adhere to any specific |
| 589 |
|
shape. Thus, gaussian fits were used to measure the width of the |
| 590 |
< |
resulting distributions. |
| 591 |
< |
% |
| 592 |
< |
%\begin{figure} |
| 593 |
< |
%\centering |
| 594 |
< |
%\includegraphics[width = \linewidth]{./gaussFit.pdf} |
| 587 |
< |
%\caption{Sample fit of the angular distribution of the force vectors |
| 588 |
< |
%accumulated using all of the studied systems. Gaussian fits were used |
| 589 |
< |
%to obtain values for the variance in force and torque vectors.} |
| 590 |
< |
%\label{fig:gaussian} |
| 591 |
< |
%\end{figure} |
| 592 |
< |
% |
| 593 |
< |
%Figure \ref{fig:gaussian} shows an example distribution with applied |
| 594 |
< |
%non-linear fits. The solid line is a Gaussian profile, while the |
| 595 |
< |
%dotted line is a Voigt profile, a convolution of a Gaussian and a |
| 596 |
< |
%Lorentzian. |
| 597 |
< |
%Since this distribution is a measure of angular error between two |
| 598 |
< |
%different electrostatic summation methods, there is no {\it a priori} |
| 599 |
< |
%reason for the profile to adhere to any specific shape. |
| 600 |
< |
%Gaussian fits was used to compare all the tested methods. |
| 601 |
< |
The variance ($\sigma^2$) was extracted from each of these fits and |
| 602 |
< |
was used to compare distribution widths. Values of $\sigma^2$ near |
| 603 |
< |
zero indicate vector directions indistinguishable from those |
| 604 |
< |
calculated when using the reference method (SPME). |
| 590 |
> |
resulting distributions. The variance ($\sigma^2$) was extracted from |
| 591 |
> |
each of these fits and was used to compare distribution widths. |
| 592 |
> |
Values of $\sigma^2$ near zero indicate vector directions |
| 593 |
> |
indistinguishable from those calculated when using the reference |
| 594 |
> |
method ({\sc spme}). |
| 595 |
|
|
| 596 |
|
\subsection{Short-time Dynamics} |
| 597 |
|
|
| 606 |
|
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
| 607 |
|
of the trajectories, |
| 608 |
|
\begin{equation} |
| 609 |
< |
C_v(t) = \frac{\langle v_i(0)\cdot v_i(t)\rangle}{\langle v_i(0)\cdot v_i(0)\rangle}. |
| 609 |
> |
C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}. |
| 610 |
|
\label{eq:vCorr} |
| 611 |
|
\end{equation} |
| 612 |
|
Velocity autocorrelation functions require detailed short time data, |
| 619 |
|
|
| 620 |
|
The effects of the same subset of alternative electrostatic methods on |
| 621 |
|
the {\it long-time} dynamics of charged systems were evaluated using |
| 622 |
< |
the same model system (NaCl crystals at 1000K). The power spectrum |
| 622 |
> |
the same model system (NaCl crystals at 1000~K). The power spectrum |
| 623 |
|
($I(\omega)$) was obtained via Fourier transform of the velocity |
| 624 |
|
autocorrelation function, \begin{equation} I(\omega) = |
| 625 |
|
\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
| 629 |
|
NaCl crystal is composed of two different atom types, the average of |
| 630 |
|
the two resulting power spectra was used for comparisons. Simulations |
| 631 |
|
were performed under the microcanonical ensemble, and velocity |
| 632 |
< |
information was saved every 5 fs over 100 ps trajectories. |
| 632 |
> |
information was saved every 5~fs over 100~ps trajectories. |
| 633 |
|
|
| 634 |
|
\subsection{Representative Simulations}\label{sec:RepSims} |
| 635 |
< |
A variety of representative simulations were analyzed to determine the |
| 636 |
< |
relative effectiveness of the pairwise summation techniques in |
| 637 |
< |
reproducing the energetics and dynamics exhibited by SPME. We wanted |
| 638 |
< |
to span the space of modern simulations (i.e. from liquids of neutral |
| 639 |
< |
molecules to ionic crystals), so the systems studied were: |
| 635 |
> |
A variety of representative molecular simulations were analyzed to |
| 636 |
> |
determine the relative effectiveness of the pairwise summation |
| 637 |
> |
techniques in reproducing the energetics and dynamics exhibited by |
| 638 |
> |
{\sc spme}. We wanted to span the space of typical molecular |
| 639 |
> |
simulations (i.e. from liquids of neutral molecules to ionic |
| 640 |
> |
crystals), so the systems studied were: |
| 641 |
|
\begin{enumerate} |
| 642 |
|
\item liquid water (SPC/E),\cite{Berendsen87} |
| 643 |
|
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), |
| 660 |
|
the crystal). The solid and liquid NaCl systems consisted of 500 |
| 661 |
|
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
| 662 |
|
these systems were selected and equilibrated in the same manner as the |
| 663 |
< |
water systems. The equilibrated temperatures were 1000~K for the NaCl |
| 664 |
< |
crystal and 7000~K for the liquid. The ionic solutions were made by |
| 665 |
< |
solvating 4 (or 40) ions in a periodic box containing 1000 SPC/E water |
| 666 |
< |
molecules. Ion and water positions were then randomly swapped, and |
| 667 |
< |
the resulting configurations were again equilibrated individually. |
| 668 |
< |
Finally, for the Argon / Water ``charge void'' systems, the identities |
| 669 |
< |
of all the SPC/E waters within 6 \AA\ of the center of the |
| 670 |
< |
equilibrated water configurations were converted to argon. |
| 671 |
< |
%(Fig. \ref{fig:argonSlice}). |
| 663 |
> |
water systems. In order to introduce measurable fluctuations in the |
| 664 |
> |
configuration energy differences, the crystalline simulations were |
| 665 |
> |
equilibrated at 1000~K, near the $T_\textrm{m}$ for NaCl. The liquid |
| 666 |
> |
NaCl configurations needed to represent a fully disordered array of |
| 667 |
> |
point charges, so the high temperature of 7000~K was selected for |
| 668 |
> |
equilibration. The ionic solutions were made by solvating 4 (or 40) |
| 669 |
> |
ions in a periodic box containing 1000 SPC/E water molecules. Ion and |
| 670 |
> |
water positions were then randomly swapped, and the resulting |
| 671 |
> |
configurations were again equilibrated individually. Finally, for the |
| 672 |
> |
Argon / Water ``charge void'' systems, the identities of all the SPC/E |
| 673 |
> |
waters within 6 \AA\ of the center of the equilibrated water |
| 674 |
> |
configurations were converted to argon. |
| 675 |
|
|
| 676 |
|
These procedures guaranteed us a set of representative configurations |
| 677 |
|
from chemically-relevant systems sampled from appropriate |
| 678 |
|
ensembles. Force field parameters for the ions and Argon were taken |
| 679 |
|
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
| 680 |
|
|
| 687 |
– |
%\begin{figure} |
| 688 |
– |
%\centering |
| 689 |
– |
%\includegraphics[width = \linewidth]{./slice.pdf} |
| 690 |
– |
%\caption{A slice from the center of a water box used in a charge void |
| 691 |
– |
%simulation. The darkened region represents the boundary sphere within |
| 692 |
– |
%which the water molecules were converted to argon atoms.} |
| 693 |
– |
%\label{fig:argonSlice} |
| 694 |
– |
%\end{figure} |
| 695 |
– |
|
| 681 |
|
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 682 |
|
We compared the following alternative summation methods with results |
| 683 |
< |
from the reference method (SPME): |
| 683 |
> |
from the reference method ({\sc spme}): |
| 684 |
|
\begin{itemize} |
| 685 |
|
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 686 |
|
and 0.3 \AA$^{-1}$, |
| 691 |
|
\end{itemize} |
| 692 |
|
Group-based cutoffs with a fifth-order polynomial switching function |
| 693 |
|
were utilized for the reaction field simulations. Additionally, we |
| 694 |
< |
investigated the use of these cutoffs with the SP, SF, and pure |
| 695 |
< |
cutoff. The SPME electrostatics were performed using the TINKER |
| 696 |
< |
implementation of SPME,\cite{Ponder87} while all other calculations |
| 694 |
> |
investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure |
| 695 |
> |
cutoff. The {\sc spme} electrostatics were performed using the {\sc tinker} |
| 696 |
> |
implementation of {\sc spme},\cite{Ponder87} while all other calculations |
| 697 |
|
were performed using the {\sc oopse} molecular mechanics |
| 698 |
|
package.\cite{Meineke05} All other portions of the energy calculation |
| 699 |
|
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 700 |
|
manner across all systems and configurations. |
| 701 |
|
|
| 702 |
< |
The althernative methods were also evaluated with three different |
| 702 |
> |
The alternative methods were also evaluated with three different |
| 703 |
|
cutoff radii (9, 12, and 15 \AA). As noted previously, the |
| 704 |
|
convergence parameter ($\alpha$) plays a role in the balance of the |
| 705 |
|
real-space and reciprocal-space portions of the Ewald calculation. |
| 708 |
|
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 709 |
|
increasing accuracy at the expense of computational time spent on the |
| 710 |
|
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 711 |
< |
The default TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 712 |
< |
in all SPME calculations, resulting in Ewald coefficients of 0.4200, |
| 711 |
> |
The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 712 |
> |
in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200, |
| 713 |
|
0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ |
| 714 |
|
respectively. |
| 715 |
|
|
| 719 |
|
In order to evaluate the performance of the pairwise electrostatic |
| 720 |
|
summation methods for Monte Carlo simulations, the energy differences |
| 721 |
|
between configurations were compared to the values obtained when using |
| 722 |
< |
SPME. The results for the subsequent regression analysis are shown in |
| 722 |
> |
{\sc spme}. The results for the subsequent regression analysis are shown in |
| 723 |
|
figure \ref{fig:delE}. |
| 724 |
|
|
| 725 |
|
\begin{figure} |
| 729 |
|
differences for a given electrostatic method compared with the |
| 730 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 731 |
|
indicate $\Delta E$ values indistinguishable from those obtained using |
| 732 |
< |
SPME. Different values of the cutoff radius are indicated with |
| 732 |
> |
{\sc spme}. Different values of the cutoff radius are indicated with |
| 733 |
|
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 734 |
|
inverted triangles).} |
| 735 |
|
\label{fig:delE} |
| 751 |
|
significant improvement using the group-switched cutoff because the |
| 752 |
|
salt and salt solution systems contain non-neutral groups. Interested |
| 753 |
|
readers can consult the accompanying supporting information for a |
| 754 |
< |
comparison where all groups are neutral. |
| 754 |
> |
comparison where all groups are neutral.\cite{EPAPSdeposit} |
| 755 |
|
|
| 756 |
|
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 757 |
|
the agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ |
| 758 |
< |
shows an excellent correlation and quality of fit with the SPME |
| 758 |
> |
shows an excellent correlation and quality of fit with the {\sc spme} |
| 759 |
|
results, particularly with a cutoff radius greater than 12 |
| 760 |
|
\AA . Use of a larger damping parameter is more helpful for the |
| 761 |
|
shortest cutoff shown, but it has a detrimental effect on simulations |
| 789 |
|
magnitudes for a given electrostatic method compared with the |
| 790 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 791 |
|
indicate force magnitude values indistinguishable from those obtained |
| 792 |
< |
using SPME. Different values of the cutoff radius are indicated with |
| 792 |
> |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 793 |
|
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 794 |
|
inverted triangles).} |
| 795 |
|
\label{fig:frcMag} |
| 796 |
|
\end{figure} |
| 797 |
|
|
| 798 |
|
Again, it is striking how well the Shifted Potential and Shifted Force |
| 799 |
< |
methods are doing at reproducing the SPME forces. The undamped and |
| 799 |
> |
methods are doing at reproducing the {\sc spme} forces. The undamped and |
| 800 |
|
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
| 801 |
|
This is perhaps expected because this method explicitly incorporates a |
| 802 |
|
smooth transition in the forces at the cutoff radius as well as the |
| 828 |
|
magnitudes for a given electrostatic method compared with the |
| 829 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 830 |
|
indicate torque magnitude values indistinguishable from those obtained |
| 831 |
< |
using SPME. Different values of the cutoff radius are indicated with |
| 831 |
> |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 832 |
|
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 833 |
|
inverted triangles).} |
| 834 |
|
\label{fig:trqMag} |
| 840 |
|
|
| 841 |
|
Torques appear to be much more sensitive to charges at a longer |
| 842 |
|
distance. The striking feature in comparing the new electrostatic |
| 843 |
< |
methods with SPME is how much the agreement improves with increasing |
| 843 |
> |
methods with {\sc spme} is how much the agreement improves with increasing |
| 844 |
|
cutoff radius. Again, the weakly damped and undamped {\sc sf} method |
| 845 |
< |
appears to be reproducing the SPME torques most accurately. |
| 845 |
> |
appears to be reproducing the {\sc spme} torques most accurately. |
| 846 |
|
|
| 847 |
|
Water molecules are dipolar, and the reaction field method reproduces |
| 848 |
|
the effect of the surrounding polarized medium on each of the |
| 857 |
|
vital in calculating dynamical quantities accurately. Force and |
| 858 |
|
torque directionalities were investigated by measuring the angles |
| 859 |
|
formed between these vectors and the same vectors calculated using |
| 860 |
< |
SPME. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 860 |
> |
{\sc spme}. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 861 |
|
variance ($\sigma^2$) of the Gaussian fits of the angle error |
| 862 |
|
distributions of the combined set over all system types. |
| 863 |
|
|
| 869 |
|
make with their counterparts obtained using the reference Ewald sum. |
| 870 |
|
Results with a variance ($\sigma^2$) equal to zero (dashed line) |
| 871 |
|
indicate force and torque directions indistinguishable from those |
| 872 |
< |
obtained using SPME. Different values of the cutoff radius are |
| 872 |
> |
obtained using {\sc spme}. Different values of the cutoff radius are |
| 873 |
|
indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, |
| 874 |
|
and 15\AA\ = inverted triangles).} |
| 875 |
|
\label{fig:frcTrqAng} |
| 878 |
|
Both the force and torque $\sigma^2$ results from the analysis of the |
| 879 |
|
total accumulated system data are tabulated in figure |
| 880 |
|
\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc |
| 881 |
< |
sp}) method would be essentially unusable for molecular dynamics until |
| 882 |
< |
the damping function is added. The Shifted Force ({\sc sf}) method, |
| 883 |
< |
however, is generating force and torque vectors which are within a few |
| 884 |
< |
degrees of the Ewald results even with weak (or no) damping. |
| 881 |
> |
sp}) method would be essentially unusable for molecular dynamics |
| 882 |
> |
unless the damping function is added. The Shifted Force ({\sc sf}) |
| 883 |
> |
method, however, is generating force and torque vectors which are |
| 884 |
> |
within a few degrees of the Ewald results even with weak (or no) |
| 885 |
> |
damping. |
| 886 |
|
|
| 887 |
|
All of the sets (aside from the over-damped case) show the improvement |
| 888 |
|
afforded by choosing a larger cutoff radius. Increasing the cutoff |
| 889 |
|
from 9 to 12 \AA\ typically results in a halving of the width of the |
| 890 |
< |
distribution, with a similar improvement going from 12 to 15 |
| 890 |
> |
distribution, with a similar improvement when going from 12 to 15 |
| 891 |
|
\AA . |
| 892 |
|
|
| 893 |
|
The undamped {\sc sf}, group-based cutoff, and reaction field methods |
| 894 |
|
all do equivalently well at capturing the direction of both the force |
| 895 |
< |
and torque vectors. Using damping improves the angular behavior |
| 896 |
< |
significantly for the {\sc sp} and moderately for the {\sc sf} |
| 897 |
< |
methods. Overdamping is detrimental to both methods. Again it is |
| 898 |
< |
important to recognize that the force vectors cover all particles in |
| 899 |
< |
the systems, while torque vectors are only available for neutral |
| 900 |
< |
molecular groups. Damping appears to have a more beneficial effect on |
| 895 |
> |
and torque vectors. Using the electrostatic damping improves the |
| 896 |
> |
angular behavior significantly for the {\sc sp} and moderately for the |
| 897 |
> |
{\sc sf} methods. Overdamping is detrimental to both methods. Again |
| 898 |
> |
it is important to recognize that the force vectors cover all |
| 899 |
> |
particles in all seven systems, while torque vectors are only |
| 900 |
> |
available for neutral molecular groups. Damping is more beneficial to |
| 901 |
|
charged bodies, and this observation is investigated further in the |
| 902 |
< |
accompanying supporting information. |
| 902 |
> |
accompanying supporting information.\cite{EPAPSdeposit} |
| 903 |
|
|
| 904 |
|
Although not discussed previously, group based cutoffs can be applied |
| 905 |
< |
to both the {\sc sp} and {\sc sf} methods. Use of a switching |
| 906 |
< |
function corrects for the discontinuities that arise when atoms of the |
| 907 |
< |
two groups exit the cutoff radius before the group centers leave each |
| 908 |
< |
other's cutoff. Though there are no significant benefits or drawbacks |
| 909 |
< |
observed in $\Delta E$ and vector magnitude results when doing this, |
| 910 |
< |
there is a measurable improvement in the vector angle results. Table |
| 911 |
< |
\ref{tab:groupAngle} shows the angular variance values obtained using |
| 912 |
< |
group based cutoffs and a switching function alongside the results |
| 913 |
< |
seen in figure \ref{fig:frcTrqAng}. The {\sc sp} shows much narrower |
| 914 |
< |
angular distributions for both the force and torque vectors when using |
| 915 |
< |
an $\alpha$ of 0.2 \AA$^{-1}$ or less, while {\sc sf} shows |
| 916 |
< |
improvements in the undamped and lightly damped cases. Thus, by |
| 931 |
< |
calculating the electrostatic interactions in terms of molecular pairs |
| 932 |
< |
rather than atomic pairs, the direction of the force and torque |
| 933 |
< |
vectors can be determined more accurately. |
| 905 |
> |
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 906 |
> |
will reintroduce small discontinuities at the cutoff radius, but the |
| 907 |
> |
effects of these can be minimized by utilizing a switching function. |
| 908 |
> |
Though there are no significant benefits or drawbacks observed in |
| 909 |
> |
$\Delta E$ and the force and torque magnitudes when doing this, there |
| 910 |
> |
is a measurable improvement in the directionality of the forces and |
| 911 |
> |
torques. Table \ref{tab:groupAngle} shows the angular variances |
| 912 |
> |
obtained using group based cutoffs along with the results seen in |
| 913 |
> |
figure \ref{fig:frcTrqAng}. The {\sc sp} (with an $\alpha$ of 0.2 |
| 914 |
> |
\AA$^{-1}$ or smaller) shows much narrower angular distributions when |
| 915 |
> |
using group-based cutoffs. The {\sc sf} method likewise shows |
| 916 |
> |
improvement in the undamped and lightly damped cases. |
| 917 |
|
|
| 918 |
|
\begin{table}[htbp] |
| 919 |
< |
\centering |
| 920 |
< |
\caption{Variance ($\sigma^2$) of the force (top set) and torque |
| 921 |
< |
(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$.} |
| 919 |
> |
\centering |
| 920 |
> |
\caption{Statistical analysis of the angular |
| 921 |
> |
distributions that the force (upper) and torque (lower) vectors |
| 922 |
> |
from a given electrostatic method make with their counterparts |
| 923 |
> |
obtained using the reference Ewald sum. Calculations were |
| 924 |
> |
performed both with (Y) and without (N) group based cutoffs and a |
| 925 |
> |
switching function. The $\alpha$ values have units of \AA$^{-1}$ |
| 926 |
> |
and the variance values have units of degrees$^2$.} |
| 927 |
> |
|
| 928 |
|
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 929 |
|
\\ |
| 930 |
|
\toprule |
| 955 |
|
\label{tab:groupAngle} |
| 956 |
|
\end{table} |
| 957 |
|
|
| 958 |
< |
One additional trend to recognize in table \ref{tab:groupAngle} is |
| 959 |
< |
that the $\sigma^2$ values for both {\sc sp} and {\sc sf} converge as |
| 960 |
< |
$\alpha$ increases, something that is easier to see when using group |
| 961 |
< |
based cutoffs. The reason for this is that the complimentary error |
| 962 |
< |
function inserted into the potential weakens the electrostatic |
| 963 |
< |
interaction as $\alpha$ increases. Thus, at larger values of |
| 964 |
< |
$\alpha$, both summation methods progress toward non-interacting |
| 976 |
< |
functions, so care is required in choosing large damping functions |
| 977 |
< |
lest one generate an undesirable loss in the pair interaction. Kast |
| 958 |
> |
One additional trend in table \ref{tab:groupAngle} is that the |
| 959 |
> |
$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$ |
| 960 |
> |
increases, something that is more obvious with group-based cutoffs. |
| 961 |
> |
The complimentary error function inserted into the potential weakens |
| 962 |
> |
the electrostatic interaction as the value of $\alpha$ is increased. |
| 963 |
> |
However, at larger values of $\alpha$, it is possible to overdamp the |
| 964 |
> |
electrostatic interaction and to remove it completely. Kast |
| 965 |
|
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 966 |
|
values for these types of electrostatic summation methods by fitting |
| 967 |
|
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
| 969 |
|
respectively.\cite{Kast03} These appear to be reasonable choices to |
| 970 |
|
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
| 971 |
|
these findings, choices this high would introduce error in the |
| 972 |
< |
molecular torques, particularly for the shorter cutoffs. Based on the |
| 973 |
< |
above findings, empirical damping up to 0.2 \AA$^{-1}$ proves to be |
| 974 |
< |
beneficial, but damping may be unnecessary when using the {\sc sf} |
| 988 |
< |
method. |
| 972 |
> |
molecular torques, particularly for the shorter cutoffs. Based on our |
| 973 |
> |
observations, empirical damping up to 0.2 \AA$^{-1}$ is beneficial, |
| 974 |
> |
but damping may be unnecessary when using the {\sc sf} method. |
| 975 |
|
|
| 976 |
|
\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals} |
| 977 |
|
|
| 984 |
|
distribution functions, diffusion constants, and velocity and |
| 985 |
|
orientational correlation functions) may not be particularly sensitive |
| 986 |
|
to the long-range and collective behavior that governs the |
| 987 |
< |
low-frequency behavior in crystalline systems. |
| 987 |
> |
low-frequency behavior in crystalline systems. Additionally, the |
| 988 |
> |
ionic crystals are the worst case scenario for the pairwise methods |
| 989 |
> |
because they lack the reciprocal space contribution contained in the |
| 990 |
> |
Ewald summation. |
| 991 |
|
|
| 992 |
|
We are using two separate measures to probe the effects of these |
| 993 |
|
alternative electrostatic methods on the dynamics in crystalline |
| 999 |
|
\begin{figure} |
| 1000 |
|
\centering |
| 1001 |
|
\includegraphics[width = \linewidth]{./vCorrPlot.pdf} |
| 1002 |
< |
\caption{Velocity auto-correlation functions of NaCl crystals at |
| 1003 |
< |
1000 K using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc |
| 1002 |
> |
\caption{Velocity autocorrelation functions of NaCl crystals at |
| 1003 |
> |
1000 K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc |
| 1004 |
|
sp} ($\alpha$ = 0.2). The inset is a magnification of the area around |
| 1005 |
|
the first minimum. The times to first collision are nearly identical, |
| 1006 |
|
but differences can be seen in the peaks and troughs, where the |
| 1007 |
|
undamped and weakly damped methods are stiffer than the moderately |
| 1008 |
< |
damped and SPME methods.} |
| 1008 |
> |
damped and {\sc spme} methods.} |
| 1009 |
|
\label{fig:vCorrPlot} |
| 1010 |
|
\end{figure} |
| 1011 |
|
|
| 1012 |
< |
The short-time decays through the first collision are nearly identical |
| 1013 |
< |
in figure \ref{fig:vCorrPlot}, but the peaks and troughs of the |
| 1014 |
< |
functions show how the methods differ. The undamped {\sc sf} method |
| 1015 |
< |
has deeper troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher |
| 1016 |
< |
peaks than any of the other methods. As the damping function is |
| 1017 |
< |
increased, these peaks are smoothed out, and approach the SPME |
| 1018 |
< |
curve. The damping acts as a distance dependent Gaussian screening of |
| 1019 |
< |
the point charges for the pairwise summation methods; thus, the |
| 1020 |
< |
collisions are more elastic in the undamped {\sc sf} potential, and the |
| 1021 |
< |
stiffness of the potential is diminished as the electrostatic |
| 1022 |
< |
interactions are softened by the damping function. With $\alpha$ |
| 1023 |
< |
values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are |
| 1024 |
< |
nearly identical and track the SPME features quite well. This is not |
| 1025 |
< |
too surprising in that the differences between the {\sc sf} and {\sc |
| 1026 |
< |
sp} potentials are mitigated with increased damping. However, this |
| 1038 |
< |
appears to indicate that once damping is utilized, the form of the |
| 1039 |
< |
potential seems to play a lesser role in the crystal dynamics. |
| 1012 |
> |
The short-time decay of the velocity autocorrelation function through |
| 1013 |
> |
the first collision are nearly identical in figure |
| 1014 |
> |
\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show |
| 1015 |
> |
how the methods differ. The undamped {\sc sf} method has deeper |
| 1016 |
> |
troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than |
| 1017 |
> |
any of the other methods. As the damping parameter ($\alpha$) is |
| 1018 |
> |
increased, these peaks are smoothed out, and the {\sc sf} method |
| 1019 |
> |
approaches the {\sc spme} results. With $\alpha$ values of 0.2 \AA$^{-1}$, |
| 1020 |
> |
the {\sc sf} and {\sc sp} functions are nearly identical and track the |
| 1021 |
> |
{\sc spme} features quite well. This is not surprising because the {\sc sf} |
| 1022 |
> |
and {\sc sp} potentials become nearly identical with increased |
| 1023 |
> |
damping. However, this appears to indicate that once damping is |
| 1024 |
> |
utilized, the details of the form of the potential (and forces) |
| 1025 |
> |
constructed out of the damped electrostatic interaction are less |
| 1026 |
> |
important. |
| 1027 |
|
|
| 1028 |
|
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 1029 |
|
|
| 1030 |
< |
The short time dynamics were extended to evaluate how the differences |
| 1031 |
< |
between the methods affect the collective long-time motion. The same |
| 1032 |
< |
electrostatic summation methods were used as in the short time |
| 1033 |
< |
velocity autocorrelation function evaluation, but the trajectories |
| 1034 |
< |
were sampled over a much longer time. The power spectra of the |
| 1035 |
< |
resulting velocity autocorrelation functions were calculated and are |
| 1036 |
< |
displayed in figure \ref{fig:methodPS}. |
| 1030 |
> |
To evaluate how the differences between the methods affect the |
| 1031 |
> |
collective long-time motion, we computed power spectra from long-time |
| 1032 |
> |
traces of the velocity autocorrelation function. The power spectra for |
| 1033 |
> |
the best-performing alternative methods are shown in |
| 1034 |
> |
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
| 1035 |
> |
a cubic switching function between 40 and 50 ps was used to reduce the |
| 1036 |
> |
ringing resulting from data truncation. This procedure had no |
| 1037 |
> |
noticeable effect on peak location or magnitude. |
| 1038 |
|
|
| 1039 |
|
\begin{figure} |
| 1040 |
|
\centering |
| 1041 |
|
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 1042 |
|
\caption{Power spectra obtained from the velocity auto-correlation |
| 1043 |
< |
functions of NaCl crystals at 1000 K while using SPME, {\sc sf} |
| 1044 |
< |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). |
| 1045 |
< |
Apodization of the correlation functions via a cubic switching |
| 1046 |
< |
function between 40 and 50 ps was used to clear up the spectral noise |
| 1059 |
< |
resulting from data truncation, and had no noticeable effect on peak |
| 1060 |
< |
location or magnitude. The inset shows the frequency region below 100 |
| 1061 |
< |
cm$^{-1}$ to highlight where the spectra begin to differ.} |
| 1043 |
> |
functions of NaCl crystals at 1000 K while using {\sc spme}, {\sc sf} |
| 1044 |
> |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset |
| 1045 |
> |
shows the frequency region below 100 cm$^{-1}$ to highlight where the |
| 1046 |
> |
spectra differ.} |
| 1047 |
|
\label{fig:methodPS} |
| 1048 |
|
\end{figure} |
| 1049 |
|
|
| 1050 |
< |
While high frequency peaks of the spectra in this figure overlap, |
| 1051 |
< |
showing the same general features, the low frequency region shows how |
| 1052 |
< |
the summation methods differ. Considering the low-frequency inset |
| 1053 |
< |
(expanded in the upper frame of figure \ref{fig:dampInc}), at |
| 1054 |
< |
frequencies below 100 cm$^{-1}$, the correlated motions are |
| 1055 |
< |
blue-shifted when using undamped or weakly damped {\sc sf}. When |
| 1056 |
< |
using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the {\sc sf} |
| 1057 |
< |
and {\sc sp} methods give near identical correlated motion behavior as |
| 1058 |
< |
the Ewald method (which has a damping value of 0.3119). This |
| 1059 |
< |
weakening of the electrostatic interaction with increased damping |
| 1060 |
< |
explains why the long-ranged correlated motions are at lower |
| 1061 |
< |
frequencies for the moderately damped methods than for undamped or |
| 1062 |
< |
weakly damped methods. To see this effect more clearly, we show how |
| 1063 |
< |
damping strength alone affects a simple real-space electrostatic |
| 1064 |
< |
potential, |
| 1065 |
< |
\begin{equation} |
| 1066 |
< |
V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r), |
| 1067 |
< |
\end{equation} |
| 1068 |
< |
where $S(r)$ is a switching function that smoothly zeroes the |
| 1069 |
< |
potential at the cutoff radius. Figure \ref{fig:dampInc} shows how |
| 1070 |
< |
the low frequency motions are dependent on the damping used in the |
| 1071 |
< |
direct electrostatic sum. As the damping increases, the peaks drop to |
| 1072 |
< |
lower frequencies. Incidentally, use of an $\alpha$ of 0.25 |
| 1073 |
< |
\AA$^{-1}$ on a simple electrostatic summation results in low |
| 1074 |
< |
frequency correlated dynamics equivalent to a simulation using SPME. |
| 1075 |
< |
When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks |
| 1076 |
< |
shift to higher frequency in exponential fashion. Though not shown, |
| 1077 |
< |
the spectrum for the simple undamped electrostatic potential is |
| 1093 |
< |
blue-shifted such that the lowest frequency peak resides near 325 |
| 1094 |
< |
cm$^{-1}$. In light of these results, the undamped {\sc sf} method |
| 1095 |
< |
producing low-lying motion peaks within 10 cm$^{-1}$ of SPME is quite |
| 1096 |
< |
respectable and shows that the shifted force procedure accounts for |
| 1097 |
< |
most of the effect afforded through use of the Ewald summation. |
| 1098 |
< |
However, it appears as though moderate damping is required for |
| 1099 |
< |
accurate reproduction of crystal dynamics. |
| 1050 |
> |
While the high frequency regions of the power spectra for the |
| 1051 |
> |
alternative methods are quantitatively identical with Ewald spectrum, |
| 1052 |
> |
the low frequency region shows how the summation methods differ. |
| 1053 |
> |
Considering the low-frequency inset (expanded in the upper frame of |
| 1054 |
> |
figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the |
| 1055 |
> |
correlated motions are blue-shifted when using undamped or weakly |
| 1056 |
> |
damped {\sc sf}. When using moderate damping ($\alpha = 0.2$ |
| 1057 |
> |
\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical |
| 1058 |
> |
correlated motion to the Ewald method (which has a convergence |
| 1059 |
> |
parameter of 0.3119 \AA$^{-1}$). This weakening of the electrostatic |
| 1060 |
> |
interaction with increased damping explains why the long-ranged |
| 1061 |
> |
correlated motions are at lower frequencies for the moderately damped |
| 1062 |
> |
methods than for undamped or weakly damped methods. |
| 1063 |
> |
|
| 1064 |
> |
To isolate the role of the damping constant, we have computed the |
| 1065 |
> |
spectra for a single method ({\sc sf}) with a range of damping |
| 1066 |
> |
constants and compared this with the {\sc spme} spectrum. |
| 1067 |
> |
Fig. \ref{fig:dampInc} shows more clearly that increasing the |
| 1068 |
> |
electrostatic damping red-shifts the lowest frequency phonon modes. |
| 1069 |
> |
However, even without any electrostatic damping, the {\sc sf} method |
| 1070 |
> |
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
| 1071 |
> |
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
| 1072 |
> |
would predict the lowest frequency peak near 325 cm$^{-1}$. {\it |
| 1073 |
> |
Most} of the collective behavior in the crystal is accurately captured |
| 1074 |
> |
using the {\sc sf} method. Quantitative agreement with Ewald can be |
| 1075 |
> |
obtained using moderate damping in addition to the shifting at the |
| 1076 |
> |
cutoff distance. |
| 1077 |
> |
|
| 1078 |
|
\begin{figure} |
| 1079 |
|
\centering |
| 1080 |
< |
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 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 |
| 1110 |
< |
$\alpha$ increases, the peaks are red-shifted toward and eventually |
| 1111 |
< |
beyond the values given by SPME. The larger $\alpha$ values weaken |
| 1112 |
< |
the real-space electrostatics, explaining this shift towards less |
| 1113 |
< |
strongly correlated motions in the crystal.} |
| 1080 |
> |
\includegraphics[width = \linewidth]{./increasedDamping.pdf} |
| 1081 |
> |
\caption{Effect of damping on the two lowest-frequency phonon modes in |
| 1082 |
> |
the NaCl crystal at 1000~K. The undamped shifted force ({\sc sf}) |
| 1083 |
> |
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1084 |
> |
electrostatic damping to 0.25 \AA$^{-1}$ gives quantitative agreement |
| 1085 |
> |
with the power spectrum obtained using the Ewald sum. Overdamping can |
| 1086 |
> |
result in underestimates of frequencies of the long-wavelength |
| 1087 |
> |
motions.} |
| 1088 |
|
\label{fig:dampInc} |
| 1089 |
|
\end{figure} |
| 1090 |
|
|
| 1091 |
|
\section{Conclusions} |
| 1092 |
|
|
| 1093 |
|
This investigation of pairwise electrostatic summation techniques |
| 1094 |
< |
shows that there are viable and more computationally efficient |
| 1095 |
< |
electrostatic summation techniques than the Ewald summation, chiefly |
| 1096 |
< |
methods derived from the damped Coulombic sum originally proposed by |
| 1097 |
< |
Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the |
| 1098 |
< |
{\sc sf} method, reformulated above as eq. (\ref{eq:DSFPot}), |
| 1099 |
< |
shows a remarkable ability to reproduce the energetic and dynamic |
| 1100 |
< |
characteristics exhibited by simulations employing lattice summation |
| 1101 |
< |
techniques. The cumulative energy difference results showed the |
| 1102 |
< |
undamped {\sc sf} and moderately damped {\sc sp} methods |
| 1103 |
< |
produced results nearly identical to SPME. Similarly for the dynamic |
| 1104 |
< |
features, the undamped or moderately damped {\sc sf} and |
| 1105 |
< |
moderately damped {\sc sp} methods produce force and torque |
| 1106 |
< |
vector magnitude and directions very similar to the expected values. |
| 1107 |
< |
These results translate into long-time dynamic behavior equivalent to |
| 1108 |
< |
that produced in simulations using SPME. |
| 1094 |
> |
shows that there are viable and computationally efficient alternatives |
| 1095 |
> |
to the Ewald summation. These methods are derived from the damped and |
| 1096 |
> |
cutoff-neutralized Coulombic sum originally proposed by Wolf |
| 1097 |
> |
\textit{et al.}\cite{Wolf99} In particular, the {\sc sf} |
| 1098 |
> |
method, reformulated above as eqs. (\ref{eq:DSFPot}) and |
| 1099 |
> |
(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the |
| 1100 |
> |
energetic and dynamic characteristics exhibited by simulations |
| 1101 |
> |
employing lattice summation techniques. The cumulative energy |
| 1102 |
> |
difference results showed the undamped {\sc sf} and moderately damped |
| 1103 |
> |
{\sc sp} methods produced results nearly identical to {\sc spme}. Similarly |
| 1104 |
> |
for the dynamic features, the undamped or moderately damped {\sc sf} |
| 1105 |
> |
and moderately damped {\sc sp} methods produce force and torque vector |
| 1106 |
> |
magnitude and directions very similar to the expected values. These |
| 1107 |
> |
results translate into long-time dynamic behavior equivalent to that |
| 1108 |
> |
produced in simulations using {\sc spme}. |
| 1109 |
|
|
| 1110 |
+ |
As in all purely-pairwise cutoff methods, these methods are expected |
| 1111 |
+ |
to scale approximately {\it linearly} with system size, and they are |
| 1112 |
+ |
easily parallelizable. This should result in substantial reductions |
| 1113 |
+ |
in the computational cost of performing large simulations. |
| 1114 |
+ |
|
| 1115 |
|
Aside from the computational cost benefit, these techniques have |
| 1116 |
|
applicability in situations where the use of the Ewald sum can prove |
| 1117 |
< |
problematic. Primary among them is their use in interfacial systems, |
| 1118 |
< |
where the unmodified lattice sum techniques artificially accentuate |
| 1119 |
< |
the periodicity of the system in an undesirable manner. There have |
| 1120 |
< |
been alterations to the standard Ewald techniques, via corrections and |
| 1121 |
< |
reformulations, to compensate for these systems; but the pairwise |
| 1122 |
< |
techniques discussed here require no modifications, making them |
| 1123 |
< |
natural tools to tackle these problems. Additionally, this |
| 1124 |
< |
transferability gives them benefits over other pairwise methods, like |
| 1125 |
< |
reaction field, because estimations of physical properties (e.g. the |
| 1126 |
< |
dielectric constant) are unnecessary. |
| 1117 |
> |
problematic. Of greatest interest is their potential use in |
| 1118 |
> |
interfacial systems, where the unmodified lattice sum techniques |
| 1119 |
> |
artificially accentuate the periodicity of the system in an |
| 1120 |
> |
undesirable manner. There have been alterations to the standard Ewald |
| 1121 |
> |
techniques, via corrections and reformulations, to compensate for |
| 1122 |
> |
these systems; but the pairwise techniques discussed here require no |
| 1123 |
> |
modifications, making them natural tools to tackle these problems. |
| 1124 |
> |
Additionally, this transferability gives them benefits over other |
| 1125 |
> |
pairwise methods, like reaction field, because estimations of physical |
| 1126 |
> |
properties (e.g. the dielectric constant) are unnecessary. |
| 1127 |
|
|
| 1128 |
< |
We are not suggesting any flaw with the Ewald sum; in fact, it is the |
| 1129 |
< |
standard by which these simple pairwise sums are judged. However, |
| 1130 |
< |
these results do suggest that in the typical simulations performed |
| 1131 |
< |
today, the Ewald summation may no longer be required to obtain the |
| 1132 |
< |
level of accuracy most researchers have come to expect |
| 1128 |
> |
If a researcher is using Monte Carlo simulations of large chemical |
| 1129 |
> |
systems containing point charges, most structural features will be |
| 1130 |
> |
accurately captured using the undamped {\sc sf} method or the {\sc sp} |
| 1131 |
> |
method with an electrostatic damping of 0.2 \AA$^{-1}$. These methods |
| 1132 |
> |
would also be appropriate for molecular dynamics simulations where the |
| 1133 |
> |
data of interest is either structural or short-time dynamical |
| 1134 |
> |
quantities. For long-time dynamics and collective motions, the safest |
| 1135 |
> |
pairwise method we have evaluated is the {\sc sf} method with an |
| 1136 |
> |
electrostatic damping between 0.2 and 0.25 |
| 1137 |
> |
\AA$^{-1}$. |
| 1138 |
|
|
| 1139 |
+ |
We are not suggesting that there is any flaw with the Ewald sum; in |
| 1140 |
+ |
fact, it is the standard by which these simple pairwise sums have been |
| 1141 |
+ |
judged. However, these results do suggest that in the typical |
| 1142 |
+ |
simulations performed today, the Ewald summation may no longer be |
| 1143 |
+ |
required to obtain the level of accuracy most researchers have come to |
| 1144 |
+ |
expect. |
| 1145 |
+ |
|
| 1146 |
|
\section{Acknowledgments} |
| 1147 |
+ |
Support for this project was provided by the National Science |
| 1148 |
+ |
Foundation under grant CHE-0134881. The authors would like to thank |
| 1149 |
+ |
Steve Corcelli and Ed Maginn for helpful discussions and comments. |
| 1150 |
+ |
|
| 1151 |
|
\newpage |
| 1152 |
|
|
| 1153 |
|
\bibliographystyle{jcp2} |
