35 |
|
%\linenumbers\relax % Commence numbering lines |
36 |
|
\usepackage{amsmath} |
37 |
|
\usepackage{times} |
38 |
< |
\usepackage{mathptm} |
38 |
> |
\usepackage{mathptmx} |
39 |
|
\usepackage{tabularx} |
40 |
|
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
41 |
|
\usepackage{url} |
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), |
66 |
> |
We report on tests of the real-space shifted potential (SP), |
67 |
|
gradient-shifted force (GSF), and Taylor-shifted force (TSF) methods |
68 |
< |
for multipoles that were developed in the first paper in this series |
69 |
< |
against a reference method. The tests were carried out in a variety |
70 |
< |
of condensed-phase environments which were designed to test all |
71 |
< |
levels of the multipole-multipole 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 sum. We have |
75 |
< |
also investigated the energy conservation properties of the new |
76 |
< |
methods in molecular dynamics simulations using all of these |
77 |
< |
methods. The SP method shows excellent agreement with |
68 |
> |
for multipole interactions developed in the first paper in this |
69 |
> |
series, using the multipolar Ewald sum as a reference method. The |
70 |
> |
tests were carried out in a variety of condensed-phase environments |
71 |
> |
designed to test up to quadrupole-quadrupole interactions. |
72 |
> |
Comparisons of the energy differences between configurations, |
73 |
> |
molecular forces, and torques were used to analyze how well the |
74 |
> |
real-space models perform relative to the more computationally |
75 |
> |
expensive Ewald treatment. We have also investigated the energy |
76 |
> |
conservation properties of the new methods in molecular dynamics |
77 |
> |
simulations. The SP method shows excellent agreement with |
78 |
|
configurational energy differences, forces, and torques, and would |
79 |
|
be suitable for use in Monte Carlo calculations. Of the two new |
80 |
|
shifted-force methods, the GSF approach shows the best agreement |
81 |
< |
with Ewald-derived energies, forces, and torques and exhibits energy |
82 |
< |
conservation properties that make it an excellent choice for |
83 |
< |
efficiently computing electrostatic interactions in molecular |
81 |
> |
with Ewald-derived energies, forces, and torques and also exhibits |
82 |
> |
energy conservation properties that make it an excellent choice for |
83 |
> |
efficient computation of electrostatic interactions in molecular |
84 |
|
dynamics simulations. |
85 |
|
\end{abstract} |
86 |
|
|
87 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
88 |
|
% Classification Scheme. |
89 |
< |
\keywords{Electrostatics, Multipoles, Real-space} |
89 |
> |
%\keywords{Electrostatics, Multipoles, Real-space} |
90 |
|
|
91 |
|
\maketitle |
92 |
|
|
97 |
– |
|
93 |
|
\section{\label{sec:intro}Introduction} |
94 |
|
Computing the interactions between electrostatic sites is one of the |
95 |
< |
most expensive aspects of molecular simulations, which is why there |
96 |
< |
have been significant efforts to develop practical, efficient and |
97 |
< |
convergent methods for handling these interactions. Ewald's method is |
98 |
< |
perhaps the best known and most accurate method for evaluating |
99 |
< |
energies, forces, and torques in explicitly-periodic simulation |
100 |
< |
cells. In this approach, the conditionally convergent electrostatic |
101 |
< |
energy is converted into two absolutely convergent contributions, one |
102 |
< |
which is carried out in real space with a cutoff radius, and one in |
103 |
< |
reciprocal space.\cite{Clarke:1986eu,Woodcock75} |
95 |
> |
most expensive aspects of molecular simulations. There have been |
96 |
> |
significant efforts to develop practical, efficient and convergent |
97 |
> |
methods for handling these interactions. Ewald's method is perhaps the |
98 |
> |
best known and most accurate method for evaluating energies, forces, |
99 |
> |
and torques in explicitly-periodic simulation cells. In this approach, |
100 |
> |
the conditionally convergent electrostatic energy is converted into |
101 |
> |
two absolutely convergent contributions, one which is carried out in |
102 |
> |
real space with a cutoff radius, and one in reciprocal |
103 |
> |
space. BETTER CITATIONS\cite{Clarke:1986eu,Woodcock75} |
104 |
|
|
105 |
|
When carried out as originally formulated, the reciprocal-space |
106 |
|
portion of the Ewald sum exhibits relatively poor computational |
107 |
|
scaling, making it prohibitive for large systems. By utilizing |
108 |
|
particle meshes and three dimensional fast Fourier transforms (FFT), |
109 |
|
the particle-mesh Ewald (PME), particle-particle particle-mesh Ewald |
110 |
< |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) methods can decrease |
111 |
< |
the computational cost from $O(N^2)$ down to $O(N \log |
110 |
> |
(P\textsuperscript{3}ME), and smooth particle mesh Ewald (SPME) |
111 |
> |
methods can decrease the computational cost from $O(N^2)$ down to $O(N |
112 |
> |
\log |
113 |
|
N)$.\cite{Takada93,Gunsteren94,Gunsteren95,Darden:1993pd,Essmann:1995pb}. |
114 |
|
|
115 |
< |
Because of the artificial periodicity required for the Ewald sum, the |
120 |
< |
method may require modification to compute interactions for |
115 |
> |
Because of the artificial periodicity required for the Ewald sum, |
116 |
|
interfacial molecular systems such as membranes and liquid-vapor |
117 |
< |
interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
118 |
< |
To simulate interfacial systems, Parry's extension of the 3D Ewald sum |
119 |
< |
is appropriate for slab geometries.\cite{Parry:1975if} The inherent |
120 |
< |
periodicity in the Ewald’s method can also be problematic for |
121 |
< |
interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald |
122 |
< |
methods that were developed to handle two-dimensional (2D) |
123 |
< |
electrostatic interactions in interfacial systems have not had similar |
124 |
< |
particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77, |
125 |
< |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} |
117 |
> |
interfaces require modifications to the |
118 |
> |
method.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
119 |
> |
Parry's extension of the three dimensional Ewald sum is appropriate |
120 |
> |
for slab geometries.\cite{Parry:1975if} Modified Ewald methods that |
121 |
> |
were developed to handle two-dimensional (2D) electrostatic |
122 |
> |
interactions in interfacial systems have not seen similar |
123 |
> |
particle-mesh treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
124 |
> |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly |
125 |
> |
with system size. The inherent periodicity in the Ewald’s method can |
126 |
> |
also be problematic for interfacial molecular |
127 |
> |
systems.\cite{Fennell:2006lq} |
128 |
|
|
129 |
|
\subsection{Real-space methods} |
130 |
|
Wolf \textit{et al.}\cite{Wolf:1999dn} proposed a real space $O(N)$ |
131 |
|
method for calculating electrostatic interactions between point |
132 |
< |
charges. They argued that the effective Coulomb interaction in |
133 |
< |
condensed systems is actually short ranged.\cite{Wolf92,Wolf95}. For |
134 |
< |
an ordered lattice (e.g., when computing the Madelung constant of an |
135 |
< |
ionic solid), the material can be considered as a set of ions |
136 |
< |
interacting with neutral dipolar or quadrupolar ``molecules'' giving |
137 |
< |
an effective distance dependence for the electrostatic interactions of |
138 |
< |
$r^{-5}$ (see figure \ref{fig:NaCl}). For this reason, careful |
139 |
< |
applications of Wolf's method are able to obtain accurate estimates of |
140 |
< |
Madelung constants using relatively short cutoff radii. Recently, |
141 |
< |
Fukuda used neutralization of the higher order moments for the |
142 |
< |
calculation of the electrostatic interaction of the point charges |
143 |
< |
system.\cite{Fukuda:2013sf} |
132 |
> |
charges. They argued that the effective Coulomb interaction in most |
133 |
> |
condensed phase systems is effectively short |
134 |
> |
ranged.\cite{Wolf92,Wolf95} For an ordered lattice (e.g., when |
135 |
> |
computing the Madelung constant of an ionic solid), the material can |
136 |
> |
be considered as a set of ions interacting with neutral dipolar or |
137 |
> |
quadrupolar ``molecules'' giving an effective distance dependence for |
138 |
> |
the electrostatic interactions of $r^{-5}$ (see figure |
139 |
> |
\ref{fig:schematic}). If one views the \ce{NaCl} crystal as a simple |
140 |
> |
cubic (SC) structure with an octupolar \ce{(NaCl)4} basis, the |
141 |
> |
electrostatic energy per ion converges more rapidly to the Madelung |
142 |
> |
energy than the dipolar approximation.\cite{Wolf92} To find the |
143 |
> |
correct Madelung constant, Lacman suggested that the NaCl structure |
144 |
> |
could be constructed in a way that the finite crystal terminates with |
145 |
> |
complete \ce{(NaCl)4} molecules.\cite{Lacman65} The central ion sees |
146 |
> |
what is effectively a set of octupoles at large distances. These facts |
147 |
> |
suggest that the Madelung constants are relatively short ranged for |
148 |
> |
perfect ionic crystals.\cite{Wolf:1999dn} For this reason, careful |
149 |
> |
application of Wolf's method can provide accurate estimates of |
150 |
> |
Madelung constants using relatively short cutoff radii. |
151 |
|
|
152 |
< |
\begin{figure}[h!] |
153 |
< |
\centering |
154 |
< |
\includegraphics[width=0.50 \textwidth]{chargesystem.pdf} |
155 |
< |
\caption{Top: NaCl crystal showing how spherical truncation can |
156 |
< |
breaking effective charge ordering, and how complete \ce{(NaCl)4} |
157 |
< |
molecules interact with the central ion. Bottom: A dipolar |
158 |
< |
crystal exhibiting similar behavior and illustrating how the |
159 |
< |
effective dipole-octupole interactions can be disrupted by |
160 |
< |
spherical truncation.} |
161 |
< |
\label{fig:NaCl} |
162 |
< |
\end{figure} |
163 |
< |
|
164 |
< |
The direct truncation of interactions at a cutoff radius creates |
165 |
< |
truncation defects. Wolf \textit{et al.} further argued that |
166 |
< |
truncation errors are due to net charge remaining inside the cutoff |
167 |
< |
sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed |
164 |
< |
placing an image charge on the surface of the cutoff sphere for every |
165 |
< |
real charge inside the cutoff. These charges are present for the |
166 |
< |
evaluation of both the pair interaction energy and the force, although |
167 |
< |
the force expression maintained a discontinuity at the cutoff sphere. |
168 |
< |
In the original Wolf formulation, the total energy for the charge and |
169 |
< |
image were not equal to the integral of their force expression, and as |
170 |
< |
a result, the total energy would not be conserved in molecular |
171 |
< |
dynamics (MD) simulations.\cite{Zahn:2002hc} Zahn \textit{et al.} and |
172 |
< |
Fennel and Gezelter later proposed shifted force variants of the Wolf |
173 |
< |
method with commensurate force and energy expressions that do not |
174 |
< |
exhibit this problem.\cite{Fennell:2006lq} Related real-space |
175 |
< |
methods were also proposed by Chen \textit{et |
152 |
> |
Direct truncation of interactions at a cutoff radius creates numerical |
153 |
> |
errors. Wolf \textit{et al.} suggest that truncation errors are due |
154 |
> |
to net charge remaining inside the cutoff sphere.\cite{Wolf:1999dn} To |
155 |
> |
neutralize this charge they proposed placing an image charge on the |
156 |
> |
surface of the cutoff sphere for every real charge inside the cutoff. |
157 |
> |
These charges are present for the evaluation of both the pair |
158 |
> |
interaction energy and the force, although the force expression |
159 |
> |
maintains a discontinuity at the cutoff sphere. In the original Wolf |
160 |
> |
formulation, the total energy for the charge and image were not equal |
161 |
> |
to the integral of the force expression, and as a result, the total |
162 |
> |
energy would not be conserved in molecular dynamics (MD) |
163 |
> |
simulations.\cite{Zahn:2002hc} Zahn \textit{et al.}, and Fennel and |
164 |
> |
Gezelter later proposed shifted force variants of the Wolf method with |
165 |
> |
commensurate force and energy expressions that do not exhibit this |
166 |
> |
problem.\cite{Zahn:2002hc,Fennell:2006lq} Related real-space methods |
167 |
> |
were also proposed by Chen \textit{et |
168 |
|
al.}\cite{Chen:2004du,Chen:2006ii,Denesyuk:2008ez,Rodgers:2006nw} |
169 |
< |
and by Wu and Brooks.\cite{Wu:044107} |
169 |
> |
and by Wu and Brooks.\cite{Wu:044107} Recently, Fukuda has successfuly |
170 |
> |
used additional neutralization of higher order moments for systems of |
171 |
> |
point charges.\cite{Fukuda:2013sf} |
172 |
|
|
173 |
< |
Considering the interaction of one central ion in an ionic crystal |
174 |
< |
with a portion of the crystal at some distance, the effective Columbic |
175 |
< |
potential is found to be decreasing as $r^{-5}$. If one views the |
176 |
< |
\ce{NaCl} crystal as simple cubic (SC) structure with an octupolar |
177 |
< |
\ce{(NaCl)4} basis, the electrostatic energy per ion converges more |
178 |
< |
rapidly to the Madelung energy than the dipolar |
179 |
< |
approximation.\cite{Wolf92} To find the correct Madelung constant, |
180 |
< |
Lacman suggested that the NaCl structure could be constructed in a way |
181 |
< |
that the finite crystal terminates with complete \ce{(NaCl)4} |
182 |
< |
molecules.\cite{Lacman65} Any charge in a NaCl crystal is surrounded |
183 |
< |
by opposite charges. Similarly for each pair of charges, there is an |
184 |
< |
opposite pair of charge adjacent to it. The central ion sees what is |
185 |
< |
effectively a set of octupoles at large distances. These facts suggest |
186 |
< |
that the Madelung constants are relatively short ranged for perfect |
187 |
< |
ionic crystals.\cite{Wolf:1999dn} |
173 |
> |
\begin{figure} |
174 |
> |
\centering |
175 |
> |
\includegraphics[width=\linewidth]{schematic.pdf} |
176 |
> |
\caption{Top: Ionic systems exhibit local clustering of dissimilar |
177 |
> |
charges (in the smaller grey circle), so interactions are |
178 |
> |
effectively charge-multipole at longer distances. With hard |
179 |
> |
cutoffs, motion of individual charges in and out of the cutoff |
180 |
> |
sphere can break the effective multipolar ordering. Bottom: |
181 |
> |
dipolar crystals and fluids have a similar effective |
182 |
> |
\textit{quadrupolar} ordering (in the smaller grey circles), and |
183 |
> |
orientational averaging helps to reduce the effective range of the |
184 |
> |
interactions in the fluid. Placement of reversed image multipoles |
185 |
> |
on the surface of the cutoff sphere recovers the effective |
186 |
> |
higher-order multipole behavior.} |
187 |
> |
\label{fig:schematic} |
188 |
> |
\end{figure} |
189 |
|
|
190 |
< |
One can make a similar argument for crystals of point multipoles. The |
191 |
< |
Luttinger and Tisza treatment of energy constants for dipolar lattices |
192 |
< |
utilizes 24 basis vectors that contain dipoles at the eight corners of |
193 |
< |
a unit cube. Only three of these basis vectors, $X_1, Y_1, |
194 |
< |
\mathrm{~and~} Z_1,$ retain net dipole moments, while the rest have |
195 |
< |
zero net dipole and retain contributions only from higher order |
196 |
< |
multipoles. The effective interaction between a dipole at the center |
190 |
> |
One can make a similar effective range argument for crystals of point |
191 |
> |
\textit{multipoles}. The Luttinger and Tisza treatment of energy |
192 |
> |
constants for dipolar lattices utilizes 24 basis vectors that contain |
193 |
> |
dipoles at the eight corners of a unit cube.\cite{LT} Only three of |
194 |
> |
these basis vectors, $X_1, Y_1, \mathrm{~and~} Z_1,$ retain net dipole |
195 |
> |
moments, while the rest have zero net dipole and retain contributions |
196 |
> |
only from higher order multipoles. The lowest-energy crystalline |
197 |
> |
structures are built out of basis vectors that have only residual |
198 |
> |
quadrupolar moments (e.g. the $Z_5$ array). In these low energy |
199 |
> |
structures, the effective interaction between a dipole at the center |
200 |
|
of a crystal and a group of eight dipoles farther away is |
201 |
|
significantly shorter ranged than the $r^{-3}$ that one would expect |
202 |
|
for raw dipole-dipole interactions. Only in crystals which retain a |
206 |
|
unstable. |
207 |
|
|
208 |
|
In ionic crystals, real-space truncation can break the effective |
209 |
< |
multipolar arrangements (see Fig. \ref{fig:NaCl}), causing significant |
210 |
< |
swings in the electrostatic energy as individual ions move back and |
211 |
< |
forth across the boundary. This is why the image charges are |
209 |
> |
multipolar arrangements (see Fig. \ref{fig:schematic}), causing |
210 |
> |
significant swings in the electrostatic energy as individual ions move |
211 |
> |
back and forth across the boundary. This is why the image charges are |
212 |
|
necessary for the Wolf sum to exhibit rapid convergence. Similarly, |
213 |
|
the real-space truncation of point multipole interactions breaks |
214 |
|
higher order multipole arrangements, and image multipoles are required |
215 |
|
for real-space treatments of electrostatic energies. |
216 |
|
|
217 |
+ |
The shorter effective range of electrostatic interactions is not |
218 |
+ |
limited to perfect crystals, but can also apply in disordered fluids. |
219 |
+ |
Even at elevated temperatures, there is local charge balance in an |
220 |
+ |
ionic liquid, where each positive ion has surroundings dominated by |
221 |
+ |
negaitve ions and vice versa. The reversed-charge images on the |
222 |
+ |
cutoff sphere that are integral to the Wolf and DSF approaches retain |
223 |
+ |
the effective multipolar interactions as the charges traverse the |
224 |
+ |
cutoff boundary. |
225 |
+ |
|
226 |
+ |
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
227 |
+ |
significant orientational averaging that additionally reduces the |
228 |
+ |
effect of long-range multipolar interactions. The image multipoles |
229 |
+ |
that are introduced in the TSF, GSF, and SP methods mimic this effect |
230 |
+ |
and reduce the effective range of the multipolar interactions as |
231 |
+ |
interacting molecules traverse each other's cutoff boundaries. |
232 |
+ |
|
233 |
|
% Because of this reason, although the nature of electrostatic |
234 |
|
% interaction short ranged, the hard cutoff sphere creates very large |
235 |
|
% fluctuation in the electrostatic energy for the perfect crystal. In |
240 |
|
% to the non-neutralized value of the higher order moments within the |
241 |
|
% cutoff sphere. |
242 |
|
|
243 |
< |
The forces and torques acting on atomic sites are the fundamental |
244 |
< |
factors driving dynamics in molecular simulations. Fennell and |
245 |
< |
Gezelter proposed the damped shifted force (DSF) energy kernel to |
246 |
< |
obtain consistent energies and forces on the atoms within the cutoff |
247 |
< |
sphere. Both the energy and the force go smoothly to zero as an atom |
248 |
< |
aproaches the cutoff radius. The comparisons of the accuracy these |
249 |
< |
quantities between the DSF kernel and SPME was surprisingly |
250 |
< |
good.\cite{Fennell:2006lq} The DSF method has seen increasing use for |
251 |
< |
calculating electrostatic interactions in molecular systems with |
252 |
< |
relatively uniform charge |
239 |
< |
densities.\cite{Shi:2013ij,Kannam:2012rr,Acevedo13,Space12,English08,Lawrence13,Vergne13} |
243 |
> |
Forces and torques acting on atomic sites are fundamental in driving |
244 |
> |
dynamics in molecular simulations, and the damped shifted force (DSF) |
245 |
> |
energy kernel provides consistent energies and forces on charged atoms |
246 |
> |
within the cutoff sphere. Both the energy and the force go smoothly to |
247 |
> |
zero as an atom aproaches the cutoff radius. The comparisons of the |
248 |
> |
accuracy these quantities between the DSF kernel and SPME was |
249 |
> |
surprisingly good.\cite{Fennell:2006lq} As a result, the DSF method |
250 |
> |
has seen increasing use in molecular systems with relatively uniform |
251 |
> |
charge |
252 |
> |
densities.\cite{English08,Kannam:2012rr,Space12,Lawrence13,Acevedo13,Shi:2013ij,Vergne13} |
253 |
|
|
254 |
|
\subsection{The damping function} |
255 |
< |
The damping function used in our research has been discussed in detail |
256 |
< |
in the first paper of this series.\cite{PaperI} The radial kernel |
257 |
< |
$1/r$ for the interactions between point charges can be replaced by |
258 |
< |
the complementary error function $\mathrm{erfc}(\alpha r)/r$ to |
259 |
< |
accelerate the rate of convergence, where $\alpha$ is a damping |
260 |
< |
parameter with units of inverse distance. Altering the value of |
261 |
< |
$\alpha$ is equivalent to changing the width of Gaussian charge |
262 |
< |
distributions that replace each point charge -- Gaussian overlap |
263 |
< |
integrals yield complementary error functions when truncated at a |
264 |
< |
finite distance. |
255 |
> |
The damping function has been discussed in detail in the first paper |
256 |
> |
of this series.\cite{PaperI} The $1/r$ Coulombic kernel for the |
257 |
> |
interactions between point charges can be replaced by the |
258 |
> |
complementary error function $\mathrm{erfc}(\alpha r)/r$ to accelerate |
259 |
> |
convergence, where $\alpha$ is a damping parameter with units of |
260 |
> |
inverse distance. Altering the value of $\alpha$ is equivalent to |
261 |
> |
changing the width of Gaussian charge distributions that replace each |
262 |
> |
point charge, as Coulomb integrals with Gaussian charge distributions |
263 |
> |
produce complementary error functions when truncated at a finite |
264 |
> |
distance. |
265 |
|
|
266 |
< |
By using suitable value of damping alpha ($\alpha \sim 0.2$) for a |
267 |
< |
cutoff radius ($r_{c}=9 A$), Fennel and Gezelter produced very good |
268 |
< |
agreement with SPME for the interaction energies, forces and torques |
269 |
< |
for charge-charge interactions.\cite{Fennell:2006lq} |
266 |
> |
With moderate damping coefficients, $\alpha \sim 0.2$, the DSF method |
267 |
> |
produced very good agreement with SPME for interaction energies, |
268 |
> |
forces and torques for charge-charge |
269 |
> |
interactions.\cite{Fennell:2006lq} |
270 |
|
|
271 |
|
\subsection{Point multipoles in molecular modeling} |
272 |
|
Coarse-graining approaches which treat entire molecular subsystems as |
273 |
|
a single rigid body are now widely used. A common feature of many |
274 |
|
coarse-graining approaches is simplification of the electrostatic |
275 |
|
interactions between bodies so that fewer site-site interactions are |
276 |
< |
required to compute configurational energies. Many coarse-grained |
277 |
< |
molecular structures would normally consist of equal positive and |
265 |
< |
negative charges, and rather than use multiple site-site interactions, |
266 |
< |
the interaction between higher order multipoles can also be used to |
267 |
< |
evaluate a single molecule-molecule |
268 |
< |
interaction.\cite{Ren06,Essex10,Essex11} |
276 |
> |
required to compute configurational |
277 |
> |
energies.\cite{Ren06,Essex10,Essex11} |
278 |
|
|
279 |
< |
Because electrons in a molecule are not localized at specific points, |
280 |
< |
the assignment of partial charges to atomic centers is a relatively |
281 |
< |
rough approximation. Atomic sites can also be assigned point |
282 |
< |
multipoles and polarizabilities to increase the accuracy of the |
283 |
< |
molecular model. Recently, water has been modeled with point |
284 |
< |
multipoles up to octupolar order using the soft sticky |
276 |
< |
dipole-quadrupole-octupole (SSDQO) |
279 |
> |
Additionally, because electrons in a molecule are not localized at |
280 |
> |
specific points, the assignment of partial charges to atomic centers |
281 |
> |
is always an approximation. For increased accuracy, atomic sites can |
282 |
> |
also be assigned point multipoles and polarizabilities. Recently, |
283 |
> |
water has been modeled with point multipoles up to octupolar order |
284 |
> |
using the soft sticky dipole-quadrupole-octupole (SSDQO) |
285 |
|
model.\cite{Ichiye10_1,Ichiye10_2,Ichiye10_3} Atom-centered point |
286 |
|
multipoles up to quadrupolar order have also been coupled with point |
287 |
|
polarizabilities in the high-quality AMOEBA and iAMOEBA water |
288 |
< |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} But |
289 |
< |
using point multipole with the real space truncation without |
290 |
< |
accounting for multipolar neutrality will create energy conservation |
291 |
< |
issues in molecular dynamics (MD) simulations. |
288 |
> |
models.\cite{Ren:2003uq,Ren:2004kx,Ponder:2010vl,Wang:2013fk} However, |
289 |
> |
truncating point multipoles without smoothing the forces and torques |
290 |
> |
can create energy conservation issues in molecular dynamics |
291 |
> |
simulations. |
292 |
|
|
293 |
|
In this paper we test a set of real-space methods that were developed |
294 |
|
for point multipolar interactions. These methods extend the damped |
295 |
|
shifted force (DSF) and Wolf methods originally developed for |
296 |
|
charge-charge interactions and generalize them for higher order |
297 |
< |
multipoles. The detailed mathematical development of these methods has |
298 |
< |
been presented in the first paper in this series, while this work |
299 |
< |
covers the testing the energies, forces, torques, and energy |
297 |
> |
multipoles. The detailed mathematical development of these methods |
298 |
> |
has been presented in the first paper in this series, while this work |
299 |
> |
covers the testing of energies, forces, torques, and energy |
300 |
|
conservation properties of the methods in realistic simulation |
301 |
|
environments. In all cases, the methods are compared with the |
302 |
< |
reference method, a full multipolar Ewald treatment. |
302 |
> |
reference method, a full multipolar Ewald treatment.\cite{Smith82,Smith98} |
303 |
|
|
304 |
|
|
305 |
|
%\subsection{Conservation of total energy } |
328 |
|
where the multipole operator for site $\bf a$, $\hat{M}_{\bf a}$, is |
329 |
|
expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf |
330 |
|
a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object |
331 |
< |
$\bf a$. |
331 |
> |
$\bf a$, etc. |
332 |
|
|
333 |
|
% Interactions between multipoles can be expressed as higher derivatives |
334 |
|
% of the bare Coulomb potential, so one way of ensuring that the forces |
356 |
|
\label{generic} |
357 |
|
\end{equation} |
358 |
|
where $f_n(r)$ is a shifted kernel that is appropriate for the order |
359 |
< |
of the interaction, with $n=0$ for charge-charge, $n=1$ for |
360 |
< |
charge-dipole, $n=2$ for charge-quadrupole and dipole-dipole, $n=3$ |
361 |
< |
for dipole-quadrupole, and $n=4$ for quadrupole-quadrupole. To ensure |
362 |
< |
smooth convergence of the energy, force, and torques, a Taylor |
363 |
< |
expansion with $n$ terms must be performed at cutoff radius ($r_c$) to |
364 |
< |
obtain $f_n(r)$. |
359 |
> |
of the interaction (see Ref. \onlinecite{PaperI}), with $n=0$ for |
360 |
> |
charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole |
361 |
> |
and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for |
362 |
> |
quadrupole-quadrupole. To ensure smooth convergence of the energy, |
363 |
> |
force, and torques, a Taylor expansion with $n$ terms must be |
364 |
> |
performed at cutoff radius ($r_c$) to obtain $f_n(r)$. |
365 |
|
|
366 |
|
% To carry out the same procedure for a damped electrostatic kernel, we |
367 |
|
% replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. |
401 |
|
connection to unmodified electrostatics as well as the smooth |
402 |
|
transition to zero in both these functions as $r\rightarrow r_c$. The |
403 |
|
electrostatic forces and torques acting on the central multipole due |
404 |
< |
to another site within cutoff sphere are derived from |
404 |
> |
to another site within the cutoff sphere are derived from |
405 |
|
Eq.~\ref{generic}, accounting for the appropriate number of |
406 |
|
derivatives. Complete energy, force, and torque expressions are |
407 |
|
presented in the first paper in this series (Reference |
415 |
|
shifted smoothly by finding the gradient for two interacting dipoles |
416 |
|
which have been projected onto the surface of the cutoff sphere |
417 |
|
without changing their relative orientation, |
418 |
< |
\begin{displaymath} |
418 |
> |
\begin{equation} |
419 |
|
U_{D_{\bf a}D_{\bf b}}(r) = U_{D_{\bf a}D_{\bf b}}(r) - |
420 |
|
U_{D_{\bf a} D_{\bf b}}(r_c) |
421 |
|
- (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot |
422 |
< |
\vec{\nabla} U_{D_{\bf a}D_{\bf b}}(r) \Big \lvert _{r_c} |
423 |
< |
\end{displaymath} |
422 |
> |
\nabla U_{D_{\bf a}D_{\bf b}}(r_c). |
423 |
> |
\end{equation} |
424 |
|
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{\bf |
425 |
|
a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance |
426 |
|
(although the signs are reversed for the dipole that has been |
443 |
|
In general, the gradient shifted potential between a central multipole |
444 |
|
and any multipolar site inside the cutoff radius is given by, |
445 |
|
\begin{equation} |
446 |
< |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
447 |
< |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} |
448 |
< |
\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} \right] |
446 |
> |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
447 |
> |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{\mathbf{r}} |
448 |
> |
\cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
449 |
|
\label{generic2} |
450 |
|
\end{equation} |
451 |
|
where the sum describes a separate force-shifting that is applied to |
452 |
< |
each orientational contribution to the energy. |
452 |
> |
each orientational contribution to the energy. In this expression, |
453 |
> |
$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles |
454 |
> |
($a$ and $b$) in space, and $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$ |
455 |
> |
represent the orientations the multipoles. |
456 |
|
|
457 |
|
The third term converges more rapidly than the first two terms as a |
458 |
|
function of radius, hence the contribution of the third term is very |
459 |
|
small for large cutoff radii. The force and torque derived from |
460 |
< |
equation \ref{generic2} are consistent with the energy expression and |
460 |
> |
Eq. \ref{generic2} are consistent with the energy expression and |
461 |
|
approach zero as $r \rightarrow r_c$. Both the GSF and TSF methods |
462 |
|
can be considered generalizations of the original DSF method for |
463 |
|
higher order multipole interactions. GSF and TSF are also identical up |
465 |
|
the energy, force and torque for higher order multipole-multipole |
466 |
|
interactions. Complete energy, force, and torque expressions for the |
467 |
|
GSF potential are presented in the first paper in this series |
468 |
< |
(Reference~\onlinecite{PaperI}) |
468 |
> |
(Reference~\onlinecite{PaperI}). |
469 |
|
|
470 |
|
|
471 |
|
\subsection{Shifted potential (SP) } |
479 |
|
effectively shifts the total potential to zero at the cutoff radius, |
480 |
|
\begin{equation} |
481 |
|
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
482 |
< |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
482 |
> |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
483 |
|
\label{eq:SP} |
484 |
|
\end{equation} |
485 |
|
where the sum describes separate potential shifting that is done for |
968 |
|
energy over time, $\delta E_1$, and the standard deviation of energy |
969 |
|
fluctuations around this drift $\delta E_0$. Both of the |
970 |
|
shifted-force methods (GSF and TSF) provide excellent energy |
971 |
< |
conservation (drift less than $10^{-6}$ kcal / mol / ns / particle), |
971 |
> |
conservation (drift less than $10^{-5}$ kcal / mol / ns / particle), |
972 |
|
while the hard cutoff is essentially unusable for molecular dynamics. |
973 |
|
SP provides some benefit over the hard cutoff because the energetic |
974 |
|
jumps that happen as particles leave and enter the cutoff sphere are |