1 |
% ****** Start of file aipsamp.tex ****** |
2 |
% |
3 |
% This file is part of the AIP files in the AIP distribution for REVTeX 4. |
4 |
% Version 4.1 of REVTeX, October 2009 |
5 |
% |
6 |
% Copyright (c) 2009 American Institute of Physics. |
7 |
% |
8 |
% See the AIP README file for restrictions and more information. |
9 |
% |
10 |
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed |
11 |
% as well as the rest of the prerequisites for REVTeX 4.1 |
12 |
% |
13 |
% It also requires running BibTeX. The commands are as follows: |
14 |
% |
15 |
% 1) latex aipsamp |
16 |
% 2) bibtex aipsamp |
17 |
% 3) latex aipsamp |
18 |
% 4) latex aipsamp |
19 |
% |
20 |
% Use this file as a source of example code for your aip document. |
21 |
% Use the file aiptemplate.tex as a template for your document. |
22 |
\documentclass[% |
23 |
aip,jcp, |
24 |
amsmath,amssymb, |
25 |
preprint, |
26 |
%reprint,% |
27 |
%author-year,% |
28 |
%author-numerical,% |
29 |
]{revtex4-1} |
30 |
|
31 |
\usepackage{graphicx}% Include figure files |
32 |
\usepackage{dcolumn}% Align table columns on decimal point |
33 |
%\usepackage{bm}% bold math |
34 |
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math |
35 |
%\linenumbers\relax % Commence numbering lines |
36 |
\usepackage{amsmath} |
37 |
\usepackage{times} |
38 |
\usepackage{mathptmx} |
39 |
\usepackage{tabularx} |
40 |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
41 |
\usepackage{url} |
42 |
\usepackage[english]{babel} |
43 |
|
44 |
\newcolumntype{Y}{>{\centering\arraybackslash}X} |
45 |
|
46 |
\begin{document} |
47 |
|
48 |
%\preprint{AIP/123-QED} |
49 |
|
50 |
\title{Real space electrostatics for multipoles. II. Comparisons with |
51 |
the Ewald Sum} |
52 |
|
53 |
\author{Madan Lamichhane} |
54 |
\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} |
55 |
|
56 |
\author{Kathie E. Newman} |
57 |
\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} |
58 |
|
59 |
\author{J. Daniel Gezelter}% |
60 |
\email{gezelter@nd.edu.} |
61 |
\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556 |
62 |
} |
63 |
|
64 |
\date{\today} |
65 |
|
66 |
\begin{abstract} |
67 |
We report on tests of the shifted potential (SP), gradient shifted |
68 |
force (GSF), and Taylor shifted force (TSF) real-space methods for |
69 |
multipole interactions developed in the first paper in this series, |
70 |
using the multipolar Ewald sum as a reference method. The tests were |
71 |
carried out in a variety of condensed-phase environments designed to |
72 |
test up to quadrupole-quadrupole interactions. Comparisons of the |
73 |
energy differences between configurations, molecular forces, and |
74 |
torques were used to analyze how well the real-space models perform |
75 |
relative to the more computationally expensive Ewald treatment. We |
76 |
have also investigated the energy conservation, structural, and |
77 |
dynamical properties of the new methods in molecular dynamics |
78 |
simulations. The SP method shows excellent agreement with |
79 |
configurational energy differences, forces, and torques, and would |
80 |
be suitable for use in Monte Carlo calculations. Of the two new |
81 |
shifted-force methods, the GSF approach shows the best agreement |
82 |
with Ewald-derived energies, forces, and torques and also exhibits |
83 |
energy conservation properties that make it an excellent choice for |
84 |
efficient computation of electrostatic interactions in molecular |
85 |
dynamics simulations. Both SP and GSF are able to reproduce |
86 |
structural and dynamical properties in the liquid models with |
87 |
excellent fidelity. |
88 |
\end{abstract} |
89 |
|
90 |
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
91 |
% Classification Scheme. |
92 |
%\keywords{Electrostatics, Multipoles, Real-space} |
93 |
|
94 |
\maketitle |
95 |
|
96 |
\section{\label{sec:intro}Introduction} |
97 |
Computing the interactions between electrostatic sites is one of the |
98 |
most expensive aspects of molecular simulations. There have been |
99 |
significant efforts to develop practical, efficient and convergent |
100 |
methods for handling these interactions. Ewald's method is perhaps the |
101 |
best known and most accurate method for evaluating energies, forces, |
102 |
and torques in explicitly-periodic simulation cells. In this approach, |
103 |
the conditionally convergent electrostatic energy is converted into |
104 |
two absolutely convergent contributions, one which is carried out in |
105 |
real space with a cutoff radius, and one in reciprocal |
106 |
space.\cite{Ewald21,deLeeuw80,Smith81,Allen87} |
107 |
|
108 |
When carried out as originally formulated, the reciprocal-space |
109 |
portion of the Ewald sum exhibits relatively poor computational |
110 |
scaling, making it prohibitive for large systems. By utilizing a |
111 |
particle mesh and three dimensional fast Fourier transforms (FFT), the |
112 |
particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
113 |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) |
114 |
methods can decrease the computational cost from $O(N^2)$ down to $O(N |
115 |
\log |
116 |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb} |
117 |
|
118 |
Because of the artificial periodicity required for the Ewald sum, |
119 |
interfacial molecular systems such as membranes and liquid-vapor |
120 |
interfaces require modifications to the method. Parry's extension of |
121 |
the three dimensional Ewald sum is appropriate for slab |
122 |
geometries.\cite{Parry:1975if} Modified Ewald methods that were |
123 |
developed to handle two-dimensional (2-D) electrostatic |
124 |
interactions.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
125 |
These methods were originally quite computationally |
126 |
expensive.\cite{Spohr97,Yeh99} There have been several successful |
127 |
efforts that reduced the computational cost of 2-D lattice summations, |
128 |
bringing them more in line with the scaling for the full 3-D |
129 |
treatments.\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} The |
130 |
inherent periodicity required by the Ewald method can also be |
131 |
problematic in a number of protein/solvent and ionic solution |
132 |
environments.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq} |
133 |
|
134 |
\subsection{Real-space methods} |
135 |
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
136 |
method for calculating electrostatic interactions between point |
137 |
charges. They argued that the effective Coulomb interaction in most |
138 |
condensed phase systems is effectively short |
139 |
ranged.\cite{Wolf92,Wolf95} For an ordered lattice (e.g., when |
140 |
computing the Madelung constant of an ionic solid), the material can |
141 |
be considered as a set of ions interacting with neutral dipolar or |
142 |
quadrupolar ``molecules'' giving an effective distance dependence for |
143 |
the electrostatic interactions of $r^{-5}$ (see figure |
144 |
\ref{fig:schematic}). If one views the \ce{NaCl} crystal as a simple |
145 |
cubic (SC) structure with an octupolar \ce{(NaCl)4} basis, the |
146 |
electrostatic energy per ion converges more rapidly to the Madelung |
147 |
energy than the dipolar approximation.\cite{Wolf92} To find the |
148 |
correct Madelung constant, Lacman suggested that the NaCl structure |
149 |
could be constructed in a way that the finite crystal terminates with |
150 |
complete \ce{(NaCl)4} molecules.\cite{Lacman65} The central ion sees |
151 |
what is effectively a set of octupoles at large distances. These facts |
152 |
suggest that the Madelung constants are relatively short ranged for |
153 |
perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful |
154 |
application of Wolf's method can provide accurate estimates of |
155 |
Madelung constants using relatively short cutoff radii. |
156 |
|
157 |
Direct truncation of interactions at a cutoff radius creates numerical |
158 |
errors. Wolf \textit{et al.} suggest that truncation errors are due |
159 |
to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To |
160 |
neutralize this charge they proposed placing an image charge on the |
161 |
surface of the cutoff sphere for every real charge inside the cutoff. |
162 |
These charges are present for the evaluation of both the pair |
163 |
interaction energy and the force, although the force expression |
164 |
maintains a discontinuity at the cutoff sphere. In the original Wolf |
165 |
formulation, the total energy for the charge and image were not equal |
166 |
to the integral of the force expression, and as a result, the total |
167 |
energy would not be conserved in molecular dynamics (MD) |
168 |
simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennell and |
169 |
Gezelter later proposed shifted force variants of the Wolf method with |
170 |
commensurate force and energy expressions that do not exhibit this |
171 |
problem.\cite{Zahn:2002hc,Fennell:2006lq} Related real-space methods |
172 |
were also proposed by Chen \textit{et |
173 |
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
174 |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfully |
175 |
used additional neutralization of higher order moments for systems of |
176 |
point charges.\cite{Fukuda:2013sf} |
177 |
|
178 |
\begin{figure} |
179 |
\centering |
180 |
\includegraphics[width=\linewidth]{schematic.eps} |
181 |
\caption{Top: Ionic systems exhibit local clustering of dissimilar |
182 |
charges (in the smaller grey circle), so interactions are |
183 |
effectively charge-multipole at longer distances. With hard |
184 |
cutoffs, motion of individual charges in and out of the cutoff |
185 |
sphere can break the effective multipolar ordering. Bottom: |
186 |
dipolar crystals and fluids have a similar effective |
187 |
\textit{quadrupolar} ordering (in the smaller grey circles), and |
188 |
orientational averaging helps to reduce the effective range of the |
189 |
interactions in the fluid. Placement of reversed image multipoles |
190 |
on the surface of the cutoff sphere recovers the effective |
191 |
higher-order multipole behavior. \label{fig:schematic}} |
192 |
\end{figure} |
193 |
|
194 |
One can make a similar effective range argument for crystals of point |
195 |
\textit{multipoles}. The Luttinger and Tisza treatment of energy |
196 |
constants for dipolar lattices utilizes 24 basis vectors that contain |
197 |
dipoles at the eight corners of a unit cube.\cite{LT} Only three of |
198 |
these basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
199 |
moments, while the rest have zero net dipole and retain contributions |
200 |
only from higher order multipoles. The lowest-energy crystalline |
201 |
structures are built out of basis vectors that have only residual |
202 |
quadrupolar moments (e.g. the $Z_5$ array). In these low energy |
203 |
structures, the effective interaction between a dipole at the center |
204 |
of a crystal and a group of eight dipoles farther away is |
205 |
significantly shorter ranged than the $r^{-3}$ that one would expect |
206 |
for raw dipole-dipole interactions. Only in crystals which retain a |
207 |
bulk dipole moment (e.g. ferroelectrics) does the analogy with the |
208 |
ionic crystal break down -- ferroelectric dipolar crystals can exist, |
209 |
while ionic crystals with net charge in each unit cell would be |
210 |
unstable. |
211 |
|
212 |
In ionic crystals, real-space truncation can break the effective |
213 |
multipolar arrangements (see Fig. \ref{fig:schematic}), causing |
214 |
significant swings in the electrostatic energy as individual ions move |
215 |
back and forth across the boundary. This is why the image charges are |
216 |
necessary for the Wolf sum to exhibit rapid convergence. Similarly, |
217 |
the real-space truncation of point multipole interactions breaks |
218 |
higher order multipole arrangements, and image multipoles are required |
219 |
for real-space treatments of electrostatic energies. |
220 |
|
221 |
The shorter effective range of electrostatic interactions is not |
222 |
limited to perfect crystals, but can also apply in disordered fluids. |
223 |
Even at elevated temperatures, there is local charge balance in an |
224 |
ionic liquid, where each positive ion has surroundings dominated by |
225 |
negative ions and vice versa. The reversed-charge images on the |
226 |
cutoff sphere that are integral to the Wolf and damped shifted force |
227 |
(DSF) approaches retain the effective multipolar interactions as the |
228 |
charges traverse the cutoff boundary. |
229 |
|
230 |
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
231 |
significant orientational averaging that additionally reduces the |
232 |
effect of long-range multipolar interactions. The image multipoles |
233 |
that are introduced in the Taylor shifted force (TSF), gradient |
234 |
shifted force (GSF), and shifted potential (SP) methods mimic this effect |
235 |
and reduce the effective range of the multipolar interactions as |
236 |
interacting molecules traverse each other's cutoff boundaries. |
237 |
|
238 |
Forces and torques acting on atomic sites are fundamental in driving |
239 |
dynamics in molecular simulations, and the DSF energy kernel provides |
240 |
consistent energies and forces on charged atoms within the cutoff |
241 |
sphere. Both the energy and the force go smoothly to zero as an atom |
242 |
approaches the cutoff radius. The comparisons of the accuracy these |
243 |
quantities between the DSF kernel and SPME was surprisingly |
244 |
good.\cite{Fennell:2006lq} As a result, the DSF method has seen |
245 |
increasing use in molecular systems with relatively uniform charge |
246 |
densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,Vergne13} |
247 |
|
248 |
\subsection{The damping function} |
249 |
The damping function has been discussed in detail in the first paper |
250 |
of this series.\cite{PaperI} The $1/r$ Coulombic kernel for the |
251 |
interactions between point charges can be replaced by the |
252 |
complementary error function $\mathrm{erfc}(\alpha r)/r$ to accelerate |
253 |
convergence, where $\alpha$ is a damping parameter with units of |
254 |
inverse distance. Altering the value of $\alpha$ is equivalent to |
255 |
changing the width of Gaussian charge distributions that replace each |
256 |
point charge, as Coulomb integrals with Gaussian charge distributions |
257 |
produce complementary error functions when truncated at a finite |
258 |
distance. |
259 |
|
260 |
With moderate damping coefficients, $\alpha \sim 0.2$ \AA$^{-1}$, the DSF method |
261 |
produced very good agreement with SPME for interaction energies, |
262 |
forces and torques for charge-charge |
263 |
interactions.\cite{Fennell:2006lq} |
264 |
|
265 |
\subsection{Point multipoles in molecular modeling} |
266 |
Coarse-graining approaches which treat entire molecular subsystems as |
267 |
a single rigid body are now widely used. A common feature of many |
268 |
coarse-graining approaches is simplification of the electrostatic |
269 |
interactions between bodies so that fewer site-site interactions are |
270 |
required to compute configurational |
271 |
energies.\cite{Ren06,Essex10,Essex11} |
272 |
|
273 |
Additionally, because electrons in a molecule are not localized at |
274 |
specific points, the assignment of partial charges to atomic centers |
275 |
is always an approximation. For increased accuracy, atomic sites can |
276 |
also be assigned point multipoles and polarizabilities. Recently, |
277 |
water has been modeled with point multipoles up to octupolar order |
278 |
using the soft sticky dipole-quadrupole-octupole (SSDQO) |
279 |
model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
280 |
multipoles up to quadrupolar order have also been coupled with point |
281 |
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
282 |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However, |
283 |
truncating point multipoles without smoothing the forces and torques |
284 |
can create energy conservation issues in molecular dynamics |
285 |
simulations. |
286 |
|
287 |
In this paper we test a set of real-space methods that were developed |
288 |
for point multipolar interactions. These methods extend the damped |
289 |
shifted force (DSF) and Wolf methods originally developed for |
290 |
charge-charge interactions and generalize them for higher order |
291 |
multipoles. The detailed mathematical development of these methods |
292 |
has been presented in the first paper in this series, while this work |
293 |
covers the testing of energies, forces, torques, and energy |
294 |
conservation properties of the methods in realistic simulation |
295 |
environments. In all cases, the methods are compared with the |
296 |
reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98} |
297 |
|
298 |
|
299 |
\section{\label{sec:method}Review of Methods} |
300 |
Any real-space electrostatic method that is suitable for MD |
301 |
simulations should have the electrostatic energy, forces and torques |
302 |
between two sites go smoothly to zero as the distance between the |
303 |
sites, $r_{ab}$ approaches the cutoff radius, $r_c$. Requiring |
304 |
this continuity at the cutoff is essential for energy conservation in |
305 |
MD simulations. The mathematical details of the shifted potential |
306 |
(SP), gradient-shifted-force (GSF) and Taylor shifted-force (TSF) |
307 |
methods have been discussed in detail in the previous paper in this |
308 |
series.\cite{PaperI} Here we briefly review the new methods and |
309 |
describe their essential features. |
310 |
|
311 |
\subsection{Taylor-shifted force (TSF)} |
312 |
|
313 |
The electrostatic potential energy between point multipoles can be |
314 |
expressed as the product of two multipole operators and a Coulombic |
315 |
kernel, |
316 |
\begin{equation} |
317 |
U_{ab}(r)= M_{a} M_{b} \frac{1}{r} \label{kernel}. |
318 |
\end{equation} |
319 |
where the multipole operator for site $a$, $M_{a}$, is |
320 |
expressed in terms of the point charge, $C_{a}$, dipole, ${\bf D}_{a}$, and quadrupole, $\mathsf{Q}_{a}$, for object |
321 |
$a$, etc. |
322 |
|
323 |
The TSF potential for any multipole-multipole interaction can be |
324 |
written |
325 |
\begin{equation} |
326 |
U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) |
327 |
\label{generic} |
328 |
\end{equation} |
329 |
where $f_n(r)$ is a shifted kernel that is appropriate for the order |
330 |
of the interaction (see Ref. \onlinecite{PaperI}), with $n=0$ for |
331 |
charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole |
332 |
and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for |
333 |
quadrupole-quadrupole. To ensure smooth convergence of the energy, |
334 |
force, and torques, a Taylor expansion with $n$ terms must be |
335 |
performed at cutoff radius ($r_c$) to obtain $f_n(r)$. |
336 |
|
337 |
For multipole-multipole interactions, following this procedure results |
338 |
in separate radial functions for each of the distinct orientational |
339 |
contributions to the potential, and ensures that the forces and |
340 |
torques from each of these contributions will vanish at the cutoff |
341 |
radius. For example, the direct dipole dot product |
342 |
($\mathbf{D}_{a} |
343 |
\cdot \mathbf{D}_{b}$) is treated differently than the dipole-distance |
344 |
dot products: |
345 |
\begin{equation} |
346 |
U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r)= -\frac{1}{4\pi \epsilon_0} \left[ \left( |
347 |
\mathbf{D}_{a} \cdot |
348 |
\mathbf{D}_{b} \right) v_{21}(r) + |
349 |
\left( \mathbf{D}_{a} \cdot \hat{\mathbf{r}} \right) |
350 |
\left( \mathbf{D}_{b} \cdot \hat{\mathbf{r}} \right) v_{22}(r) \right] |
351 |
\end{equation} |
352 |
|
353 |
For the Taylor shifted (TSF) method with the undamped kernel, |
354 |
$v_{21}(r) = -\frac{1}{r^3} + \frac{3 r}{r_c^4} - \frac{8}{r_c^3} + |
355 |
\frac{6}{r r_c^2}$ and $v_{22}(r) = \frac{3}{r^3} + \frac{3 r}{r_c^4} |
356 |
- \frac{6}{r r_c^2}$. In these functions, one can easily see the |
357 |
connection to unmodified electrostatics as well as the smooth |
358 |
transition to zero in both these functions as $r\rightarrow r_c$. The |
359 |
electrostatic forces and torques acting on the central multipole due |
360 |
to another site within the cutoff sphere are derived from |
361 |
Eq.~\ref{generic}, accounting for the appropriate number of |
362 |
derivatives. Complete energy, force, and torque expressions are |
363 |
presented in the first paper in this series (Reference |
364 |
\onlinecite{PaperI}). |
365 |
|
366 |
\subsection{Gradient-shifted force (GSF)} |
367 |
|
368 |
A second (and conceptually simpler) method involves shifting the |
369 |
gradient of the raw Coulomb potential for each particular multipole |
370 |
order. For example, the raw dipole-dipole potential energy may be |
371 |
shifted smoothly by finding the gradient for two interacting dipoles |
372 |
which have been projected onto the surface of the cutoff sphere |
373 |
without changing their relative orientation, |
374 |
\begin{equation} |
375 |
U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r) = U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r) - |
376 |
U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r_c) |
377 |
- (r_{ab}-r_c) ~~~\hat{\mathbf{r}}_{ab} \cdot |
378 |
\nabla U_{\mathbf{D}_{a}\mathbf{D}_{b}}(r_c). |
379 |
\end{equation} |
380 |
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{a}$ and $\mathbf{D}_{b}$, are retained at the cutoff distance |
381 |
(although the signs are reversed for the dipole that has been |
382 |
projected onto the cutoff sphere). In many ways, this simpler |
383 |
approach is closer in spirit to the original shifted force method, in |
384 |
that it projects a neutralizing multipole (and the resulting forces |
385 |
from this multipole) onto a cutoff sphere. The resulting functional |
386 |
forms for the potentials, forces, and torques turn out to be quite |
387 |
similar in form to the Taylor-shifted approach, although the radial |
388 |
contributions are significantly less perturbed by the gradient-shifted |
389 |
approach than they are in the Taylor-shifted method. |
390 |
|
391 |
For the gradient shifted (GSF) method with the undamped kernel, |
392 |
$v_{21}(r) = -\frac{1}{r^3} - \frac{3 r}{r_c^4} + \frac{4}{r_c^3}$ and |
393 |
$v_{22}(r) = \frac{3}{r^3} + \frac{9 r}{r_c^4} - \frac{12}{r_c^3}$. |
394 |
Again, these functions go smoothly to zero as $r\rightarrow r_c$, and |
395 |
because the Taylor expansion retains only one term, they are |
396 |
significantly less perturbed than the TSF functions. |
397 |
|
398 |
In general, the gradient shifted potential between a central multipole |
399 |
and any multipolar site inside the cutoff radius is given by, |
400 |
\begin{equation} |
401 |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) - |
402 |
U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) - (r-r_c) |
403 |
\hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right] |
404 |
\label{generic2} |
405 |
\end{equation} |
406 |
where the sum describes a separate force-shifting that is applied to |
407 |
each orientational contribution to the energy. In this expression, |
408 |
$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles |
409 |
($a$ and $b$) in space, and $\mathsf{A}$ and $\mathsf{B}$ |
410 |
represent the orientations the multipoles. |
411 |
|
412 |
The third term converges more rapidly than the first two terms as a |
413 |
function of radius, hence the contribution of the third term is very |
414 |
small for large cutoff radii. The force and torque derived from |
415 |
Eq. \ref{generic2} are consistent with the energy expression and |
416 |
approach zero as $r \rightarrow r_c$. Both the GSF and TSF methods |
417 |
can be considered generalizations of the original DSF method for |
418 |
higher order multipole interactions. GSF and TSF are also identical up |
419 |
to the charge-dipole interaction but generate different expressions in |
420 |
the energy, force and torque for higher order multipole-multipole |
421 |
interactions. Complete energy, force, and torque expressions for the |
422 |
GSF potential are presented in the first paper in this series |
423 |
(Reference~\onlinecite{PaperI}). |
424 |
|
425 |
|
426 |
\subsection{Shifted potential (SP) } |
427 |
A discontinuous truncation of the electrostatic potential at the |
428 |
cutoff sphere introduces a severe artifact (oscillation in the |
429 |
electrostatic energy) even for molecules with the higher-order |
430 |
multipoles.\cite{PaperI} We have also formulated an extension of the |
431 |
Wolf approach for point multipoles by simply projecting the image |
432 |
multipole onto the surface of the cutoff sphere, and including the |
433 |
interactions with the central multipole and the image. This |
434 |
effectively shifts the total potential to zero at the cutoff radius, |
435 |
\begin{equation} |
436 |
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \mathsf{A}, \mathsf{B}) - |
437 |
U(r_c \hat{\mathbf{r}},\mathsf{A}, \mathsf{B}) \right] |
438 |
\label{eq:SP} |
439 |
\end{equation} |
440 |
where the sum describes separate potential shifting that is done for |
441 |
each orientational contribution to the energy (e.g. the direct dipole |
442 |
product contribution is shifted {\it separately} from the |
443 |
dipole-distance terms in dipole-dipole interactions). Note that this |
444 |
is not a simple shifting of the total potential at $r_c$. Each radial |
445 |
contribution is shifted separately. One consequence of this is that |
446 |
multipoles that reorient after leaving the cutoff sphere can re-enter |
447 |
the cutoff sphere without perturbing the total energy. |
448 |
|
449 |
For the shifted potential (SP) method with the undamped kernel, |
450 |
$v_{21}(r) = -\frac{1}{r^3} + \frac{1}{r_c^3}$ and $v_{22}(r) = |
451 |
\frac{3}{r^3} - \frac{3}{r_c^3}$. The potential energy between a |
452 |
central multipole and other multipolar sites goes smoothly to zero as |
453 |
$r \rightarrow r_c$. However, the force and torque obtained from the |
454 |
shifted potential (SP) are discontinuous at $r_c$. MD simulations |
455 |
will still experience energy drift while operating under the SP |
456 |
potential, but it may be suitable for Monte Carlo approaches where the |
457 |
configurational energy differences are the primary quantity of |
458 |
interest. |
459 |
|
460 |
\subsection{The Self Term} |
461 |
In the TSF, GSF, and SP methods, a self-interaction is retained for |
462 |
the central multipole interacting with its own image on the surface of |
463 |
the cutoff sphere. This self interaction is nearly identical with the |
464 |
self-terms that arise in the Ewald sum for multipoles. Complete |
465 |
expressions for the self terms are presented in the first paper in |
466 |
this series (Reference \onlinecite{PaperI}). |
467 |
|
468 |
|
469 |
\section{\label{sec:methodology}Methodology} |
470 |
|
471 |
To understand how the real-space multipole methods behave in computer |
472 |
simulations, it is vital to test against established methods for |
473 |
computing electrostatic interactions in periodic systems, and to |
474 |
evaluate the size and sources of any errors that arise from the |
475 |
real-space cutoffs. In the first paper of this series, we compared |
476 |
the dipolar and quadrupolar energy expressions against analytic |
477 |
expressions for ordered dipolar and quadrupolar |
478 |
arrays.\cite{Sauer,LT,Nagai01081960,Nagai01091963} In this work, we |
479 |
used the multipolar Ewald sum as a reference method for comparing |
480 |
energies, forces, and torques for molecular models that mimic |
481 |
disordered and ordered condensed-phase systems. The parameters used |
482 |
in the test cases are given in table~\ref{tab:pars}. |
483 |
|
484 |
\begin{table} |
485 |
\caption{The parameters used in the systems used to evaluate the new |
486 |
real-space methods. The most comprehensive test was a liquid |
487 |
composed of 2000 soft DQ liquid molecules with 48 dissolved ions (24 \ce{Na+} and 24 \ce{Cl-} |
488 |
ions). This test exercises all orders of the multipolar |
489 |
interactions developed in the first paper.\label{tab:pars}} |
490 |
\begin{tabularx}{\textwidth}{r|cc|YYccc|Yccc} \hline |
491 |
& \multicolumn{2}{c|}{LJ parameters} & |
492 |
\multicolumn{5}{c|}{Electrostatic moments} & & & & \\ |
493 |
Test system & $\sigma$& $\epsilon$ & $C$ & $D$ & |
494 |
$Q_{xx}$ & $Q_{yy}$ & $Q_{zz}$ & mass & $I_{xx}$ & $I_{yy}$ & |
495 |
$I_{zz}$ \\ \cline{6-8}\cline{10-12} |
496 |
& (\AA) & (kcal/mol) & (e) & (debye) & \multicolumn{3}{c|}{(debye \AA)} & (amu) & \multicolumn{3}{c}{(amu |
497 |
\AA\textsuperscript{2})} \\ \hline |
498 |
Soft Dipolar fluid & 3.051 & 0.152 & & 2.35 & & & & 18.0153 & 1.77&0.6145& 1.155 \\ |
499 |
Soft Dipolar solid & 2.837 & 1.0 & & 2.35 & & & & $10^4$ & 17.6 &17.6 & 0 \\ |
500 |
Soft Quadrupolar fluid & 3.051 & 0.152 & & & -1&-1&-2.5 & 18.0153 & 1.77&0.6145&1.155 \\ |
501 |
Soft Quadrupolar solid & 2.837 & 1.0 & & & -1&-1&-2.5 & $10^4$ & 17.6&17.6&0 \\ |
502 |
Soft DQ liquid & 3.051 & 0.152 & & 2.35 & -1.35&0&-0.68 & 18.0153 & 1.77&0.6145&1.155 \\ |
503 |
\ce{Na+} & 2.579 & 0.118 & +1& & & & & 22.99 & & &\\ |
504 |
\ce{Cl-} & 4.445 & 0.1 & -1& & & & & 35.4527& & & \\ \hline |
505 |
\end{tabularx} |
506 |
\end{table} |
507 |
The systems consist of pure multipolar solids (both dipole and |
508 |
quadrupole), pure multipolar liquids (both dipole and quadrupole), a |
509 |
fluid composed of sites containing both dipoles and quadrupoles |
510 |
simultaneously, and a final test case that includes ions with point |
511 |
charges in addition to the multipolar fluid. The solid-phase |
512 |
parameters were chosen so that the systems can explore some |
513 |
orientational freedom for the multipolar sites, while maintaining |
514 |
relatively strict translational order. The soft DQ liquid model used |
515 |
here based loosely on the SSDQO water |
516 |
model,\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} but is not itself a |
517 |
particularly accurate water model. However, the soft DQ model does |
518 |
test dipole-dipole, dipole-quadrupole, and quadrupole-quadrupole |
519 |
interactions at roughly the same magnitudes. The last test case, a |
520 |
soft DQ liquid with dissolved ions, exercises \textit{all} levels of |
521 |
the multipole-multipole interactions we have derived so far and |
522 |
represents the most complete test of the new methods. |
523 |
|
524 |
In the following section, we present results for the total |
525 |
electrostatic energy, as well as the electrostatic contributions to |
526 |
the force and torque on each molecule. These quantities have been |
527 |
computed using the SP, TSF, and GSF methods, as well as a hard cutoff, |
528 |
and have been compared with the values obtained from the multipolar |
529 |
Ewald sum. In Monte Carlo (MC) simulations, the energy differences |
530 |
between two configurations is the primary quantity that governs how |
531 |
the simulation proceeds. These differences are the most important |
532 |
indicators of the reliability of a method even if the absolute |
533 |
energies are not exact. For each of the multipolar systems listed |
534 |
above, we have compared the change in electrostatic potential energy |
535 |
($\Delta E$) between 250 statistically-independent configurations. In |
536 |
molecular dynamics (MD) simulations, the forces and torques govern the |
537 |
behavior of the simulation, so we also compute the electrostatic |
538 |
contributions to the forces and torques. |
539 |
|
540 |
\subsection{Implementation} |
541 |
The real-space methods developed in the first paper in this series |
542 |
have been implemented in our group's open source molecular simulation |
543 |
program, OpenMD,\cite{Meineke05,openmd} which was used for all calculations in |
544 |
this work. The complementary error function can be a relatively slow |
545 |
function on some processors, so all of the radial functions are |
546 |
precomputed on a fine grid and are spline-interpolated to provide |
547 |
values when required. |
548 |
|
549 |
Using the same simulation code, we compare to a multipolar Ewald sum |
550 |
with a reciprocal space cutoff, $k_\mathrm{max} = 7$. Our version of |
551 |
the Ewald sum is a re-implementation of the algorithm originally |
552 |
proposed by Smith that does not use the particle mesh or smoothing |
553 |
approximations.\cite{Smith82,Smith98} This implementation was tested |
554 |
extensively against the analytic energy constants for the multipolar |
555 |
lattices that are discussed in reference \onlinecite{PaperI}. In all |
556 |
cases discussed below, the quantities being compared are the |
557 |
electrostatic contributions to energies, force, and torques. All |
558 |
other contributions to these quantities (i.e. from Lennard-Jones |
559 |
interactions) are removed prior to the comparisons. |
560 |
|
561 |
The convergence parameter ($\alpha$) also plays a role in the balance |
562 |
of the real-space and reciprocal-space portions of the Ewald |
563 |
calculation. Typical molecular mechanics packages set this to a value |
564 |
that depends on the cutoff radius and a tolerance (typically less than |
565 |
$1 \times 10^{-4}$ kcal/mol). Smaller tolerances are typically |
566 |
associated with increasing accuracy at the expense of computational |
567 |
time spent on the reciprocal-space portion of the |
568 |
summation.\cite{Perram88,Essmann:1995pb} A default tolerance of $1 \times |
569 |
10^{-8}$ kcal/mol was used in all Ewald calculations, resulting in |
570 |
Ewald coefficient 0.3119 \AA$^{-1}$ for a cutoff radius of 12 \AA. |
571 |
|
572 |
The real-space models have self-interactions that provide |
573 |
contributions to the energies only. Although the self interaction is |
574 |
a rapid calculation, we note that in systems with fluctuating charges |
575 |
or point polarizabilities, the self-term is not static and must be |
576 |
recomputed at each time step. |
577 |
|
578 |
\subsection{Model systems} |
579 |
To sample independent configurations of the multipolar crystals, body |
580 |
centered cubic (bcc) crystals, which exhibit the minimum energy |
581 |
structures for point dipoles, were generated using 3,456 molecules. |
582 |
The multipoles were translationally locked in their respective crystal |
583 |
sites for equilibration at a relatively low temperature (50K) so that |
584 |
dipoles or quadrupoles could freely explore all accessible |
585 |
orientations. The translational constraints were then removed, the |
586 |
systems were re-equilibrated, and the crystals were simulated for an |
587 |
additional 10 ps in the microcanonical (NVE) ensemble with an average |
588 |
temperature of 50 K. The balance between moments of inertia and |
589 |
particle mass were chosen to allow orientational sampling without |
590 |
significant translational motion. Configurations were sampled at |
591 |
equal time intervals in order to compare configurational energy |
592 |
differences. The crystals were simulated far from the melting point |
593 |
in order to avoid translational deformation away of the ideal lattice |
594 |
geometry. |
595 |
|
596 |
For dipolar, quadrupolar, and mixed-multipole \textit{liquid} |
597 |
simulations, each system was created with 2,048 randomly-oriented |
598 |
molecules. These were equilibrated at a temperature of 300K for 1 ns. |
599 |
Each system was then simulated for 1 ns in the microcanonical (NVE) |
600 |
ensemble with the Dullweber, Leimkuhler, and McLachlan (DLM) |
601 |
symplectic splitting integrator using 1 fs |
602 |
timesteps.\cite{Dullweber1997} We collected 250 different |
603 |
configurations at equal time intervals. For the liquid system that |
604 |
included ionic species, we converted 48 randomly-distributed molecules |
605 |
into 24 \ce{Na+} and 24 \ce{Cl-} ions and re-equilibrated. After |
606 |
equilibration, the system was run under the same conditions for 1 |
607 |
ns. A total of 250 configurations were collected. In the following |
608 |
comparisons of energies, forces, and torques, the Lennard-Jones |
609 |
potentials were turned off and only the purely electrostatic |
610 |
quantities were compared with the same values obtained via the Ewald |
611 |
sum. |
612 |
|
613 |
\subsection{Accuracy of Energy Differences, Forces and Torques} |
614 |
The pairwise summation techniques (outlined above) were evaluated for |
615 |
use in MC simulations by studying the energy differences between |
616 |
different configurations. We took the Ewald-computed energy |
617 |
difference between two conformations to be the correct behavior. An |
618 |
ideal performance by one of the new methods would reproduce these |
619 |
energy differences exactly. The configurational energies being used |
620 |
here contain only contributions from electrostatic interactions. |
621 |
Lennard-Jones interactions were omitted from the comparison as they |
622 |
should be identical for all methods. |
623 |
|
624 |
Since none of the real-space methods provide exact energy differences, |
625 |
we used least square regressions analysis for the six different |
626 |
molecular systems to compare $\Delta E$ from Hard, SP, GSF, and TSF |
627 |
with the multipolar Ewald reference method. A result of unity for |
628 |
both the correlation (slope) and coefficient of determination ($R^2$) |
629 |
for these regressions would indicate perfect agreement between the |
630 |
real-space method and the multipolar Ewald sum. |
631 |
|
632 |
Molecular systems were run long enough to explore independent |
633 |
configurations and 250 configurations were recorded for comparison. |
634 |
Each system provided 31,125 energy differences for a total of 186,750 |
635 |
data points. Similarly, the magnitudes of the forces and torques have |
636 |
also been compared using least squares regression analysis. In the |
637 |
forces and torques comparison, the magnitudes of the forces acting in |
638 |
each molecule for each configuration were evaluated. For example, our |
639 |
dipolar liquid simulation contains 2048 molecules and there are 250 |
640 |
different configurations for each system resulting in 3,072,000 data |
641 |
points for comparison of forces and torques. |
642 |
|
643 |
\subsection{Analysis of vector quantities} |
644 |
Getting the magnitudes of the force and torque vectors correct is only |
645 |
part of the issue for carrying out accurate molecular dynamics |
646 |
simulations. Because the real space methods reweight the different |
647 |
orientational contributions to the energies, it is also important to |
648 |
understand how the methods impact the \textit{directionality} of the |
649 |
force and torque vectors. Fisher developed a probability density |
650 |
function to analyse directional data sets, |
651 |
\begin{equation} |
652 |
p_f(\theta) = \frac{\kappa}{2 \sinh\kappa}\sin\theta e^{\kappa \cos\theta} |
653 |
\label{eq:pdf} |
654 |
\end{equation} |
655 |
where $\kappa$ measures directional dispersion of the data around the |
656 |
mean direction.\cite{fisher53} This quantity $(\kappa)$ can be |
657 |
estimated as a reciprocal of the circular variance.\cite{Allen91} To |
658 |
quantify the directional error, forces obtained from the Ewald sum |
659 |
were taken as the mean (or correct) direction and the angle between |
660 |
the forces obtained via the Ewald sum and the real-space methods were |
661 |
evaluated, |
662 |
\begin{equation} |
663 |
\cos\theta_i = \frac{\mathbf{f}_i^\mathrm{~Ewald} \cdot |
664 |
\mathbf{f}_i^\mathrm{~GSF}}{\left|\mathbf{f}_i^\mathrm{~Ewald}\right| \left|\mathbf{f}_i^\mathrm{~GSF}\right|} |
665 |
\end{equation} |
666 |
The total angular displacement of the vectors was calculated as, |
667 |
\begin{equation} |
668 |
R = \sqrt{\left(\sum\limits_{i=1}^N \cos\theta_i\right)^2 + \left(\sum\limits_{i=1}^N \sin\theta_i\right)^2} |
669 |
\label{eq:displacement} |
670 |
\end{equation} |
671 |
where $N$ is number of force vectors. The circular variance is |
672 |
defined as |
673 |
\begin{equation} |
674 |
\mathrm{Var}(\theta) \approx 1/\kappa = 1 - R/N |
675 |
\end{equation} |
676 |
The circular variance takes on values between from 0 to 1, with 0 |
677 |
indicating a perfect directional match between the Ewald force vectors |
678 |
and the real-space forces. Lower values of $\mathrm{Var}(\theta)$ |
679 |
correspond to higher values of $\kappa$, which indicates tighter |
680 |
clustering of the real-space force vectors around the Ewald forces. |
681 |
|
682 |
A similar analysis was carried out for the electrostatic contribution |
683 |
to the molecular torques as well as forces. |
684 |
|
685 |
\subsection{Energy conservation} |
686 |
To test conservation the energy for the methods, the mixed molecular |
687 |
system of 2000 soft DQ liquid molecules with 24 \ce{Na+} and 24 |
688 |
\ce{Cl-} ions was run for 1 ns in the microcanonical ensemble at an |
689 |
average temperature of 300K. Each of the different electrostatic |
690 |
methods (Ewald, Hard, SP, GSF, and TSF) was tested for a range of |
691 |
different damping values. The molecular system was started with same |
692 |
initial positions and velocities for all cutoff methods. The energy |
693 |
drift ($\delta E_1$) and standard deviation of the energy about the |
694 |
slope ($\delta E_0$) were evaluated from the total energy of the |
695 |
system as a function of time. Although both measures are valuable at |
696 |
investigating new methods for molecular dynamics, a useful interaction |
697 |
model must allow for long simulation times with minimal energy drift. |
698 |
|
699 |
\section{\label{sec:result}RESULTS} |
700 |
\subsection{Configurational energy differences} |
701 |
|
702 |
\begin{figure} |
703 |
\centering |
704 |
\includegraphics[width=0.6\linewidth]{energyPlot_slopeCorrelation_combined.eps} |
705 |
\caption{Statistical analysis of the quality of configurational |
706 |
energy differences for the real-space electrostatic methods |
707 |
compared with the reference Ewald sum. Results with a value equal |
708 |
to 1 (dashed line) indicate $\Delta E$ values indistinguishable |
709 |
from those obtained using the multipolar Ewald sum. Different |
710 |
values of the cutoff radius are indicated with different symbols |
711 |
(9~\AA\ = circles, 12~\AA\ = squares, and 15~\AA\ = inverted |
712 |
triangles).\label{fig:slopeCorr_energy}} |
713 |
\end{figure} |
714 |
|
715 |
The combined coefficient of determination and slope for all six |
716 |
systems is shown in Figure ~\ref{fig:slopeCorr_energy}. Most of the |
717 |
methods reproduce the Ewald configurational energy differences with |
718 |
remarkable fidelity. Undamped hard cutoffs introduce a significant |
719 |
amount of random scatter in the energy differences which is apparent |
720 |
in the reduced value of $R^2$ for this method. This can be easily |
721 |
understood as configurations which exhibit small traversals of a few |
722 |
dipoles or quadrupoles out of the cutoff sphere will see large energy |
723 |
jumps when hard cutoffs are used. The orientations of the multipoles |
724 |
(particularly in the ordered crystals) mean that these energy jumps |
725 |
can go in either direction, producing a significant amount of random |
726 |
scatter, but no systematic error. |
727 |
|
728 |
The TSF method produces energy differences that are highly correlated |
729 |
with the Ewald results, but it also introduces a significant |
730 |
systematic bias in the values of the energies, particularly for |
731 |
smaller cutoff values. The TSF method alters the distance dependence |
732 |
of different orientational contributions to the energy in a |
733 |
non-uniform way, so the size of the cutoff sphere can have a large |
734 |
effect, particularly for the crystalline systems. |
735 |
|
736 |
Both the SP and GSF methods appear to reproduce the Ewald results with |
737 |
excellent fidelity, particularly for moderate damping ($\alpha \approx |
738 |
0.2$~\AA$^{-1}$) and with a commonly-used cutoff value ($r_c = |
739 |
12$~\AA). With the exception of the undamped hard cutoff, and the TSF |
740 |
method with short cutoffs, all of the methods would be appropriate for |
741 |
use in Monte Carlo simulations. |
742 |
|
743 |
\subsection{Magnitude of the force and torque vectors} |
744 |
|
745 |
The comparisons of the magnitudes of the forces and torques for the |
746 |
data accumulated from all six systems are shown in Figures |
747 |
~\ref{fig:slopeCorr_force} and \ref{fig:slopeCorr_torque}. The |
748 |
correlation and slope for the forces agree well with the Ewald sum |
749 |
even for the hard cutoffs. |
750 |
|
751 |
For systems of molecules with only multipolar interactions, the pair |
752 |
energy contributions are quite short ranged. Moreover, the force |
753 |
decays more rapidly than the electrostatic energy, hence the hard |
754 |
cutoff method can also produce reasonable agreement for this quantity. |
755 |
Although the pure cutoff gives reasonably good electrostatic forces |
756 |
for pairs of molecules included within each other's cutoff spheres, |
757 |
the discontinuity in the force at the cutoff radius can potentially |
758 |
cause energy conservation problems as molecules enter and leave the |
759 |
cutoff spheres. This is discussed in detail in section |
760 |
\ref{sec:conservation}. |
761 |
|
762 |
The two shifted-force methods (GSF and TSF) exhibit a small amount of |
763 |
systematic variation and scatter compared with the Ewald forces. The |
764 |
shifted-force models intentionally perturb the forces between pairs of |
765 |
molecules inside each other's cutoff spheres in order to correct the |
766 |
energy conservation issues, and this perturbation is evident in the |
767 |
statistics accumulated for the molecular forces. The GSF |
768 |
perturbations are minimal, particularly for moderate damping and |
769 |
commonly-used cutoff values ($r_c = 12$~\AA). The TSF method shows |
770 |
reasonable agreement in $R^2$, but again the systematic error in the |
771 |
forces is concerning if replication of Ewald forces is desired. |
772 |
|
773 |
It is important to note that the forces and torques from the SP and |
774 |
the Hard cutoffs are not identical. The SP method shifts each |
775 |
orientational contribution separately (e.g. the dipole-dipole dot |
776 |
product is shifted by a different function than the dipole-distance |
777 |
products), while the hard cutoff contains no orientation-dependent |
778 |
shifting. The forces and torques for these methods therefore diverge |
779 |
for multipoles even though the forces for point charges are identical. |
780 |
|
781 |
\begin{figure} |
782 |
\centering |
783 |
\includegraphics[width=0.6\linewidth]{forcePlot_slopeCorrelation_combined.eps} |
784 |
\caption{Statistical analysis of the quality of the force vector |
785 |
magnitudes for the real-space electrostatic methods compared with |
786 |
the reference Ewald sum. Results with a value equal to 1 (dashed |
787 |
line) indicate force magnitude values indistinguishable from those |
788 |
obtained using the multipolar Ewald sum. Different values of the |
789 |
cutoff radius are indicated with different symbols (9~\AA\ = |
790 |
circles, 12~\AA\ = squares, and 15~\AA\ = inverted |
791 |
triangles).\label{fig:slopeCorr_force}} |
792 |
\end{figure} |
793 |
|
794 |
|
795 |
\begin{figure} |
796 |
\centering |
797 |
\includegraphics[width=0.6\linewidth]{torquePlot_slopeCorrelation_combined.eps} |
798 |
\caption{Statistical analysis of the quality of the torque vector |
799 |
magnitudes for the real-space electrostatic methods compared with |
800 |
the reference Ewald sum. Results with a value equal to 1 (dashed |
801 |
line) indicate force magnitude values indistinguishable from those |
802 |
obtained using the multipolar Ewald sum. Different values of the |
803 |
cutoff radius are indicated with different symbols (9~\AA\ = |
804 |
circles, 12~\AA\ = squares, and 15~\AA\ = inverted |
805 |
triangles).\label{fig:slopeCorr_torque}} |
806 |
\end{figure} |
807 |
|
808 |
The torques (Fig. \ref{fig:slopeCorr_torque}) appear to be |
809 |
significantly influenced by the choice of real-space method. The |
810 |
torque expressions have the same distance dependence as the energies, |
811 |
which are naturally longer-ranged expressions than the inter-site |
812 |
forces. Torques are also quite sensitive to orientations of |
813 |
neighboring molecules, even those that are near the cutoff distance. |
814 |
|
815 |
The results shows that the torque from the hard cutoff method |
816 |
reproduces the torques in quite good agreement with the Ewald sum. |
817 |
The other real-space methods can cause some deviations, but excellent |
818 |
agreement with the Ewald sum torques is recovered at moderate values |
819 |
of the damping coefficient ($\alpha \approx 0.2$~\AA$^{-1}$) and cutoff |
820 |
radius ($r_c \ge 12$~\AA). The TSF method exhibits only fair agreement |
821 |
in the slope when compared with the Ewald torques even for larger |
822 |
cutoff radii. It appears that the severity of the perturbations in |
823 |
the TSF method are most in evidence for the torques. |
824 |
|
825 |
\subsection{Directionality of the force and torque vectors} |
826 |
|
827 |
The accurate evaluation of force and torque directions is just as |
828 |
important for molecular dynamics simulations as the magnitudes of |
829 |
these quantities. Force and torque vectors for all six systems were |
830 |
analyzed using Fisher statistics, and the quality of the vector |
831 |
directionality is shown in terms of circular variance |
832 |
($\mathrm{Var}(\theta)$) in |
833 |
Fig. \ref{fig:slopeCorr_circularVariance}. The force and torque |
834 |
vectors from the new real-space methods exhibit nearly-ideal Fisher |
835 |
probability distributions (Eq.~\ref{eq:pdf}). Both the hard and SP |
836 |
cutoff methods exhibit the best vectorial agreement with the Ewald |
837 |
sum. The force and torque vectors from GSF method also show good |
838 |
agreement with the Ewald method, which can also be systematically |
839 |
improved by using moderate damping and a reasonable cutoff radius. For |
840 |
$\alpha = 0.2$~\AA$^{-1}$ and $r_c = 12$~\AA, we observe |
841 |
$\mathrm{Var}(\theta) = 0.00206$, which corresponds to a distribution |
842 |
with 95\% of force vectors within $6.37^\circ$ of the corresponding |
843 |
Ewald forces. The TSF method produces the poorest agreement with the |
844 |
Ewald force directions. |
845 |
|
846 |
Torques are again more perturbed than the forces by the new real-space |
847 |
methods, but even here the variance is reasonably small. For the same |
848 |
method (GSF) with the same parameters ($\alpha = 0.2$~\AA$^{-1}$, $r_c |
849 |
= 12$~\AA), the circular variance was 0.01415, corresponds to a |
850 |
distribution which has 95\% of torque vectors are within $16.75^\circ$ |
851 |
of the Ewald results. Again, the direction of the force and torque |
852 |
vectors can be systematically improved by varying $\alpha$ and $r_c$. |
853 |
|
854 |
\begin{figure} |
855 |
\centering |
856 |
\includegraphics[width=0.65\linewidth]{Variance_forceNtorque_modified.eps} |
857 |
\caption{The circular variance of the direction of the force and |
858 |
torque vectors obtained from the real-space methods around the |
859 |
reference Ewald vectors. A variance equal to 0 (dashed line) |
860 |
indicates direction of the force or torque vectors are |
861 |
indistinguishable from those obtained from the Ewald sum. Here |
862 |
different symbols represent different values of the cutoff radius |
863 |
(9~\AA\ = circle, 12~\AA\ = square, 15~\AA\ = inverted triangle)\label{fig:slopeCorr_circularVariance}} |
864 |
\end{figure} |
865 |
|
866 |
\subsection{Energy conservation\label{sec:conservation}} |
867 |
|
868 |
We have tested the conservation of energy one can expect to see with |
869 |
the new real-space methods using the soft DQ liquid model with a small |
870 |
fraction of solvated ions. This is a test system which exercises all |
871 |
orders of multipole-multipole interactions derived in the first paper |
872 |
in this series and provides the most comprehensive test of the new |
873 |
methods. A liquid-phase system was created with 2000 liquid-phase |
874 |
molecules and 48 dissolved ions at a density of 0.98 g cm$^{-3}$ and a |
875 |
temperature of 300K. After equilibration in the canonical (NVT) |
876 |
ensemble using a Nos\'e-Hoover thermostat, this liquid-phase system |
877 |
was run for 1 ns in the microcanonical (NVE) ensemble under the Ewald, |
878 |
Hard, SP, GSF, and TSF methods with a cutoff radius of 12~\AA. The |
879 |
value of the damping coefficient was also varied from the undamped |
880 |
case ($\alpha = 0$) to a heavily damped case ($\alpha = |
881 |
0.3$~\AA$^{-1}$) for all of the real space methods. A sample was also |
882 |
run using the multipolar Ewald sum with the same real-space cutoff. |
883 |
|
884 |
In figure~\ref{fig:energyDrift} we show the both the linear drift in |
885 |
energy over time, $\delta E_1$, and the standard deviation of energy |
886 |
fluctuations around this drift $\delta E_0$. Both of the |
887 |
shifted-force methods (GSF and TSF) provide excellent energy |
888 |
conservation (drift less than $10^{-5}$ kcal / mol / ns / particle), |
889 |
while the hard cutoff is essentially unusable for molecular dynamics. |
890 |
SP provides some benefit over the hard cutoff because the energetic |
891 |
jumps that happen as particles leave and enter the cutoff sphere are |
892 |
somewhat reduced, but like the Wolf method for charges, the SP method |
893 |
would not be as useful for molecular dynamics as either of the |
894 |
shifted-force methods. |
895 |
|
896 |
We note that for all tested values of the cutoff radius, the new |
897 |
real-space methods can provide better energy conservation behavior |
898 |
than the multipolar Ewald sum, even when relatively large $k$-space |
899 |
cutoff values are utilized. |
900 |
|
901 |
\begin{figure} |
902 |
\centering |
903 |
\includegraphics[width=\textwidth]{newDrift_12.eps} |
904 |
\caption{Energy conservation of the real-space methods for the soft |
905 |
DQ liauid / ion system. $\delta \mathrm{E}_1$ is the linear drift |
906 |
in energy over time (in kcal/mol/particle/ns) and $\delta |
907 |
\mathrm{E}_0$ is the standard deviation of energy fluctuations |
908 |
around this drift (in kcal/mol/particle). Points that appear in |
909 |
the green region at the bottom exhibit better energy conservation |
910 |
than would be obtained using common parameters for Ewald-based |
911 |
electrostatics.\label{fig:energyDrift}} |
912 |
\end{figure} |
913 |
|
914 |
\subsection{Reproduction of Structural \& Dynamical Features\label{sec:structure}} |
915 |
The most important test of the modified interaction potentials is the |
916 |
fidelity with which they can reproduce structural features and |
917 |
dynamical properties in a liquid. One commonly-utilized measure of |
918 |
structural ordering is the pair distribution function, $g(r)$, which |
919 |
measures local density deviations in relation to the bulk density. In |
920 |
the electrostatic approaches studied here, the short-range repulsion |
921 |
from the Lennard-Jones potential is identical for the various |
922 |
electrostatic methods, and since short range repulsion determines much |
923 |
of the local liquid ordering, one would not expect to see many |
924 |
differences in $g(r)$. Indeed, the pair distributions are essentially |
925 |
identical for all of the electrostatic methods studied (for each of |
926 |
the different systems under investigation). |
927 |
|
928 |
% An example of this agreement for the soft DQ liquid/ion system is |
929 |
% shown in Fig. \ref{fig:gofr}. |
930 |
|
931 |
% \begin{figure} |
932 |
% \centering |
933 |
% \includegraphics[width=\textwidth]{gofr_ssdqc.eps} |
934 |
% \caption{The pair distribution functions, $g(r)$, for the SSDQ |
935 |
% water/ion system obtained using the different real-space methods are |
936 |
% essentially identical with the result from the Ewald |
937 |
% treatment.\label{fig:gofr}} |
938 |
% \end{figure} |
939 |
|
940 |
There is a minor over-structuring of the first solvation shell when |
941 |
using TSF or when overdamping with any of the real-space methods. |
942 |
With moderate damping, GSF and SP produce pair distributions that are |
943 |
identical (within numerical noise) to their Ewald counterparts. The |
944 |
degree of over-structuring can be measured most easily using the |
945 |
coordination number, |
946 |
\begin{equation} |
947 |
n_C = 4\pi\rho \int_{0}^{a}r^2\text{g}(r)dr, |
948 |
\end{equation} |
949 |
where $\rho$ is the number density of the site-site pair interactions, |
950 |
and $a$ is the radial location of the minima following the first peak |
951 |
in $g(r)$ ($a = 4.2$~\AA\ for the soft DQ liquid / ion system). The |
952 |
coordination number is shown as a function of the damping coefficient |
953 |
for all of the real space methods in Fig. \ref{fig:Props}. |
954 |
|
955 |
A more demanding test of modified electrostatics is the average value |
956 |
of the electrostatic energy $\langle U_\mathrm{elect} \rangle / N$ |
957 |
which is obtained by sampling the liquid-state configurations |
958 |
experienced by a liquid evolving entirely under the influence of each |
959 |
of the methods. In Fig. \ref{fig:Props} we demonstrate how $\langle |
960 |
U_\mathrm{elect} \rangle / N$ varies with the damping parameter, |
961 |
$\alpha$, for each of the methods. |
962 |
|
963 |
As in the crystals studied in the first paper, damping is important |
964 |
for converging the mean electrostatic energy values, particularly for |
965 |
the two shifted force methods (GSF and TSF). A value of $\alpha |
966 |
\approx 0.2$~\AA$^{-1}$ is sufficient to converge the SP and GSF |
967 |
energies with a cutoff of 12 \AA, while shorter cutoffs require more |
968 |
dramatic damping ($\alpha \approx 0.28$~\AA$^{-1}$ for $r_c = 9$~\AA). |
969 |
Overdamping the real-space electrostatic methods occurs with $\alpha > |
970 |
0.3$~\AA$^{-1}$, causing the estimate of the electrostatic energy to |
971 |
drop below the Ewald results. |
972 |
|
973 |
These ``optimal'' values of the damping coefficient for structural |
974 |
features are similar to those observed for DSF electrostatics for |
975 |
purely point-charge systems, and the range $\alpha= 0.175 \rightarrow |
976 |
0.225$~\AA$^{-1}$ for $r_c = 12$~\AA\ appears to be an excellent |
977 |
compromise for mixed charge/multipolar systems. |
978 |
|
979 |
To test the fidelity of the electrostatic methods at reproducing |
980 |
\textit{dynamics} in a multipolar liquid, it is also useful to look at |
981 |
transport properties, particularly the diffusion constant, |
982 |
\begin{equation} |
983 |
D = \lim_{t \rightarrow \infty} \frac{1}{6 t} \langle \left| |
984 |
\mathbf{r}(t) -\mathbf{r}(0) \right|^2 \rangle |
985 |
\label{eq:diff} |
986 |
\end{equation} |
987 |
which measures long-time behavior and is sensitive to the forces on |
988 |
the multipoles. The self-diffusion constants (D) were calculated from |
989 |
linear fits to the long-time portion of the mean square displacement, |
990 |
$\langle r^{2}(t) \rangle$.\cite{Allen87} In Fig. \ref{fig:Props} we |
991 |
demonstrate how the diffusion constant depends on the choice of |
992 |
real-space methods and the damping coefficient. Both the SP and GSF |
993 |
methods can obtain excellent agreement with Ewald again using moderate |
994 |
damping. |
995 |
|
996 |
In addition to translational diffusion, orientational relaxation times |
997 |
were calculated for comparisons with the Ewald simulations and with |
998 |
experiments. These values were determined by calculating the |
999 |
orientational time correlation function, |
1000 |
\begin{equation} |
1001 |
C_l^\gamma(t) = \left\langle P_l\left[\hat{\mathbf{A}}_\gamma(t) |
1002 |
\cdot\hat{\mathbf{A}}_\gamma(0)\right]\right\rangle, |
1003 |
\label{eq:OrientCorr} |
1004 |
\end{equation} |
1005 |
from the same 350 ps microcanonical trajectories that were used for |
1006 |
translational diffusion. Here, $P_l$ is the Legendre polynomial of |
1007 |
order $l$ and $\hat{\mathbf{A}}_\gamma$ is the unit vector for body |
1008 |
axis $\gamma$. The reference frame used for our sample dipolar |
1009 |
systems has the $z$-axis running along the dipoles, and for the soft |
1010 |
DQ liquid model, the $y$-axis connects the two implied hydrogen-like |
1011 |
positions. From the orientation autocorrelation functions, we can |
1012 |
obtain time constants for rotational relaxation either by fitting to a |
1013 |
multi-exponential model for the orientational relaxation, or by |
1014 |
integrating the correlation functions. |
1015 |
|
1016 |
In a good model for water, the orientational decay times would be |
1017 |
comparable to water orientational relaxation times from nuclear |
1018 |
magnetic resonance (NMR). The relaxation constant obtained from |
1019 |
$C_2^y(t)$ is normally of experimental interest because it describes |
1020 |
the relaxation of the principle axis connecting the hydrogen |
1021 |
atoms. Thus, $C_2^y(t)$ can be compared to the intermolecular portion |
1022 |
of the dipole-dipole relaxation from a proton NMR signal and can |
1023 |
provide an estimate of the NMR relaxation time constant.\cite{Impey82} |
1024 |
In Fig. \ref{fig:Props} we compare the $\tau_2^y$ and $\tau_2^z$ |
1025 |
values for the various real-space methods over a range of different |
1026 |
damping coefficients. The rotational relaxation for the $z$ axis |
1027 |
primarily probes the torques on the dipoles, while the relaxation for |
1028 |
the $y$ axis is sensitive primarily to the quadrupolar torques. |
1029 |
|
1030 |
\begin{figure} |
1031 |
\includegraphics[width=\textwidth]{properties.eps} |
1032 |
\caption{Comparison of the structural and dynamic properties for the |
1033 |
combined multipolar liquid (soft DQ liquid + ions) for all of the |
1034 |
real-space methods with $r_c = 12$~\AA. Electrostatic energies, |
1035 |
$\langle U_\mathrm{elect} \rangle / N$ (in kcal / mol), |
1036 |
coordination numbers, $n_C$, diffusion constants (in $10^{-5} |
1037 |
\mathrm{cm}^2\mathrm{s}^{-1}$), and rotational correlation times |
1038 |
(in ps) all show excellent agreement with Ewald results for |
1039 |
damping coefficients in the range $\alpha= 0.175 \rightarrow |
1040 |
0.225$~\AA$^{-1}$. \label{fig:Props}} |
1041 |
\end{figure} |
1042 |
|
1043 |
In Fig. \ref{fig:Props} it appears that values for $D$, $\tau_2^y$, |
1044 |
and $\tau_2^z$ using the Ewald sum are reproduced with excellent |
1045 |
fidelity by the GSF and SP methods. All of the real space methods can |
1046 |
be \textit{overdamped}, which reduces the effective range of multipole |
1047 |
interactions, causing structural and dynamical changes from the |
1048 |
correct behavior. Because overdamping weakens orientational |
1049 |
preferences between adjacent molecules, it manifests as too-rapid |
1050 |
orientational decay coupled with faster diffusion and |
1051 |
over-coordination of the liquid. Underdamping is less problematic for |
1052 |
the SP and GSF methods, as their structural and dynamical properties |
1053 |
still reproduce the Ewald results even in the completely undamped |
1054 |
($\alpha = 0$) case. An optimal range for the electrostatic damping |
1055 |
parameter appears to be $\alpha= 0.175 \rightarrow 0.225$~\AA$^{-1}$ |
1056 |
for $r_c = 12$~\AA, which similar to the optimal range found for the |
1057 |
damped shifted force potential for point charges.\cite{Fennell:2006lq} |
1058 |
|
1059 |
\section{CONCLUSION} |
1060 |
In the first paper in this series, we generalized the |
1061 |
charge-neutralized electrostatic energy originally developed by Wolf |
1062 |
\textit{et al.}\cite{Wolf:1999dn} to multipole-multipole interactions |
1063 |
up to quadrupolar order. The SP method is essentially a |
1064 |
multipole-capable version of the Wolf model. The SP method for |
1065 |
multipoles provides excellent agreement with Ewald-derived energies, |
1066 |
forces and torques, and is suitable for Monte Carlo simulations, |
1067 |
although the forces and torques retain discontinuities at the cutoff |
1068 |
distance that prevents its use in molecular dynamics. |
1069 |
|
1070 |
We also developed two natural extensions of the damped shifted-force |
1071 |
(DSF) model originally proposed by Zahn {\it et al.} and extended by |
1072 |
Fennell and Gezelter.\cite{Zahn:2002hc,Fennell:2006lq} The GSF and TSF |
1073 |
approaches provide smooth truncation of energies, forces, and torques |
1074 |
at the real-space cutoff, and both converge to DSF electrostatics for |
1075 |
point-charge interactions. The TSF model is based on a high-order |
1076 |
truncated Taylor expansion which can be relatively perturbative inside |
1077 |
the cutoff sphere. The GSF model takes the gradient from an images of |
1078 |
the interacting multipole that has been projected onto the cutoff |
1079 |
sphere to derive shifted force and torque expressions, and is a |
1080 |
significantly more gentle approach. |
1081 |
|
1082 |
The GSF method produces quantitative agreement with Ewald energies, |
1083 |
forces, and torques. It also performs well in conserving energy in MD |
1084 |
simulations. The Taylor-shifted (TSF) model provides smooth dynamics, |
1085 |
but these take place on a potential energy surface that is |
1086 |
significantly perturbed from Ewald-based electrostatics. Because it |
1087 |
performs relatively poorly compared with GSF, it may seem odd that |
1088 |
that the TSF model was included in this work. However, the functional |
1089 |
forms derived for the SP and GSF methods depend on the separation of |
1090 |
orientational contributions that were made visible by the Taylor |
1091 |
series of the electrostatic kernel at the cutoff radius. The TSF |
1092 |
method also has the unique property that a large number of derivatives |
1093 |
can be made to vanish at the cutoff radius. This property has proven |
1094 |
useful in past treatments of the corrections to the Clausius-Mossotti |
1095 |
fluctuation formula for dielectric constants.\cite{Izvekov:2008wo} |
1096 |
|
1097 |
Reproduction of both structural and dynamical features in the liquid |
1098 |
systems is remarkably good for both the SP and GSF models. Pair |
1099 |
distribution functions are essentially equivalent to the same |
1100 |
functions produced using Ewald-based electrostatics, and with moderate |
1101 |
damping, a structural feature that directly probes the electrostatic |
1102 |
interaction (e.g. the mean electrostatic potential energy) can also be |
1103 |
made quantitative. Dynamical features are sensitive probes of the |
1104 |
forces and torques produced by these methods, and even though the |
1105 |
smooth behavior of forces is produced by perturbing the overall |
1106 |
potential, the diffusion constants and orientational correlation times |
1107 |
are quite close to the Ewald-based results. |
1108 |
|
1109 |
The only cases we have found where the new GSF and SP real-space |
1110 |
methods can be problematic are those which retain a bulk dipole moment |
1111 |
at large distances (e.g. the $Z_1$ dipolar lattice). In ferroelectric |
1112 |
materials, uniform weighting of the orientational contributions can be |
1113 |
important for converging the total energy. In these cases, the |
1114 |
damping function which causes the non-uniform weighting can be |
1115 |
replaced by the bare electrostatic kernel, and the energies return to |
1116 |
the expected converged values. |
1117 |
|
1118 |
Based on the results of this work, we can conclude that the GSF method |
1119 |
is a suitable and efficient replacement for the Ewald sum for |
1120 |
evaluating electrostatic interactions in modern MD simulations, and |
1121 |
the SP method would be an excellent choice for Monte Carlo |
1122 |
simulations where smooth forces and energy conservation are not |
1123 |
important. Both the SP and GSF methods retain excellent fidelity to |
1124 |
the Ewald energies, forces and torques. Additionally, the energy |
1125 |
drift and fluctuations from the GSF electrostatics are significantly |
1126 |
better than a multipolar Ewald sum for finite-sized reciprocal spaces, |
1127 |
and physical properties are reproduced accurately. |
1128 |
|
1129 |
As in all purely pairwise cutoff methods, the SP, GSF and TSF methods |
1130 |
are expected to scale approximately {\it linearly} with system size, |
1131 |
and are easily parallelizable. This should result in substantial |
1132 |
reductions in the computational cost of performing large simulations. |
1133 |
With the proper use of pre-computation and spline interpolation of the |
1134 |
radial functions, the real-space methods are essentially the same cost |
1135 |
as a simple real-space cutoff. They require no Fourier transforms or |
1136 |
$k$-space sums, and guarantee the smooth handling of energies, forces, |
1137 |
and torques as multipoles cross the real-space cutoff boundary. |
1138 |
|
1139 |
We are not suggesting that there is any flaw with the Ewald sum; in |
1140 |
fact, it is the standard by which the SP, GSF, and TSF methods have |
1141 |
been judged in this work. However, these results provide evidence |
1142 |
that in the typical simulations performed today, the Ewald summation |
1143 |
may no longer be required to obtain the level of accuracy most |
1144 |
researchers have come to expect. |
1145 |
|
1146 |
\begin{acknowledgments} |
1147 |
JDG acknowledges helpful discussions with Christopher |
1148 |
Fennell. Support for this project was provided by the National |
1149 |
Science Foundation under grant CHE-1362211. Computational time was |
1150 |
provided by the Center for Research Computing (CRC) at the |
1151 |
University of Notre Dame. |
1152 |
\end{acknowledgments} |
1153 |
|
1154 |
%\bibliographystyle{aip} |
1155 |
\newpage |
1156 |
\bibliography{references} |
1157 |
\end{document} |
1158 |
|
1159 |
% |
1160 |
% ****** End of file aipsamp.tex ****** |