| 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 an |
| 138 |
< |
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:schematic}). 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} |
| 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} |
| 150 |
|
\centering |
| 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 in order at longer distances. With |
| 155 |
< |
hard cutoffs, motion of individual charges in and out of the |
| 156 |
< |
cutoff sphere can break the effective multipolar ordering. |
| 157 |
< |
Bottom: dipolar crystals and fluids have a similar effective |
| 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 |
| 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 |
| 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 |
| 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 |
| 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} |
| 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 |
| 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 |
