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