1 |
chrisfen |
2575 |
%\documentclass[prb,aps,twocolumn,tabularx]{revtex4} |
2 |
gezelter |
2617 |
%\documentclass[aps,prb,preprint]{revtex4} |
3 |
chrisfen |
2741 |
\documentclass[11pt]{article} |
4 |
|
|
\usepackage{endfloat} |
5 |
chrisfen |
2636 |
\usepackage{amsmath,bm} |
6 |
chrisfen |
2594 |
\usepackage{amssymb} |
7 |
chrisfen |
2575 |
\usepackage{epsf} |
8 |
|
|
\usepackage{times} |
9 |
gezelter |
2617 |
\usepackage{mathptmx} |
10 |
chrisfen |
2575 |
\usepackage{setspace} |
11 |
|
|
\usepackage{tabularx} |
12 |
|
|
\usepackage{graphicx} |
13 |
chrisfen |
2595 |
\usepackage{booktabs} |
14 |
chrisfen |
2605 |
\usepackage{bibentry} |
15 |
|
|
\usepackage{mathrsfs} |
16 |
chrisfen |
2575 |
\usepackage[ref]{overcite} |
17 |
|
|
\pagestyle{plain} |
18 |
|
|
\pagenumbering{arabic} |
19 |
|
|
\oddsidemargin 0.0cm \evensidemargin 0.0cm |
20 |
|
|
\topmargin -21pt \headsep 10pt |
21 |
|
|
\textheight 9.0in \textwidth 6.5in |
22 |
|
|
\brokenpenalty=10000 |
23 |
chrisfen |
2742 |
%\renewcommand{\baselinestretch}{1.2} |
24 |
|
|
\renewcommand{\baselinestretch}{2} |
25 |
chrisfen |
2575 |
\renewcommand\citemid{\ } % no comma in optional reference note |
26 |
chrisfen |
2742 |
\AtBeginDelayedFloats{\renewcommand{\baselinestretch}{2}} %doublespace captions |
27 |
|
|
\let\Caption\caption |
28 |
|
|
\renewcommand\caption[1]{% |
29 |
|
|
\Caption[#1]{}% |
30 |
|
|
} |
31 |
gezelter |
2744 |
\setlength{\abovecaptionskip}{1.2in} |
32 |
|
|
\setlength{\belowcaptionskip}{1.2in} |
33 |
chrisfen |
2575 |
|
34 |
|
|
\begin{document} |
35 |
|
|
|
36 |
chrisfen |
2667 |
\title{Is the Ewald summation still necessary? \\ |
37 |
gezelter |
2744 |
Pairwise alternatives to the accepted standard \\ |
38 |
|
|
for long-range electrostatics} |
39 |
chrisfen |
2575 |
|
40 |
gezelter |
2617 |
\author{Christopher J. Fennell and J. Daniel Gezelter\footnote{Corresponding author. \ Electronic mail: |
41 |
|
|
gezelter@nd.edu} \\ |
42 |
chrisfen |
2575 |
Department of Chemistry and Biochemistry\\ |
43 |
|
|
University of Notre Dame\\ |
44 |
|
|
Notre Dame, Indiana 46556} |
45 |
|
|
|
46 |
|
|
\date{\today} |
47 |
|
|
|
48 |
|
|
\maketitle |
49 |
chrisfen |
2742 |
%\doublespacing |
50 |
gezelter |
2617 |
|
51 |
chrisfen |
2575 |
\begin{abstract} |
52 |
gezelter |
2656 |
We investigate pairwise electrostatic interaction methods and show |
53 |
|
|
that there are viable and computationally efficient $(\mathscr{O}(N))$ |
54 |
|
|
alternatives to the Ewald summation for typical modern molecular |
55 |
|
|
simulations. These methods are extended from the damped and |
56 |
chrisfen |
2667 |
cutoff-neutralized Coulombic sum originally proposed by |
57 |
|
|
[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 |
58 |
gezelter |
2656 |
a remarkable ability to reproduce the energetic and dynamic |
59 |
|
|
characteristics exhibited by simulations employing lattice summation |
60 |
|
|
techniques. Comparisons were performed with this and other pairwise |
61 |
chrisfen |
2667 |
methods against the smooth particle mesh Ewald ({\sc spme}) summation |
62 |
|
|
to see how well they reproduce the energetics and dynamics of a |
63 |
gezelter |
2744 |
variety of molecular simulations. |
64 |
chrisfen |
2575 |
\end{abstract} |
65 |
|
|
|
66 |
gezelter |
2617 |
\newpage |
67 |
|
|
|
68 |
chrisfen |
2575 |
%\narrowtext |
69 |
|
|
|
70 |
chrisfen |
2620 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
71 |
chrisfen |
2575 |
% BODY OF TEXT |
72 |
chrisfen |
2620 |
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
73 |
chrisfen |
2575 |
|
74 |
|
|
\section{Introduction} |
75 |
|
|
|
76 |
gezelter |
2617 |
In molecular simulations, proper accumulation of the electrostatic |
77 |
gezelter |
2643 |
interactions is essential and is one of the most |
78 |
|
|
computationally-demanding tasks. The common molecular mechanics force |
79 |
|
|
fields represent atomic sites with full or partial charges protected |
80 |
|
|
by Lennard-Jones (short range) interactions. This means that nearly |
81 |
|
|
every pair interaction involves a calculation of charge-charge forces. |
82 |
|
|
Coupled with relatively long-ranged $r^{-1}$ decay, the monopole |
83 |
|
|
interactions quickly become the most expensive part of molecular |
84 |
|
|
simulations. Historically, the electrostatic pair interaction would |
85 |
|
|
not have decayed appreciably within the typical box lengths that could |
86 |
|
|
be feasibly simulated. In the larger systems that are more typical of |
87 |
|
|
modern simulations, large cutoffs should be used to incorporate |
88 |
|
|
electrostatics correctly. |
89 |
chrisfen |
2604 |
|
90 |
gezelter |
2643 |
There have been many efforts to address the proper and practical |
91 |
|
|
handling of electrostatic interactions, and these have resulted in a |
92 |
|
|
variety of techniques.\cite{Roux99,Sagui99,Tobias01} These are |
93 |
|
|
typically classified as implicit methods (i.e., continuum dielectrics, |
94 |
|
|
static dipolar fields),\cite{Born20,Grossfield00} explicit methods |
95 |
|
|
(i.e., Ewald summations, interaction shifting or |
96 |
chrisfen |
2640 |
truncation),\cite{Ewald21,Steinbach94} or a mixture of the two (i.e., |
97 |
chrisfen |
2639 |
reaction field type methods, fast multipole |
98 |
|
|
methods).\cite{Onsager36,Rokhlin85} The explicit or mixed methods are |
99 |
gezelter |
2643 |
often preferred because they physically incorporate solvent molecules |
100 |
|
|
in the system of interest, but these methods are sometimes difficult |
101 |
|
|
to utilize because of their high computational cost.\cite{Roux99} In |
102 |
|
|
addition to the computational cost, there have been some questions |
103 |
|
|
regarding possible artifacts caused by the inherent periodicity of the |
104 |
|
|
explicit Ewald summation.\cite{Tobias01} |
105 |
chrisfen |
2639 |
|
106 |
chrisfen |
2667 |
In this paper, we focus on a new set of pairwise methods devised by |
107 |
gezelter |
2643 |
Wolf {\it et al.},\cite{Wolf99} which we further extend. These |
108 |
|
|
methods along with a few other mixed methods (i.e. reaction field) are |
109 |
|
|
compared with the smooth particle mesh Ewald |
110 |
|
|
sum,\cite{Onsager36,Essmann99} which is our reference method for |
111 |
|
|
handling long-range electrostatic interactions. The new methods for |
112 |
|
|
handling electrostatics have the potential to scale linearly with |
113 |
|
|
increasing system size since they involve only a simple modification |
114 |
|
|
to the direct pairwise sum. They also lack the added periodicity of |
115 |
|
|
the Ewald sum, so they can be used for systems which are non-periodic |
116 |
|
|
or which have one- or two-dimensional periodicity. Below, these |
117 |
chrisfen |
2740 |
methods are evaluated using a variety of model systems to |
118 |
|
|
establish their usability in molecular simulations. |
119 |
chrisfen |
2639 |
|
120 |
chrisfen |
2608 |
\subsection{The Ewald Sum} |
121 |
chrisfen |
2667 |
The complete accumulation of the electrostatic interactions in a system with |
122 |
chrisfen |
2639 |
periodic boundary conditions (PBC) requires the consideration of the |
123 |
gezelter |
2643 |
effect of all charges within a (cubic) simulation box as well as those |
124 |
|
|
in the periodic replicas, |
125 |
chrisfen |
2636 |
\begin{equation} |
126 |
|
|
V_\textrm{elec} = \frac{1}{2} {\sum_{\mathbf{n}}}^\prime \left[ \sum_{i=1}^N\sum_{j=1}^N \phi\left( \mathbf{r}_{ij} + L\mathbf{n},\bm{\Omega}_i,\bm{\Omega}_j\right) \right], |
127 |
|
|
\label{eq:PBCSum} |
128 |
|
|
\end{equation} |
129 |
chrisfen |
2639 |
where the sum over $\mathbf{n}$ is a sum over all periodic box |
130 |
|
|
replicas with integer coordinates $\mathbf{n} = (l,m,n)$, and the |
131 |
|
|
prime indicates $i = j$ are neglected for $\mathbf{n} = |
132 |
|
|
0$.\cite{deLeeuw80} Within the sum, $N$ is the number of electrostatic |
133 |
|
|
particles, $\mathbf{r}_{ij}$ is $\mathbf{r}_j - \mathbf{r}_i$, $L$ is |
134 |
|
|
the cell length, $\bm{\Omega}_{i,j}$ are the Euler angles for $i$ and |
135 |
gezelter |
2643 |
$j$, and $\phi$ is the solution to Poisson's equation |
136 |
|
|
($\phi(\mathbf{r}_{ij}) = q_i q_j |\mathbf{r}_{ij}|^{-1}$ for |
137 |
|
|
charge-charge interactions). In the case of monopole electrostatics, |
138 |
|
|
eq. (\ref{eq:PBCSum}) is conditionally convergent and is divergent for |
139 |
|
|
non-neutral systems. |
140 |
chrisfen |
2604 |
|
141 |
gezelter |
2643 |
The electrostatic summation problem was originally studied by Ewald |
142 |
chrisfen |
2636 |
for the case of an infinite crystal.\cite{Ewald21}. The approach he |
143 |
|
|
took was to convert this conditionally convergent sum into two |
144 |
|
|
absolutely convergent summations: a short-ranged real-space summation |
145 |
|
|
and a long-ranged reciprocal-space summation, |
146 |
|
|
\begin{equation} |
147 |
|
|
\begin{split} |
148 |
chrisfen |
2637 |
V_\textrm{elec} = \frac{1}{2}& \sum_{i=1}^N\sum_{j=1}^N \Biggr[ q_i q_j\Biggr( {\sum_{\mathbf{n}}}^\prime \frac{\textrm{erfc}\left( \alpha|\mathbf{r}_{ij}+\mathbf{n}|\right)}{|\mathbf{r}_{ij}+\mathbf{n}|} \\ &+ \frac{1}{\pi L^3}\sum_{\mathbf{k}\ne 0}\frac{4\pi^2}{|\mathbf{k}|^2} \exp{\left(-\frac{\pi^2|\mathbf{k}|^2}{\alpha^2}\right) \cos\left(\mathbf{k}\cdot\mathbf{r}_{ij}\right)}\Biggr)\Biggr] \\ &- \frac{\alpha}{\pi^{1/2}}\sum_{i=1}^N q_i^2 + \frac{2\pi}{(2\epsilon_\textrm{S}+1)L^3}\left|\sum_{i=1}^N q_i\mathbf{r}_i\right|^2, |
149 |
chrisfen |
2636 |
\end{split} |
150 |
|
|
\label{eq:EwaldSum} |
151 |
|
|
\end{equation} |
152 |
chrisfen |
2649 |
where $\alpha$ is the damping or convergence parameter with units of |
153 |
|
|
\AA$^{-1}$, $\mathbf{k}$ are the reciprocal vectors and are equal to |
154 |
|
|
$2\pi\mathbf{n}/L^2$, and $\epsilon_\textrm{S}$ is the dielectric |
155 |
|
|
constant of the surrounding medium. The final two terms of |
156 |
chrisfen |
2636 |
eq. (\ref{eq:EwaldSum}) are a particle-self term and a dipolar term |
157 |
|
|
for interacting with a surrounding dielectric.\cite{Allen87} This |
158 |
|
|
dipolar term was neglected in early applications in molecular |
159 |
|
|
simulations,\cite{Brush66,Woodcock71} until it was introduced by de |
160 |
|
|
Leeuw {\it et al.} to address situations where the unit cell has a |
161 |
gezelter |
2643 |
dipole moment which is magnified through replication of the periodic |
162 |
|
|
images.\cite{deLeeuw80,Smith81} If this term is taken to be zero, the |
163 |
|
|
system is said to be using conducting (or ``tin-foil'') boundary |
164 |
chrisfen |
2637 |
conditions, $\epsilon_{\rm S} = \infty$. Figure |
165 |
chrisfen |
2636 |
\ref{fig:ewaldTime} shows how the Ewald sum has been applied over |
166 |
gezelter |
2653 |
time. Initially, due to the small system sizes that could be |
167 |
|
|
simulated feasibly, the entire simulation box was replicated to |
168 |
|
|
convergence. In more modern simulations, the systems have grown large |
169 |
|
|
enough that a real-space cutoff could potentially give convergent |
170 |
|
|
behavior. Indeed, it has been observed that with the choice of a |
171 |
|
|
small $\alpha$, the reciprocal-space portion of the Ewald sum can be |
172 |
|
|
rapidly convergent and small relative to the real-space |
173 |
|
|
portion.\cite{Karasawa89,Kolafa92} |
174 |
gezelter |
2643 |
|
175 |
chrisfen |
2610 |
\begin{figure} |
176 |
|
|
\centering |
177 |
gezelter |
2656 |
\includegraphics[width = \linewidth]{./ewaldProgression.pdf} |
178 |
gezelter |
2669 |
\caption{The change in the need for the Ewald sum with |
179 |
|
|
increasing computational power. A:~Initially, only small systems |
180 |
|
|
could be studied, and the Ewald sum replicated the simulation box to |
181 |
|
|
convergence. B:~Now, radial cutoff methods should be able to reach |
182 |
|
|
convergence for the larger systems of charges that are common today.} |
183 |
chrisfen |
2610 |
\label{fig:ewaldTime} |
184 |
|
|
\end{figure} |
185 |
|
|
|
186 |
gezelter |
2643 |
The original Ewald summation is an $\mathscr{O}(N^2)$ algorithm. The |
187 |
chrisfen |
2649 |
convergence parameter $(\alpha)$ plays an important role in balancing |
188 |
gezelter |
2643 |
the computational cost between the direct and reciprocal-space |
189 |
|
|
portions of the summation. The choice of this value allows one to |
190 |
|
|
select whether the real-space or reciprocal space portion of the |
191 |
|
|
summation is an $\mathscr{O}(N^2)$ calculation (with the other being |
192 |
|
|
$\mathscr{O}(N)$).\cite{Sagui99} With the appropriate choice of |
193 |
|
|
$\alpha$ and thoughtful algorithm development, this cost can be |
194 |
|
|
reduced to $\mathscr{O}(N^{3/2})$.\cite{Perram88} The typical route |
195 |
|
|
taken to reduce the cost of the Ewald summation even further is to set |
196 |
|
|
$\alpha$ such that the real-space interactions decay rapidly, allowing |
197 |
|
|
for a short spherical cutoff. Then the reciprocal space summation is |
198 |
|
|
optimized. These optimizations usually involve utilization of the |
199 |
|
|
fast Fourier transform (FFT),\cite{Hockney81} leading to the |
200 |
chrisfen |
2637 |
particle-particle particle-mesh (P3M) and particle mesh Ewald (PME) |
201 |
|
|
methods.\cite{Shimada93,Luty94,Luty95,Darden93,Essmann95} In these |
202 |
|
|
methods, the cost of the reciprocal-space portion of the Ewald |
203 |
gezelter |
2643 |
summation is reduced from $\mathscr{O}(N^2)$ down to $\mathscr{O}(N |
204 |
|
|
\log N)$. |
205 |
chrisfen |
2636 |
|
206 |
gezelter |
2643 |
These developments and optimizations have made the use of the Ewald |
207 |
|
|
summation routine in simulations with periodic boundary |
208 |
|
|
conditions. However, in certain systems, such as vapor-liquid |
209 |
|
|
interfaces and membranes, the intrinsic three-dimensional periodicity |
210 |
|
|
can prove problematic. The Ewald sum has been reformulated to handle |
211 |
gezelter |
2744 |
2-D systems,\cite{Parry75,Parry76,Heyes77,deLeeuw79,Rhee89} but these |
212 |
|
|
methods are computationally expensive.\cite{Spohr97,Yeh99} More |
213 |
chrisfen |
2740 |
recently, there have been several successful efforts toward reducing |
214 |
gezelter |
2744 |
the computational cost of 2-D lattice |
215 |
|
|
summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} |
216 |
|
|
bringing them more in line with the cost of the full 3-D summation. |
217 |
chrisfen |
2637 |
|
218 |
gezelter |
2744 |
|
219 |
chrisfen |
2637 |
Several studies have recognized that the inherent periodicity in the |
220 |
gezelter |
2643 |
Ewald sum can also have an effect on three-dimensional |
221 |
|
|
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00} |
222 |
|
|
Solvated proteins are essentially kept at high concentration due to |
223 |
|
|
the periodicity of the electrostatic summation method. In these |
224 |
chrisfen |
2637 |
systems, the more compact folded states of a protein can be |
225 |
|
|
artificially stabilized by the periodic replicas introduced by the |
226 |
gezelter |
2643 |
Ewald summation.\cite{Weber00} Thus, care must be taken when |
227 |
|
|
considering the use of the Ewald summation where the assumed |
228 |
|
|
periodicity would introduce spurious effects in the system dynamics. |
229 |
chrisfen |
2637 |
|
230 |
chrisfen |
2608 |
\subsection{The Wolf and Zahn Methods} |
231 |
gezelter |
2617 |
In a recent paper by Wolf \textit{et al.}, a procedure was outlined |
232 |
gezelter |
2624 |
for the accurate accumulation of electrostatic interactions in an |
233 |
gezelter |
2643 |
efficient pairwise fashion. This procedure lacks the inherent |
234 |
|
|
periodicity of the Ewald summation.\cite{Wolf99} Wolf \textit{et al.} |
235 |
|
|
observed that the electrostatic interaction is effectively |
236 |
|
|
short-ranged in condensed phase systems and that neutralization of the |
237 |
|
|
charge contained within the cutoff radius is crucial for potential |
238 |
|
|
stability. They devised a pairwise summation method that ensures |
239 |
|
|
charge neutrality and gives results similar to those obtained with the |
240 |
chrisfen |
2667 |
Ewald summation. The resulting shifted Coulomb potential includes |
241 |
|
|
image-charges subtracted out through placement on the cutoff sphere |
242 |
|
|
and a distance-dependent damping function (identical to that seen in |
243 |
|
|
the real-space portion of the Ewald sum) to aid convergence |
244 |
chrisfen |
2601 |
\begin{equation} |
245 |
chrisfen |
2640 |
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\}. |
246 |
chrisfen |
2601 |
\label{eq:WolfPot} |
247 |
|
|
\end{equation} |
248 |
gezelter |
2617 |
Eq. (\ref{eq:WolfPot}) is essentially the common form of a shifted |
249 |
|
|
potential. However, neutralizing the charge contained within each |
250 |
|
|
cutoff sphere requires the placement of a self-image charge on the |
251 |
|
|
surface of the cutoff sphere. This additional self-term in the total |
252 |
gezelter |
2624 |
potential enabled Wolf {\it et al.} to obtain excellent estimates of |
253 |
gezelter |
2617 |
Madelung energies for many crystals. |
254 |
|
|
|
255 |
|
|
In order to use their charge-neutralized potential in molecular |
256 |
|
|
dynamics simulations, Wolf \textit{et al.} suggested taking the |
257 |
|
|
derivative of this potential prior to evaluation of the limit. This |
258 |
|
|
procedure gives an expression for the forces, |
259 |
chrisfen |
2601 |
\begin{equation} |
260 |
chrisfen |
2636 |
F_{\textrm{Wolf}}(r_{ij}) = q_i q_j\left\{\left[\frac{\textrm{erfc}\left(\alpha r_{ij}\right)}{r^2_{ij}}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r_{ij}^2\right)}}{r_{ij}}\right]-\left[\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right]\right\}, |
261 |
chrisfen |
2601 |
\label{eq:WolfForces} |
262 |
|
|
\end{equation} |
263 |
gezelter |
2617 |
that incorporates both image charges and damping of the electrostatic |
264 |
|
|
interaction. |
265 |
|
|
|
266 |
|
|
More recently, Zahn \textit{et al.} investigated these potential and |
267 |
|
|
force expressions for use in simulations involving water.\cite{Zahn02} |
268 |
gezelter |
2624 |
In their work, they pointed out that the forces and derivative of |
269 |
|
|
the potential are not commensurate. Attempts to use both |
270 |
gezelter |
2643 |
eqs. (\ref{eq:WolfPot}) and (\ref{eq:WolfForces}) together will lead |
271 |
gezelter |
2624 |
to poor energy conservation. They correctly observed that taking the |
272 |
|
|
limit shown in equation (\ref{eq:WolfPot}) {\it after} calculating the |
273 |
|
|
derivatives gives forces for a different potential energy function |
274 |
gezelter |
2643 |
than the one shown in eq. (\ref{eq:WolfPot}). |
275 |
gezelter |
2617 |
|
276 |
gezelter |
2643 |
Zahn \textit{et al.} introduced a modified form of this summation |
277 |
|
|
method as a way to use the technique in Molecular Dynamics |
278 |
|
|
simulations. They proposed a new damped Coulomb potential, |
279 |
chrisfen |
2601 |
\begin{equation} |
280 |
gezelter |
2643 |
V_{\textrm{Zahn}}(r_{ij}) = q_iq_j\left\{\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}-\left[\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right]\left(r_{ij}-R_\mathrm{c}\right)\right\}, |
281 |
chrisfen |
2601 |
\label{eq:ZahnPot} |
282 |
|
|
\end{equation} |
283 |
gezelter |
2643 |
and showed that this potential does fairly well at capturing the |
284 |
gezelter |
2617 |
structural and dynamic properties of water compared the same |
285 |
|
|
properties obtained using the Ewald sum. |
286 |
chrisfen |
2601 |
|
287 |
chrisfen |
2608 |
\subsection{Simple Forms for Pairwise Electrostatics} |
288 |
|
|
|
289 |
gezelter |
2617 |
The potentials proposed by Wolf \textit{et al.} and Zahn \textit{et |
290 |
|
|
al.} are constructed using two different (and separable) computational |
291 |
gezelter |
2624 |
tricks: \begin{enumerate} |
292 |
gezelter |
2617 |
\item shifting through the use of image charges, and |
293 |
|
|
\item damping the electrostatic interaction. |
294 |
gezelter |
2624 |
\end{enumerate} Wolf \textit{et al.} treated the |
295 |
gezelter |
2617 |
development of their summation method as a progressive application of |
296 |
|
|
these techniques,\cite{Wolf99} while Zahn \textit{et al.} founded |
297 |
|
|
their damped Coulomb modification (Eq. (\ref{eq:ZahnPot})) on the |
298 |
|
|
post-limit forces (Eq. (\ref{eq:WolfForces})) which were derived using |
299 |
|
|
both techniques. It is possible, however, to separate these |
300 |
|
|
tricks and study their effects independently. |
301 |
|
|
|
302 |
|
|
Starting with the original observation that the effective range of the |
303 |
|
|
electrostatic interaction in condensed phases is considerably less |
304 |
|
|
than $r^{-1}$, either the cutoff sphere neutralization or the |
305 |
|
|
distance-dependent damping technique could be used as a foundation for |
306 |
|
|
a new pairwise summation method. Wolf \textit{et al.} made the |
307 |
|
|
observation that charge neutralization within the cutoff sphere plays |
308 |
|
|
a significant role in energy convergence; therefore we will begin our |
309 |
|
|
analysis with the various shifted forms that maintain this charge |
310 |
|
|
neutralization. We can evaluate the methods of Wolf |
311 |
|
|
\textit{et al.} and Zahn \textit{et al.} by considering the standard |
312 |
|
|
shifted potential, |
313 |
chrisfen |
2601 |
\begin{equation} |
314 |
gezelter |
2643 |
V_\textrm{SP}(r) = \begin{cases} |
315 |
gezelter |
2617 |
v(r)-v_\textrm{c} &\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > |
316 |
|
|
R_\textrm{c} |
317 |
|
|
\end{cases}, |
318 |
|
|
\label{eq:shiftingPotForm} |
319 |
|
|
\end{equation} |
320 |
|
|
and shifted force, |
321 |
|
|
\begin{equation} |
322 |
gezelter |
2643 |
V_\textrm{SF}(r) = \begin{cases} |
323 |
gezelter |
2624 |
v(r)-v_\textrm{c}-\left(\frac{d v(r)}{d r}\right)_{r=R_\textrm{c}}(r-R_\textrm{c}) |
324 |
gezelter |
2617 |
&\quad r\leqslant R_\textrm{c} \\ 0 &\quad r > R_\textrm{c} |
325 |
chrisfen |
2601 |
\end{cases}, |
326 |
chrisfen |
2612 |
\label{eq:shiftingForm} |
327 |
chrisfen |
2601 |
\end{equation} |
328 |
gezelter |
2617 |
functions where $v(r)$ is the unshifted form of the potential, and |
329 |
|
|
$v_c$ is $v(R_\textrm{c})$. The Shifted Force ({\sc sf}) form ensures |
330 |
|
|
that both the potential and the forces goes to zero at the cutoff |
331 |
|
|
radius, while the Shifted Potential ({\sc sp}) form only ensures the |
332 |
|
|
potential is smooth at the cutoff radius |
333 |
|
|
($R_\textrm{c}$).\cite{Allen87} |
334 |
|
|
|
335 |
gezelter |
2624 |
The forces associated with the shifted potential are simply the forces |
336 |
|
|
of the unshifted potential itself (when inside the cutoff sphere), |
337 |
chrisfen |
2601 |
\begin{equation} |
338 |
gezelter |
2643 |
F_{\textrm{SP}} = -\left( \frac{d v(r)}{dr} \right), |
339 |
chrisfen |
2612 |
\end{equation} |
340 |
gezelter |
2624 |
and are zero outside. Inside the cutoff sphere, the forces associated |
341 |
|
|
with the shifted force form can be written, |
342 |
chrisfen |
2612 |
\begin{equation} |
343 |
gezelter |
2643 |
F_{\textrm{SF}} = -\left( \frac{d v(r)}{dr} \right) + \left(\frac{d |
344 |
gezelter |
2624 |
v(r)}{dr} \right)_{r=R_\textrm{c}}. |
345 |
|
|
\end{equation} |
346 |
|
|
|
347 |
gezelter |
2643 |
If the potential, $v(r)$, is taken to be the normal Coulomb potential, |
348 |
gezelter |
2624 |
\begin{equation} |
349 |
|
|
v(r) = \frac{q_i q_j}{r}, |
350 |
|
|
\label{eq:Coulomb} |
351 |
|
|
\end{equation} |
352 |
|
|
then the Shifted Potential ({\sc sp}) forms will give Wolf {\it et |
353 |
|
|
al.}'s undamped prescription: |
354 |
|
|
\begin{equation} |
355 |
gezelter |
2643 |
V_\textrm{SP}(r) = |
356 |
gezelter |
2624 |
q_iq_j\left(\frac{1}{r}-\frac{1}{R_\textrm{c}}\right) \quad |
357 |
|
|
r\leqslant R_\textrm{c}, |
358 |
chrisfen |
2636 |
\label{eq:SPPot} |
359 |
gezelter |
2624 |
\end{equation} |
360 |
|
|
with associated forces, |
361 |
|
|
\begin{equation} |
362 |
gezelter |
2643 |
F_\textrm{SP}(r) = q_iq_j\left(\frac{1}{r^2}\right) \quad r\leqslant R_\textrm{c}. |
363 |
chrisfen |
2636 |
\label{eq:SPForces} |
364 |
chrisfen |
2612 |
\end{equation} |
365 |
gezelter |
2624 |
These forces are identical to the forces of the standard Coulomb |
366 |
|
|
interaction, and cutting these off at $R_c$ was addressed by Wolf |
367 |
|
|
\textit{et al.} as undesirable. They pointed out that the effect of |
368 |
|
|
the image charges is neglected in the forces when this form is |
369 |
|
|
used,\cite{Wolf99} thereby eliminating any benefit from the method in |
370 |
|
|
molecular dynamics. Additionally, there is a discontinuity in the |
371 |
|
|
forces at the cutoff radius which results in energy drift during MD |
372 |
|
|
simulations. |
373 |
chrisfen |
2612 |
|
374 |
gezelter |
2624 |
The shifted force ({\sc sf}) form using the normal Coulomb potential |
375 |
|
|
will give, |
376 |
chrisfen |
2612 |
\begin{equation} |
377 |
gezelter |
2643 |
V_\textrm{SF}(r) = q_iq_j\left[\frac{1}{r}-\frac{1}{R_\textrm{c}}+\left(\frac{1}{R_\textrm{c}^2}\right)(r-R_\textrm{c})\right] \quad r\leqslant R_\textrm{c}. |
378 |
chrisfen |
2612 |
\label{eq:SFPot} |
379 |
|
|
\end{equation} |
380 |
gezelter |
2624 |
with associated forces, |
381 |
chrisfen |
2612 |
\begin{equation} |
382 |
gezelter |
2643 |
F_\textrm{SF}(r) = q_iq_j\left(\frac{1}{r^2}-\frac{1}{R_\textrm{c}^2}\right) \quad r\leqslant R_\textrm{c}. |
383 |
chrisfen |
2612 |
\label{eq:SFForces} |
384 |
|
|
\end{equation} |
385 |
gezelter |
2624 |
This formulation has the benefits that there are no discontinuities at |
386 |
gezelter |
2643 |
the cutoff radius, while the neutralizing image charges are present in |
387 |
|
|
both the energy and force expressions. It would be simple to add the |
388 |
|
|
self-neutralizing term back when computing the total energy of the |
389 |
gezelter |
2624 |
system, thereby maintaining the agreement with the Madelung energies. |
390 |
|
|
A side effect of this treatment is the alteration in the shape of the |
391 |
|
|
potential that comes from the derivative term. Thus, a degree of |
392 |
|
|
clarity about agreement with the empirical potential is lost in order |
393 |
|
|
to gain functionality in dynamics simulations. |
394 |
chrisfen |
2612 |
|
395 |
chrisfen |
2620 |
Wolf \textit{et al.} originally discussed the energetics of the |
396 |
gezelter |
2643 |
shifted Coulomb potential (Eq. \ref{eq:SPPot}) and found that it was |
397 |
|
|
insufficient for accurate determination of the energy with reasonable |
398 |
|
|
cutoff distances. The calculated Madelung energies fluctuated around |
399 |
|
|
the expected value as the cutoff radius was increased, but the |
400 |
|
|
oscillations converged toward the correct value.\cite{Wolf99} A |
401 |
gezelter |
2624 |
damping function was incorporated to accelerate the convergence; and |
402 |
gezelter |
2643 |
though alternative forms for the damping function could be |
403 |
gezelter |
2624 |
used,\cite{Jones56,Heyes81} the complimentary error function was |
404 |
|
|
chosen to mirror the effective screening used in the Ewald summation. |
405 |
|
|
Incorporating this error function damping into the simple Coulomb |
406 |
|
|
potential, |
407 |
chrisfen |
2612 |
\begin{equation} |
408 |
gezelter |
2624 |
v(r) = \frac{\mathrm{erfc}\left(\alpha r\right)}{r}, |
409 |
chrisfen |
2601 |
\label{eq:dampCoulomb} |
410 |
|
|
\end{equation} |
411 |
gezelter |
2643 |
the shifted potential (eq. (\ref{eq:SPPot})) becomes |
412 |
chrisfen |
2601 |
\begin{equation} |
413 |
gezelter |
2643 |
V_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r}-\frac{\textrm{erfc}\left(\alpha R_\textrm{c}\right)}{R_\textrm{c}}\right) \quad r\leqslant R_\textrm{c}, |
414 |
chrisfen |
2612 |
\label{eq:DSPPot} |
415 |
chrisfen |
2629 |
\end{equation} |
416 |
gezelter |
2624 |
with associated forces, |
417 |
chrisfen |
2612 |
\begin{equation} |
418 |
gezelter |
2643 |
F_{\textrm{DSP}}(r) = q_iq_j\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \quad r\leqslant R_\textrm{c}. |
419 |
chrisfen |
2612 |
\label{eq:DSPForces} |
420 |
|
|
\end{equation} |
421 |
gezelter |
2643 |
Again, this damped shifted potential suffers from a |
422 |
|
|
force-discontinuity at the cutoff radius, and the image charges play |
423 |
|
|
no role in the forces. To remedy these concerns, one may derive a |
424 |
|
|
{\sc sf} variant by including the derivative term in |
425 |
|
|
eq. (\ref{eq:shiftingForm}), |
426 |
chrisfen |
2612 |
\begin{equation} |
427 |
chrisfen |
2620 |
\begin{split} |
428 |
gezelter |
2643 |
V_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\frac{\mathrm{erfc}\left(\alpha r\right)}{r}-\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} \\ &\left.+\left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}\right)\left(r-R_\mathrm{c}\right)\ \right] \quad r\leqslant R_\textrm{c}. |
429 |
chrisfen |
2612 |
\label{eq:DSFPot} |
430 |
chrisfen |
2620 |
\end{split} |
431 |
chrisfen |
2612 |
\end{equation} |
432 |
chrisfen |
2636 |
The derivative of the above potential will lead to the following forces, |
433 |
chrisfen |
2612 |
\begin{equation} |
434 |
chrisfen |
2620 |
\begin{split} |
435 |
gezelter |
2643 |
F_\mathrm{DSF}(r) = q_iq_j\Biggr{[}&\left(\frac{\textrm{erfc}\left(\alpha r\right)}{r^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2r^2\right)}}{r}\right) \\ &\left.-\left(\frac{\textrm{erfc}\left(\alpha R_{\textrm{c}}\right)}{R_{\textrm{c}}^2}+\frac{2\alpha}{\pi^{1/2}}\frac{\exp{\left(-\alpha^2R_{\textrm{c}}^2\right)}}{R_{\textrm{c}}}\right)\right] \quad r\leqslant R_\textrm{c}. |
436 |
chrisfen |
2612 |
\label{eq:DSFForces} |
437 |
chrisfen |
2620 |
\end{split} |
438 |
chrisfen |
2612 |
\end{equation} |
439 |
gezelter |
2643 |
If the damping parameter $(\alpha)$ is set to zero, the undamped case, |
440 |
|
|
eqs. (\ref{eq:SPPot} through \ref{eq:SFForces}) are correctly |
441 |
|
|
recovered from eqs. (\ref{eq:DSPPot} through \ref{eq:DSFForces}). |
442 |
chrisfen |
2601 |
|
443 |
chrisfen |
2636 |
This new {\sc sf} potential is similar to equation \ref{eq:ZahnPot} |
444 |
|
|
derived by Zahn \textit{et al.}; however, there are two important |
445 |
|
|
differences.\cite{Zahn02} First, the $v_\textrm{c}$ term from |
446 |
|
|
eq. (\ref{eq:shiftingForm}) is equal to eq. (\ref{eq:dampCoulomb}) |
447 |
|
|
with $r$ replaced by $R_\textrm{c}$. This term is {\it not} present |
448 |
|
|
in the Zahn potential, resulting in a potential discontinuity as |
449 |
|
|
particles cross $R_\textrm{c}$. Second, the sign of the derivative |
450 |
|
|
portion is different. The missing $v_\textrm{c}$ term would not |
451 |
|
|
affect molecular dynamics simulations (although the computed energy |
452 |
|
|
would be expected to have sudden jumps as particle distances crossed |
453 |
gezelter |
2643 |
$R_c$). The sign problem is a potential source of errors, however. |
454 |
|
|
In fact, it introduces a discontinuity in the forces at the cutoff, |
455 |
|
|
because the force function is shifted in the wrong direction and |
456 |
|
|
doesn't cross zero at $R_\textrm{c}$. |
457 |
chrisfen |
2602 |
|
458 |
gezelter |
2624 |
Eqs. (\ref{eq:DSFPot}) and (\ref{eq:DSFForces}) result in an |
459 |
gezelter |
2643 |
electrostatic summation method in which the potential and forces are |
460 |
|
|
continuous at the cutoff radius and which incorporates the damping |
461 |
|
|
function proposed by Wolf \textit{et al.}\cite{Wolf99} In the rest of |
462 |
|
|
this paper, we will evaluate exactly how good these methods ({\sc sp}, |
463 |
|
|
{\sc sf}, damping) are at reproducing the correct electrostatic |
464 |
|
|
summation performed by the Ewald sum. |
465 |
gezelter |
2624 |
|
466 |
|
|
\subsection{Other alternatives} |
467 |
gezelter |
2643 |
In addition to the methods described above, we considered some other |
468 |
|
|
techniques that are commonly used in molecular simulations. The |
469 |
chrisfen |
2629 |
simplest of these is group-based cutoffs. Though of little use for |
470 |
gezelter |
2643 |
charged molecules, collecting atoms into neutral groups takes |
471 |
chrisfen |
2629 |
advantage of the observation that the electrostatic interactions decay |
472 |
|
|
faster than those for monopolar pairs.\cite{Steinbach94} When |
473 |
gezelter |
2643 |
considering these molecules as neutral groups, the relative |
474 |
|
|
orientations of the molecules control the strength of the interactions |
475 |
|
|
at the cutoff radius. Consequently, as these molecular particles move |
476 |
|
|
through $R_\textrm{c}$, the energy will drift upward due to the |
477 |
|
|
anisotropy of the net molecular dipole interactions.\cite{Rahman71} To |
478 |
|
|
maintain good energy conservation, both the potential and derivative |
479 |
|
|
need to be smoothly switched to zero at $R_\textrm{c}$.\cite{Adams79} |
480 |
|
|
This is accomplished using a standard switching function. If a smooth |
481 |
|
|
second derivative is desired, a fifth (or higher) order polynomial can |
482 |
|
|
be used.\cite{Andrea83} |
483 |
gezelter |
2624 |
|
484 |
chrisfen |
2629 |
Group-based cutoffs neglect the surroundings beyond $R_\textrm{c}$, |
485 |
gezelter |
2643 |
and to incorporate the effects of the surroundings, a method like |
486 |
|
|
Reaction Field ({\sc rf}) can be used. The original theory for {\sc |
487 |
|
|
rf} was originally developed by Onsager,\cite{Onsager36} and it was |
488 |
|
|
applied in simulations for the study of water by Barker and |
489 |
|
|
Watts.\cite{Barker73} In modern simulation codes, {\sc rf} is simply |
490 |
|
|
an extension of the group-based cutoff method where the net dipole |
491 |
|
|
within the cutoff sphere polarizes an external dielectric, which |
492 |
|
|
reacts back on the central dipole. The same switching function |
493 |
|
|
considerations for group-based cutoffs need to made for {\sc rf}, with |
494 |
|
|
the additional pre-specification of a dielectric constant. |
495 |
gezelter |
2624 |
|
496 |
chrisfen |
2608 |
\section{Methods} |
497 |
|
|
|
498 |
chrisfen |
2620 |
In classical molecular mechanics simulations, there are two primary |
499 |
|
|
techniques utilized to obtain information about the system of |
500 |
|
|
interest: Monte Carlo (MC) and Molecular Dynamics (MD). Both of these |
501 |
|
|
techniques utilize pairwise summations of interactions between |
502 |
|
|
particle sites, but they use these summations in different ways. |
503 |
chrisfen |
2608 |
|
504 |
gezelter |
2645 |
In MC, the potential energy difference between configurations dictates |
505 |
|
|
the progression of MC sampling. Going back to the origins of this |
506 |
|
|
method, the acceptance criterion for the canonical ensemble laid out |
507 |
|
|
by Metropolis \textit{et al.} states that a subsequent configuration |
508 |
|
|
is accepted if $\Delta E < 0$ or if $\xi < \exp(-\Delta E/kT)$, where |
509 |
|
|
$\xi$ is a random number between 0 and 1.\cite{Metropolis53} |
510 |
|
|
Maintaining the correct $\Delta E$ when using an alternate method for |
511 |
|
|
handling the long-range electrostatics will ensure proper sampling |
512 |
|
|
from the ensemble. |
513 |
chrisfen |
2608 |
|
514 |
gezelter |
2624 |
In MD, the derivative of the potential governs how the system will |
515 |
chrisfen |
2620 |
progress in time. Consequently, the force and torque vectors on each |
516 |
gezelter |
2624 |
body in the system dictate how the system evolves. If the magnitude |
517 |
|
|
and direction of these vectors are similar when using alternate |
518 |
|
|
electrostatic summation techniques, the dynamics in the short term |
519 |
|
|
will be indistinguishable. Because error in MD calculations is |
520 |
|
|
cumulative, one should expect greater deviation at longer times, |
521 |
|
|
although methods which have large differences in the force and torque |
522 |
|
|
vectors will diverge from each other more rapidly. |
523 |
chrisfen |
2608 |
|
524 |
chrisfen |
2609 |
\subsection{Monte Carlo and the Energy Gap}\label{sec:MCMethods} |
525 |
gezelter |
2645 |
|
526 |
gezelter |
2624 |
The pairwise summation techniques (outlined in section |
527 |
|
|
\ref{sec:ESMethods}) were evaluated for use in MC simulations by |
528 |
|
|
studying the energy differences between conformations. We took the |
529 |
gezelter |
2656 |
{\sc spme}-computed energy difference between two conformations to be the |
530 |
gezelter |
2624 |
correct behavior. An ideal performance by an alternative method would |
531 |
gezelter |
2645 |
reproduce these energy differences exactly (even if the absolute |
532 |
|
|
energies calculated by the methods are different). Since none of the |
533 |
|
|
methods provide exact energy differences, we used linear least squares |
534 |
|
|
regressions of energy gap data to evaluate how closely the methods |
535 |
|
|
mimicked the Ewald energy gaps. Unitary results for both the |
536 |
|
|
correlation (slope) and correlation coefficient for these regressions |
537 |
gezelter |
2656 |
indicate perfect agreement between the alternative method and {\sc spme}. |
538 |
gezelter |
2645 |
Sample correlation plots for two alternate methods are shown in |
539 |
|
|
Fig. \ref{fig:linearFit}. |
540 |
chrisfen |
2608 |
|
541 |
chrisfen |
2609 |
\begin{figure} |
542 |
|
|
\centering |
543 |
chrisfen |
2741 |
\includegraphics[width = \linewidth]{./dualLinear.pdf} |
544 |
gezelter |
2645 |
\caption{Example least squares regressions of the configuration energy |
545 |
|
|
differences for SPC/E water systems. The upper plot shows a data set |
546 |
|
|
with a poor correlation coefficient ($R^2$), while the lower plot |
547 |
|
|
shows a data set with a good correlation coefficient.} |
548 |
|
|
\label{fig:linearFit} |
549 |
chrisfen |
2609 |
\end{figure} |
550 |
|
|
|
551 |
chrisfen |
2740 |
Each of the seven system types (detailed in section \ref{sec:RepSims}) |
552 |
|
|
were represented using 500 independent configurations. Thus, each of |
553 |
|
|
the alternative (non-Ewald) electrostatic summation methods was |
554 |
|
|
evaluated using an accumulated 873,250 configurational energy |
555 |
|
|
differences. |
556 |
chrisfen |
2609 |
|
557 |
gezelter |
2624 |
Results and discussion for the individual analysis of each of the |
558 |
chrisfen |
2742 |
system types appear in the supporting information,\cite{EPAPSdeposit} |
559 |
|
|
while the cumulative results over all the investigated systems appears |
560 |
|
|
below in section \ref{sec:EnergyResults}. |
561 |
gezelter |
2624 |
|
562 |
chrisfen |
2609 |
\subsection{Molecular Dynamics and the Force and Torque Vectors}\label{sec:MDMethods} |
563 |
gezelter |
2624 |
We evaluated the pairwise methods (outlined in section |
564 |
|
|
\ref{sec:ESMethods}) for use in MD simulations by |
565 |
|
|
comparing the force and torque vectors with those obtained using the |
566 |
gezelter |
2656 |
reference Ewald summation ({\sc spme}). Both the magnitude and the |
567 |
gezelter |
2624 |
direction of these vectors on each of the bodies in the system were |
568 |
|
|
analyzed. For the magnitude of these vectors, linear least squares |
569 |
|
|
regression analyses were performed as described previously for |
570 |
|
|
comparing $\Delta E$ values. Instead of a single energy difference |
571 |
|
|
between two system configurations, we compared the magnitudes of the |
572 |
|
|
forces (and torques) on each molecule in each configuration. For a |
573 |
|
|
system of 1000 water molecules and 40 ions, there are 1040 force |
574 |
|
|
vectors and 1000 torque vectors. With 500 configurations, this |
575 |
|
|
results in 520,000 force and 500,000 torque vector comparisons. |
576 |
|
|
Additionally, data from seven different system types was aggregated |
577 |
|
|
before the comparison was made. |
578 |
chrisfen |
2609 |
|
579 |
gezelter |
2624 |
The {\it directionality} of the force and torque vectors was |
580 |
|
|
investigated through measurement of the angle ($\theta$) formed |
581 |
gezelter |
2656 |
between those computed from the particular method and those from {\sc spme}, |
582 |
chrisfen |
2610 |
\begin{equation} |
583 |
gezelter |
2645 |
\theta_f = \cos^{-1} \left(\hat{F}_\textrm{SPME} \cdot \hat{F}_\textrm{M}\right), |
584 |
chrisfen |
2610 |
\end{equation} |
585 |
gezelter |
2656 |
where $\hat{F}_\textrm{M}$ is the unit vector pointing along the force |
586 |
chrisfen |
2651 |
vector computed using method M. Each of these $\theta$ values was |
587 |
|
|
accumulated in a distribution function and weighted by the area on the |
588 |
chrisfen |
2652 |
unit sphere. Since this distribution is a measure of angular error |
589 |
|
|
between two different electrostatic summation methods, there is no |
590 |
|
|
{\it a priori} reason for the profile to adhere to any specific |
591 |
|
|
shape. Thus, gaussian fits were used to measure the width of the |
592 |
chrisfen |
2667 |
resulting distributions. The variance ($\sigma^2$) was extracted from |
593 |
|
|
each of these fits and was used to compare distribution widths. |
594 |
|
|
Values of $\sigma^2$ near zero indicate vector directions |
595 |
|
|
indistinguishable from those calculated when using the reference |
596 |
|
|
method ({\sc spme}). |
597 |
gezelter |
2624 |
|
598 |
|
|
\subsection{Short-time Dynamics} |
599 |
gezelter |
2645 |
|
600 |
|
|
The effects of the alternative electrostatic summation methods on the |
601 |
|
|
short-time dynamics of charged systems were evaluated by considering a |
602 |
|
|
NaCl crystal at a temperature of 1000 K. A subset of the best |
603 |
|
|
performing pairwise methods was used in this comparison. The NaCl |
604 |
|
|
crystal was chosen to avoid possible complications from the treatment |
605 |
|
|
of orientational motion in molecular systems. All systems were |
606 |
|
|
started with the same initial positions and velocities. Simulations |
607 |
|
|
were performed under the microcanonical ensemble, and velocity |
608 |
chrisfen |
2638 |
autocorrelation functions (Eq. \ref{eq:vCorr}) were computed for each |
609 |
|
|
of the trajectories, |
610 |
chrisfen |
2609 |
\begin{equation} |
611 |
gezelter |
2656 |
C_v(t) = \frac{\langle v(0)\cdot v(t)\rangle}{\langle v^2\rangle}. |
612 |
chrisfen |
2609 |
\label{eq:vCorr} |
613 |
|
|
\end{equation} |
614 |
chrisfen |
2638 |
Velocity autocorrelation functions require detailed short time data, |
615 |
|
|
thus velocity information was saved every 2 fs over 10 ps |
616 |
|
|
trajectories. Because the NaCl crystal is composed of two different |
617 |
|
|
atom types, the average of the two resulting velocity autocorrelation |
618 |
|
|
functions was used for comparisons. |
619 |
|
|
|
620 |
|
|
\subsection{Long-Time and Collective Motion}\label{sec:LongTimeMethods} |
621 |
gezelter |
2645 |
|
622 |
|
|
The effects of the same subset of alternative electrostatic methods on |
623 |
|
|
the {\it long-time} dynamics of charged systems were evaluated using |
624 |
chrisfen |
2667 |
the same model system (NaCl crystals at 1000~K). The power spectrum |
625 |
gezelter |
2645 |
($I(\omega)$) was obtained via Fourier transform of the velocity |
626 |
|
|
autocorrelation function, \begin{equation} I(\omega) = |
627 |
|
|
\frac{1}{2\pi}\int^{\infty}_{-\infty}C_v(t)e^{-i\omega t}dt, |
628 |
chrisfen |
2609 |
\label{eq:powerSpec} |
629 |
|
|
\end{equation} |
630 |
chrisfen |
2638 |
where the frequency, $\omega=0,\ 1,\ ...,\ N-1$. Again, because the |
631 |
|
|
NaCl crystal is composed of two different atom types, the average of |
632 |
gezelter |
2645 |
the two resulting power spectra was used for comparisons. Simulations |
633 |
|
|
were performed under the microcanonical ensemble, and velocity |
634 |
chrisfen |
2740 |
information was saved every 5~fs over 100~ps trajectories. |
635 |
chrisfen |
2609 |
|
636 |
|
|
\subsection{Representative Simulations}\label{sec:RepSims} |
637 |
chrisfen |
2740 |
A variety of representative molecular simulations were analyzed to |
638 |
|
|
determine the relative effectiveness of the pairwise summation |
639 |
|
|
techniques in reproducing the energetics and dynamics exhibited by |
640 |
|
|
{\sc spme}. We wanted to span the space of typical molecular |
641 |
|
|
simulations (i.e. from liquids of neutral molecules to ionic |
642 |
|
|
crystals), so the systems studied were: |
643 |
chrisfen |
2599 |
\begin{enumerate} |
644 |
gezelter |
2645 |
\item liquid water (SPC/E),\cite{Berendsen87} |
645 |
|
|
\item crystalline water (Ice I$_\textrm{c}$ crystals of SPC/E), |
646 |
|
|
\item NaCl crystals, |
647 |
|
|
\item NaCl melts, |
648 |
|
|
\item a low ionic strength solution of NaCl in water (0.11 M), |
649 |
|
|
\item a high ionic strength solution of NaCl in water (1.1 M), and |
650 |
|
|
\item a 6 \AA\ radius sphere of Argon in water. |
651 |
chrisfen |
2599 |
\end{enumerate} |
652 |
chrisfen |
2620 |
By utilizing the pairwise techniques (outlined in section |
653 |
|
|
\ref{sec:ESMethods}) in systems composed entirely of neutral groups, |
654 |
gezelter |
2645 |
charged particles, and mixtures of the two, we hope to discern under |
655 |
|
|
which conditions it will be possible to use one of the alternative |
656 |
|
|
summation methodologies instead of the Ewald sum. |
657 |
chrisfen |
2586 |
|
658 |
gezelter |
2645 |
For the solid and liquid water configurations, configurations were |
659 |
|
|
taken at regular intervals from high temperature trajectories of 1000 |
660 |
|
|
SPC/E water molecules. Each configuration was equilibrated |
661 |
|
|
independently at a lower temperature (300~K for the liquid, 200~K for |
662 |
|
|
the crystal). The solid and liquid NaCl systems consisted of 500 |
663 |
|
|
$\textrm{Na}^{+}$ and 500 $\textrm{Cl}^{-}$ ions. Configurations for |
664 |
|
|
these systems were selected and equilibrated in the same manner as the |
665 |
chrisfen |
2667 |
water systems. In order to introduce measurable fluctuations in the |
666 |
|
|
configuration energy differences, the crystalline simulations were |
667 |
|
|
equilibrated at 1000~K, near the $T_\textrm{m}$ for NaCl. The liquid |
668 |
|
|
NaCl configurations needed to represent a fully disordered array of |
669 |
|
|
point charges, so the high temperature of 7000~K was selected for |
670 |
|
|
equilibration. The ionic solutions were made by solvating 4 (or 40) |
671 |
|
|
ions in a periodic box containing 1000 SPC/E water molecules. Ion and |
672 |
|
|
water positions were then randomly swapped, and the resulting |
673 |
|
|
configurations were again equilibrated individually. Finally, for the |
674 |
|
|
Argon / Water ``charge void'' systems, the identities of all the SPC/E |
675 |
|
|
waters within 6 \AA\ of the center of the equilibrated water |
676 |
|
|
configurations were converted to argon. |
677 |
chrisfen |
2586 |
|
678 |
gezelter |
2645 |
These procedures guaranteed us a set of representative configurations |
679 |
gezelter |
2653 |
from chemically-relevant systems sampled from appropriate |
680 |
|
|
ensembles. Force field parameters for the ions and Argon were taken |
681 |
gezelter |
2645 |
from the force field utilized by {\sc oopse}.\cite{Meineke05} |
682 |
|
|
|
683 |
|
|
\subsection{Comparison of Summation Methods}\label{sec:ESMethods} |
684 |
|
|
We compared the following alternative summation methods with results |
685 |
gezelter |
2656 |
from the reference method ({\sc spme}): |
686 |
gezelter |
2645 |
\begin{itemize} |
687 |
|
|
\item {\sc sp} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
688 |
|
|
and 0.3 \AA$^{-1}$, |
689 |
|
|
\item {\sc sf} with damping parameters ($\alpha$) of 0.0, 0.1, 0.2, |
690 |
|
|
and 0.3 \AA$^{-1}$, |
691 |
|
|
\item reaction field with an infinite dielectric constant, and |
692 |
|
|
\item an unmodified cutoff. |
693 |
|
|
\end{itemize} |
694 |
|
|
Group-based cutoffs with a fifth-order polynomial switching function |
695 |
|
|
were utilized for the reaction field simulations. Additionally, we |
696 |
gezelter |
2656 |
investigated the use of these cutoffs with the {\sc sp}, {\sc sf}, and pure |
697 |
|
|
cutoff. The {\sc spme} electrostatics were performed using the {\sc tinker} |
698 |
|
|
implementation of {\sc spme},\cite{Ponder87} while all other calculations |
699 |
gezelter |
2653 |
were performed using the {\sc oopse} molecular mechanics |
700 |
gezelter |
2645 |
package.\cite{Meineke05} All other portions of the energy calculation |
701 |
|
|
(i.e. Lennard-Jones interactions) were handled in exactly the same |
702 |
|
|
manner across all systems and configurations. |
703 |
chrisfen |
2586 |
|
704 |
chrisfen |
2667 |
The alternative methods were also evaluated with three different |
705 |
chrisfen |
2649 |
cutoff radii (9, 12, and 15 \AA). As noted previously, the |
706 |
|
|
convergence parameter ($\alpha$) plays a role in the balance of the |
707 |
|
|
real-space and reciprocal-space portions of the Ewald calculation. |
708 |
|
|
Typical molecular mechanics packages set this to a value dependent on |
709 |
|
|
the cutoff radius and a tolerance (typically less than $1 \times |
710 |
|
|
10^{-4}$ kcal/mol). Smaller tolerances are typically associated with |
711 |
gezelter |
2653 |
increasing accuracy at the expense of computational time spent on the |
712 |
|
|
reciprocal-space portion of the summation.\cite{Perram88,Essmann95} |
713 |
gezelter |
2656 |
The default {\sc tinker} tolerance of $1 \times 10^{-8}$ kcal/mol was used |
714 |
|
|
in all {\sc spme} calculations, resulting in Ewald coefficients of 0.4200, |
715 |
gezelter |
2653 |
0.3119, and 0.2476 \AA$^{-1}$ for cutoff radii of 9, 12, and 15 \AA\ |
716 |
|
|
respectively. |
717 |
chrisfen |
2609 |
|
718 |
chrisfen |
2575 |
\section{Results and Discussion} |
719 |
|
|
|
720 |
chrisfen |
2609 |
\subsection{Configuration Energy Differences}\label{sec:EnergyResults} |
721 |
chrisfen |
2620 |
In order to evaluate the performance of the pairwise electrostatic |
722 |
|
|
summation methods for Monte Carlo simulations, the energy differences |
723 |
|
|
between configurations were compared to the values obtained when using |
724 |
gezelter |
2656 |
{\sc spme}. The results for the subsequent regression analysis are shown in |
725 |
chrisfen |
2620 |
figure \ref{fig:delE}. |
726 |
chrisfen |
2590 |
|
727 |
|
|
\begin{figure} |
728 |
|
|
\centering |
729 |
chrisfen |
2741 |
\includegraphics[width=5.5in]{./delEplot.pdf} |
730 |
gezelter |
2645 |
\caption{Statistical analysis of the quality of configurational energy |
731 |
|
|
differences for a given electrostatic method compared with the |
732 |
|
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
733 |
|
|
indicate $\Delta E$ values indistinguishable from those obtained using |
734 |
gezelter |
2656 |
{\sc spme}. Different values of the cutoff radius are indicated with |
735 |
gezelter |
2645 |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
736 |
|
|
inverted triangles).} |
737 |
chrisfen |
2601 |
\label{fig:delE} |
738 |
chrisfen |
2594 |
\end{figure} |
739 |
|
|
|
740 |
gezelter |
2645 |
The most striking feature of this plot is how well the Shifted Force |
741 |
|
|
({\sc sf}) and Shifted Potential ({\sc sp}) methods capture the energy |
742 |
|
|
differences. For the undamped {\sc sf} method, and the |
743 |
|
|
moderately-damped {\sc sp} methods, the results are nearly |
744 |
|
|
indistinguishable from the Ewald results. The other common methods do |
745 |
|
|
significantly less well. |
746 |
chrisfen |
2594 |
|
747 |
gezelter |
2645 |
The unmodified cutoff method is essentially unusable. This is not |
748 |
|
|
surprising since hard cutoffs give large energy fluctuations as atoms |
749 |
|
|
or molecules move in and out of the cutoff |
750 |
|
|
radius.\cite{Rahman71,Adams79} These fluctuations can be alleviated to |
751 |
|
|
some degree by using group based cutoffs with a switching |
752 |
|
|
function.\cite{Adams79,Steinbach94,Leach01} However, we do not see |
753 |
|
|
significant improvement using the group-switched cutoff because the |
754 |
|
|
salt and salt solution systems contain non-neutral groups. Interested |
755 |
|
|
readers can consult the accompanying supporting information for a |
756 |
chrisfen |
2742 |
comparison where all groups are neutral.\cite{EPAPSdeposit} |
757 |
gezelter |
2645 |
|
758 |
gezelter |
2653 |
For the {\sc sp} method, inclusion of electrostatic damping improves |
759 |
|
|
the agreement with Ewald, and using an $\alpha$ of 0.2 \AA $^{-1}$ |
760 |
gezelter |
2656 |
shows an excellent correlation and quality of fit with the {\sc spme} |
761 |
gezelter |
2653 |
results, particularly with a cutoff radius greater than 12 |
762 |
gezelter |
2645 |
\AA . Use of a larger damping parameter is more helpful for the |
763 |
|
|
shortest cutoff shown, but it has a detrimental effect on simulations |
764 |
|
|
with larger cutoffs. |
765 |
chrisfen |
2609 |
|
766 |
gezelter |
2653 |
In the {\sc sf} sets, increasing damping results in progressively {\it |
767 |
|
|
worse} correlation with Ewald. Overall, the undamped case is the best |
768 |
gezelter |
2645 |
performing set, as the correlation and quality of fits are |
769 |
|
|
consistently superior regardless of the cutoff distance. The undamped |
770 |
|
|
case is also less computationally demanding (because no evaluation of |
771 |
|
|
the complementary error function is required). |
772 |
|
|
|
773 |
|
|
The reaction field results illustrates some of that method's |
774 |
|
|
limitations, primarily that it was developed for use in homogenous |
775 |
|
|
systems; although it does provide results that are an improvement over |
776 |
|
|
those from an unmodified cutoff. |
777 |
|
|
|
778 |
chrisfen |
2608 |
\subsection{Magnitudes of the Force and Torque Vectors} |
779 |
chrisfen |
2599 |
|
780 |
chrisfen |
2620 |
Evaluation of pairwise methods for use in Molecular Dynamics |
781 |
|
|
simulations requires consideration of effects on the forces and |
782 |
gezelter |
2653 |
torques. Figures \ref{fig:frcMag} and \ref{fig:trqMag} show the |
783 |
|
|
regression results for the force and torque vector magnitudes, |
784 |
|
|
respectively. The data in these figures was generated from an |
785 |
|
|
accumulation of the statistics from all of the system types. |
786 |
chrisfen |
2594 |
|
787 |
|
|
\begin{figure} |
788 |
|
|
\centering |
789 |
chrisfen |
2741 |
\includegraphics[width=5.5in]{./frcMagplot.pdf} |
790 |
chrisfen |
2651 |
\caption{Statistical analysis of the quality of the force vector |
791 |
|
|
magnitudes for a given electrostatic method compared with the |
792 |
|
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
793 |
|
|
indicate force magnitude values indistinguishable from those obtained |
794 |
gezelter |
2656 |
using {\sc spme}. Different values of the cutoff radius are indicated with |
795 |
chrisfen |
2651 |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
796 |
|
|
inverted triangles).} |
797 |
chrisfen |
2601 |
\label{fig:frcMag} |
798 |
chrisfen |
2594 |
\end{figure} |
799 |
|
|
|
800 |
gezelter |
2653 |
Again, it is striking how well the Shifted Potential and Shifted Force |
801 |
gezelter |
2656 |
methods are doing at reproducing the {\sc spme} forces. The undamped and |
802 |
gezelter |
2653 |
weakly-damped {\sc sf} method gives the best agreement with Ewald. |
803 |
|
|
This is perhaps expected because this method explicitly incorporates a |
804 |
|
|
smooth transition in the forces at the cutoff radius as well as the |
805 |
|
|
neutralizing image charges. |
806 |
|
|
|
807 |
chrisfen |
2620 |
Figure \ref{fig:frcMag}, for the most part, parallels the results seen |
808 |
|
|
in the previous $\Delta E$ section. The unmodified cutoff results are |
809 |
|
|
poor, but using group based cutoffs and a switching function provides |
810 |
gezelter |
2653 |
an improvement much more significant than what was seen with $\Delta |
811 |
|
|
E$. |
812 |
|
|
|
813 |
|
|
With moderate damping and a large enough cutoff radius, the {\sc sp} |
814 |
|
|
method is generating usable forces. Further increases in damping, |
815 |
chrisfen |
2620 |
while beneficial for simulations with a cutoff radius of 9 \AA\ , is |
816 |
gezelter |
2653 |
detrimental to simulations with larger cutoff radii. |
817 |
|
|
|
818 |
|
|
The reaction field results are surprisingly good, considering the poor |
819 |
chrisfen |
2620 |
quality of the fits for the $\Delta E$ results. There is still a |
820 |
gezelter |
2653 |
considerable degree of scatter in the data, but the forces correlate |
821 |
|
|
well with the Ewald forces in general. We note that the reaction |
822 |
|
|
field calculations do not include the pure NaCl systems, so these |
823 |
chrisfen |
2620 |
results are partly biased towards conditions in which the method |
824 |
|
|
performs more favorably. |
825 |
chrisfen |
2594 |
|
826 |
|
|
\begin{figure} |
827 |
|
|
\centering |
828 |
chrisfen |
2741 |
\includegraphics[width=5.5in]{./trqMagplot.pdf} |
829 |
chrisfen |
2651 |
\caption{Statistical analysis of the quality of the torque vector |
830 |
|
|
magnitudes for a given electrostatic method compared with the |
831 |
|
|
reference Ewald sum. Results with a value equal to 1 (dashed line) |
832 |
|
|
indicate torque magnitude values indistinguishable from those obtained |
833 |
gezelter |
2656 |
using {\sc spme}. Different values of the cutoff radius are indicated with |
834 |
chrisfen |
2651 |
different symbols (9\AA\ = circles, 12\AA\ = squares, and 15\AA\ = |
835 |
|
|
inverted triangles).} |
836 |
chrisfen |
2601 |
\label{fig:trqMag} |
837 |
chrisfen |
2594 |
\end{figure} |
838 |
|
|
|
839 |
gezelter |
2653 |
Molecular torques were only available from the systems which contained |
840 |
|
|
rigid molecules (i.e. the systems containing water). The data in |
841 |
|
|
fig. \ref{fig:trqMag} is taken from this smaller sampling pool. |
842 |
chrisfen |
2594 |
|
843 |
gezelter |
2653 |
Torques appear to be much more sensitive to charges at a longer |
844 |
|
|
distance. The striking feature in comparing the new electrostatic |
845 |
gezelter |
2656 |
methods with {\sc spme} is how much the agreement improves with increasing |
846 |
gezelter |
2653 |
cutoff radius. Again, the weakly damped and undamped {\sc sf} method |
847 |
gezelter |
2656 |
appears to be reproducing the {\sc spme} torques most accurately. |
848 |
gezelter |
2653 |
|
849 |
|
|
Water molecules are dipolar, and the reaction field method reproduces |
850 |
|
|
the effect of the surrounding polarized medium on each of the |
851 |
|
|
molecular bodies. Therefore it is not surprising that reaction field |
852 |
|
|
performs best of all of the methods on molecular torques. |
853 |
|
|
|
854 |
chrisfen |
2608 |
\subsection{Directionality of the Force and Torque Vectors} |
855 |
chrisfen |
2599 |
|
856 |
gezelter |
2653 |
It is clearly important that a new electrostatic method can reproduce |
857 |
|
|
the magnitudes of the force and torque vectors obtained via the Ewald |
858 |
|
|
sum. However, the {\it directionality} of these vectors will also be |
859 |
|
|
vital in calculating dynamical quantities accurately. Force and |
860 |
|
|
torque directionalities were investigated by measuring the angles |
861 |
|
|
formed between these vectors and the same vectors calculated using |
862 |
gezelter |
2656 |
{\sc spme}. The results (Fig. \ref{fig:frcTrqAng}) are compared through the |
863 |
gezelter |
2653 |
variance ($\sigma^2$) of the Gaussian fits of the angle error |
864 |
|
|
distributions of the combined set over all system types. |
865 |
chrisfen |
2594 |
|
866 |
|
|
\begin{figure} |
867 |
|
|
\centering |
868 |
chrisfen |
2741 |
\includegraphics[width=5.5in]{./frcTrqAngplot.pdf} |
869 |
gezelter |
2653 |
\caption{Statistical analysis of the width of the angular distribution |
870 |
|
|
that the force and torque vectors from a given electrostatic method |
871 |
|
|
make with their counterparts obtained using the reference Ewald sum. |
872 |
|
|
Results with a variance ($\sigma^2$) equal to zero (dashed line) |
873 |
|
|
indicate force and torque directions indistinguishable from those |
874 |
gezelter |
2656 |
obtained using {\sc spme}. Different values of the cutoff radius are |
875 |
gezelter |
2653 |
indicated with different symbols (9\AA\ = circles, 12\AA\ = squares, |
876 |
|
|
and 15\AA\ = inverted triangles).} |
877 |
chrisfen |
2601 |
\label{fig:frcTrqAng} |
878 |
chrisfen |
2594 |
\end{figure} |
879 |
|
|
|
880 |
chrisfen |
2620 |
Both the force and torque $\sigma^2$ results from the analysis of the |
881 |
|
|
total accumulated system data are tabulated in figure |
882 |
gezelter |
2653 |
\ref{fig:frcTrqAng}. Here it is clear that the Shifted Potential ({\sc |
883 |
gezelter |
2656 |
sp}) method would be essentially unusable for molecular dynamics |
884 |
|
|
unless the damping function is added. The Shifted Force ({\sc sf}) |
885 |
|
|
method, however, is generating force and torque vectors which are |
886 |
|
|
within a few degrees of the Ewald results even with weak (or no) |
887 |
|
|
damping. |
888 |
chrisfen |
2594 |
|
889 |
gezelter |
2653 |
All of the sets (aside from the over-damped case) show the improvement |
890 |
|
|
afforded by choosing a larger cutoff radius. Increasing the cutoff |
891 |
|
|
from 9 to 12 \AA\ typically results in a halving of the width of the |
892 |
gezelter |
2656 |
distribution, with a similar improvement when going from 12 to 15 |
893 |
gezelter |
2653 |
\AA . |
894 |
|
|
|
895 |
|
|
The undamped {\sc sf}, group-based cutoff, and reaction field methods |
896 |
|
|
all do equivalently well at capturing the direction of both the force |
897 |
gezelter |
2656 |
and torque vectors. Using the electrostatic damping improves the |
898 |
|
|
angular behavior significantly for the {\sc sp} and moderately for the |
899 |
|
|
{\sc sf} methods. Overdamping is detrimental to both methods. Again |
900 |
|
|
it is important to recognize that the force vectors cover all |
901 |
|
|
particles in all seven systems, while torque vectors are only |
902 |
|
|
available for neutral molecular groups. Damping is more beneficial to |
903 |
gezelter |
2653 |
charged bodies, and this observation is investigated further in the |
904 |
chrisfen |
2742 |
accompanying supporting information.\cite{EPAPSdeposit} |
905 |
gezelter |
2653 |
|
906 |
|
|
Although not discussed previously, group based cutoffs can be applied |
907 |
gezelter |
2656 |
to both the {\sc sp} and {\sc sf} methods. The group-based cutoffs |
908 |
|
|
will reintroduce small discontinuities at the cutoff radius, but the |
909 |
|
|
effects of these can be minimized by utilizing a switching function. |
910 |
|
|
Though there are no significant benefits or drawbacks observed in |
911 |
|
|
$\Delta E$ and the force and torque magnitudes when doing this, there |
912 |
|
|
is a measurable improvement in the directionality of the forces and |
913 |
|
|
torques. Table \ref{tab:groupAngle} shows the angular variances |
914 |
|
|
obtained using group based cutoffs along with the results seen in |
915 |
|
|
figure \ref{fig:frcTrqAng}. The {\sc sp} (with an $\alpha$ of 0.2 |
916 |
|
|
\AA$^{-1}$ or smaller) shows much narrower angular distributions when |
917 |
|
|
using group-based cutoffs. The {\sc sf} method likewise shows |
918 |
|
|
improvement in the undamped and lightly damped cases. |
919 |
gezelter |
2653 |
|
920 |
chrisfen |
2595 |
\begin{table}[htbp] |
921 |
gezelter |
2656 |
\centering |
922 |
|
|
\caption{Statistical analysis of the angular |
923 |
|
|
distributions that the force (upper) and torque (lower) vectors |
924 |
|
|
from a given electrostatic method make with their counterparts |
925 |
|
|
obtained using the reference Ewald sum. Calculations were |
926 |
|
|
performed both with (Y) and without (N) group based cutoffs and a |
927 |
|
|
switching function. The $\alpha$ values have units of \AA$^{-1}$ |
928 |
|
|
and the variance values have units of degrees$^2$.} |
929 |
|
|
|
930 |
chrisfen |
2599 |
\begin{tabular}{@{} ccrrrrrrrr @{}} |
931 |
chrisfen |
2595 |
\\ |
932 |
|
|
\toprule |
933 |
|
|
& & \multicolumn{4}{c}{Shifted Potential} & \multicolumn{4}{c}{Shifted Force} \\ |
934 |
|
|
\cmidrule(lr){3-6} |
935 |
|
|
\cmidrule(l){7-10} |
936 |
chrisfen |
2599 |
$R_\textrm{c}$ & Groups & $\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$ & $\alpha = 0$ & $\alpha = 0.1$ & $\alpha = 0.2$ & $\alpha = 0.3$\\ |
937 |
chrisfen |
2595 |
\midrule |
938 |
chrisfen |
2599 |
|
939 |
|
|
9 \AA & N & 29.545 & 12.003 & 5.489 & 0.610 & 2.323 & 2.321 & 0.429 & 0.603 \\ |
940 |
|
|
& \textbf{Y} & \textbf{2.486} & \textbf{2.160} & \textbf{0.667} & \textbf{0.608} & \textbf{1.768} & \textbf{1.766} & \textbf{0.676} & \textbf{0.609} \\ |
941 |
|
|
12 \AA & N & 19.381 & 3.097 & 0.190 & 0.608 & 0.920 & 0.736 & 0.133 & 0.612 \\ |
942 |
|
|
& \textbf{Y} & \textbf{0.515} & \textbf{0.288} & \textbf{0.127} & \textbf{0.586} & \textbf{0.308} & \textbf{0.249} & \textbf{0.127} & \textbf{0.586} \\ |
943 |
|
|
15 \AA & N & 12.700 & 1.196 & 0.123 & 0.601 & 0.339 & 0.160 & 0.123 & 0.601 \\ |
944 |
|
|
& \textbf{Y} & \textbf{0.228} & \textbf{0.099} & \textbf{0.121} & \textbf{0.598} & \textbf{0.144} & \textbf{0.090} & \textbf{0.121} & \textbf{0.598} \\ |
945 |
chrisfen |
2594 |
|
946 |
chrisfen |
2595 |
\midrule |
947 |
|
|
|
948 |
chrisfen |
2599 |
9 \AA & N & 262.716 & 116.585 & 5.234 & 5.103 & 2.392 & 2.350 & 1.770 & 5.122 \\ |
949 |
|
|
& \textbf{Y} & \textbf{2.115} & \textbf{1.914} & \textbf{1.878} & \textbf{5.142} & \textbf{2.076} & \textbf{2.039} & \textbf{1.972} & \textbf{5.146} \\ |
950 |
|
|
12 \AA & N & 129.576 & 25.560 & 1.369 & 5.080 & 0.913 & 0.790 & 1.362 & 5.124 \\ |
951 |
|
|
& \textbf{Y} & \textbf{0.810} & \textbf{0.685} & \textbf{1.352} & \textbf{5.082} & \textbf{0.765} & \textbf{0.714} & \textbf{1.360} & \textbf{5.082} \\ |
952 |
|
|
15 \AA & N & 87.275 & 4.473 & 1.271 & 5.000 & 0.372 & 0.312 & 1.271 & 5.000 \\ |
953 |
|
|
& \textbf{Y} & \textbf{0.282} & \textbf{0.294} & \textbf{1.272} & \textbf{4.999} & \textbf{0.324} & \textbf{0.318} & \textbf{1.272} & \textbf{4.999} \\ |
954 |
chrisfen |
2595 |
|
955 |
|
|
\bottomrule |
956 |
|
|
\end{tabular} |
957 |
chrisfen |
2601 |
\label{tab:groupAngle} |
958 |
chrisfen |
2595 |
\end{table} |
959 |
|
|
|
960 |
gezelter |
2656 |
One additional trend in table \ref{tab:groupAngle} is that the |
961 |
|
|
$\sigma^2$ values for both {\sc sp} and {\sc sf} converge as $\alpha$ |
962 |
|
|
increases, something that is more obvious with group-based cutoffs. |
963 |
|
|
The complimentary error function inserted into the potential weakens |
964 |
|
|
the electrostatic interaction as the value of $\alpha$ is increased. |
965 |
|
|
However, at larger values of $\alpha$, it is possible to overdamp the |
966 |
|
|
electrostatic interaction and to remove it completely. Kast |
967 |
gezelter |
2653 |
\textit{et al.} developed a method for choosing appropriate $\alpha$ |
968 |
|
|
values for these types of electrostatic summation methods by fitting |
969 |
|
|
to $g(r)$ data, and their methods indicate optimal values of 0.34, |
970 |
|
|
0.25, and 0.16 \AA$^{-1}$ for cutoff values of 9, 12, and 15 \AA\ |
971 |
|
|
respectively.\cite{Kast03} These appear to be reasonable choices to |
972 |
|
|
obtain proper MC behavior (Fig. \ref{fig:delE}); however, based on |
973 |
|
|
these findings, choices this high would introduce error in the |
974 |
gezelter |
2656 |
molecular torques, particularly for the shorter cutoffs. Based on our |
975 |
|
|
observations, empirical damping up to 0.2 \AA$^{-1}$ is beneficial, |
976 |
|
|
but damping may be unnecessary when using the {\sc sf} method. |
977 |
chrisfen |
2595 |
|
978 |
chrisfen |
2638 |
\subsection{Short-Time Dynamics: Velocity Autocorrelation Functions of NaCl Crystals} |
979 |
chrisfen |
2601 |
|
980 |
gezelter |
2653 |
Zahn {\it et al.} investigated the structure and dynamics of water |
981 |
|
|
using eqs. (\ref{eq:ZahnPot}) and |
982 |
|
|
(\ref{eq:WolfForces}).\cite{Zahn02,Kast03} Their results indicated |
983 |
|
|
that a method similar (but not identical with) the damped {\sc sf} |
984 |
|
|
method resulted in properties very similar to those obtained when |
985 |
|
|
using the Ewald summation. The properties they studied (pair |
986 |
|
|
distribution functions, diffusion constants, and velocity and |
987 |
|
|
orientational correlation functions) may not be particularly sensitive |
988 |
|
|
to the long-range and collective behavior that governs the |
989 |
gezelter |
2656 |
low-frequency behavior in crystalline systems. Additionally, the |
990 |
|
|
ionic crystals are the worst case scenario for the pairwise methods |
991 |
|
|
because they lack the reciprocal space contribution contained in the |
992 |
|
|
Ewald summation. |
993 |
chrisfen |
2601 |
|
994 |
gezelter |
2653 |
We are using two separate measures to probe the effects of these |
995 |
|
|
alternative electrostatic methods on the dynamics in crystalline |
996 |
|
|
materials. For short- and intermediate-time dynamics, we are |
997 |
|
|
computing the velocity autocorrelation function, and for long-time |
998 |
|
|
and large length-scale collective motions, we are looking at the |
999 |
|
|
low-frequency portion of the power spectrum. |
1000 |
|
|
|
1001 |
chrisfen |
2601 |
\begin{figure} |
1002 |
|
|
\centering |
1003 |
chrisfen |
2741 |
\includegraphics[width = \linewidth]{./vCorrPlot.pdf} |
1004 |
gezelter |
2656 |
\caption{Velocity autocorrelation functions of NaCl crystals at |
1005 |
|
|
1000 K using {\sc spme}, {\sc sf} ($\alpha$ = 0.0, 0.1, \& 0.2), and {\sc |
1006 |
gezelter |
2653 |
sp} ($\alpha$ = 0.2). The inset is a magnification of the area around |
1007 |
|
|
the first minimum. The times to first collision are nearly identical, |
1008 |
|
|
but differences can be seen in the peaks and troughs, where the |
1009 |
|
|
undamped and weakly damped methods are stiffer than the moderately |
1010 |
gezelter |
2656 |
damped and {\sc spme} methods.} |
1011 |
chrisfen |
2638 |
\label{fig:vCorrPlot} |
1012 |
|
|
\end{figure} |
1013 |
|
|
|
1014 |
gezelter |
2656 |
The short-time decay of the velocity autocorrelation function through |
1015 |
|
|
the first collision are nearly identical in figure |
1016 |
|
|
\ref{fig:vCorrPlot}, but the peaks and troughs of the functions show |
1017 |
|
|
how the methods differ. The undamped {\sc sf} method has deeper |
1018 |
|
|
troughs (see inset in Fig. \ref{fig:vCorrPlot}) and higher peaks than |
1019 |
|
|
any of the other methods. As the damping parameter ($\alpha$) is |
1020 |
|
|
increased, these peaks are smoothed out, and the {\sc sf} method |
1021 |
|
|
approaches the {\sc spme} results. With $\alpha$ values of 0.2 \AA$^{-1}$, |
1022 |
|
|
the {\sc sf} and {\sc sp} functions are nearly identical and track the |
1023 |
|
|
{\sc spme} features quite well. This is not surprising because the {\sc sf} |
1024 |
|
|
and {\sc sp} potentials become nearly identical with increased |
1025 |
|
|
damping. However, this appears to indicate that once damping is |
1026 |
|
|
utilized, the details of the form of the potential (and forces) |
1027 |
|
|
constructed out of the damped electrostatic interaction are less |
1028 |
|
|
important. |
1029 |
chrisfen |
2638 |
|
1030 |
|
|
\subsection{Collective Motion: Power Spectra of NaCl Crystals} |
1031 |
|
|
|
1032 |
gezelter |
2656 |
To evaluate how the differences between the methods affect the |
1033 |
|
|
collective long-time motion, we computed power spectra from long-time |
1034 |
|
|
traces of the velocity autocorrelation function. The power spectra for |
1035 |
|
|
the best-performing alternative methods are shown in |
1036 |
|
|
fig. \ref{fig:methodPS}. Apodization of the correlation functions via |
1037 |
|
|
a cubic switching function between 40 and 50 ps was used to reduce the |
1038 |
|
|
ringing resulting from data truncation. This procedure had no |
1039 |
|
|
noticeable effect on peak location or magnitude. |
1040 |
chrisfen |
2638 |
|
1041 |
|
|
\begin{figure} |
1042 |
|
|
\centering |
1043 |
chrisfen |
2741 |
\includegraphics[width = \linewidth]{./spectraSquare.pdf} |
1044 |
chrisfen |
2651 |
\caption{Power spectra obtained from the velocity auto-correlation |
1045 |
gezelter |
2656 |
functions of NaCl crystals at 1000 K while using {\sc spme}, {\sc sf} |
1046 |
|
|
($\alpha$ = 0, 0.1, \& 0.2), and {\sc sp} ($\alpha$ = 0.2). The inset |
1047 |
|
|
shows the frequency region below 100 cm$^{-1}$ to highlight where the |
1048 |
|
|
spectra differ.} |
1049 |
chrisfen |
2610 |
\label{fig:methodPS} |
1050 |
chrisfen |
2601 |
\end{figure} |
1051 |
|
|
|
1052 |
gezelter |
2656 |
While the high frequency regions of the power spectra for the |
1053 |
|
|
alternative methods are quantitatively identical with Ewald spectrum, |
1054 |
|
|
the low frequency region shows how the summation methods differ. |
1055 |
|
|
Considering the low-frequency inset (expanded in the upper frame of |
1056 |
|
|
figure \ref{fig:dampInc}), at frequencies below 100 cm$^{-1}$, the |
1057 |
|
|
correlated motions are blue-shifted when using undamped or weakly |
1058 |
|
|
damped {\sc sf}. When using moderate damping ($\alpha = 0.2$ |
1059 |
|
|
\AA$^{-1}$) both the {\sc sf} and {\sc sp} methods give nearly identical |
1060 |
|
|
correlated motion to the Ewald method (which has a convergence |
1061 |
|
|
parameter of 0.3119 \AA$^{-1}$). This weakening of the electrostatic |
1062 |
|
|
interaction with increased damping explains why the long-ranged |
1063 |
|
|
correlated motions are at lower frequencies for the moderately damped |
1064 |
|
|
methods than for undamped or weakly damped methods. |
1065 |
|
|
|
1066 |
|
|
To isolate the role of the damping constant, we have computed the |
1067 |
|
|
spectra for a single method ({\sc sf}) with a range of damping |
1068 |
|
|
constants and compared this with the {\sc spme} spectrum. |
1069 |
|
|
Fig. \ref{fig:dampInc} shows more clearly that increasing the |
1070 |
|
|
electrostatic damping red-shifts the lowest frequency phonon modes. |
1071 |
|
|
However, even without any electrostatic damping, the {\sc sf} method |
1072 |
|
|
has at most a 10 cm$^{-1}$ error in the lowest frequency phonon mode. |
1073 |
|
|
Without the {\sc sf} modifications, an undamped (pure cutoff) method |
1074 |
|
|
would predict the lowest frequency peak near 325 cm$^{-1}$. {\it |
1075 |
|
|
Most} of the collective behavior in the crystal is accurately captured |
1076 |
|
|
using the {\sc sf} method. Quantitative agreement with Ewald can be |
1077 |
|
|
obtained using moderate damping in addition to the shifting at the |
1078 |
|
|
cutoff distance. |
1079 |
|
|
|
1080 |
chrisfen |
2601 |
\begin{figure} |
1081 |
|
|
\centering |
1082 |
chrisfen |
2741 |
\includegraphics[width = \linewidth]{./increasedDamping.pdf} |
1083 |
gezelter |
2656 |
\caption{Effect of damping on the two lowest-frequency phonon modes in |
1084 |
chrisfen |
2667 |
the NaCl crystal at 1000~K. The undamped shifted force ({\sc sf}) |
1085 |
gezelter |
2656 |
method is off by less than 10 cm$^{-1}$, and increasing the |
1086 |
|
|
electrostatic damping to 0.25 \AA$^{-1}$ gives quantitative agreement |
1087 |
|
|
with the power spectrum obtained using the Ewald sum. Overdamping can |
1088 |
|
|
result in underestimates of frequencies of the long-wavelength |
1089 |
|
|
motions.} |
1090 |
chrisfen |
2601 |
\label{fig:dampInc} |
1091 |
|
|
\end{figure} |
1092 |
|
|
|
1093 |
chrisfen |
2575 |
\section{Conclusions} |
1094 |
|
|
|
1095 |
chrisfen |
2620 |
This investigation of pairwise electrostatic summation techniques |
1096 |
gezelter |
2656 |
shows that there are viable and computationally efficient alternatives |
1097 |
|
|
to the Ewald summation. These methods are derived from the damped and |
1098 |
|
|
cutoff-neutralized Coulombic sum originally proposed by Wolf |
1099 |
|
|
\textit{et al.}\cite{Wolf99} In particular, the {\sc sf} |
1100 |
|
|
method, reformulated above as eqs. (\ref{eq:DSFPot}) and |
1101 |
|
|
(\ref{eq:DSFForces}), shows a remarkable ability to reproduce the |
1102 |
|
|
energetic and dynamic characteristics exhibited by simulations |
1103 |
|
|
employing lattice summation techniques. The cumulative energy |
1104 |
|
|
difference results showed the undamped {\sc sf} and moderately damped |
1105 |
|
|
{\sc sp} methods produced results nearly identical to {\sc spme}. Similarly |
1106 |
|
|
for the dynamic features, the undamped or moderately damped {\sc sf} |
1107 |
|
|
and moderately damped {\sc sp} methods produce force and torque vector |
1108 |
|
|
magnitude and directions very similar to the expected values. These |
1109 |
|
|
results translate into long-time dynamic behavior equivalent to that |
1110 |
|
|
produced in simulations using {\sc spme}. |
1111 |
chrisfen |
2604 |
|
1112 |
gezelter |
2656 |
As in all purely-pairwise cutoff methods, these methods are expected |
1113 |
|
|
to scale approximately {\it linearly} with system size, and they are |
1114 |
|
|
easily parallelizable. This should result in substantial reductions |
1115 |
|
|
in the computational cost of performing large simulations. |
1116 |
|
|
|
1117 |
chrisfen |
2620 |
Aside from the computational cost benefit, these techniques have |
1118 |
|
|
applicability in situations where the use of the Ewald sum can prove |
1119 |
gezelter |
2656 |
problematic. Of greatest interest is their potential use in |
1120 |
|
|
interfacial systems, where the unmodified lattice sum techniques |
1121 |
|
|
artificially accentuate the periodicity of the system in an |
1122 |
|
|
undesirable manner. There have been alterations to the standard Ewald |
1123 |
|
|
techniques, via corrections and reformulations, to compensate for |
1124 |
|
|
these systems; but the pairwise techniques discussed here require no |
1125 |
|
|
modifications, making them natural tools to tackle these problems. |
1126 |
|
|
Additionally, this transferability gives them benefits over other |
1127 |
|
|
pairwise methods, like reaction field, because estimations of physical |
1128 |
|
|
properties (e.g. the dielectric constant) are unnecessary. |
1129 |
chrisfen |
2605 |
|
1130 |
gezelter |
2656 |
If a researcher is using Monte Carlo simulations of large chemical |
1131 |
|
|
systems containing point charges, most structural features will be |
1132 |
|
|
accurately captured using the undamped {\sc sf} method or the {\sc sp} |
1133 |
|
|
method with an electrostatic damping of 0.2 \AA$^{-1}$. These methods |
1134 |
|
|
would also be appropriate for molecular dynamics simulations where the |
1135 |
|
|
data of interest is either structural or short-time dynamical |
1136 |
|
|
quantities. For long-time dynamics and collective motions, the safest |
1137 |
|
|
pairwise method we have evaluated is the {\sc sf} method with an |
1138 |
|
|
electrostatic damping between 0.2 and 0.25 |
1139 |
|
|
\AA$^{-1}$. |
1140 |
chrisfen |
2605 |
|
1141 |
gezelter |
2656 |
We are not suggesting that there is any flaw with the Ewald sum; in |
1142 |
|
|
fact, it is the standard by which these simple pairwise sums have been |
1143 |
|
|
judged. However, these results do suggest that in the typical |
1144 |
|
|
simulations performed today, the Ewald summation may no longer be |
1145 |
|
|
required to obtain the level of accuracy most researchers have come to |
1146 |
|
|
expect. |
1147 |
|
|
|
1148 |
chrisfen |
2575 |
\section{Acknowledgments} |
1149 |
gezelter |
2656 |
Support for this project was provided by the National Science |
1150 |
|
|
Foundation under grant CHE-0134881. The authors would like to thank |
1151 |
|
|
Steve Corcelli and Ed Maginn for helpful discussions and comments. |
1152 |
|
|
|
1153 |
chrisfen |
2594 |
\newpage |
1154 |
|
|
|
1155 |
gezelter |
2617 |
\bibliographystyle{jcp2} |
1156 |
chrisfen |
2575 |
\bibliography{electrostaticMethods} |
1157 |
|
|
|
1158 |
|
|
|
1159 |
|
|
\end{document} |