| 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 necessary? \\ |
| 29 |
> |
Pairwise alternatives to the accepted standard for \\ |
| 30 |
> |
long-range electrostatics} |
| 31 |
|
|
| 32 |
|
\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
| 33 |
|
gezelter@nd.edu} \\ |
| 42 |
|
|
| 43 |
|
\nobibliography{} |
| 44 |
|
\begin{abstract} |
| 45 |
< |
A new method for accumulating electrostatic interactions was derived |
| 46 |
< |
from the previous efforts described in \bibentry{Wolf99} and |
| 47 |
< |
\bibentry{Zahn02} as a possible replacement for lattice sum methods in |
| 48 |
< |
molecular simulations. Comparisons were performed with this and other |
| 49 |
< |
pairwise electrostatic summation techniques against the smooth |
| 50 |
< |
particle mesh Ewald (SPME) summation to see how well they reproduce |
| 51 |
< |
the energetics and dynamics of a variety of simulation types. The |
| 52 |
< |
newly derived Shifted-Force technique shows a remarkable ability to |
| 53 |
< |
reproduce the behavior exhibited in simulations using SPME with an |
| 54 |
< |
$\mathscr{O}(N)$ computational cost, equivalent to merely the |
| 55 |
< |
real-space portion of the lattice summation. |
| 56 |
< |
|
| 45 |
> |
We investigate pairwise electrostatic interaction methods and show |
| 46 |
> |
that there are viable and computationally efficient $(\mathscr{O}(N))$ |
| 47 |
> |
alternatives to the Ewald summation for typical modern molecular |
| 48 |
> |
simulations. These methods are extended from the damped and |
| 49 |
> |
cutoff-neutralized Coulombic sum originally proposed by Wolf |
| 50 |
> |
\textit{et al.} One of these, the damped shifted force method, shows |
| 51 |
> |
a remarkable ability to reproduce the energetic and dynamic |
| 52 |
> |
characteristics exhibited by simulations employing lattice summation |
| 53 |
> |
techniques. Comparisons were performed with this and other pairwise |
| 54 |
> |
methods against the smooth particle mesh Ewald ({\sc spme}) summation to see |
| 55 |
> |
how well they reproduce the energetics and dynamics of a variety of |
| 56 |
> |
simulation types. |
| 57 |
|
\end{abstract} |
| 58 |
|
|
| 59 |
|
\newpage |
| 167 |
|
|
| 168 |
|
\begin{figure} |
| 169 |
|
\centering |
| 170 |
< |
\includegraphics[width = \linewidth]{./ewaldProgression2.pdf} |
| 170 |
> |
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
| 171 |
|
\caption{The change in the application of the Ewald sum with |
| 172 |
|
increasing computational power. Initially, only small systems could |
| 173 |
|
be studied, and the Ewald sum replicated the simulation box to |
| 518 |
|
The pairwise summation techniques (outlined in section |
| 519 |
|
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
| 520 |
|
studying the energy differences between conformations. We took the |
| 521 |
< |
SPME-computed energy difference between two conformations to be the |
| 521 |
> |
{\sc spme}-computed energy difference between two conformations to be the |
| 522 |
|
correct behavior. An ideal performance by an alternative method would |
| 523 |
|
reproduce these energy differences exactly (even if the absolute |
| 524 |
|
energies calculated by the methods are different). Since none of the |
| 526 |
|
regressions of energy gap data to evaluate how closely the methods |
| 527 |
|
mimicked the Ewald energy gaps. Unitary results for both the |
| 528 |
|
correlation (slope) and correlation coefficient for these regressions |
| 529 |
< |
indicate perfect agreement between the alternative method and SPME. |
| 529 |
> |
indicate perfect agreement between the alternative method and {\sc spme}. |
| 530 |
|
Sample correlation plots for two alternate methods are shown in |
| 531 |
|
Fig. \ref{fig:linearFit}. |
| 532 |
|
|
| 555 |
|
We evaluated the pairwise methods (outlined in section |
| 556 |
|
\ref{sec:ESMethods}) for use in MD simulations by |
| 557 |
|
comparing the force and torque vectors with those obtained using the |
| 558 |
< |
reference Ewald summation (SPME). Both the magnitude and the |
| 558 |
> |
reference Ewald summation ({\sc spme}). Both the magnitude and the |
| 559 |
|
direction of these vectors on each of the bodies in the system were |
| 560 |
|
analyzed. For the magnitude of these vectors, linear least squares |
| 561 |
|
regression analyses were performed as described previously for |
| 570 |
|
|
| 571 |
|
The {\it directionality} of the force and torque vectors was |
| 572 |
|
investigated through measurement of the angle ($\theta$) formed |
| 573 |
< |
between those computed from the particular method and those from SPME, |
| 573 |
> |
between those computed from the particular method and those from {\sc spme}, |
| 574 |
|
\begin{equation} |
| 575 |
|
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right), |
| 576 |
|
\end{equation} |
| 577 |
< |
where $\hat{f}_\textrm{M}$ is the unit vector pointing along the force |
| 577 |
> |
where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force |
| 578 |
|
vector computed using method M. Each of these $\theta$ values was |
| 579 |
|
accumulated in a distribution function and weighted by the area on the |
| 580 |
|
unit sphere. Since this distribution is a measure of angular error |
| 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). |
| 606 |
> |
calculated when using the reference method ({\sc spme}). |
| 607 |
|
|
| 608 |
|
\subsection{Short-time Dynamics} |
| 609 |
|
|
| 618 |
|
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
| 619 |
|
of the trajectories, |
| 620 |
|
\begin{equation} |
| 621 |
< |
C_v(t) = \frac{\langle v_i(0)\cdot v_i(t)\rangle}{\langle v_i(0)\cdot v_i(0)\rangle}. |
| 621 |
> |
C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}. |
| 622 |
|
\label{eq:vCorr} |
| 623 |
|
\end{equation} |
| 624 |
|
Velocity autocorrelation functions require detailed short time data, |
| 646 |
|
\subsection{Representative Simulations}\label{sec:RepSims} |
| 647 |
|
A variety of representative simulations were analyzed to determine the |
| 648 |
|
relative effectiveness of the pairwise summation techniques in |
| 649 |
< |
reproducing the energetics and dynamics exhibited by SPME. We wanted |
| 649 |
> |
reproducing the energetics and dynamics exhibited by {\sc spme}. We wanted |
| 650 |
|
to span the space of modern simulations (i.e. from liquids of neutral |
| 651 |
|
molecules to ionic crystals), so the systems studied were: |
| 652 |
|
\begin{enumerate} |
| 697 |
|
|
| 698 |
|
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
| 699 |
|
We compared the following alternative summation methods with results |
| 700 |
< |
from the reference method (SPME): |
| 700 |
> |
from the reference method ({\sc spme}): |
| 701 |
|
\begin{itemize} |
| 702 |
|
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
| 703 |
|
and 0.3 \AA$^{-1}$, |
| 708 |
|
\end{itemize} |
| 709 |
|
Group-based cutoffs with a fifth-order polynomial switching function |
| 710 |
|
were utilized for the reaction field simulations. Additionally, we |
| 711 |
< |
investigated the use of these cutoffs with the SP, SF, and pure |
| 712 |
< |
cutoff. The SPME electrostatics were performed using the TINKER |
| 713 |
< |
implementation of SPME,\cite{Ponder87} while all other calculations |
| 711 |
> |
investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure |
| 712 |
> |
cutoff. The {\sc spme} electrostatics were performed using the {\sc tinker} |
| 713 |
> |
implementation of {\sc spme},\cite{Ponder87} while all other calculations |
| 714 |
|
were performed using the {\sc oopse} molecular mechanics |
| 715 |
|
package.\cite{Meineke05} All other portions of the energy calculation |
| 716 |
|
(i.e. Lennard-Jones interactions) were handled in exactly the same |
| 725 |
|
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
| 726 |
|
increasing accuracy at the expense of computational time spent on the |
| 727 |
|
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
| 728 |
< |
The default TINKER tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 729 |
< |
in all SPME calculations, resulting in Ewald coefficients of 0.4200, |
| 728 |
> |
The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used |
| 729 |
> |
in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200, |
| 730 |
|
0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ |
| 731 |
|
respectively. |
| 732 |
|
|
| 736 |
|
In order to evaluate the performance of the pairwise electrostatic |
| 737 |
|
summation methods for Monte Carlo simulations, the energy differences |
| 738 |
|
between configurations were compared to the values obtained when using |
| 739 |
< |
SPME. The results for the subsequent regression analysis are shown in |
| 739 |
> |
{\sc spme}. The results for the subsequent regression analysis are shown in |
| 740 |
|
figure \ref{fig:delE}. |
| 741 |
|
|
| 742 |
|
\begin{figure} |
| 746 |
|
differences for a given electrostatic method compared with the |
| 747 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 748 |
|
indicate $\Delta E$ values indistinguishable from those obtained using |
| 749 |
< |
SPME. Different values of the cutoff radius are indicated with |
| 749 |
> |
{\sc spme}. Different values of the cutoff radius are indicated with |
| 750 |
|
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 751 |
|
inverted triangles).} |
| 752 |
|
\label{fig:delE} |
| 772 |
|
|
| 773 |
|
For the {\sc sp} method, inclusion of electrostatic damping improves |
| 774 |
|
the agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ |
| 775 |
< |
shows an excellent correlation and quality of fit with the SPME |
| 775 |
> |
shows an excellent correlation and quality of fit with the {\sc spme} |
| 776 |
|
results, particularly with a cutoff radius greater than 12 |
| 777 |
|
\AA . Use of a larger damping parameter is more helpful for the |
| 778 |
|
shortest cutoff shown, but it has a detrimental effect on simulations |
| 806 |
|
magnitudes for a given electrostatic method compared with the |
| 807 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 808 |
|
indicate force magnitude values indistinguishable from those obtained |
| 809 |
< |
using SPME. Different values of the cutoff radius are indicated with |
| 809 |
> |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 810 |
|
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 811 |
|
inverted triangles).} |
| 812 |
|
\label{fig:frcMag} |
| 813 |
|
\end{figure} |
| 814 |
|
|
| 815 |
|
Again, it is striking how well the Shifted Potential and Shifted Force |
| 816 |
< |
methods are doing at reproducing the SPME forces. The undamped and |
| 816 |
> |
methods are doing at reproducing the {\sc spme} forces. The undamped and |
| 817 |
|
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
| 818 |
|
This is perhaps expected because this method explicitly incorporates a |
| 819 |
|
smooth transition in the forces at the cutoff radius as well as the |
| 845 |
|
magnitudes for a given electrostatic method compared with the |
| 846 |
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
| 847 |
|
indicate torque magnitude values indistinguishable from those obtained |
| 848 |
< |
using SPME. Different values of the cutoff radius are indicated with |
| 848 |
> |
using {\sc spme}. Different values of the cutoff radius are indicated with |
| 849 |
|
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
| 850 |
|
inverted triangles).} |
| 851 |
|
\label{fig:trqMag} |
| 857 |
|
|
| 858 |
|
Torques appear to be much more sensitive to charges at a longer |
| 859 |
|
distance. The striking feature in comparing the new electrostatic |
| 860 |
< |
methods with SPME is how much the agreement improves with increasing |
| 860 |
> |
methods with {\sc spme} is how much the agreement improves with increasing |
| 861 |
|
cutoff radius. Again, the weakly damped and undamped {\sc sf} method |
| 862 |
< |
appears to be reproducing the SPME torques most accurately. |
| 862 |
> |
appears to be reproducing the {\sc spme} torques most accurately. |
| 863 |
|
|
| 864 |
|
Water molecules are dipolar, and the reaction field method reproduces |
| 865 |
|
the effect of the surrounding polarized medium on each of the |
| 874 |
|
vital in calculating dynamical quantities accurately. Force and |
| 875 |
|
torque directionalities were investigated by measuring the angles |
| 876 |
|
formed between these vectors and the same vectors calculated using |
| 877 |
< |
SPME. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 877 |
> |
{\sc spme}. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
| 878 |
|
variance ($\sigma^2$) of the Gaussian fits of the angle error |
| 879 |
|
distributions of the combined set over all system types. |
| 880 |
|
|
| 886 |
|
make with their counterparts obtained using the reference Ewald sum. |
| 887 |
|
Results with a variance ($\sigma^2$) equal to zero (dashed line) |
| 888 |
|
indicate force and torque directions indistinguishable from those |
| 889 |
< |
obtained using SPME. Different values of the cutoff radius are |
| 889 |
> |
obtained using {\sc spme}. Different values of the cutoff radius are |
| 890 |
|
indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, |
| 891 |
|
and 15\AA\ = inverted triangles).} |
| 892 |
|
\label{fig:frcTrqAng} |
| 895 |
|
Both the force and torque $\sigma^2$ results from the analysis of the |
| 896 |
|
total accumulated system data are tabulated in figure |
| 897 |
|
\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc |
| 898 |
< |
sp}) method would be essentially unusable for molecular dynamics until |
| 899 |
< |
the damping function is added. The Shifted Force ({\sc sf}) method, |
| 900 |
< |
however, is generating force and torque vectors which are within a few |
| 901 |
< |
degrees of the Ewald results even with weak (or no) damping. |
| 898 |
> |
sp}) method would be essentially unusable for molecular dynamics |
| 899 |
> |
unless the damping function is added. The Shifted Force ({\sc sf}) |
| 900 |
> |
method, however, is generating force and torque vectors which are |
| 901 |
> |
within a few degrees of the Ewald results even with weak (or no) |
| 902 |
> |
damping. |
| 903 |
|
|
| 904 |
|
All of the sets (aside from the over-damped case) show the improvement |
| 905 |
|
afforded by choosing a larger cutoff radius. Increasing the cutoff |
| 906 |
|
from 9 to 12 \AA\ typically results in a halving of the width of the |
| 907 |
< |
distribution, with a similar improvement going from 12 to 15 |
| 907 |
> |
distribution, with a similar improvement when going from 12 to 15 |
| 908 |
|
\AA . |
| 909 |
|
|
| 910 |
|
The undamped {\sc sf}, group-based cutoff, and reaction field methods |
| 911 |
|
all do equivalently well at capturing the direction of both the force |
| 912 |
< |
and torque vectors. Using damping improves the angular behavior |
| 913 |
< |
significantly for the {\sc sp} and moderately for the {\sc sf} |
| 914 |
< |
methods. Overdamping is detrimental to both methods. Again it is |
| 915 |
< |
important to recognize that the force vectors cover all particles in |
| 916 |
< |
the systems, while torque vectors are only available for neutral |
| 917 |
< |
molecular groups. Damping appears to have a more beneficial effect on |
| 912 |
> |
and torque vectors. Using the electrostatic damping improves the |
| 913 |
> |
angular behavior significantly for the {\sc sp} and moderately for the |
| 914 |
> |
{\sc sf} methods. Overdamping is detrimental to both methods. Again |
| 915 |
> |
it is important to recognize that the force vectors cover all |
| 916 |
> |
particles in all seven systems, while torque vectors are only |
| 917 |
> |
available for neutral molecular groups. Damping is more beneficial to |
| 918 |
|
charged bodies, and this observation is investigated further in the |
| 919 |
|
accompanying supporting information. |
| 920 |
|
|
| 921 |
|
Although not discussed previously, group based cutoffs can be applied |
| 922 |
< |
to both the {\sc sp} and {\sc sf} methods. Use of a switching |
| 923 |
< |
function corrects for the discontinuities that arise when atoms of the |
| 924 |
< |
two groups exit the cutoff radius before the group centers leave each |
| 925 |
< |
other's cutoff. Though there are no significant benefits or drawbacks |
| 926 |
< |
observed in $\Delta E$ and vector magnitude results when doing this, |
| 927 |
< |
there is a measurable improvement in the vector angle results. Table |
| 928 |
< |
\ref{tab:groupAngle} shows the angular variance values obtained using |
| 929 |
< |
group based cutoffs and a switching function alongside the results |
| 930 |
< |
seen in figure \ref{fig:frcTrqAng}. The {\sc sp} shows much narrower |
| 931 |
< |
angular distributions for both the force and torque vectors when using |
| 932 |
< |
an $\alpha$ of 0.2 \AA$^{-1}$ or less, while {\sc sf} shows |
| 933 |
< |
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. |
| 922 |
> |
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
| 923 |
> |
will reintroduce small discontinuities at the cutoff radius, but the |
| 924 |
> |
effects of these can be minimized by utilizing a switching function. |
| 925 |
> |
Though there are no significant benefits or drawbacks observed in |
| 926 |
> |
$\Delta E$ and the force and torque magnitudes when doing this, there |
| 927 |
> |
is a measurable improvement in the directionality of the forces and |
| 928 |
> |
torques. Table \ref{tab:groupAngle} shows the angular variances |
| 929 |
> |
obtained using group based cutoffs along with the results seen in |
| 930 |
> |
figure \ref{fig:frcTrqAng}. The {\sc sp} (with an $\alpha$ of 0.2 |
| 931 |
> |
\AA$^{-1}$ or smaller) shows much narrower angular distributions when |
| 932 |
> |
using group-based cutoffs. The {\sc sf} method likewise shows |
| 933 |
> |
improvement in the undamped and lightly damped cases. |
| 934 |
|
|
| 935 |
|
\begin{table}[htbp] |
| 936 |
< |
\centering |
| 937 |
< |
\caption{Variance ($\sigma^2$) of the force (top set) and torque |
| 938 |
< |
(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$.} |
| 936 |
> |
\centering |
| 937 |
> |
\caption{Statistical analysis of the angular |
| 938 |
> |
distributions that the force (upper) and torque (lower) vectors |
| 939 |
> |
from a given electrostatic method make with their counterparts |
| 940 |
> |
obtained using the reference Ewald sum. Calculations were |
| 941 |
> |
performed both with (Y) and without (N) group based cutoffs and a |
| 942 |
> |
switching function. The $\alpha$ values have units of \AA$^{-1}$ |
| 943 |
> |
and the variance values have units of degrees$^2$.} |
| 944 |
> |
|
| 945 |
|
\begin{tabular}{@{} ccrrrrrrrr @{}} |
| 946 |
|
\\ |
| 947 |
|
\toprule |
| 972 |
|
\label{tab:groupAngle} |
| 973 |
|
\end{table} |
| 974 |
|
|
| 975 |
< |
One additional trend to recognize in table \ref{tab:groupAngle} is |
| 976 |
< |
that the $\sigma^2$ values for both {\sc sp} and {\sc sf} converge as |
| 977 |
< |
$\alpha$ increases, something that is easier to see when using group |
| 978 |
< |
based cutoffs. The reason for this is that the complimentary error |
| 979 |
< |
function inserted into the potential weakens the electrostatic |
| 980 |
< |
interaction as $\alpha$ increases. Thus, at larger values of |
| 981 |
< |
$\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 |
| 975 |
> |
One additional trend in table \ref{tab:groupAngle} is that the |
| 976 |
> |
$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$ |
| 977 |
> |
increases, something that is more obvious with group-based cutoffs. |
| 978 |
> |
The complimentary error function inserted into the potential weakens |
| 979 |
> |
the electrostatic interaction as the value of $\alpha$ is increased. |
| 980 |
> |
However, at larger values of $\alpha$, it is possible to overdamp the |
| 981 |
> |
electrostatic interaction and to remove it completely. Kast |
| 982 |
|
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
| 983 |
|
values for these types of electrostatic summation methods by fitting |
| 984 |
|
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
| 986 |
|
respectively.\cite{Kast03} These appear to be reasonable choices to |
| 987 |
|
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
| 988 |
|
these findings, choices this high would introduce error in the |
| 989 |
< |
molecular torques, particularly for the shorter cutoffs. Based on the |
| 990 |
< |
above findings, empirical damping up to 0.2 \AA$^{-1}$ proves to be |
| 991 |
< |
beneficial, but damping may be unnecessary when using the {\sc sf} |
| 988 |
< |
method. |
| 989 |
> |
molecular torques, particularly for the shorter cutoffs. Based on our |
| 990 |
> |
observations, empirical damping up to 0.2 \AA$^{-1}$ is beneficial, |
| 991 |
> |
but damping may be unnecessary when using the {\sc sf} method. |
| 992 |
|
|
| 993 |
|
\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals} |
| 994 |
|
|
| 1001 |
|
distribution functions, diffusion constants, and velocity and |
| 1002 |
|
orientational correlation functions) may not be particularly sensitive |
| 1003 |
|
to the long-range and collective behavior that governs the |
| 1004 |
< |
low-frequency behavior in crystalline systems. |
| 1004 |
> |
low-frequency behavior in crystalline systems. Additionally, the |
| 1005 |
> |
ionic crystals are the worst case scenario for the pairwise methods |
| 1006 |
> |
because they lack the reciprocal space contribution contained in the |
| 1007 |
> |
Ewald summation. |
| 1008 |
|
|
| 1009 |
|
We are using two separate measures to probe the effects of these |
| 1010 |
|
alternative electrostatic methods on the dynamics in crystalline |
| 1016 |
|
\begin{figure} |
| 1017 |
|
\centering |
| 1018 |
|
\includegraphics[width = \linewidth]{./vCorrPlot.pdf} |
| 1019 |
< |
\caption{Velocity auto-correlation functions of NaCl crystals at |
| 1020 |
< |
1000 K using SPME, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc |
| 1019 |
> |
\caption{Velocity autocorrelation functions of NaCl crystals at |
| 1020 |
> |
1000 K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc |
| 1021 |
|
sp} ($\alpha$ = 0.2). The inset is a magnification of the area around |
| 1022 |
|
the first minimum. The times to first collision are nearly identical, |
| 1023 |
|
but differences can be seen in the peaks and troughs, where the |
| 1024 |
|
undamped and weakly damped methods are stiffer than the moderately |
| 1025 |
< |
damped and SPME methods.} |
| 1025 |
> |
damped and {\sc spme} methods.} |
| 1026 |
|
\label{fig:vCorrPlot} |
| 1027 |
|
\end{figure} |
| 1028 |
|
|
| 1029 |
< |
The short-time decays through the first collision are nearly identical |
| 1030 |
< |
in figure \ref{fig:vCorrPlot}, but the peaks and troughs of the |
| 1031 |
< |
functions show how the methods differ. The undamped {\sc sf} method |
| 1032 |
< |
has deeper troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher |
| 1033 |
< |
peaks than any of the other methods. As the damping function is |
| 1034 |
< |
increased, these peaks are smoothed out, and approach the SPME |
| 1035 |
< |
curve. The damping acts as a distance dependent Gaussian screening of |
| 1036 |
< |
the point charges for the pairwise summation methods; thus, the |
| 1037 |
< |
collisions are more elastic in the undamped {\sc sf} potential, and the |
| 1038 |
< |
stiffness of the potential is diminished as the electrostatic |
| 1039 |
< |
interactions are softened by the damping function. With $\alpha$ |
| 1040 |
< |
values of 0.2 \AA$^{-1}$, the {\sc sf} and {\sc sp} functions are |
| 1041 |
< |
nearly identical and track the SPME features quite well. This is not |
| 1042 |
< |
too surprising in that the differences between the {\sc sf} and {\sc |
| 1043 |
< |
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. |
| 1029 |
> |
The short-time decay of the velocity autocorrelation function through |
| 1030 |
> |
the first collision are nearly identical in figure |
| 1031 |
> |
\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show |
| 1032 |
> |
how the methods differ. The undamped {\sc sf} method has deeper |
| 1033 |
> |
troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than |
| 1034 |
> |
any of the other methods. As the damping parameter ($\alpha$) is |
| 1035 |
> |
increased, these peaks are smoothed out, and the {\sc sf} method |
| 1036 |
> |
approaches the {\sc spme} results. With $\alpha$ values of 0.2 \AA$^{-1}$, |
| 1037 |
> |
the {\sc sf} and {\sc sp} functions are nearly identical and track the |
| 1038 |
> |
{\sc spme} features quite well. This is not surprising because the {\sc sf} |
| 1039 |
> |
and {\sc sp} potentials become nearly identical with increased |
| 1040 |
> |
damping. However, this appears to indicate that once damping is |
| 1041 |
> |
utilized, the details of the form of the potential (and forces) |
| 1042 |
> |
constructed out of the damped electrostatic interaction are less |
| 1043 |
> |
important. |
| 1044 |
|
|
| 1045 |
|
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
| 1046 |
|
|
| 1047 |
< |
The short time dynamics were extended to evaluate how the differences |
| 1048 |
< |
between the methods affect the collective long-time motion. The same |
| 1049 |
< |
electrostatic summation methods were used as in the short time |
| 1050 |
< |
velocity autocorrelation function evaluation, but the trajectories |
| 1051 |
< |
were sampled over a much longer time. The power spectra of the |
| 1052 |
< |
resulting velocity autocorrelation functions were calculated and are |
| 1053 |
< |
displayed in figure \ref{fig:methodPS}. |
| 1047 |
> |
To evaluate how the differences between the methods affect the |
| 1048 |
> |
collective long-time motion, we computed power spectra from long-time |
| 1049 |
> |
traces of the velocity autocorrelation function. The power spectra for |
| 1050 |
> |
the best-performing alternative methods are shown in |
| 1051 |
> |
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
| 1052 |
> |
a cubic switching function between 40 and 50 ps was used to reduce the |
| 1053 |
> |
ringing resulting from data truncation. This procedure had no |
| 1054 |
> |
noticeable effect on peak location or magnitude. |
| 1055 |
|
|
| 1056 |
|
\begin{figure} |
| 1057 |
|
\centering |
| 1058 |
|
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
| 1059 |
|
\caption{Power spectra obtained from the velocity auto-correlation |
| 1060 |
< |
functions of NaCl crystals at 1000 K while using SPME, {\sc sf} |
| 1061 |
< |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). |
| 1062 |
< |
Apodization of the correlation functions via a cubic switching |
| 1063 |
< |
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.} |
| 1060 |
> |
functions of NaCl crystals at 1000 K while using {\sc spme}, {\sc sf} |
| 1061 |
> |
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset |
| 1062 |
> |
shows the frequency region below 100 cm$^{-1}$ to highlight where the |
| 1063 |
> |
spectra differ.} |
| 1064 |
|
\label{fig:methodPS} |
| 1065 |
|
\end{figure} |
| 1066 |
|
|
| 1067 |
< |
While high frequency peaks of the spectra in this figure overlap, |
| 1068 |
< |
showing the same general features, the low frequency region shows how |
| 1069 |
< |
the summation methods differ. Considering the low-frequency inset |
| 1070 |
< |
(expanded in the upper frame of figure \ref{fig:dampInc}), at |
| 1071 |
< |
frequencies below 100 cm$^{-1}$, the correlated motions are |
| 1072 |
< |
blue-shifted when using undamped or weakly damped {\sc sf}. When |
| 1073 |
< |
using moderate damping ($\alpha = 0.2$ \AA$^{-1}$) both the {\sc sf} |
| 1074 |
< |
and {\sc sp} methods give near identical correlated motion behavior as |
| 1075 |
< |
the Ewald method (which has a damping value of 0.3119). This |
| 1076 |
< |
weakening of the electrostatic interaction with increased damping |
| 1077 |
< |
explains why the long-ranged correlated motions are at lower |
| 1078 |
< |
frequencies for the moderately damped methods than for undamped or |
| 1079 |
< |
weakly damped methods. To see this effect more clearly, we show how |
| 1080 |
< |
damping strength alone affects a simple real-space electrostatic |
| 1081 |
< |
potential, |
| 1082 |
< |
\begin{equation} |
| 1083 |
< |
V_\textrm{damped}=\sum^N_i\sum^N_{j\ne i}q_iq_j\left[\frac{\textrm{erfc}({\alpha r})}{r}\right]S(r), |
| 1084 |
< |
\end{equation} |
| 1085 |
< |
where $S(r)$ is a switching function that smoothly zeroes the |
| 1086 |
< |
potential at the cutoff radius. Figure \ref{fig:dampInc} shows how |
| 1087 |
< |
the low frequency motions are dependent on the damping used in the |
| 1088 |
< |
direct electrostatic sum. As the damping increases, the peaks drop to |
| 1089 |
< |
lower frequencies. Incidentally, use of an $\alpha$ of 0.25 |
| 1090 |
< |
\AA$^{-1}$ on a simple electrostatic summation results in low |
| 1091 |
< |
frequency correlated dynamics equivalent to a simulation using SPME. |
| 1092 |
< |
When the coefficient lowers to 0.15 \AA$^{-1}$ and below, these peaks |
| 1093 |
< |
shift to higher frequency in exponential fashion. Though not shown, |
| 1094 |
< |
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. |
| 1067 |
> |
While the high frequency regions of the power spectra for the |
| 1068 |
> |
alternative methods are quantitatively identical with Ewald spectrum, |
| 1069 |
> |
the low frequency region shows how the summation methods differ. |
| 1070 |
> |
Considering the low-frequency inset (expanded in the upper frame of |
| 1071 |
> |
figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the |
| 1072 |
> |
correlated motions are blue-shifted when using undamped or weakly |
| 1073 |
> |
damped {\sc sf}. When using moderate damping ($\alpha = 0.2$ |
| 1074 |
> |
\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical |
| 1075 |
> |
correlated motion to the Ewald method (which has a convergence |
| 1076 |
> |
parameter of 0.3119 \AA$^{-1}$). This weakening of the electrostatic |
| 1077 |
> |
interaction with increased damping explains why the long-ranged |
| 1078 |
> |
correlated motions are at lower frequencies for the moderately damped |
| 1079 |
> |
methods than for undamped or weakly damped methods. |
| 1080 |
> |
|
| 1081 |
> |
To isolate the role of the damping constant, we have computed the |
| 1082 |
> |
spectra for a single method ({\sc sf}) with a range of damping |
| 1083 |
> |
constants and compared this with the {\sc spme} spectrum. |
| 1084 |
> |
Fig. \ref{fig:dampInc} shows more clearly that increasing the |
| 1085 |
> |
electrostatic damping red-shifts the lowest frequency phonon modes. |
| 1086 |
> |
However, even without any electrostatic damping, the {\sc sf} method |
| 1087 |
> |
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
| 1088 |
> |
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
| 1089 |
> |
would predict the lowest frequency peak near 325 cm$^{-1}$. {\it |
| 1090 |
> |
Most} of the collective behavior in the crystal is accurately captured |
| 1091 |
> |
using the {\sc sf} method. Quantitative agreement with Ewald can be |
| 1092 |
> |
obtained using moderate damping in addition to the shifting at the |
| 1093 |
> |
cutoff distance. |
| 1094 |
> |
|
| 1095 |
|
\begin{figure} |
| 1096 |
|
\centering |
| 1097 |
|
\includegraphics[width = \linewidth]{./comboSquare.pdf} |
| 1098 |
< |
\caption{Regions of spectra showing the low-frequency correlated |
| 1099 |
< |
motions for NaCl crystals at 1000 K using various electrostatic |
| 1100 |
< |
summation methods. The upper plot is a zoomed inset from figure |
| 1101 |
< |
\ref{fig:methodPS}. As the damping value for the {\sc sf} potential |
| 1102 |
< |
increases, the low-frequency peaks red-shift. The lower plot is of |
| 1103 |
< |
spectra when using SPME and a simple damped Coulombic sum with damping |
| 1104 |
< |
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.} |
| 1098 |
> |
\caption{Effect of damping on the two lowest-frequency phonon modes in |
| 1099 |
> |
the NaCl crystal at 1000K. The undamped shifted force ({\sc sf}) |
| 1100 |
> |
method is off by less than 10 cm$^{-1}$, and increasing the |
| 1101 |
> |
electrostatic damping to 0.25 \AA$^{-1}$ gives quantitative agreement |
| 1102 |
> |
with the power spectrum obtained using the Ewald sum. Overdamping can |
| 1103 |
> |
result in underestimates of frequencies of the long-wavelength |
| 1104 |
> |
motions.} |
| 1105 |
|
\label{fig:dampInc} |
| 1106 |
|
\end{figure} |
| 1107 |
|
|
| 1108 |
|
\section{Conclusions} |
| 1109 |
|
|
| 1110 |
|
This investigation of pairwise electrostatic summation techniques |
| 1111 |
< |
shows that there are viable and more computationally efficient |
| 1112 |
< |
electrostatic summation techniques than the Ewald summation, chiefly |
| 1113 |
< |
methods derived from the damped Coulombic sum originally proposed by |
| 1114 |
< |
Wolf \textit{et al.}\cite{Wolf99,Zahn02} In particular, the |
| 1115 |
< |
{\sc sf} method, reformulated above as eq. (\ref{eq:DSFPot}), |
| 1116 |
< |
shows a remarkable ability to reproduce the energetic and dynamic |
| 1117 |
< |
characteristics exhibited by simulations employing lattice summation |
| 1118 |
< |
techniques. The cumulative energy difference results showed the |
| 1119 |
< |
undamped {\sc sf} and moderately damped {\sc sp} methods |
| 1120 |
< |
produced results nearly identical to SPME. Similarly for the dynamic |
| 1121 |
< |
features, the undamped or moderately damped {\sc sf} and |
| 1122 |
< |
moderately damped {\sc sp} methods produce force and torque |
| 1123 |
< |
vector magnitude and directions very similar to the expected values. |
| 1124 |
< |
These results translate into long-time dynamic behavior equivalent to |
| 1125 |
< |
that produced in simulations using SPME. |
| 1111 |
> |
shows that there are viable and computationally efficient alternatives |
| 1112 |
> |
to the Ewald summation. These methods are derived from the damped and |
| 1113 |
> |
cutoff-neutralized Coulombic sum originally proposed by Wolf |
| 1114 |
> |
\textit{et al.}\cite{Wolf99} In particular, the {\sc sf} |
| 1115 |
> |
method, reformulated above as eqs. (\ref{eq:DSFPot}) and |
| 1116 |
> |
(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the |
| 1117 |
> |
energetic and dynamic characteristics exhibited by simulations |
| 1118 |
> |
employing lattice summation techniques. The cumulative energy |
| 1119 |
> |
difference results showed the undamped {\sc sf} and moderately damped |
| 1120 |
> |
{\sc sp} methods produced results nearly identical to {\sc spme}. Similarly |
| 1121 |
> |
for the dynamic features, the undamped or moderately damped {\sc sf} |
| 1122 |
> |
and moderately damped {\sc sp} methods produce force and torque vector |
| 1123 |
> |
magnitude and directions very similar to the expected values. These |
| 1124 |
> |
results translate into long-time dynamic behavior equivalent to that |
| 1125 |
> |
produced in simulations using {\sc spme}. |
| 1126 |
|
|
| 1127 |
+ |
As in all purely-pairwise cutoff methods, these methods are expected |
| 1128 |
+ |
to scale approximately {\it linearly} with system size, and they are |
| 1129 |
+ |
easily parallelizable. This should result in substantial reductions |
| 1130 |
+ |
in the computational cost of performing large simulations. |
| 1131 |
+ |
|
| 1132 |
|
Aside from the computational cost benefit, these techniques have |
| 1133 |
|
applicability in situations where the use of the Ewald sum can prove |
| 1134 |
< |
problematic. Primary among them is their use in interfacial systems, |
| 1135 |
< |
where the unmodified lattice sum techniques artificially accentuate |
| 1136 |
< |
the periodicity of the system in an undesirable manner. There have |
| 1137 |
< |
been alterations to the standard Ewald techniques, via corrections and |
| 1138 |
< |
reformulations, to compensate for these systems; but the pairwise |
| 1139 |
< |
techniques discussed here require no modifications, making them |
| 1140 |
< |
natural tools to tackle these problems. Additionally, this |
| 1141 |
< |
transferability gives them benefits over other pairwise methods, like |
| 1142 |
< |
reaction field, because estimations of physical properties (e.g. the |
| 1143 |
< |
dielectric constant) are unnecessary. |
| 1134 |
> |
problematic. Of greatest interest is their potential use in |
| 1135 |
> |
interfacial systems, where the unmodified lattice sum techniques |
| 1136 |
> |
artificially accentuate the periodicity of the system in an |
| 1137 |
> |
undesirable manner. There have been alterations to the standard Ewald |
| 1138 |
> |
techniques, via corrections and reformulations, to compensate for |
| 1139 |
> |
these systems; but the pairwise techniques discussed here require no |
| 1140 |
> |
modifications, making them natural tools to tackle these problems. |
| 1141 |
> |
Additionally, this transferability gives them benefits over other |
| 1142 |
> |
pairwise methods, like reaction field, because estimations of physical |
| 1143 |
> |
properties (e.g. the dielectric constant) are unnecessary. |
| 1144 |
|
|
| 1145 |
< |
We are not suggesting any flaw with the Ewald sum; in fact, it is the |
| 1146 |
< |
standard by which these simple pairwise sums are judged. However, |
| 1147 |
< |
these results do suggest that in the typical simulations performed |
| 1148 |
< |
today, the Ewald summation may no longer be required to obtain the |
| 1149 |
< |
level of accuracy most researchers have come to expect |
| 1145 |
> |
If a researcher is using Monte Carlo simulations of large chemical |
| 1146 |
> |
systems containing point charges, most structural features will be |
| 1147 |
> |
accurately captured using the undamped {\sc sf} method or the {\sc sp} |
| 1148 |
> |
method with an electrostatic damping of 0.2 \AA$^{-1}$. These methods |
| 1149 |
> |
would also be appropriate for molecular dynamics simulations where the |
| 1150 |
> |
data of interest is either structural or short-time dynamical |
| 1151 |
> |
quantities. For long-time dynamics and collective motions, the safest |
| 1152 |
> |
pairwise method we have evaluated is the {\sc sf} method with an |
| 1153 |
> |
electrostatic damping between 0.2 and 0.25 |
| 1154 |
> |
\AA$^{-1}$. |
| 1155 |
|
|
| 1156 |
+ |
We are not suggesting that there is any flaw with the Ewald sum; in |
| 1157 |
+ |
fact, it is the standard by which these simple pairwise sums have been |
| 1158 |
+ |
judged. However, these results do suggest that in the typical |
| 1159 |
+ |
simulations performed today, the Ewald summation may no longer be |
| 1160 |
+ |
required to obtain the level of accuracy most researchers have come to |
| 1161 |
+ |
expect. |
| 1162 |
+ |
|
| 1163 |
|
\section{Acknowledgments} |
| 1164 |
+ |
Support for this project was provided by the National Science |
| 1165 |
+ |
Foundation under grant CHE-0134881. The authors would like to thank |
| 1166 |
+ |
Steve Corcelli and Ed Maginn for helpful discussions and comments. |
| 1167 |
+ |
|
| 1168 |
|
\newpage |
| 1169 |
|
|
| 1170 |
|
\bibliographystyle{jcp2} |
