| 20 |
|
% Use this file as a source of example code for your aip document. |
| 21 |
|
% Use the file aiptemplate.tex as a template for your document. |
| 22 |
|
\documentclass[% |
| 23 |
< |
aip, |
| 24 |
< |
jmp, |
| 23 |
> |
aip,jcp, |
| 24 |
|
amsmath,amssymb, |
| 25 |
|
preprint,% |
| 26 |
|
% reprint,% |
| 30 |
|
|
| 31 |
|
\usepackage{graphicx}% Include figure files |
| 32 |
|
\usepackage{dcolumn}% Align table columns on decimal point |
| 33 |
< |
\usepackage{bm}% bold math |
| 33 |
> |
%\usepackage{bm}% bold math |
| 34 |
> |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
| 35 |
> |
\usepackage{url} |
| 36 |
> |
|
| 37 |
|
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math |
| 38 |
|
%\linenumbers\relax % Commence numbering lines |
| 39 |
|
|
| 41 |
|
|
| 42 |
|
\preprint{AIP/123-QED} |
| 43 |
|
|
| 44 |
< |
\title[Generalization of the Shifted-Force Potential to Higher-Order Potentials] |
| 45 |
< |
{Generalization of the Shifted-Force Potential to Higher-Order Potentials} |
| 44 |
> |
\title[Taylor-shifted and Gradient-shifted electrostatics for multipoles] |
| 45 |
> |
{Real space alternatives to the Ewald |
| 46 |
> |
Sum. I. Taylor-shifted and Gradient-shifted electrostatics for multipoles} |
| 47 |
|
|
| 48 |
|
\author{Madan Lamichhane} |
| 49 |
|
\affiliation{Department of Physics, University |
| 63 |
|
% but any date may be explicitly specified |
| 64 |
|
|
| 65 |
|
\begin{abstract} |
| 66 |
< |
Over the past several years, there has been increasing interest |
| 67 |
< |
in real space methods for calculating electrostatic interactions |
| 68 |
< |
in computer simulations of condensed molecular systems. We |
| 69 |
< |
have extended our original damped-shifted force (DSF) |
| 70 |
< |
electrostatic kernel and have been able to derive a set of |
| 71 |
< |
interaction models for higher-order multipoles based on |
| 72 |
< |
truncated Taylor expansions around the cutoff. For multipolemultipole |
| 73 |
< |
interactions, we find that each of the distinct |
| 74 |
< |
orientational contributions has a separate radial function to |
| 75 |
< |
ensure that the overall forces and torques vanish at the cutoff |
| 73 |
< |
radius. |
| 66 |
> |
We have extended the original damped-shifted force (DSF) |
| 67 |
> |
electrostatic kernel and have been able to derive two new |
| 68 |
> |
electrostatic potentials for higher-order multipoles that are based |
| 69 |
> |
on truncated Taylor expansions around the cutoff radius. For |
| 70 |
> |
multipole-multipole interactions, we find that each of the distinct |
| 71 |
> |
orientational contributions has a separate radial function to ensure |
| 72 |
> |
that the overall forces and torques vanish at the cutoff radius. In |
| 73 |
> |
this paper, we present energy, force, and torque expressions for the |
| 74 |
> |
new models, and compare these real-space interaction models to exact |
| 75 |
> |
results for ordered arrays of multipoles. |
| 76 |
|
\end{abstract} |
| 77 |
|
|
| 78 |
|
\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 83 |
|
|
| 84 |
|
\section{Introduction} |
| 85 |
|
|
| 86 |
< |
The Coulomb electrostatic interaction is of importance in a number of physical chemistry problems |
| 87 |
< |
[background needed, do we mention gases, liquids, solids?]. |
| 86 |
> |
There is increasing interest in real-space methods for calculating |
| 87 |
> |
electrostatic interactions in computer simulations of condensed |
| 88 |
> |
molecular |
| 89 |
> |
systems.\cite{Wolf99,Zahn02,Kast03,BeckD.A.C._bi0486381,Ma05,Fennell:2006zl,Chen:2004du,Chen:2006ii,Rodgers:2006nw,Denesyuk:2008ez,Izvekov:2008wo} |
| 90 |
> |
The simplest of these techniques was developed by Wolf {\it et al.} |
| 91 |
> |
in their work towards an $\mathcal{O}(N)$ Coulombic sum.\cite{Wolf99} |
| 92 |
> |
Fennell and Gezelter showed that a simple damped shifted force (DSF) |
| 93 |
> |
modification to Wolf's method could give nearly quantitative agreement |
| 94 |
> |
with smooth particle mesh Ewald (SPME)\cite{Essmann95} for quantities |
| 95 |
> |
of interest in Monte Carlo (i.e., configurational energy differences) |
| 96 |
> |
and Molecular Dynamics (i.e., atomic force and molecular torque |
| 97 |
> |
vectors).\cite{Fennell:2006zl} |
| 98 |
|
|
| 99 |
< |
[...] |
| 99 |
> |
The computational efficiency and the accuracy of the DSF method are |
| 100 |
> |
surprisingly good, particularly for systems with uniform charge |
| 101 |
> |
density. Additionally, dielectric constants obtained using DSF and |
| 102 |
> |
similar methods where the force vanishes at $R_\textrm{c}$ are |
| 103 |
> |
essentially quantitative.\cite{Izvekov:2008wo} The DSF and other |
| 104 |
> |
related methods have now been widely investigated,\cite{Hansen:2012uq} |
| 105 |
> |
and DSF is now used routinely in simulations of ionic |
| 106 |
> |
liquids,\cite{doi:10.1021/la400226g,McCann:2013fk} flow in carbon |
| 107 |
> |
nanotubes,\cite{kannam:094701} gas sorption in metal-organic framework |
| 108 |
> |
materials,\cite{Forrest:2012ly} thermal conductivity of methane |
| 109 |
> |
hydrates,\cite{English:2008kx} condensation coefficients of |
| 110 |
> |
water,\cite{Louden:2013ve} and in tribology at solid-liquid-solid |
| 111 |
> |
interfaces.\cite{Tokumasu:2013zr} DSF electrostatics provides a |
| 112 |
> |
compromise between the computational speed of real-space cutoffs and |
| 113 |
> |
the accuracy of fully-periodic Ewald treatments. |
| 114 |
|
|
| 115 |
< |
The methods that we develop in this paper are meant specifically for problems involving interacting rigid molecules which will be described |
| 116 |
< |
in terms of classical mechanics and electrodynamics. From mechanics, the molecule's mass distribution |
| 117 |
< |
determines its total mass and moment of inertia tensor. |
| 118 |
< |
From electrostatics, its charge distribution is conveniently described using multipoles. |
| 119 |
< |
Our goal is to advance methods for handling inter-molecular interactions in molecular dynamics simulations. |
| 120 |
< |
To do this, we must develop consistent approximate equations for |
| 121 |
< |
interaction energies, forces, and torques. |
| 115 |
> |
One common feature of many coarse-graining approaches, which treat |
| 116 |
> |
entire molecular subsystems as a single rigid body, is simplification |
| 117 |
> |
of the electrostatic interactions between these bodies so that fewer |
| 118 |
> |
site-site interactions are required to compute configurational |
| 119 |
> |
energies. Notably, the force matching approaches of Voth and coworkers |
| 120 |
> |
are an exciting development in their ability to represent realistic |
| 121 |
> |
(and {\it reactive}) chemical systems for very large length scales and |
| 122 |
> |
long times. This approach utilized a coarse-graining in interaction |
| 123 |
> |
space (CGIS) which fits an effective force for the coarse grained |
| 124 |
> |
system using a variational force-matching method to a fine-grained |
| 125 |
> |
simulation.\cite{Izvekov:2008wo} |
| 126 |
|
|
| 127 |
< |
[...] |
| 127 |
> |
The coarse-graining approaches of Ren \& coworkers,\cite{Golubkov06} |
| 128 |
> |
and Essex \& |
| 129 |
> |
coworkers,\cite{ISI:000276097500009,ISI:000298664400012} |
| 130 |
> |
both utilize Gay-Berne |
| 131 |
> |
ellipsoids~\cite{Berne72,Gay81,Luckhurst90,Cleaver96,Berardi98,Ravichandran:1999fk,Berardi99,Pasterny00} |
| 132 |
> |
to model dispersive interactions and point multipoles to model |
| 133 |
> |
electrostatics for entire molecules or functional groups. |
| 134 |
|
|
| 135 |
< |
This paper extends the shifted-force potential method |
| 136 |
< |
of Fennel and Gezelter to higher-order multipole interactions. [describe?] |
| 137 |
< |
Extending an idea from Wolf, multipole images are placed on the surface of a |
| 138 |
< |
``cutoff'' sphere of radius $r_c$. These images are used to make all interaction energies, forces, and torques |
| 139 |
< |
be zero for $r < r_c$. |
| 140 |
< |
Two such methods have been developed, both based on Taylor-series expansions. |
| 141 |
< |
The first is applied to the Coulomb kernel of the multipole expansion. The second is |
| 142 |
< |
applied to individual terms for interaction energies in the multipole expansion. |
| 143 |
< |
Because of differences in the initial assumptions, the two methods yield different results. |
| 135 |
> |
Ichiye and coworkers have recently introduced a number of very fast |
| 136 |
> |
water models based on a ``sticky'' multipole model which are |
| 137 |
> |
qualitatively better at reproducing the behavior of real water than |
| 138 |
> |
the more common point-charge models (SPC/E, TIPnP). The point charge |
| 139 |
> |
models are also substantially more computationally demanding than the |
| 140 |
> |
sticky multipole |
| 141 |
> |
approach.\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} The |
| 142 |
> |
SSDQO model requires the use of an approximate multipole expansion |
| 143 |
> |
(AME) as the exact multipole expansion is quite expensive |
| 144 |
> |
(particularly when handled via the Ewald sum).\cite{Ichiye:2006qy} |
| 145 |
> |
Another particularly important use of point multipoles (and multipole |
| 146 |
> |
polarizability) is in the very high-quality AMOEBA water model and |
| 147 |
> |
related force fields.\cite{Ponder:2010fk,schnieders:124114,Ren:2011uq} |
| 148 |
|
|
| 149 |
< |
Also explored here is the effect of replacing the bare Coulomb kernel with that of a smeared |
| 150 |
< |
charge distribution. Thus four methods are compared in this paper: |
| 151 |
< |
(1) Shifted force, Coulomb, method 1; |
| 152 |
< |
(2) Sihfted force, Coulomb, method 2; |
| 153 |
< |
(3) Shifted force, smeared charge, method 1; and |
| 154 |
< |
(4) Shifted force, smeared charge, method 2. |
| 115 |
< |
The last of these methods is our preferred method and is called the Extended Shifted Force Method. |
| 116 |
< |
Subsequent papers will apply this method to various problems of physical and chemical interest. |
| 149 |
> |
Higher-order multipoles present a peculiar issue for molecular |
| 150 |
> |
dynamics. Multipolar interactions are inherently short-ranged, and |
| 151 |
> |
should not need the relatively expensive Ewald treatment. However, |
| 152 |
> |
real-space cutoff methods are normally applied in an orientation-blind |
| 153 |
> |
fashion so multipoles which leave and then re-enter a cutoff sphere in |
| 154 |
> |
a different orientation can cause energy discontinuities. |
| 155 |
|
|
| 156 |
< |
[...] |
| 156 |
> |
This paper outlines an extension of the original DSF electrostatic |
| 157 |
> |
kernel to point multipoles. We present in this paper two distinct |
| 158 |
> |
real-space interaction models for higher-order multipoles based on two |
| 159 |
> |
truncated Taylor expansions that are carried out at the cutoff radius. |
| 160 |
> |
We are calling these models {\bf Taylor-shifted} and {\bf |
| 161 |
> |
Gradient-shifted} electrostatics. Because of differences in the |
| 162 |
> |
initial assumptions, the two methods yield related, but different |
| 163 |
> |
expressions for energies, forces, and torques. |
| 164 |
|
|
| 165 |
+ |
\section{Methodology} |
| 166 |
|
|
| 167 |
< |
\section{Development of the Methods} |
| 167 |
> |
\subsection{Self-neutralization, damping, and force-shifting} |
| 168 |
|
|
| 169 |
+ |
The DSF and Wolf methods operate by neutralizing the total charge |
| 170 |
+ |
contained within the cutoff sphere surrounding each particle. This is |
| 171 |
+ |
accomplished by shifting the potential functions to generate image |
| 172 |
+ |
charges on the surface of the cutoff sphere for each pair interaction |
| 173 |
+ |
computed within $R_\textrm{c}$. An Ewald-like complementary error |
| 174 |
+ |
function damping is applied to the potential to accelerate |
| 175 |
+ |
convergence. The potential for the DSF method, where $\alpha$ is the |
| 176 |
+ |
adjustable damping parameter, is given by |
| 177 |
+ |
\begin{equation*} |
| 178 |
+ |
V_\mathrm{DSF}(r) = C_a C_b \Biggr{[} |
| 179 |
+ |
\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} |
| 180 |
+ |
- \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} + \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} |
| 181 |
+ |
+ \frac{2\alpha}{\pi^{1/2}} |
| 182 |
+ |
\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}} |
| 183 |
+ |
\right)\left(r_{ij}-R_\mathrm{c}\right)\ \Biggr{]} |
| 184 |
+ |
\label{eq:DSFPot} |
| 185 |
+ |
\end{equation*} |
| 186 |
+ |
|
| 187 |
+ |
To insure net charge neutrality within each cutoff sphere, an |
| 188 |
+ |
additional ``self'' term is added to the potential. This term is |
| 189 |
+ |
constant (as long as the charges and cutoff radius do not change), and |
| 190 |
+ |
exists outside the normal pair-loop for molecular simulations. It can |
| 191 |
+ |
be thought of as a contribution from a charge opposite in sign, but |
| 192 |
+ |
equal in magnitude, to the central charge, which has been spread out |
| 193 |
+ |
over the surface of the cutoff sphere. A portion of the self term is |
| 194 |
+ |
identical to the self term in the Ewald summation, and comes from the |
| 195 |
+ |
utilization of the complimentary error function for electrostatic |
| 196 |
+ |
damping.\cite{deLeeuw80,Wolf99} |
| 197 |
+ |
|
| 198 |
+ |
\begin{figure} |
| 199 |
+ |
\includegraphics[width=\linewidth]{SM} |
| 200 |
+ |
\caption{Reversed multipoles are projected onto the surface of the |
| 201 |
+ |
cutoff sphere. The forces, torques, and potential are then smoothly |
| 202 |
+ |
shifted to zero as the sites leave the cutoff region.} |
| 203 |
+ |
\label{fig:shiftedMultipoles} |
| 204 |
+ |
\end{figure} |
| 205 |
+ |
|
| 206 |
|
\subsection{Multipole Expansion} |
| 207 |
|
|
| 208 |
< |
Consider two discrete rigid collections of atoms and ions, denoted as objects $\bf a$ and $\bf b$. |
| 209 |
< |
In the following, we assume |
| 210 |
< |
that the two objects only interact via electrostatics and describe those interactions in terms of |
| 211 |
< |
a multipole expansion. Putting the origin of the coordinate system at the center of mass of $\bf a$, we use |
| 212 |
< |
vectors $\mathbf{r}_k$ to denote the positions of all charges $q_k$ in $\bf a$. |
| 213 |
< |
Then the electrostatic potential of object $\bf a$ at $\mathbf{r}$ is given by |
| 208 |
> |
Consider two discrete rigid collections of point charges, denoted as |
| 209 |
> |
objects $\bf a$ and $\bf b$. In the following, we assume that the two |
| 210 |
> |
objects only interact via electrostatics and describe those |
| 211 |
> |
interactions in terms of a multipole expansion. Putting the origin of |
| 212 |
> |
the coordinate system at the center of mass of $\bf a$, we use vectors |
| 213 |
> |
$\mathbf{r}_k$ to denote the positions of all charges $q_k$ in $\bf |
| 214 |
> |
a$. Then the electrostatic potential of object $\bf a$ at |
| 215 |
> |
$\mathbf{r}$ is given by |
| 216 |
|
\begin{equation} |
| 217 |
|
V_a(\mathbf r) = \frac{1}{4\pi\epsilon_0} |
| 218 |
|
\sum_{k \, \text{in \bf a}} \frac{q_k}{\lvert \mathbf{r} - \mathbf{r}_k \rvert}. |
| 219 |
|
\end{equation} |
| 220 |
< |
We write the Taylor series expansion in $r$ using an implied summation notation, |
| 221 |
< |
Greek indices are used to indicate space coordinates $x$, $y$, $z$ and the subscripts |
| 222 |
< |
$k$ and $j$ are reserved for labelling specific charges in $\bf a$ and $\bf b$ respectively, and find: |
| 220 |
> |
We write the Taylor series expansion in $r$ using an implied summation |
| 221 |
> |
notation, Greek indices are used to indicate space coordinates ($x$, |
| 222 |
> |
$y$, $z$) and the subscripts $k$ and $j$ are reserved for labelling |
| 223 |
> |
specific charges in $\bf a$ and $\bf b$ respectively, and find: |
| 224 |
|
\begin{equation} |
| 225 |
|
\frac{1}{\lvert \mathbf{r} - \mathbf{r}_k \rvert} = |
| 226 |
|
\left( 1 |
| 229 |
|
\right) |
| 230 |
|
\frac{1}{r} . |
| 231 |
|
\end{equation} |
| 232 |
< |
We then follow Smith in defining an operator for the expansion: |
| 232 |
> |
The electrostatic potential on $\bf a$ can then be expressed in terms |
| 233 |
> |
of multipole operators, |
| 234 |
|
\begin{equation} |
| 235 |
|
V_{\bf a}(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\hat{M}_{\bf a} \frac{1}{r} |
| 236 |
|
\end{equation} |
| 240 |
|
+ Q_{{\bf a}\alpha\beta} |
| 241 |
|
\frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots |
| 242 |
|
\end{equation} |
| 243 |
< |
and the charge $C_{\bf a}$, dipole $D_{{\bf a}\alpha}$, |
| 244 |
< |
and quadrupole $Q_{{\bf a}\alpha\beta}$ are defined by |
| 243 |
> |
Here, the point charge, dipole, and quadrupole for object $\bf a$ are |
| 244 |
> |
given by $C_{\bf a}$, $D_{{\bf a}\alpha}$, and $Q_{{\bf |
| 245 |
> |
a}\alpha\beta}$, respectively. These can be tied to distribution |
| 246 |
> |
of charges |
| 247 |
|
\begin{equation} |
| 248 |
|
C_{\bf a}=\sum_{k \, \text{in \bf a}} q_k , |
| 249 |
|
\end{equation} |
| 492 |
|
We find the first form to be useful in writing equations prior to converting to computer code. The second form is helpful |
| 493 |
|
in derivations of the interaction energy expressions. The third one is specifically helpful when deriving forces and torques, as will |
| 494 |
|
be discussed below. |
| 495 |
+ |
|
| 496 |
+ |
|
| 497 |
+ |
\section{The Self-Interaction} |
| 498 |
+ |
The Wolf summation~\cite{Wolf99} and the later damped shifted force |
| 499 |
+ |
(DSF) extension~\cite{Fennell06} included self-interactions that are |
| 500 |
+ |
handled separately from the pairwise interactions between sites. The |
| 501 |
+ |
self-term is normally calculated via a single loop over all sites in |
| 502 |
+ |
the system, and is relatively cheap to evaluate. The self-interaction |
| 503 |
+ |
has contributions from two sources: |
| 504 |
+ |
\begin{itemize} |
| 505 |
+ |
\item The neutralization procedure within the cutoff radius requires a |
| 506 |
+ |
contribution from a charge opposite in sign, but equal in magnitude, |
| 507 |
+ |
to the central charge, which has been spread out over the surface of |
| 508 |
+ |
the cutoff sphere. For a system of undamped charges, the total |
| 509 |
+ |
self-term is |
| 510 |
+ |
\begin{equation} |
| 511 |
+ |
V_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2} |
| 512 |
+ |
\label{eq:selfTerm} |
| 513 |
+ |
\end{equation} |
| 514 |
+ |
Note that in this potential and in all electrostatic quantities that |
| 515 |
+ |
follow, the standard $4 \pi \epsilon_{0}$ has been omitted for |
| 516 |
+ |
clarity. |
| 517 |
+ |
\item Charge damping with the complementary error function is a |
| 518 |
+ |
partial analogy to the Ewald procedure which splits the interaction |
| 519 |
+ |
into real and reciprocal space sums. The real space sum is retained |
| 520 |
+ |
in the Wolf and DSF methods. The reciprocal space sum is first |
| 521 |
+ |
minimized by folding the largest contribution (the self-interaction) |
| 522 |
+ |
into the self-interaction from charge neutralization of the damped |
| 523 |
+ |
potential. The remainder of the reciprocal space portion is then |
| 524 |
+ |
discarded (as this contributes the largest computational cost and |
| 525 |
+ |
complexity to the Ewald sum). For a system containing only damped |
| 526 |
+ |
charges, the complete self-interaction can be written as |
| 527 |
+ |
\begin{equation} |
| 528 |
+ |
V_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} + |
| 529 |
+ |
\frac{\alpha}{\sqrt{\pi}} \right) \sum_{{\bf a}=1}^N |
| 530 |
+ |
C_{\bf a}^{2}. |
| 531 |
+ |
\label{eq:dampSelfTerm} |
| 532 |
+ |
\end{equation} |
| 533 |
+ |
\end{itemize} |
| 534 |
|
|
| 535 |
+ |
The extension of DSF electrostatics to point multipoles requires |
| 536 |
+ |
treatment of {\it both} the self-neutralization and reciprocal |
| 537 |
+ |
contributions to the self-interaction for higher order multipoles. In |
| 538 |
+ |
this section we give formulae for these interactions up to quadrupolar |
| 539 |
+ |
order. |
| 540 |
+ |
|
| 541 |
+ |
The self-neutralization term is computed by taking the {\it |
| 542 |
+ |
non-shifted} kernel for each interaction, placing a multipole of |
| 543 |
+ |
equal magnitude (but opposite in polarization) on the surface of the |
| 544 |
+ |
cutoff sphere, and averaging over the surface of the cutoff sphere. |
| 545 |
+ |
Because the self term is carried out as a single sum over sites, the |
| 546 |
+ |
reciprocal-space portion is identical to half of the self-term |
| 547 |
+ |
obtained by Smith and Aguado and Madden for the application of the |
| 548 |
+ |
Ewald sum to multipoles.\cite{Smith82,Smith98,Aguado03} For a given |
| 549 |
+ |
site which posesses a charge, dipole, and multipole, both types of |
| 550 |
+ |
contribution are given in table \ref{tab:tableSelf}. |
| 551 |
+ |
|
| 552 |
+ |
\begin{table*} |
| 553 |
+ |
\caption{\label{tab:tableSelf} Self-interaction contributions for |
| 554 |
+ |
site ({\bf a}) that has a charge $(C_{\bf a})$, dipole |
| 555 |
+ |
$(\mathbf{D}_{\bf a})$, and quadrupole $(\mathbf{Q}_{\bf a})$} |
| 556 |
+ |
\begin{ruledtabular} |
| 557 |
+ |
\begin{tabular}{lccc} |
| 558 |
+ |
Multipole order & Summed Quantity & Self-neutralization & Reciprocal \\ \hline |
| 559 |
+ |
Charge & $C_{\bf a}^2$ & $-f(r_c)$ & $-\frac{\alpha}{\sqrt{\pi}}$ \\ |
| 560 |
+ |
Dipole & $|\mathbf{D}_{\bf a}|^2$ & $\frac{1}{3} \left( h(r_c) + |
| 561 |
+ |
\frac{2 g(r_c)}{r_c} \right)$ & $-\frac{2 \alpha^3}{3 \sqrt{\pi}}$\\ |
| 562 |
+ |
Quadrupole & $2 \text{Tr}(\mathbf{Q}_{\bf a}^2) + \text{Tr}(\mathbf{Q}_{\bf a})^2$ & |
| 563 |
+ |
$- \frac{1}{15} \left( t(r_c)+ \frac{4 s(r_c)}{r_c} \right)$ & |
| 564 |
+ |
$-\frac{4 \alpha^5}{5 \sqrt{\pi}}$ \\ |
| 565 |
+ |
Charge-Quadrupole & $-2 C_{\bf a} \text{Tr}(\mathbf{Q}_{\bf a})$ & $\frac{1}{3} \left( |
| 566 |
+ |
h(r_c) + \frac{2 g(r_c)}{r_c} \right)$& $-\frac{2 \alpha^3}{3 \sqrt{\pi}}$ \\ |
| 567 |
+ |
\end{tabular} |
| 568 |
+ |
\end{ruledtabular} |
| 569 |
+ |
\end{table*} |
| 570 |
+ |
|
| 571 |
+ |
For sites which simultaneously contain charges and quadrupoles, the |
| 572 |
+ |
self-interaction includes a cross-interaction between these two |
| 573 |
+ |
multipole orders. Symmetry prevents the charge-dipole and |
| 574 |
+ |
dipole-quadrupole interactions from contributing to the |
| 575 |
+ |
self-interaction. The functions that go into the self-neutralization |
| 576 |
+ |
terms, $f(r), g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive |
| 577 |
+ |
derivatives of the electrostatic kernel (either the undamped $1/r$ or |
| 578 |
+ |
the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that are |
| 579 |
+ |
evaluated at the cutoff distance. For undamped interactions, $f(r_c) |
| 580 |
+ |
= 1/r_c$, $g(r_c) = -1/r_c^{2}$, and so on. For damped interactions, |
| 581 |
+ |
$f(r_c) = B_0(r_c)$, $g(r_c) = B_0'(r_c)$, and so on. Appendix XX |
| 582 |
+ |
contains recursion relations that allow rapid evaluation of these |
| 583 |
+ |
derivatives. |
| 584 |
+ |
|
| 585 |
|
\section{Energies, forces, and torques} |
| 586 |
|
\subsection{Interaction energies} |
| 587 |
|
|
| 1261 |
|
\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r) |
| 1262 |
|
\end{split} |
| 1263 |
|
\end{equation} |
| 1086 |
– |
% |
| 1087 |
– |
% |
| 1264 |
|
% |
| 1265 |
+ |
|
| 1266 |
+ |
\section{Comparison to known multipolar energies} |
| 1267 |
+ |
|
| 1268 |
+ |
To understand how these new real-space multipole methods behave in |
| 1269 |
+ |
computer simulations, it is vital to test against established methods |
| 1270 |
+ |
for computing electrostatic interactions in periodic systems, and to |
| 1271 |
+ |
evaluate the size and sources of any errors that arise from the |
| 1272 |
+ |
real-space cutoffs. In this paper we test Taylor-shifted and |
| 1273 |
+ |
Gradient-shifted electrostatics against analytical methods for |
| 1274 |
+ |
computing the energies of ordered multipolar arrays. In the following |
| 1275 |
+ |
paper, we test the new methods against the multipolar Ewald sum for |
| 1276 |
+ |
computing the energies, forces and torques for a wide range of typical |
| 1277 |
+ |
condensed-phase (disordered) systems. |
| 1278 |
+ |
|
| 1279 |
+ |
Because long-range electrostatic effects can be significant in |
| 1280 |
+ |
crystalline materials, ordered multipolar arrays present one of the |
| 1281 |
+ |
biggest challenges for real-space cutoff methods. The dipolar |
| 1282 |
+ |
analogues to the Madelung constants were first worked out by Sauer, |
| 1283 |
+ |
who computed the energies of ordered dipole arrays of zero |
| 1284 |
+ |
magnetization and obtained a number of these constants.\cite{Sauer} |
| 1285 |
+ |
This theory was developed more completely by Luttinger and |
| 1286 |
+ |
Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays and |
| 1287 |
+ |
other periodic structures. We have repeated the Luttinger \& Tisza |
| 1288 |
+ |
series summations to much higher order and obtained the following |
| 1289 |
+ |
energy constants (converged to one part in $10^9$): |
| 1290 |
+ |
\begin{table*} |
| 1291 |
+ |
\centering{ |
| 1292 |
+ |
\caption{Luttinger \& Tisza arrays and their associated |
| 1293 |
+ |
energy constants. Type "A" arrays have nearest neighbor strings of |
| 1294 |
+ |
antiparallel dipoles. Type "B" arrays have nearest neighbor |
| 1295 |
+ |
strings of antiparallel dipoles if the dipoles are contained in a |
| 1296 |
+ |
plane perpendicular to the dipole direction that passes through |
| 1297 |
+ |
the dipole.} |
| 1298 |
+ |
} |
| 1299 |
+ |
\label{tab:LT} |
| 1300 |
+ |
\begin{ruledtabular} |
| 1301 |
+ |
\begin{tabular}{cccc} |
| 1302 |
+ |
Array Type & Lattice & Dipole Direction & Energy constants \\ \hline |
| 1303 |
+ |
A & SC & 001 & -2.676788684 \\ |
| 1304 |
+ |
A & BCC & 001 & 0 \\ |
| 1305 |
+ |
A & BCC & 111 & -1.770078733 \\ |
| 1306 |
+ |
A & FCC & 001 & 2.166932835 \\ |
| 1307 |
+ |
A & FCC & 011 & -1.083466417 \\ |
| 1308 |
+ |
|
| 1309 |
+ |
* & BCC & minimum & -1.985920929 \\ |
| 1310 |
+ |
|
| 1311 |
+ |
B & SC & 001 & -2.676788684 \\ |
| 1312 |
+ |
B & BCC & 001 & -1.338394342 \\ |
| 1313 |
+ |
B & BCC & 111 & -1.770078733 \\ |
| 1314 |
+ |
B & FCC & 001 & -1.083466417 \\ |
| 1315 |
+ |
B & FCC & 011 & -1.807573634 |
| 1316 |
+ |
\end{tabular} |
| 1317 |
+ |
\end{ruledtabular} |
| 1318 |
+ |
\end{table*} |
| 1319 |
+ |
|
| 1320 |
+ |
In addition to the A and B arrays, there is an additional minimum |
| 1321 |
+ |
energy structure for the BCC lattice that was found by Luttinger \& |
| 1322 |
+ |
Tisza. The total electrostatic energy for an array is given by: |
| 1323 |
+ |
\begin{equation} |
| 1324 |
+ |
E = C N^2 \mu^2 |
| 1325 |
+ |
\end{equation} |
| 1326 |
+ |
where $C$ is the energy constant given above, $N$ is the number of |
| 1327 |
+ |
dipoles per unit volume, and $\mu$ is the strength of the dipole. |
| 1328 |
+ |
|
| 1329 |
+ |
{\it Quadrupolar} analogues to the Madelung constants were first worked out by Nagai and Nakamura who |
| 1330 |
+ |
computed the energies of selected quadrupole arrays based on |
| 1331 |
+ |
extensions to the Luttinger and Tisza |
| 1332 |
+ |
approach.\cite{Nagai01081960,Nagai01091963} We have compared the |
| 1333 |
+ |
energy constants for the lowest energy configurations for linear |
| 1334 |
+ |
quadrupoles shown in table \ref{tab:NNQ} |
| 1335 |
+ |
|
| 1336 |
+ |
\begin{table*} |
| 1337 |
+ |
\centering{ |
| 1338 |
+ |
\caption{Nagai and Nakamura Quadurpolar arrays}} |
| 1339 |
+ |
\label{tab:NNQ} |
| 1340 |
+ |
\begin{ruledtabular} |
| 1341 |
+ |
\begin{tabular}{ccc} |
| 1342 |
+ |
Lattice & Quadrupole Direction & Energy constants \\ \hline |
| 1343 |
+ |
SC & 111 & -8.3 \\ |
| 1344 |
+ |
BCC & 011 & -21.7 \\ |
| 1345 |
+ |
FCC & 111 & -80.5 |
| 1346 |
+ |
\end{tabular} |
| 1347 |
+ |
\end{ruledtabular} |
| 1348 |
+ |
\end{table*} |
| 1349 |
+ |
|
| 1350 |
+ |
In analogy to the dipolar arrays, the total electrostatic energy for |
| 1351 |
+ |
the quadrupolar arrays is: |
| 1352 |
+ |
\begin{equation} |
| 1353 |
+ |
E = C \frac{3}{4} N^2 Q^2 |
| 1354 |
+ |
\end{equation} |
| 1355 |
+ |
where $Q$ is the quadrupole moment. |
| 1356 |
+ |
|
| 1357 |
+ |
|
| 1358 |
+ |
|
| 1359 |
+ |
|
| 1360 |
+ |
|
| 1361 |
+ |
|
| 1362 |
+ |
|
| 1363 |
+ |
|
| 1364 |
|
\begin{acknowledgments} |
| 1365 |
< |
We wish to acknowledge the support of the author community in using |
| 1366 |
< |
REV\TeX{}, offering suggestions and encouragement, testing new versions, |
| 1367 |
< |
\dots. |
| 1365 |
> |
Support for this project was provided by the National Science |
| 1366 |
> |
Foundation under grant CHE-0848243. Computational time was provided by |
| 1367 |
> |
the Center for Research Computing (CRC) at the University of Notre |
| 1368 |
> |
Dame. |
| 1369 |
|
\end{acknowledgments} |
| 1370 |
|
|
| 1371 |
|
\appendix |
| 2217 |
|
\end{tabular} |
| 2218 |
|
\end{ruledtabular} |
| 2219 |
|
\end{table*} |
| 2220 |
+ |
|
| 2221 |
+ |
\newpage |
| 2222 |
+ |
|
| 2223 |
+ |
\bibliography{multipole} |
| 2224 |
+ |
|
| 2225 |
|
\end{document} |
| 2226 |
|
% |
| 2227 |
|
% ****** End of file multipole.tex ****** |
