69 |
|
\begin{abstract} |
70 |
|
We have tested the real-space shifted potential (SP), |
71 |
|
gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods |
72 |
< |
for multipoles that were developed in the first paper in this series |
73 |
< |
against the multipolar Ewald sum as a reference method. The tests |
74 |
< |
were carried out in a variety of condensed-phase environments which |
75 |
< |
were designed to test all levels of the multipole-multipole |
76 |
< |
interactions. Comparisons of the energy differences between |
77 |
< |
configurations, molecular forces, and torques were used to analyze |
78 |
< |
how well the real-space models perform relative to the more |
79 |
< |
computationally expensive Ewald treatment. We have also investigated the |
80 |
< |
energy conservation properties of the new methods in molecular |
81 |
< |
dynamics simulations using all of these methods. The SP method shows |
82 |
< |
excellent agreement with configurational energy differences, forces, |
83 |
< |
and torques, and would be suitable for use in Monte Carlo |
84 |
< |
calculations. Of the two new shifted-force methods, the GSF |
85 |
< |
approach shows the best agreement with Ewald-derived energies, |
86 |
< |
forces, and torques and exhibits energy conservation properties that |
87 |
< |
make it an excellent choice for efficiently computing electrostatic |
88 |
< |
interactions in molecular dynamics simulations. |
72 |
> |
for multipole interactions that were developed in the first paper in |
73 |
> |
this series, using the multipolar Ewald sum as a reference |
74 |
> |
method. The tests were carried out in a variety of condensed-phase |
75 |
> |
environments which were designed to test all levels of the |
76 |
> |
multipole-multipole interactions. Comparisons of the energy |
77 |
> |
differences between configurations, molecular forces, and torques |
78 |
> |
were used to analyze how well the real-space models perform relative |
79 |
> |
to the more computationally expensive Ewald treatment. We have also |
80 |
> |
investigated the energy conservation properties of the new methods |
81 |
> |
in molecular dynamics simulations. The SP method shows excellent |
82 |
> |
agreement with configurational energy differences, forces, and |
83 |
> |
torques, and would be suitable for use in Monte Carlo calculations. |
84 |
> |
Of the two new shifted-force methods, the GSF approach shows the |
85 |
> |
best agreement with Ewald-derived energies, forces, and torques and |
86 |
> |
exhibits energy conservation properties that make it an excellent |
87 |
> |
choice for efficient computation of electrostatic interactions in |
88 |
> |
molecular dynamics simulations. |
89 |
|
\end{abstract} |
90 |
|
|
91 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
97 |
|
|
98 |
|
\section{\label{sec:intro}Introduction} |
99 |
|
Computing the interactions between electrostatic sites is one of the |
100 |
< |
most expensive aspects of molecular simulations, which is why there |
101 |
< |
have been significant efforts to develop practical, efficient and |
102 |
< |
convergent methods for handling these interactions. Ewald's method is |
103 |
< |
perhaps the best known and most accurate method for evaluating |
104 |
< |
energies, forces, and torques in explicitly-periodic simulation |
105 |
< |
cells. In this approach, the conditionally convergent electrostatic |
106 |
< |
energy is converted into two absolutely convergent contributions, one |
107 |
< |
which is carried out in real space with a cutoff radius, and one in |
108 |
< |
reciprocal space.\cite{Clarke:1986eu,Woodcock75} |
100 |
> |
most expensive aspects of molecular simulations. There have been |
101 |
> |
significant efforts to develop practical, efficient and convergent |
102 |
> |
methods for handling these interactions. Ewald's method is perhaps the |
103 |
> |
best known and most accurate method for evaluating energies, forces, |
104 |
> |
and torques in explicitly-periodic simulation cells. In this approach, |
105 |
> |
the conditionally convergent electrostatic energy is converted into |
106 |
> |
two absolutely convergent contributions, one which is carried out in |
107 |
> |
real space with a cutoff radius, and one in reciprocal |
108 |
> |
space. BETTER CITATIONS\cite{Clarke:1986eu,Woodcock75} |
109 |
|
|
110 |
|
When carried out as originally formulated, the reciprocal-space |
111 |
|
portion of the Ewald sum exhibits relatively poor computational |
116 |
|
the computational cost from $O(N^2)$ down to $O(N \log |
117 |
|
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
118 |
|
|
119 |
< |
Because of the artificial periodicity required for the Ewald sum, the |
120 |
< |
method may require modification to compute interactions for |
119 |
> |
Because of the artificial periodicity required for the Ewald sum, |
120 |
|
interfacial molecular systems such as membranes and liquid-vapor |
121 |
< |
interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
122 |
< |
To simulate interfacial systems, Parry's extension of the 3D Ewald sum |
123 |
< |
is appropriate for slab geometries.\cite{Parry:1975if} The inherent |
124 |
< |
periodicity in the Ewald method can also be problematic for molecular |
125 |
< |
interfaces.\cite{Fennell:2006lq} Modified Ewald methods that were |
126 |
< |
developed to handle two-dimensional (2D) electrostatic interactions in |
127 |
< |
interfacial systems have not seen similar particle-mesh |
129 |
< |
treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
121 |
> |
interfaces require modifications to the |
122 |
> |
method.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
123 |
> |
Parry's extension of the three dimensional Ewald sum is appropriate |
124 |
> |
for slab geometries.\cite{Parry:1975if} Modified Ewald methods that |
125 |
> |
were developed to handle two-dimensional (2D) electrostatic |
126 |
> |
interactions in interfacial systems have not seen similar |
127 |
> |
particle-mesh treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
128 |
|
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly |
129 |
< |
with system size. |
129 |
> |
with system size. The inherent periodicity in the Ewald’s method can |
130 |
> |
also be problematic for interfacial molecular |
131 |
> |
systems.\cite{Fennell:2006lq} |
132 |
|
|
133 |
|
\subsection{Real-space methods} |
134 |
|
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
135 |
|
method for calculating electrostatic interactions between point |
136 |
< |
charges. They argued that the effective Coulomb interaction in |
137 |
< |
condensed phase systems is actually short ranged.\cite{Wolf92,Wolf95} |
138 |
< |
For an ordered lattice (e.g., when computing the Madelung constant of |
139 |
< |
an ionic solid), the material can be considered as a set of ions |
140 |
< |
interacting with neutral dipolar or quadrupolar ``molecules'' giving |
141 |
< |
an effective distance dependence for the electrostatic interactions of |
142 |
< |
$r^{-5}$ (see figure \ref{fig:schematic}). For this reason, careful |
143 |
< |
application of Wolf's method can obtain accurate estimates of Madelung |
144 |
< |
constants using relatively short cutoff radii. Recently, Fukuda used |
145 |
< |
neutralization of the higher order moments for calculation of the |
146 |
< |
electrostatic interactions in point charge |
136 |
> |
charges. They argued that the effective Coulomb interaction in most |
137 |
> |
condensed phase systems is effectively short |
138 |
> |
ranged.\cite{Wolf92,Wolf95} For an ordered lattice (e.g., when |
139 |
> |
computing the Madelung constant of an ionic solid), the material can |
140 |
> |
be considered as a set of ions interacting with neutral dipolar or |
141 |
> |
quadrupolar ``molecules'' giving an effective distance dependence for |
142 |
> |
the electrostatic interactions of $r^{-5}$ (see figure |
143 |
> |
\ref{fig:schematic}). If one views the \ce{NaCl} crystal as a simple |
144 |
> |
cubic (SC) structure with an octupolar \ce{(NaCl)4} basis, the |
145 |
> |
electrostatic energy per ion converges more rapidly to the Madelung |
146 |
> |
energy than the dipolar approximation.\cite{Wolf92} To find the |
147 |
> |
correct Madelung constant, Lacman suggested that the NaCl structure |
148 |
> |
could be constructed in a way that the finite crystal terminates with |
149 |
> |
complete \ce{(NaCl)4} molecules.\cite{Lacman65} The central ion sees |
150 |
> |
what is effectively a set of octupoles at large distances. These facts |
151 |
> |
suggest that the Madelung constants are relatively short ranged for |
152 |
> |
perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful |
153 |
> |
application of Wolf's method are able to obtain accurate estimates of |
154 |
> |
Madelung constants using relatively short cutoff radii. |
155 |
> |
|
156 |
> |
Direct truncation of interactions at a cutoff radius creates numerical |
157 |
> |
errors. Wolf \textit{et al.} argued that truncation errors are due |
158 |
> |
to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To |
159 |
> |
neutralize this charge they proposed placing an image charge on the |
160 |
> |
surface of the cutoff sphere for every real charge inside the cutoff. |
161 |
> |
These charges are present for the evaluation of both the pair |
162 |
> |
interaction energy and the force, although the force expression |
163 |
> |
maintained a discontinuity at the cutoff sphere. In the original Wolf |
164 |
> |
formulation, the total energy for the charge and image were not equal |
165 |
> |
to the integral of their force expression, and as a result, the total |
166 |
> |
energy would not be conserved in molecular dynamics (MD) |
167 |
> |
simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and |
168 |
> |
Gezelter later proposed shifted force variants of the Wolf method with |
169 |
> |
commensurate force and energy expressions that do not exhibit this |
170 |
> |
problem.\cite{Fennell:2006lq} Related real-space methods were also |
171 |
> |
proposed by Chen \textit{et |
172 |
> |
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
173 |
> |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has used |
174 |
> |
neutralization of the higher order moments for the calculation of the |
175 |
> |
electrostatic interaction of the point charge |
176 |
|
systems.\cite{Fukuda:2013sf} |
177 |
|
|
178 |
|
\begin{figure} |
192 |
|
\label{fig:schematic} |
193 |
|
\end{figure} |
194 |
|
|
195 |
< |
The direct truncation of interactions at a cutoff radius creates |
196 |
< |
truncation defects. Wolf \textit{et al.} argued that |
197 |
< |
truncation errors are due to net charge remaining inside the cutoff |
198 |
< |
sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed |
199 |
< |
placing an image charge on the surface of the cutoff sphere for every |
200 |
< |
real charge inside the cutoff. These charges are present for the |
201 |
< |
evaluation of both the pair interaction energy and the force, although |
202 |
< |
the force expression maintained a discontinuity at the cutoff sphere. |
203 |
< |
In the original Wolf formulation, the total energy for the charge and |
204 |
< |
image were not equal to the integral of their force expression, and as |
176 |
< |
a result, the total energy would not be conserved in molecular |
177 |
< |
dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and |
178 |
< |
Fennel and Gezelter later proposed shifted force variants of the Wolf |
179 |
< |
method with commensurate force and energy expressions that do not |
180 |
< |
exhibit this problem.\cite{Fennell:2006lq} Related real-space |
181 |
< |
methods were also proposed by Chen \textit{et |
182 |
< |
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
183 |
< |
and by Wu and Brooks.\cite{Wu:044107} |
184 |
< |
|
185 |
< |
Considering the interaction of one central ion in an ionic crystal |
186 |
< |
with a portion of the crystal at some distance, the effective Columbic |
187 |
< |
potential is found to decrease as $r^{-5}$. If one views the \ce{NaCl} |
188 |
< |
crystal as a simple cubic (SC) structure with an octupolar |
189 |
< |
\ce{(NaCl)4} basis, the electrostatic energy per ion converges more |
190 |
< |
rapidly to the Madelung energy than the dipolar |
191 |
< |
approximation.\cite{Wolf92} To find the correct Madelung constant, |
192 |
< |
Lacman suggested that the NaCl structure could be constructed in a way |
193 |
< |
that the finite crystal terminates with complete \ce{(NaCl)4} |
194 |
< |
molecules.\cite{Lacman65} The central ion sees what is effectively a |
195 |
< |
set of octupoles at large distances. These facts suggest that the |
196 |
< |
Madelung constants are relatively short ranged for perfect ionic |
197 |
< |
crystals.\cite{Wolf:1999dn} |
198 |
< |
|
199 |
< |
One can make a similar argument for crystals of point multipoles. The |
200 |
< |
Luttinger and Tisza treatment of energy constants for dipolar lattices |
201 |
< |
utilizes 24 basis vectors that contain dipoles at the eight corners of |
202 |
< |
a unit cube. Only three of these basis vectors, $X_1, Y_1, |
203 |
< |
\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have |
204 |
< |
zero net dipole and retain contributions only from higher order |
205 |
< |
multipoles. The effective interaction between a dipole at the center |
195 |
> |
One can make a similar effective range argument for crystals of point |
196 |
> |
\textit{multipoles}. The Luttinger and Tisza treatment of energy |
197 |
> |
constants for dipolar lattices utilizes 24 basis vectors that contain |
198 |
> |
dipoles at the eight corners of a unit cube. Only three of these |
199 |
> |
basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
200 |
> |
moments, while the rest have zero net dipole and retain contributions |
201 |
> |
only from higher order multipoles. The lowest energy crystalline |
202 |
> |
structures are built out of basis vectors that have only residual |
203 |
> |
quadrupolar moments (e.g. the $Z_5$ array). In these low energy |
204 |
> |
structures, the effective interaction between a dipole at the center |
205 |
|
of a crystal and a group of eight dipoles farther away is |
206 |
|
significantly shorter ranged than the $r^{-3}$ that one would expect |
207 |
|
for raw dipole-dipole interactions. Only in crystals which retain a |
258 |
|
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
259 |
|
|
260 |
|
\subsection{The damping function} |
261 |
< |
The damping function used in our research has been discussed in detail |
262 |
< |
in the first paper of this series.\cite{PaperI} The radial kernel |
263 |
< |
$1/r$ for the interactions between point charges can be replaced by |
264 |
< |
the complementary error function $\mathrm{erfc}(\alpha r)/r$ to |
265 |
< |
accelerate the rate of convergence, where $\alpha$ is a damping |
266 |
< |
parameter with units of inverse distance. Altering the value of |
267 |
< |
$\alpha$ is equivalent to changing the width of Gaussian charge |
268 |
< |
distributions that replace each point charge -- Gaussian overlap |
269 |
< |
integrals yield complementary error functions when truncated at a |
270 |
< |
finite distance. |
261 |
> |
The damping function has been discussed in detail in the first paper |
262 |
> |
of this series.\cite{PaperI} The $1/r$ Coulombic kernel for the |
263 |
> |
interactions between point charges can be replaced by the |
264 |
> |
complementary error function $\mathrm{erfc}(\alpha r)/r$ to accelerate |
265 |
> |
convergence, where $\alpha$ is a damping parameter with units of |
266 |
> |
inverse distance. Altering the value of $\alpha$ is equivalent to |
267 |
> |
changing the width of Gaussian charge distributions that replace each |
268 |
> |
point charge, as Coulomb integrals with Gaussian charge distributions |
269 |
> |
produce complementary error functions when truncated at a finite |
270 |
> |
distance. |
271 |
|
|
272 |
< |
By using suitable value of damping alpha ($\alpha \sim 0.2$) for a |
273 |
< |
cutoff radius ($r_{c}=9 A$), Fennel and Gezelter produced very good |
274 |
< |
agreement with SPME for the interaction energies, forces and torques |
275 |
< |
for charge-charge interactions.\cite{Fennell:2006lq} |
272 |
> |
With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method |
273 |
> |
produced very good agreement with SPME for interaction energies, |
274 |
> |
forces and torques for charge-charge |
275 |
> |
interactions.\cite{Fennell:2006lq} |
276 |
|
|
277 |
|
\subsection{Point multipoles in molecular modeling} |
278 |
|
Coarse-graining approaches which treat entire molecular subsystems as |
279 |
|
a single rigid body are now widely used. A common feature of many |
280 |
|
coarse-graining approaches is simplification of the electrostatic |
281 |
|
interactions between bodies so that fewer site-site interactions are |
282 |
< |
required to compute configurational energies. Many coarse-grained |
283 |
< |
molecular structures would normally consist of equal positive and |
285 |
< |
negative charges, and rather than use multiple site-site interactions, |
286 |
< |
the interaction between higher order multipoles can also be used to |
287 |
< |
evaluate a single molecule-molecule |
288 |
< |
interaction.\cite{Ren06,Essex10,Essex11} |
282 |
> |
required to compute configurational |
283 |
> |
energies.\cite{Ren06,Essex10,Essex11} |
284 |
|
|
285 |
|
Because electrons in a molecule are not localized at specific points, |
286 |
< |
the assignment of partial charges to atomic centers is a relatively |
287 |
< |
rough approximation. Atomic sites can also be assigned point |
288 |
< |
multipoles and polarizabilities to increase the accuracy of the |
289 |
< |
molecular model. Recently, water has been modeled with point |
290 |
< |
multipoles up to octupolar order using the soft sticky |
296 |
< |
dipole-quadrupole-octupole (SSDQO) |
286 |
> |
the assignment of partial charges to atomic centers is always an |
287 |
> |
approximation. Atomic sites can also be assigned point multipoles and |
288 |
> |
polarizabilities to increase the accuracy of the molecular model. |
289 |
> |
Recently, water has been modeled with point multipoles up to octupolar |
290 |
> |
order using the soft sticky dipole-quadrupole-octupole (SSDQO) |
291 |
|
model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
292 |
|
multipoles up to quadrupolar order have also been coupled with point |
293 |
|
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
294 |
< |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} But |
295 |
< |
using point multipole with the real space truncation without |
296 |
< |
accounting for multipolar neutrality will create energy conservation |
303 |
< |
issues in molecular dynamics (MD) simulations. |
294 |
> |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However, |
295 |
> |
truncating point multipoles without smoothing the forces and torques |
296 |
> |
will create energy conservation issues in molecular dynamics simulations. |
297 |
|
|
298 |
|
In this paper we test a set of real-space methods that were developed |
299 |
|
for point multipolar interactions. These methods extend the damped |