47 |
|
|
48 |
|
%\preprint{AIP/123-QED} |
49 |
|
|
50 |
< |
\title{Real space alternatives to the Ewald |
51 |
< |
Sum. II. Comparison of Methods} % Force line breaks with \\ |
50 |
> |
\title{Real space alternatives to the Ewald Sum. II. Comparison of Methods} |
51 |
|
|
52 |
|
\author{Madan Lamichhane} |
53 |
< |
\affiliation{Department of Physics, University |
55 |
< |
of Notre Dame, Notre Dame, IN 46556}%Lines break automatically or can be forced with \\ |
53 |
> |
\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} |
54 |
|
|
55 |
|
\author{Kathie E. Newman} |
56 |
< |
\affiliation{Department of Physics, University |
59 |
< |
of Notre Dame, Notre Dame, IN 46556} |
56 |
> |
\affiliation{Department of Physics, University of Notre Dame, Notre Dame, IN 46556} |
57 |
|
|
58 |
|
\author{J. Daniel Gezelter}% |
59 |
|
\email{gezelter@nd.edu.} |
60 |
< |
\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556%\\This line break forced with \textbackslash\textbackslash |
61 |
< |
}% |
60 |
> |
\affiliation{Department of Chemistry and Biochemistry, University of Notre Dame, Notre Dame, IN 46556 |
61 |
> |
} |
62 |
|
|
63 |
< |
\date{\today}% It is always \today, today, |
67 |
< |
% but any date may be explicitly specified |
63 |
> |
\date{\today} |
64 |
|
|
65 |
|
\begin{abstract} |
66 |
< |
We have tested the real-space shifted potential (SP), |
67 |
< |
gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods |
68 |
< |
for multipole interactions that were developed in the first paper in |
69 |
< |
this series, using the multipolar Ewald sum as a reference |
70 |
< |
method. The tests were carried out in a variety of condensed-phase |
71 |
< |
environments which were designed to test all levels of the |
72 |
< |
multipole-multipole interactions. Comparisons of the energy |
73 |
< |
differences between configurations, molecular forces, and torques |
74 |
< |
were used to analyze how well the real-space models perform relative |
75 |
< |
to the more computationally expensive Ewald treatment. We have also |
76 |
< |
investigated the energy conservation properties of the new methods |
77 |
< |
in molecular dynamics simulations. The SP method shows excellent |
78 |
< |
agreement with configurational energy differences, forces, and |
79 |
< |
torques, and would be suitable for use in Monte Carlo calculations. |
80 |
< |
Of the two new shifted-force methods, the GSF approach shows the |
81 |
< |
best agreement with Ewald-derived energies, forces, and torques and |
82 |
< |
exhibits energy conservation properties that make it an excellent |
83 |
< |
choice for efficient computation of electrostatic interactions in |
88 |
< |
molecular dynamics simulations. |
66 |
> |
We report on tests of the shifted potential (SP), gradient shifted |
67 |
> |
force (GSF), and Taylor shifted force (TSF) real-space methods for |
68 |
> |
multipole interactions developed in the first paper in this series, |
69 |
> |
using the multipolar Ewald sum as a reference method. The tests were |
70 |
> |
carried out in a variety of condensed-phase environments designed to |
71 |
> |
test up to quadrupole-quadrupole interactions. Comparisons of the |
72 |
> |
energy differences between configurations, molecular forces, and |
73 |
> |
torques were used to analyze how well the real-space models perform |
74 |
> |
relative to the more computationally expensive Ewald treatment. We |
75 |
> |
have also investigated the energy conservation properties of the new |
76 |
> |
methods in molecular dynamics simulations. The SP method shows |
77 |
> |
excellent agreement with configurational energy differences, forces, |
78 |
> |
and torques, and would be suitable for use in Monte Carlo |
79 |
> |
calculations. Of the two new shifted-force methods, the GSF |
80 |
> |
approach shows the best agreement with Ewald-derived energies, |
81 |
> |
forces, and torques and also exhibits energy conservation properties |
82 |
> |
that make it an excellent choice for efficient computation of |
83 |
> |
electrostatic interactions in molecular dynamics simulations. |
84 |
|
\end{abstract} |
85 |
|
|
86 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
89 |
|
|
90 |
|
\maketitle |
91 |
|
|
97 |
– |
|
92 |
|
\section{\label{sec:intro}Introduction} |
93 |
|
Computing the interactions between electrostatic sites is one of the |
94 |
|
most expensive aspects of molecular simulations. There have been |
99 |
|
the conditionally convergent electrostatic energy is converted into |
100 |
|
two absolutely convergent contributions, one which is carried out in |
101 |
|
real space with a cutoff radius, and one in reciprocal |
102 |
< |
space. BETTER CITATIONS\cite{Clarke:1986eu,Woodcock75} |
102 |
> |
space.\cite{Ewald21,deLeeuw80,Smith81,Allen87} |
103 |
|
|
104 |
|
When carried out as originally formulated, the reciprocal-space |
105 |
|
portion of the Ewald sum exhibits relatively poor computational |
106 |
< |
scaling, making it prohibitive for large systems. By utilizing |
107 |
< |
particle meshes and three dimensional fast Fourier transforms (FFT), |
108 |
< |
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
109 |
< |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease |
110 |
< |
the computational cost from $O(N^2)$ down to $O(N \log |
111 |
< |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
106 |
> |
scaling, making it prohibitive for large systems. By utilizing a |
107 |
> |
particle mesh and three dimensional fast Fourier transforms (FFT), the |
108 |
> |
particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
109 |
> |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) |
110 |
> |
methods can decrease the computational cost from $O(N^2)$ down to $O(N |
111 |
> |
\log |
112 |
> |
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb} |
113 |
|
|
114 |
|
Because of the artificial periodicity required for the Ewald sum, |
115 |
|
interfacial molecular systems such as membranes and liquid-vapor |
116 |
< |
interfaces require modifications to the |
117 |
< |
method.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
118 |
< |
Parry's extension of the three dimensional Ewald sum is appropriate |
119 |
< |
for slab geometries.\cite{Parry:1975if} Modified Ewald methods that |
120 |
< |
were developed to handle two-dimensional (2D) electrostatic |
121 |
< |
interactions in interfacial systems have not seen similar |
122 |
< |
particle-mesh treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
123 |
< |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly |
124 |
< |
with system size. The inherent periodicity in the Ewald’s method can |
125 |
< |
also be problematic for interfacial molecular |
126 |
< |
systems.\cite{Fennell:2006lq} |
116 |
> |
interfaces require modifications to the method. Parry's extension of |
117 |
> |
the three dimensional Ewald sum is appropriate for slab |
118 |
> |
geometries.\cite{Parry:1975if} Modified Ewald methods that were |
119 |
> |
developed to handle two-dimensional (2-D) electrostatic |
120 |
> |
interactions,\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
121 |
> |
but these methods were originally quite computationally |
122 |
> |
expensive.\cite{Spohr97,Yeh99} There have been several successful |
123 |
> |
efforts that reduced the computational cost of 2-D lattice |
124 |
> |
summations,\cite{Yeh99,Kawata01,Arnold02,deJoannis02,Brodka04} |
125 |
> |
bringing them more in line with the scaling for the full 3-D |
126 |
> |
treatments. The inherent periodicity in the Ewald method can also |
127 |
> |
be problematic for interfacial molecular |
128 |
> |
systems.\cite{Roberts94,Roberts95,Luty96,Hunenberger99a,Hunenberger99b,Weber00,Fennell:2006lq} |
129 |
|
|
130 |
|
\subsection{Real-space methods} |
131 |
|
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
147 |
|
what is effectively a set of octupoles at large distances. These facts |
148 |
|
suggest that the Madelung constants are relatively short ranged for |
149 |
|
perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful |
150 |
< |
application of Wolf's method are able to obtain accurate estimates of |
150 |
> |
application of Wolf's method can provide accurate estimates of |
151 |
|
Madelung constants using relatively short cutoff radii. |
152 |
|
|
153 |
|
Direct truncation of interactions at a cutoff radius creates numerical |
154 |
< |
errors. Wolf \textit{et al.} argued that truncation errors are due |
154 |
> |
errors. Wolf \textit{et al.} suggest that truncation errors are due |
155 |
|
to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To |
156 |
|
neutralize this charge they proposed placing an image charge on the |
157 |
|
surface of the cutoff sphere for every real charge inside the cutoff. |
158 |
|
These charges are present for the evaluation of both the pair |
159 |
|
interaction energy and the force, although the force expression |
160 |
< |
maintained a discontinuity at the cutoff sphere. In the original Wolf |
160 |
> |
maintains a discontinuity at the cutoff sphere. In the original Wolf |
161 |
|
formulation, the total energy for the charge and image were not equal |
162 |
< |
to the integral of their force expression, and as a result, the total |
162 |
> |
to the integral of the force expression, and as a result, the total |
163 |
|
energy would not be conserved in molecular dynamics (MD) |
164 |
|
simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and |
165 |
|
Gezelter later proposed shifted force variants of the Wolf method with |
166 |
|
commensurate force and energy expressions that do not exhibit this |
167 |
< |
problem.\cite{Fennell:2006lq} Related real-space methods were also |
168 |
< |
proposed by Chen \textit{et |
167 |
> |
problem.\cite{Zahn:2002hc,Fennell:2006lq} Related real-space methods |
168 |
> |
were also proposed by Chen \textit{et |
169 |
|
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
170 |
< |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has used |
171 |
< |
neutralization of the higher order moments for the calculation of the |
172 |
< |
electrostatic interaction of the point charge |
176 |
< |
systems.\cite{Fukuda:2013sf} |
170 |
> |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfuly |
171 |
> |
used additional neutralization of higher order moments for systems of |
172 |
> |
point charges.\cite{Fukuda:2013sf} |
173 |
|
|
174 |
|
\begin{figure} |
175 |
|
\centering |
176 |
< |
\includegraphics[width=\linewidth]{schematic.pdf} |
176 |
> |
\includegraphics[width=\linewidth]{schematic.eps} |
177 |
|
\caption{Top: Ionic systems exhibit local clustering of dissimilar |
178 |
|
charges (in the smaller grey circle), so interactions are |
179 |
|
effectively charge-multipole at longer distances. With hard |
191 |
|
One can make a similar effective range argument for crystals of point |
192 |
|
\textit{multipoles}. The Luttinger and Tisza treatment of energy |
193 |
|
constants for dipolar lattices utilizes 24 basis vectors that contain |
194 |
< |
dipoles at the eight corners of a unit cube. Only three of these |
195 |
< |
basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
194 |
> |
dipoles at the eight corners of a unit cube.\cite{LT} Only three of |
195 |
> |
these basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
196 |
|
moments, while the rest have zero net dipole and retain contributions |
197 |
< |
only from higher order multipoles. The lowest energy crystalline |
197 |
> |
only from higher order multipoles. The lowest-energy crystalline |
198 |
|
structures are built out of basis vectors that have only residual |
199 |
|
quadrupolar moments (e.g. the $Z_5$ array). In these low energy |
200 |
|
structures, the effective interaction between a dipole at the center |
217 |
|
|
218 |
|
The shorter effective range of electrostatic interactions is not |
219 |
|
limited to perfect crystals, but can also apply in disordered fluids. |
220 |
< |
Even at elevated temperatures, there is, on average, local charge |
221 |
< |
balance in an ionic liquid, where each positive ion has surroundings |
222 |
< |
dominated by negaitve ions and vice versa. The reversed-charge images |
223 |
< |
on the cutoff sphere that are integral to the Wolf and DSF approaches |
224 |
< |
retain the effective multipolar interactions as the charges traverse |
225 |
< |
the cutoff boundary. |
220 |
> |
Even at elevated temperatures, there is local charge balance in an |
221 |
> |
ionic liquid, where each positive ion has surroundings dominated by |
222 |
> |
negaitve ions and vice versa. The reversed-charge images on the |
223 |
> |
cutoff sphere that are integral to the Wolf and DSF approaches retain |
224 |
> |
the effective multipolar interactions as the charges traverse the |
225 |
> |
cutoff boundary. |
226 |
|
|
227 |
|
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
228 |
|
significant orientational averaging that additionally reduces the |
241 |
|
% to the non-neutralized value of the higher order moments within the |
242 |
|
% cutoff sphere. |
243 |
|
|
244 |
< |
The forces and torques acting on atomic sites are the fundamental |
245 |
< |
factors driving dynamics in molecular simulations. Fennell and |
246 |
< |
Gezelter proposed the damped shifted force (DSF) energy kernel to |
247 |
< |
obtain consistent energies and forces on the atoms within the cutoff |
248 |
< |
sphere. Both the energy and the force go smoothly to zero as an atom |
249 |
< |
aproaches the cutoff radius. The comparisons of the accuracy these |
250 |
< |
quantities between the DSF kernel and SPME was surprisingly |
251 |
< |
good.\cite{Fennell:2006lq} The DSF method has seen increasing use for |
252 |
< |
calculating electrostatic interactions in molecular systems with |
253 |
< |
relatively uniform charge |
258 |
< |
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
244 |
> |
Forces and torques acting on atomic sites are fundamental in driving |
245 |
> |
dynamics in molecular simulations, and the damped shifted force (DSF) |
246 |
> |
energy kernel provides consistent energies and forces on charged atoms |
247 |
> |
within the cutoff sphere. Both the energy and the force go smoothly to |
248 |
> |
zero as an atom aproaches the cutoff radius. The comparisons of the |
249 |
> |
accuracy these quantities between the DSF kernel and SPME was |
250 |
> |
surprisingly good.\cite{Fennell:2006lq} As a result, the DSF method |
251 |
> |
has seen increasing use in molecular systems with relatively uniform |
252 |
> |
charge |
253 |
> |
densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,Vergne13} |
254 |
|
|
255 |
|
\subsection{The damping function} |
256 |
|
The damping function has been discussed in detail in the first paper |
277 |
|
required to compute configurational |
278 |
|
energies.\cite{Ren06,Essex10,Essex11} |
279 |
|
|
280 |
< |
Because electrons in a molecule are not localized at specific points, |
281 |
< |
the assignment of partial charges to atomic centers is always an |
282 |
< |
approximation. Atomic sites can also be assigned point multipoles and |
283 |
< |
polarizabilities to increase the accuracy of the molecular model. |
284 |
< |
Recently, water has been modeled with point multipoles up to octupolar |
285 |
< |
order using the soft sticky dipole-quadrupole-octupole (SSDQO) |
280 |
> |
Additionally, because electrons in a molecule are not localized at |
281 |
> |
specific points, the assignment of partial charges to atomic centers |
282 |
> |
is always an approximation. For increased accuracy, atomic sites can |
283 |
> |
also be assigned point multipoles and polarizabilities. Recently, |
284 |
> |
water has been modeled with point multipoles up to octupolar order |
285 |
> |
using the soft sticky dipole-quadrupole-octupole (SSDQO) |
286 |
|
model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
287 |
|
multipoles up to quadrupolar order have also been coupled with point |
288 |
|
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
289 |
|
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However, |
290 |
|
truncating point multipoles without smoothing the forces and torques |
291 |
< |
will create energy conservation issues in molecular dynamics simulations. |
291 |
> |
can create energy conservation issues in molecular dynamics |
292 |
> |
simulations. |
293 |
|
|
294 |
|
In this paper we test a set of real-space methods that were developed |
295 |
|
for point multipolar interactions. These methods extend the damped |
296 |
|
shifted force (DSF) and Wolf methods originally developed for |
297 |
|
charge-charge interactions and generalize them for higher order |
298 |
< |
multipoles. The detailed mathematical development of these methods has |
299 |
< |
been presented in the first paper in this series, while this work |
300 |
< |
covers the testing the energies, forces, torques, and energy |
298 |
> |
multipoles. The detailed mathematical development of these methods |
299 |
> |
has been presented in the first paper in this series, while this work |
300 |
> |
covers the testing of energies, forces, torques, and energy |
301 |
|
conservation properties of the methods in realistic simulation |
302 |
|
environments. In all cases, the methods are compared with the |
303 |
< |
reference method, a full multipolar Ewald treatment. |
303 |
> |
reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98} |
304 |
|
|
305 |
|
|
306 |
|
%\subsection{Conservation of total energy } |
585 |
|
\subsection{Implementation} |
586 |
|
The real-space methods developed in the first paper in this series |
587 |
|
have been implemented in our group's open source molecular simulation |
588 |
< |
program, OpenMD,\cite{openmd} which was used for all calculations in |
588 |
> |
program, OpenMD,\cite{Meineke05,openmd} which was used for all calculations in |
589 |
|
this work. The complementary error function can be a relatively slow |
590 |
|
function on some processors, so all of the radial functions are |
591 |
|
precomputed on a fine grid and are spline-interpolated to provide |
984 |
|
|
985 |
|
\begin{figure} |
986 |
|
\centering |
987 |
< |
\includegraphics[width=\textwidth]{newDrift_12.pdf} |
987 |
> |
\includegraphics[width=\textwidth]{newDrift_12.eps} |
988 |
|
\label{fig:energyDrift} |
989 |
|
\caption{Analysis of the energy conservation of the real-space |
990 |
|
electrostatic methods. $\delta \mathrm{E}_1$ is the linear drift in |