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} |
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 a reference method. The tests were carried out in a variety |
74 |
< |
of condensed-phase environments which were designed to test all |
75 |
< |
levels of the multipole-multipole interactions. Comparisons of the |
76 |
< |
energy differences between configurations, molecular forces, and |
77 |
< |
torques were used to analyze how well the real-space models perform |
78 |
< |
relative to the more computationally expensive Ewald sum. We have |
79 |
< |
also investigated the energy conservation properties of the new |
80 |
< |
methods in molecular dynamics simulations using all of these |
81 |
< |
methods. The SP method shows excellent agreement with |
82 |
< |
configurational energy differences, forces, and torques, and would |
83 |
< |
be suitable for use in Monte Carlo calculations. Of the two new |
84 |
< |
shifted-force methods, the GSF approach shows the best agreement |
85 |
< |
with Ewald-derived energies, forces, and torques and exhibits energy |
86 |
< |
conservation properties that make it an excellent choice for |
87 |
< |
efficiently computing electrostatic interactions in molecular |
88 |
< |
dynamics simulations. |
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. |
89 |
|
\end{abstract} |
90 |
|
|
91 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
92 |
|
% Classification Scheme. |
93 |
< |
\keywords{Electrostatics, Multipoles, Real-space} |
93 |
> |
%\keywords{Electrostatics, Multipoles, Real-space} |
94 |
|
|
95 |
|
\maketitle |
96 |
|
|
122 |
|
interfaces.\cite{Parry:1975if,Parry:1976fq,Clarke77,Perram79,Rhee:1989kl} |
123 |
|
To simulate interfacial systems, Parry's extension of the 3D Ewald sum |
124 |
|
is appropriate for slab geometries.\cite{Parry:1975if} The inherent |
125 |
< |
periodicity in the Ewald’s method can also be problematic for |
126 |
< |
interfacial molecular systems.\cite{Fennell:2006lq} Modified Ewald |
127 |
< |
methods that were developed to handle two-dimensional (2D) |
128 |
< |
electrostatic interactions in interfacial systems have not had similar |
129 |
< |
particle-mesh treatments.\cite{Parry:1975if, Parry:1976fq, Clarke77, |
130 |
< |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} |
125 |
> |
periodicity in the Ewald method can also be problematic for molecular |
126 |
> |
interfaces.\cite{Fennell:2006lq} Modified Ewald methods that were |
127 |
> |
developed to handle two-dimensional (2D) electrostatic interactions in |
128 |
> |
interfacial systems have not seen similar particle-mesh |
129 |
> |
treatments,\cite{Parry:1975if, Parry:1976fq, Clarke77, |
130 |
> |
Perram79,Rhee:1989kl,Spohr:1997sf,Yeh:1999oq} and still scale poorly |
131 |
> |
with system size. |
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 systems is actually short ranged.\cite{Wolf92,Wolf95}. For |
138 |
< |
an ordered lattice (e.g., when computing the Madelung constant of an |
139 |
< |
ionic solid), the material can be considered as a set of ions |
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:NaCl}). For this reason, careful |
143 |
< |
applications of Wolf's method are able to obtain accurate estimates of |
144 |
< |
Madelung constants using relatively short cutoff radii. Recently, |
145 |
< |
Fukuda used neutralization of the higher order moments for the |
146 |
< |
calculation of the electrostatic interaction of the point charges |
147 |
< |
system.\cite{Fukuda:2013sf} |
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 |
147 |
> |
systems.\cite{Fukuda:2013sf} |
148 |
|
|
149 |
< |
\begin{figure}[h!] |
149 |
> |
\begin{figure} |
150 |
|
\centering |
151 |
< |
\includegraphics[width=0.50 \textwidth]{chargesystem.pdf} |
152 |
< |
\caption{Top: NaCl crystal showing how spherical truncation can |
153 |
< |
breaking effective charge ordering, and how complete \ce{(NaCl)4} |
154 |
< |
molecules interact with the central ion. Bottom: A dipolar |
155 |
< |
crystal exhibiting similar behavior and illustrating how the |
156 |
< |
effective dipole-octupole interactions can be disrupted by |
157 |
< |
spherical truncation.} |
158 |
< |
\label{fig:NaCl} |
151 |
> |
\includegraphics[width=\linewidth]{schematic.pdf} |
152 |
> |
\caption{Top: Ionic systems exhibit local clustering of dissimilar |
153 |
> |
charges (in the smaller grey circle), so interactions are |
154 |
> |
effectively charge-multipole at longer distances. With hard |
155 |
> |
cutoffs, motion of individual charges in and out of the cutoff |
156 |
> |
sphere can break the effective multipolar ordering. Bottom: |
157 |
> |
dipolar crystals and fluids have a similar effective |
158 |
> |
\textit{quadrupolar} ordering (in the smaller grey circles), and |
159 |
> |
orientational averaging helps to reduce the effective range of the |
160 |
> |
interactions in the fluid. Placement of reversed image multipoles |
161 |
> |
on the surface of the cutoff sphere recovers the effective |
162 |
> |
higher-order multipole behavior.} |
163 |
> |
\label{fig:schematic} |
164 |
|
\end{figure} |
165 |
|
|
166 |
|
The direct truncation of interactions at a cutoff radius creates |
167 |
< |
truncation defects. Wolf \textit{et al.} further argued that |
167 |
> |
truncation defects. Wolf \textit{et al.} argued that |
168 |
|
truncation errors are due to net charge remaining inside the cutoff |
169 |
|
sphere.\cite{Wolf:1999dn} To neutralize this charge they proposed |
170 |
|
placing an image charge on the surface of the cutoff sphere for every |
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 be decreasing as $r^{-5}$. If one views the |
188 |
< |
\ce{NaCl} crystal as simple cubic (SC) structure with an octupolar |
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} Any charge in a NaCl crystal is surrounded |
195 |
< |
by opposite charges. Similarly for each pair of charges, there is an |
196 |
< |
opposite pair of charge adjacent to it. The central ion sees what is |
197 |
< |
effectively a set of octupoles at large distances. These facts suggest |
192 |
< |
that the Madelung constants are relatively short ranged for perfect |
193 |
< |
ionic crystals.\cite{Wolf:1999dn} |
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 |
212 |
|
unstable. |
213 |
|
|
214 |
|
In ionic crystals, real-space truncation can break the effective |
215 |
< |
multipolar arrangements (see Fig. \ref{fig:NaCl}), causing significant |
216 |
< |
swings in the electrostatic energy as individual ions move back and |
217 |
< |
forth across the boundary. This is why the image charges are |
215 |
> |
multipolar arrangements (see Fig. \ref{fig:schematic}), causing |
216 |
> |
significant swings in the electrostatic energy as individual ions move |
217 |
> |
back and forth across the boundary. This is why the image charges are |
218 |
|
necessary for the Wolf sum to exhibit rapid convergence. Similarly, |
219 |
|
the real-space truncation of point multipole interactions breaks |
220 |
|
higher order multipole arrangements, and image multipoles are required |
221 |
|
for real-space treatments of electrostatic energies. |
222 |
+ |
|
223 |
+ |
The shorter effective range of electrostatic interactions is not |
224 |
+ |
limited to perfect crystals, but can also apply in disordered fluids. |
225 |
+ |
Even at elevated temperatures, there is, on average, local charge |
226 |
+ |
balance in an ionic liquid, where each positive ion has surroundings |
227 |
+ |
dominated by negaitve ions and vice versa. The reversed-charge images |
228 |
+ |
on the cutoff sphere that are integral to the Wolf and DSF approaches |
229 |
+ |
retain the effective multipolar interactions as the charges traverse |
230 |
+ |
the cutoff boundary. |
231 |
+ |
|
232 |
+ |
In multipolar fluids (see Fig. \ref{fig:schematic}) there is |
233 |
+ |
significant orientational averaging that additionally reduces the |
234 |
+ |
effect of long-range multipolar interactions. The image multipoles |
235 |
+ |
that are introduced in the TSF, GSF, and SP methods mimic this effect |
236 |
+ |
and reduce the effective range of the multipolar interactions as |
237 |
+ |
interacting molecules traverse each other's cutoff boundaries. |
238 |
|
|
239 |
|
% Because of this reason, although the nature of electrostatic |
240 |
|
% interaction short ranged, the hard cutoff sphere creates very large |
340 |
|
where the multipole operator for site $\bf a$, $\hat{M}_{\bf a}$, is |
341 |
|
expressed in terms of the point charge, $C_{\bf a}$, dipole, ${\bf D}_{\bf |
342 |
|
a}$, and quadrupole, $\mathbf{Q}_{\bf a}$, for object |
343 |
< |
$\bf a$. |
343 |
> |
$\bf a$, etc. |
344 |
|
|
345 |
|
% Interactions between multipoles can be expressed as higher derivatives |
346 |
|
% of the bare Coulomb potential, so one way of ensuring that the forces |
368 |
|
\label{generic} |
369 |
|
\end{equation} |
370 |
|
where $f_n(r)$ is a shifted kernel that is appropriate for the order |
371 |
< |
of the interaction, with $n=0$ for charge-charge, $n=1$ for |
372 |
< |
charge-dipole, $n=2$ for charge-quadrupole and dipole-dipole, $n=3$ |
373 |
< |
for dipole-quadrupole, and $n=4$ for quadrupole-quadrupole. To ensure |
374 |
< |
smooth convergence of the energy, force, and torques, a Taylor |
375 |
< |
expansion with $n$ terms must be performed at cutoff radius ($r_c$) to |
376 |
< |
obtain $f_n(r)$. |
371 |
> |
of the interaction (see Ref. \onlinecite{PaperI}), with $n=0$ for |
372 |
> |
charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole |
373 |
> |
and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for |
374 |
> |
quadrupole-quadrupole. To ensure smooth convergence of the energy, |
375 |
> |
force, and torques, a Taylor expansion with $n$ terms must be |
376 |
> |
performed at cutoff radius ($r_c$) to obtain $f_n(r)$. |
377 |
|
|
378 |
|
% To carry out the same procedure for a damped electrostatic kernel, we |
379 |
|
% replace $1/r$ in the Coulomb potential with $\text{erfc}(\alpha r)/r$. |
413 |
|
connection to unmodified electrostatics as well as the smooth |
414 |
|
transition to zero in both these functions as $r\rightarrow r_c$. The |
415 |
|
electrostatic forces and torques acting on the central multipole due |
416 |
< |
to another site within cutoff sphere are derived from |
416 |
> |
to another site within the cutoff sphere are derived from |
417 |
|
Eq.~\ref{generic}, accounting for the appropriate number of |
418 |
|
derivatives. Complete energy, force, and torque expressions are |
419 |
|
presented in the first paper in this series (Reference |
427 |
|
shifted smoothly by finding the gradient for two interacting dipoles |
428 |
|
which have been projected onto the surface of the cutoff sphere |
429 |
|
without changing their relative orientation, |
430 |
< |
\begin{displaymath} |
430 |
> |
\begin{equation} |
431 |
|
U_{D_{\bf a}D_{\bf b}}(r) = U_{D_{\bf a}D_{\bf b}}(r) - |
432 |
|
U_{D_{\bf a} D_{\bf b}}(r_c) |
433 |
|
- (r_{ab}-r_c) ~~~\hat{r}_{ab} \cdot |
434 |
< |
\vec{\nabla} U_{D_{\bf a}D_{\bf b}}(r) \Big \lvert _{r_c} |
435 |
< |
\end{displaymath} |
434 |
> |
\nabla U_{D_{\bf a}D_{\bf b}}(r_c). |
435 |
> |
\end{equation} |
436 |
|
Here the lab-frame orientations of the two dipoles, $\mathbf{D}_{\bf |
437 |
|
a}$ and $\mathbf{D}_{\bf b}$, are retained at the cutoff distance |
438 |
|
(although the signs are reversed for the dipole that has been |
455 |
|
In general, the gradient shifted potential between a central multipole |
456 |
|
and any multipolar site inside the cutoff radius is given by, |
457 |
|
\begin{equation} |
458 |
< |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
459 |
< |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r} |
460 |
< |
\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert _{r_c} \right] |
458 |
> |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
459 |
> |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{\mathbf{r}} |
460 |
> |
\cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
461 |
|
\label{generic2} |
462 |
|
\end{equation} |
463 |
|
where the sum describes a separate force-shifting that is applied to |
464 |
< |
each orientational contribution to the energy. |
464 |
> |
each orientational contribution to the energy. In this expression, |
465 |
> |
$\hat{\mathbf{r}}$ is the unit vector connecting the two multipoles |
466 |
> |
($a$ and $b$) in space, and $\hat{\mathbf{a}}$ and $\hat{\mathbf{b}}$ |
467 |
> |
represent the orientations the multipoles. |
468 |
|
|
469 |
|
The third term converges more rapidly than the first two terms as a |
470 |
|
function of radius, hence the contribution of the third term is very |
471 |
|
small for large cutoff radii. The force and torque derived from |
472 |
< |
equation \ref{generic2} are consistent with the energy expression and |
472 |
> |
Eq. \ref{generic2} are consistent with the energy expression and |
473 |
|
approach zero as $r \rightarrow r_c$. Both the GSF and TSF methods |
474 |
|
can be considered generalizations of the original DSF method for |
475 |
|
higher order multipole interactions. GSF and TSF are also identical up |
477 |
|
the energy, force and torque for higher order multipole-multipole |
478 |
|
interactions. Complete energy, force, and torque expressions for the |
479 |
|
GSF potential are presented in the first paper in this series |
480 |
< |
(Reference~\onlinecite{PaperI}) |
480 |
> |
(Reference~\onlinecite{PaperI}). |
481 |
|
|
482 |
|
|
483 |
|
\subsection{Shifted potential (SP) } |
491 |
|
effectively shifts the total potential to zero at the cutoff radius, |
492 |
|
\begin{equation} |
493 |
|
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
494 |
< |
U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
494 |
> |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
495 |
|
\label{eq:SP} |
496 |
|
\end{equation} |
497 |
|
where the sum describes separate potential shifting that is done for |
980 |
|
energy over time, $\delta E_1$, and the standard deviation of energy |
981 |
|
fluctuations around this drift $\delta E_0$. Both of the |
982 |
|
shifted-force methods (GSF and TSF) provide excellent energy |
983 |
< |
conservation (drift less than $10^{-6}$ kcal / mol / ns / particle), |
983 |
> |
conservation (drift less than $10^{-5}$ kcal / mol / ns / particle), |
984 |
|
while the hard cutoff is essentially unusable for molecular dynamics. |
985 |
|
SP provides some benefit over the hard cutoff because the energetic |
986 |
|
jumps that happen as particles leave and enter the cutoff sphere are |