1 |
gezelter |
3906 |
% ****** Start of file aipsamp.tex ****** |
2 |
|
|
% |
3 |
|
|
% This file is part of the AIP files in the AIP distribution for REVTeX 4. |
4 |
|
|
% Version 4.1 of REVTeX, October 2009 |
5 |
|
|
% |
6 |
|
|
% Copyright (c) 2009 American Institute of Physics. |
7 |
|
|
% |
8 |
|
|
% See the AIP README file for restrictions and more information. |
9 |
|
|
% |
10 |
|
|
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed |
11 |
|
|
% as well as the rest of the prerequisites for REVTeX 4.1 |
12 |
|
|
% |
13 |
|
|
% It also requires running BibTeX. The commands are as follows: |
14 |
|
|
% |
15 |
|
|
% 1) latex aipsamp |
16 |
|
|
% 2) bibtex aipsamp |
17 |
|
|
% 3) latex aipsamp |
18 |
|
|
% 4) latex aipsamp |
19 |
|
|
% |
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 |
gezelter |
3980 |
aip,jcp, |
24 |
gezelter |
3906 |
amsmath,amssymb, |
25 |
|
|
preprint,% |
26 |
gezelter |
4172 |
% reprint,% |
27 |
gezelter |
3906 |
%author-year,% |
28 |
|
|
%author-numerical,% |
29 |
gezelter |
3984 |
jcp]{revtex4-1} |
30 |
gezelter |
3906 |
|
31 |
|
|
\usepackage{graphicx}% Include figure files |
32 |
|
|
\usepackage{dcolumn}% Align table columns on decimal point |
33 |
gezelter |
3980 |
%\usepackage{bm}% bold math |
34 |
gezelter |
3982 |
\usepackage{times} |
35 |
gezelter |
3980 |
\usepackage[version=3]{mhchem} % this is a great package for formatting chemical reactions |
36 |
|
|
\usepackage{url} |
37 |
gezelter |
3985 |
\usepackage{rotating} |
38 |
gezelter |
4172 |
|
39 |
gezelter |
3906 |
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math |
40 |
|
|
%\linenumbers\relax % Commence numbering lines |
41 |
|
|
|
42 |
|
|
\begin{document} |
43 |
|
|
|
44 |
gezelter |
3988 |
%\preprint{AIP/123-QED} |
45 |
gezelter |
3906 |
|
46 |
gezelter |
3988 |
\title{Real space alternatives to the Ewald |
47 |
gezelter |
4174 |
Sum. I. Shifted electrostatics for multipoles} |
48 |
gezelter |
3906 |
|
49 |
|
|
\author{Madan Lamichhane} |
50 |
|
|
\affiliation{Department of Physics, University |
51 |
|
|
of Notre Dame, Notre Dame, IN 46556} |
52 |
|
|
|
53 |
|
|
\author{J. Daniel Gezelter} |
54 |
|
|
\email{gezelter@nd.edu.} |
55 |
|
|
\affiliation{Department of Chemistry and Biochemistry, University |
56 |
|
|
of Notre Dame, Notre Dame, IN 46556} |
57 |
|
|
|
58 |
|
|
\author{Kathie E. Newman} |
59 |
|
|
\affiliation{Department of Physics, University |
60 |
|
|
of Notre Dame, Notre Dame, IN 46556} |
61 |
|
|
|
62 |
|
|
|
63 |
|
|
\date{\today}% It is always \today, today, |
64 |
|
|
% but any date may be explicitly specified |
65 |
|
|
|
66 |
|
|
\begin{abstract} |
67 |
gezelter |
3980 |
We have extended the original damped-shifted force (DSF) |
68 |
gezelter |
4179 |
electrostatic kernel and have been able to derive three new |
69 |
gezelter |
3980 |
electrostatic potentials for higher-order multipoles that are based |
70 |
gezelter |
4179 |
on truncated Taylor expansions around the cutoff radius. These |
71 |
|
|
include a shifted potential (SP) that generalizes the Wolf method |
72 |
|
|
for point multipoles, and Taylor-shifted force (TSF) and |
73 |
|
|
gradient-shifted force (GSF) potentials that are both |
74 |
|
|
generalizations of DSF electrostatics for multipoles. We find that |
75 |
|
|
each of the distinct orientational contributions requires a separate |
76 |
|
|
radial function to ensure that pairwise energies, forces and torques |
77 |
|
|
all vanish at the cutoff radius. In this paper, we present energy, |
78 |
|
|
force, and torque expressions for the new models, and compare these |
79 |
|
|
real-space interaction models to exact results for ordered arrays of |
80 |
|
|
multipoles. We find that the GSF and SP methods converge rapidly to |
81 |
|
|
the correct lattice energies for ordered dipolar and quadrupolar |
82 |
|
|
arrays, while the Taylor-Shifted Force (TSF) is too severe an |
83 |
|
|
approximation to provide accurate convergence to lattice energies. |
84 |
|
|
Because real-space methods can be made to scale linearly with system |
85 |
gezelter |
4183 |
size, SP and GSF are attractive options for large Monte |
86 |
|
|
Carlo and molecular dynamics simulations, respectively. |
87 |
gezelter |
3906 |
\end{abstract} |
88 |
|
|
|
89 |
gezelter |
3988 |
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
90 |
gezelter |
3906 |
% Classification Scheme. |
91 |
gezelter |
3988 |
%\keywords{Suggested keywords}%Use showkeys class option if keyword |
92 |
gezelter |
3906 |
%display desired |
93 |
|
|
\maketitle |
94 |
|
|
|
95 |
|
|
\section{Introduction} |
96 |
gezelter |
3982 |
There has been increasing interest in real-space methods for |
97 |
|
|
calculating electrostatic interactions in computer simulations of |
98 |
|
|
condensed molecular |
99 |
gezelter |
3980 |
systems.\cite{Wolf99,Zahn02,Kast03,BeckD.A.C._bi0486381,Ma05,Fennell:2006zl,Chen:2004du,Chen:2006ii,Rodgers:2006nw,Denesyuk:2008ez,Izvekov:2008wo} |
100 |
|
|
The simplest of these techniques was developed by Wolf {\it et al.} |
101 |
|
|
in their work towards an $\mathcal{O}(N)$ Coulombic sum.\cite{Wolf99} |
102 |
gezelter |
3982 |
For systems of point charges, Fennell and Gezelter showed that a |
103 |
|
|
simple damped shifted force (DSF) modification to Wolf's method could |
104 |
|
|
give nearly quantitative agreement with smooth particle mesh Ewald |
105 |
|
|
(SPME)\cite{Essmann95} configurational energy differences as well as |
106 |
|
|
atomic force and molecular torque vectors.\cite{Fennell:2006zl} |
107 |
gezelter |
3906 |
|
108 |
gezelter |
3980 |
The computational efficiency and the accuracy of the DSF method are |
109 |
|
|
surprisingly good, particularly for systems with uniform charge |
110 |
|
|
density. Additionally, dielectric constants obtained using DSF and |
111 |
gezelter |
3986 |
similar methods where the force vanishes at $r_{c}$ are |
112 |
gezelter |
3980 |
essentially quantitative.\cite{Izvekov:2008wo} The DSF and other |
113 |
|
|
related methods have now been widely investigated,\cite{Hansen:2012uq} |
114 |
gezelter |
3985 |
and DSF is now used routinely in a diverse set of chemical |
115 |
|
|
environments.\cite{doi:10.1021/la400226g,McCann:2013fk,kannam:094701,Forrest:2012ly,English:2008kx,Louden:2013ve,Tokumasu:2013zr} |
116 |
|
|
DSF electrostatics provides a compromise between the computational |
117 |
|
|
speed of real-space cutoffs and the accuracy of fully-periodic Ewald |
118 |
|
|
treatments. |
119 |
gezelter |
3906 |
|
120 |
gezelter |
3980 |
One common feature of many coarse-graining approaches, which treat |
121 |
|
|
entire molecular subsystems as a single rigid body, is simplification |
122 |
|
|
of the electrostatic interactions between these bodies so that fewer |
123 |
|
|
site-site interactions are required to compute configurational |
124 |
gezelter |
3986 |
energies. To do this, the interactions between coarse-grained sites |
125 |
|
|
are typically taken to be point |
126 |
|
|
multipoles.\cite{Golubkov06,ISI:000276097500009,ISI:000298664400012} |
127 |
gezelter |
3906 |
|
128 |
gezelter |
3986 |
Water, in particular, has been modeled recently with point multipoles |
129 |
|
|
up to octupolar |
130 |
|
|
order.\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} For |
131 |
|
|
maximum efficiency, these models require the use of an approximate |
132 |
|
|
multipole expansion as the exact multipole expansion can become quite |
133 |
|
|
expensive (particularly when handled via the Ewald |
134 |
|
|
sum).\cite{Ichiye:2006qy} Point multipoles and multipole |
135 |
|
|
polarizability have also been utilized in the AMOEBA water model and |
136 |
gezelter |
3980 |
related force fields.\cite{Ponder:2010fk,schnieders:124114,Ren:2011uq} |
137 |
gezelter |
3906 |
|
138 |
gezelter |
3980 |
Higher-order multipoles present a peculiar issue for molecular |
139 |
|
|
dynamics. Multipolar interactions are inherently short-ranged, and |
140 |
|
|
should not need the relatively expensive Ewald treatment. However, |
141 |
|
|
real-space cutoff methods are normally applied in an orientation-blind |
142 |
|
|
fashion so multipoles which leave and then re-enter a cutoff sphere in |
143 |
|
|
a different orientation can cause energy discontinuities. |
144 |
gezelter |
3906 |
|
145 |
gezelter |
3980 |
This paper outlines an extension of the original DSF electrostatic |
146 |
gezelter |
4183 |
kernel to point multipoles. We describe three distinct real-space |
147 |
|
|
interaction models for higher-order multipoles based on truncated |
148 |
gezelter |
3982 |
Taylor expansions that are carried out at the cutoff radius. We are |
149 |
gezelter |
4183 |
calling these models {\bf Taylor-shifted} (TSF), {\bf |
150 |
|
|
gradient-shifted} (GSF) and {\bf shifted potential} (SP) |
151 |
gezelter |
3982 |
electrostatics. Because of differences in the initial assumptions, |
152 |
gezelter |
4183 |
the two methods yield related, but distinct expressions for energies, |
153 |
|
|
forces, and torques. |
154 |
gezelter |
3906 |
|
155 |
gezelter |
3982 |
In this paper we outline the new methodology and give functional forms |
156 |
|
|
for the energies, forces, and torques up to quadrupole-quadrupole |
157 |
|
|
order. We also compare the new methods to analytic energy constants |
158 |
gezelter |
3986 |
for periodic arrays of point multipoles. In the following paper, we |
159 |
gezelter |
3985 |
provide numerical comparisons to Ewald-based electrostatics in common |
160 |
|
|
simulation enviornments. |
161 |
gezelter |
3982 |
|
162 |
gezelter |
3980 |
\section{Methodology} |
163 |
gezelter |
3986 |
An efficient real-space electrostatic method involves the use of a |
164 |
|
|
pair-wise functional form, |
165 |
|
|
\begin{equation} |
166 |
gezelter |
4175 |
U = \sum_i \sum_{j>i} U_\mathrm{pair}(\mathbf{r}_{ij}, \Omega_i, \Omega_j) + |
167 |
|
|
\sum_i U_i^\mathrm{self} |
168 |
gezelter |
3986 |
\end{equation} |
169 |
|
|
that is short-ranged and easily truncated at a cutoff radius, |
170 |
|
|
\begin{equation} |
171 |
gezelter |
4175 |
U_\mathrm{pair}(\mathbf{r}_{ij},\Omega_i, \Omega_j) = \left\{ |
172 |
gezelter |
3986 |
\begin{array}{ll} |
173 |
gezelter |
4175 |
U_\mathrm{approx} (\mathbf{r}_{ij}, \Omega_i, \Omega_j) & \quad \left| \mathbf{r}_{ij} \right| \le r_c \\ |
174 |
gezelter |
4098 |
0 & \quad \left| \mathbf{r}_{ij} \right| > r_c , |
175 |
gezelter |
3986 |
\end{array} |
176 |
|
|
\right. |
177 |
|
|
\end{equation} |
178 |
gezelter |
3989 |
along with an easily computed self-interaction term ($\sum_i |
179 |
gezelter |
4175 |
U_i^\mathrm{self}$) which scales linearly with the number of |
180 |
gezelter |
3986 |
particles. Here $\Omega_i$ and $\Omega_j$ represent orientational |
181 |
gezelter |
4098 |
coordinates of the two sites, and $\mathbf{r}_{ij}$ is the vector |
182 |
|
|
between the two sites. The computational efficiency, energy |
183 |
gezelter |
3986 |
conservation, and even some physical properties of a simulation can |
184 |
gezelter |
4175 |
depend dramatically on how the $U_\mathrm{approx}$ function behaves at |
185 |
gezelter |
3986 |
the cutoff radius. The goal of any approximation method should be to |
186 |
|
|
mimic the real behavior of the electrostatic interactions as closely |
187 |
|
|
as possible without sacrificing the near-linear scaling of a cutoff |
188 |
|
|
method. |
189 |
gezelter |
3906 |
|
190 |
gezelter |
3980 |
\subsection{Self-neutralization, damping, and force-shifting} |
191 |
|
|
The DSF and Wolf methods operate by neutralizing the total charge |
192 |
|
|
contained within the cutoff sphere surrounding each particle. This is |
193 |
|
|
accomplished by shifting the potential functions to generate image |
194 |
|
|
charges on the surface of the cutoff sphere for each pair interaction |
195 |
gezelter |
4098 |
computed within $r_c$. Damping using a complementary error function is |
196 |
|
|
applied to the potential to accelerate convergence. The interaction |
197 |
|
|
for a pair of charges ($C_i$ and $C_j$) in the DSF method, |
198 |
gezelter |
3980 |
\begin{equation*} |
199 |
gezelter |
4175 |
U_\mathrm{DSF}(r_{ij}) = C_i C_j \Biggr{[} |
200 |
gezelter |
3980 |
\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} |
201 |
gezelter |
3986 |
- \frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c} + \left(\frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c^2} |
202 |
gezelter |
3980 |
+ \frac{2\alpha}{\pi^{1/2}} |
203 |
gezelter |
3986 |
\frac{\exp\left(-\alpha^2r_c^2\right)}{r_c} |
204 |
|
|
\right)\left(r_{ij}-r_c\right)\ \Biggr{]} |
205 |
gezelter |
3980 |
\label{eq:DSFPot} |
206 |
|
|
\end{equation*} |
207 |
gezelter |
4098 |
where $\alpha$ is the adjustable damping parameter. Note that in this |
208 |
|
|
potential and in all electrostatic quantities that follow, the |
209 |
|
|
standard $1/4 \pi \epsilon_{0}$ has been omitted for clarity. |
210 |
gezelter |
3980 |
|
211 |
|
|
To insure net charge neutrality within each cutoff sphere, an |
212 |
|
|
additional ``self'' term is added to the potential. This term is |
213 |
|
|
constant (as long as the charges and cutoff radius do not change), and |
214 |
|
|
exists outside the normal pair-loop for molecular simulations. It can |
215 |
|
|
be thought of as a contribution from a charge opposite in sign, but |
216 |
|
|
equal in magnitude, to the central charge, which has been spread out |
217 |
|
|
over the surface of the cutoff sphere. A portion of the self term is |
218 |
|
|
identical to the self term in the Ewald summation, and comes from the |
219 |
|
|
utilization of the complimentary error function for electrostatic |
220 |
gezelter |
3986 |
damping.\cite{deLeeuw80,Wolf99} There have also been recent efforts to |
221 |
|
|
extend the Wolf self-neutralization method to zero out the dipole and |
222 |
|
|
higher order multipoles contained within the cutoff |
223 |
gezelter |
3985 |
sphere.\cite{Fukuda:2011jk,Fukuda:2012yu,Fukuda:2013qv} |
224 |
gezelter |
3982 |
|
225 |
gezelter |
3985 |
In this work, we extend the idea of self-neutralization for the point |
226 |
|
|
multipoles by insuring net charge-neutrality and net-zero moments |
227 |
|
|
within each cutoff sphere. In Figure \ref{fig:shiftedMultipoles}, the |
228 |
|
|
central dipolar site $\mathbf{D}_i$ is interacting with point dipole |
229 |
|
|
$\mathbf{D}_j$ and point quadrupole, $\mathbf{Q}_k$. The |
230 |
|
|
self-neutralization scheme for point multipoles involves projecting |
231 |
|
|
opposing multipoles for sites $j$ and $k$ on the surface of the cutoff |
232 |
|
|
sphere. There are also significant modifications made to make the |
233 |
|
|
forces and torques go smoothly to zero at the cutoff distance. |
234 |
gezelter |
3982 |
|
235 |
gezelter |
4172 |
\begin{figure} |
236 |
|
|
\includegraphics[width=3in]{SM} |
237 |
|
|
\caption{Reversed multipoles are projected onto the surface of the |
238 |
|
|
cutoff sphere. The forces, torques, and potential are then smoothly |
239 |
|
|
shifted to zero as the sites leave the cutoff region.} |
240 |
|
|
\label{fig:shiftedMultipoles} |
241 |
|
|
\end{figure} |
242 |
gezelter |
3980 |
|
243 |
gezelter |
3986 |
As in the point-charge approach, there is an additional contribution |
244 |
|
|
from self-neutralization of site $i$. The self term for multipoles is |
245 |
gezelter |
3982 |
described in section \ref{sec:selfTerm}. |
246 |
gezelter |
3906 |
|
247 |
gezelter |
3982 |
\subsection{The multipole expansion} |
248 |
|
|
|
249 |
gezelter |
3980 |
Consider two discrete rigid collections of point charges, denoted as |
250 |
gezelter |
3982 |
$\bf a$ and $\bf b$. In the following, we assume that the two objects |
251 |
|
|
interact via electrostatics only and describe those interactions in |
252 |
|
|
terms of a standard multipole expansion. Putting the origin of the |
253 |
|
|
coordinate system at the center of mass of $\bf a$, we use vectors |
254 |
gezelter |
3980 |
$\mathbf{r}_k$ to denote the positions of all charges $q_k$ in $\bf |
255 |
|
|
a$. Then the electrostatic potential of object $\bf a$ at |
256 |
|
|
$\mathbf{r}$ is given by |
257 |
gezelter |
3906 |
\begin{equation} |
258 |
gezelter |
4175 |
\phi_a(\mathbf r) = |
259 |
gezelter |
3906 |
\sum_{k \, \text{in \bf a}} \frac{q_k}{\lvert \mathbf{r} - \mathbf{r}_k \rvert}. |
260 |
|
|
\end{equation} |
261 |
gezelter |
3982 |
The Taylor expansion in $r$ can be written using an implied summation |
262 |
|
|
notation. Here Greek indices are used to indicate space coordinates |
263 |
|
|
($x$, $y$, $z$) and the subscripts $k$ and $j$ are reserved for |
264 |
gezelter |
4098 |
labeling specific charges in $\bf a$ and $\bf b$ respectively. The |
265 |
gezelter |
3982 |
Taylor expansion, |
266 |
gezelter |
3906 |
\begin{equation} |
267 |
|
|
\frac{1}{\lvert \mathbf{r} - \mathbf{r}_k \rvert} = |
268 |
|
|
\left( 1 |
269 |
|
|
- r_{k\alpha} \frac{\partial}{\partial r_{\alpha}} |
270 |
|
|
+ \frac{1}{2} r_{k\alpha} r_{k\beta} \frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} +\dots |
271 |
|
|
\right) |
272 |
gezelter |
3982 |
\frac{1}{r} , |
273 |
gezelter |
3906 |
\end{equation} |
274 |
gezelter |
3982 |
can then be used to express the electrostatic potential on $\bf a$ in |
275 |
|
|
terms of multipole operators, |
276 |
gezelter |
3906 |
\begin{equation} |
277 |
gezelter |
4175 |
\phi_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r} |
278 |
gezelter |
3906 |
\end{equation} |
279 |
|
|
where |
280 |
|
|
\begin{equation} |
281 |
|
|
\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}} |
282 |
|
|
+ Q_{{\bf a}\alpha\beta} |
283 |
|
|
\frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots |
284 |
|
|
\end{equation} |
285 |
gezelter |
3980 |
Here, the point charge, dipole, and quadrupole for object $\bf a$ are |
286 |
|
|
given by $C_{\bf a}$, $D_{{\bf a}\alpha}$, and $Q_{{\bf |
287 |
gezelter |
3982 |
a}\alpha\beta}$, respectively. These are the primitive multipoles |
288 |
|
|
which can be expressed as a distribution of charges, |
289 |
|
|
\begin{align} |
290 |
gezelter |
3996 |
C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \label{eq:charge} \\ |
291 |
|
|
D_{{\bf a}\alpha} =&\sum_{k \, \text{in \bf a}} q_k r_{k\alpha}, \label{eq:dipole}\\ |
292 |
|
|
Q_{{\bf a}\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k |
293 |
|
|
r_{k\alpha} r_{k\beta} . \label{eq:quadrupole} |
294 |
gezelter |
3982 |
\end{align} |
295 |
|
|
Note that the definition of the primitive quadrupole here differs from |
296 |
|
|
the standard traceless form, and contains an additional Taylor-series |
297 |
gezelter |
3996 |
based factor of $1/2$. We are essentially treating the mass |
298 |
|
|
distribution with higher priority; the moment of inertia tensor, |
299 |
|
|
$\overleftrightarrow{\mathsf I}$, is diagonalized to obtain body-fixed |
300 |
|
|
axes, and the charge distribution may result in a quadrupole tensor |
301 |
|
|
that is not necessarily diagonal in the body frame. Additional |
302 |
|
|
reasons for utilizing the primitive quadrupole are discussed in |
303 |
|
|
section \ref{sec:damped}. |
304 |
gezelter |
3906 |
|
305 |
|
|
It is convenient to locate charges $q_j$ relative to the center of mass of $\bf b$. Then with $\bf{r}$ pointing from |
306 |
gezelter |
4179 |
$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $), the interaction energy is given by |
307 |
gezelter |
3906 |
\begin{equation} |
308 |
gezelter |
3982 |
U_{\bf{ab}}(r) |
309 |
gezelter |
3985 |
= \hat{M}_a \sum_{j \, \text{in \bf b}} \frac {q_j}{\vert \bf{r}+\bf{r}_j \vert} . |
310 |
gezelter |
3982 |
\end{equation} |
311 |
|
|
This can also be expanded as a Taylor series in $r$. Using a notation |
312 |
|
|
similar to before to define the multipoles on object {\bf b}, |
313 |
|
|
\begin{equation} |
314 |
gezelter |
3906 |
\hat{M}_{\bf b} = C_{\bf b} + D_{{\bf b}\alpha} \frac{\partial}{\partial r_{\alpha}} |
315 |
|
|
+ Q_{{\bf b}\alpha\beta} |
316 |
|
|
\frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots |
317 |
|
|
\end{equation} |
318 |
gezelter |
3982 |
we arrive at the multipole expression for the total interaction energy. |
319 |
gezelter |
3906 |
\begin{equation} |
320 |
gezelter |
3985 |
U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r} \label{kernel}. |
321 |
gezelter |
3906 |
\end{equation} |
322 |
gezelter |
3982 |
This form has the benefit of separating out the energies of |
323 |
|
|
interaction into contributions from the charge, dipole, and quadrupole |
324 |
gezelter |
4179 |
of $\bf a$ interacting with the same multipoles in $\bf b$. |
325 |
gezelter |
3906 |
|
326 |
gezelter |
3982 |
\subsection{Damped Coulomb interactions} |
327 |
gezelter |
3996 |
\label{sec:damped} |
328 |
gezelter |
3982 |
In the standard multipole expansion, one typically uses the bare |
329 |
|
|
Coulomb potential, with radial dependence $1/r$, as shown in |
330 |
|
|
Eq.~(\ref{kernel}). It is also quite common to use a damped Coulomb |
331 |
|
|
interaction, which results from replacing point charges with Gaussian |
332 |
|
|
distributions of charge with width $\alpha$. In damped multipole |
333 |
|
|
electrostatics, the kernel ($1/r$) of the expansion is replaced with |
334 |
|
|
the function: |
335 |
gezelter |
3906 |
\begin{equation} |
336 |
|
|
B_0(r)=\frac{\text{erfc}(\alpha r)}{r} = \frac{2}{\sqrt{\pi}r} |
337 |
|
|
\int_{\alpha r}^{\infty} \text{e}^{-s^2} ds . |
338 |
|
|
\end{equation} |
339 |
gezelter |
3982 |
We develop equations below using the function $f(r)$ to represent |
340 |
gezelter |
3986 |
either $1/r$ or $B_0(r)$, and all of the techniques can be applied to |
341 |
|
|
bare or damped Coulomb kernels (or any other function) as long as |
342 |
|
|
derivatives of these functions are known. Smith's convenient |
343 |
gezelter |
4179 |
functions $B_l(r)$, which are used for derivatives of the damped |
344 |
|
|
kernel, are summarized in Appendix A. (N.B. there is one important |
345 |
|
|
distinction between the two kernels, which is the behavior of |
346 |
|
|
$\nabla^2 \frac{1}{r}$ compared with $\nabla^2 B_0(r)$. The former is |
347 |
|
|
zero everywhere except for a delta function evaluated at the origin. |
348 |
|
|
The latter also has delta function behavior, but is non-zero for $r |
349 |
|
|
\neq 0$. Thus the standard justification for using a traceless |
350 |
gezelter |
3996 |
quadrupole tensor fails for the damped case.) |
351 |
gezelter |
3906 |
|
352 |
gezelter |
3982 |
The main goal of this work is to smoothly cut off the interaction |
353 |
|
|
energy as well as forces and torques as $r\rightarrow r_c$. To |
354 |
|
|
describe how this goal may be met, we use two examples, charge-charge |
355 |
gezelter |
3986 |
and charge-dipole, using the bare Coulomb kernel, $f(r)=1/r$, to |
356 |
|
|
explain the idea. |
357 |
gezelter |
3906 |
|
358 |
gezelter |
3984 |
\subsection{Shifted-force methods} |
359 |
gezelter |
3982 |
In the shifted-force approximation, the interaction energy for two |
360 |
|
|
charges $C_{\bf a}$ and $C_{\bf b}$ separated by a distance $r$ is |
361 |
|
|
written: |
362 |
gezelter |
3906 |
\begin{equation} |
363 |
gezelter |
3985 |
U_{C_{\bf a}C_{\bf b}}(r)= C_{\bf a} C_{\bf b} |
364 |
gezelter |
3906 |
\left({ \frac{1}{r} - \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} } |
365 |
|
|
\right) . |
366 |
|
|
\end{equation} |
367 |
gezelter |
3982 |
Two shifting terms appear in this equations, one from the |
368 |
gezelter |
3984 |
neutralization procedure ($-1/r_c$), and one that causes the first |
369 |
|
|
derivative to vanish at the cutoff radius. |
370 |
gezelter |
3982 |
|
371 |
|
|
Since one derivative of the interaction energy is needed for the |
372 |
|
|
force, the minimal perturbation is a term linear in $(r-r_c)$ in the |
373 |
|
|
interaction energy, that is: |
374 |
gezelter |
3906 |
\begin{equation} |
375 |
|
|
\frac{d\,}{dr} |
376 |
|
|
\left( {\frac{1}{r} - \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} } |
377 |
|
|
\right) = \left(- \frac{1}{r^2} + \frac{1}{r_c^2} |
378 |
|
|
\right) . |
379 |
|
|
\end{equation} |
380 |
gezelter |
3985 |
which clearly vanishes as the $r$ approaches the cutoff radius. There |
381 |
|
|
are a number of ways to generalize this derivative shift for |
382 |
gezelter |
3984 |
higher-order multipoles. Below, we present two methods, one based on |
383 |
|
|
higher-order Taylor series for $r$ near $r_c$, and the other based on |
384 |
|
|
linear shift of the kernel gradients at the cutoff itself. |
385 |
gezelter |
3906 |
|
386 |
gezelter |
3984 |
\subsection{Taylor-shifted force (TSF) electrostatics} |
387 |
gezelter |
3982 |
In the Taylor-shifted force (TSF) method, the procedure that we follow |
388 |
|
|
is based on a Taylor expansion containing the same number of |
389 |
|
|
derivatives required for each force term to vanish at the cutoff. For |
390 |
|
|
example, the quadrupole-quadrupole interaction energy requires four |
391 |
|
|
derivatives of the kernel, and the force requires one additional |
392 |
gezelter |
3986 |
derivative. For quadrupole-quadrupole interactions, we therefore |
393 |
|
|
require shifted energy expressions that include up to $(r-r_c)^5$ so |
394 |
|
|
that all energies, forces, and torques are zero as $r \rightarrow |
395 |
|
|
r_c$. In each case, we subtract off a function $f_n^{\text{shift}}(r)$ |
396 |
|
|
from the kernel $f(r)=1/r$. The subscript $n$ indicates the number of |
397 |
|
|
derivatives to be taken when deriving a given multipole energy. We |
398 |
|
|
choose a function with guaranteed smooth derivatives -- a truncated |
399 |
|
|
Taylor series of the function $f(r)$, e.g., |
400 |
gezelter |
3906 |
% |
401 |
|
|
\begin{equation} |
402 |
gezelter |
3984 |
f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)}(r_c) . |
403 |
gezelter |
3906 |
\end{equation} |
404 |
|
|
% |
405 |
|
|
The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$. |
406 |
|
|
Thus, for $f(r)=1/r$, we find |
407 |
|
|
% |
408 |
|
|
\begin{equation} |
409 |
|
|
f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} . |
410 |
|
|
\end{equation} |
411 |
|
|
% |
412 |
gezelter |
3982 |
Continuing with the example of a charge $\bf a$ interacting with a |
413 |
|
|
dipole $\bf b$, we write |
414 |
gezelter |
3906 |
% |
415 |
|
|
\begin{equation} |
416 |
|
|
U_{C_{\bf a}D_{\bf b}}(r)= |
417 |
gezelter |
3985 |
C_{\bf a} D_{{\bf b}\alpha} \frac {\partial f_1(r) }{\partial r_\alpha} |
418 |
|
|
= C_{\bf a} D_{{\bf b}\alpha} |
419 |
gezelter |
3906 |
\frac {r_\alpha}{r} \frac {\partial f_1(r)}{\partial r} . |
420 |
|
|
\end{equation} |
421 |
|
|
% |
422 |
gezelter |
3984 |
The force that dipole $\bf b$ exerts on charge $\bf a$ is |
423 |
gezelter |
3906 |
% |
424 |
|
|
\begin{equation} |
425 |
gezelter |
3985 |
F_{C_{\bf a}D_{\bf b}\beta} = C_{\bf a} D_{{\bf b}\alpha} |
426 |
gezelter |
3906 |
\left[ \frac{\delta_{\alpha\beta}}{r} \frac {\partial}{\partial r} + |
427 |
|
|
\frac{r_\alpha r_\beta}{r^2} |
428 |
|
|
\left( -\frac{1}{r} \frac {\partial} {\partial r} |
429 |
|
|
+ \frac {\partial ^2} {\partial r^2} \right) \right] f_1(r) . |
430 |
|
|
\end{equation} |
431 |
|
|
% |
432 |
gezelter |
3984 |
For undamped coulombic interactions, $f(r)=1/r$, we find |
433 |
gezelter |
3906 |
% |
434 |
|
|
\begin{equation} |
435 |
|
|
F_{C_{\bf a}D_{\bf b}\beta} = |
436 |
gezelter |
3985 |
\frac{C_{\bf a} D_{{\bf b}\beta}}{r} |
437 |
gezelter |
3906 |
\left[ -\frac{1}{r^2}+\frac{1}{r_c^2}-\frac{2(r-r_c)}{r_c^3} \right] |
438 |
gezelter |
3985 |
+C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta |
439 |
gezelter |
3906 |
\left[ \frac{3}{r^5}-\frac{3}{r^3r_c^2} \right] . |
440 |
|
|
\end{equation} |
441 |
|
|
% |
442 |
|
|
This expansion shows the expected $1/r^3$ dependence of the force. |
443 |
|
|
|
444 |
gezelter |
3984 |
In general, we can write |
445 |
gezelter |
3906 |
% |
446 |
|
|
\begin{equation} |
447 |
gezelter |
4098 |
U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) |
448 |
gezelter |
3906 |
\label{generic} |
449 |
|
|
\end{equation} |
450 |
|
|
% |
451 |
gezelter |
3985 |
with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for |
452 |
|
|
charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and |
453 |
|
|
$n=4$ for quadrupole-quadrupole. For example, in |
454 |
|
|
quadrupole-quadrupole interactions for which the $\text{prefactor}$ is |
455 |
|
|
$Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$, the derivatives are |
456 |
|
|
$\partial^4/\partial r_\alpha \partial r_\beta \partial |
457 |
|
|
r_\gamma \partial r_\delta$, with implied summation combining the |
458 |
gezelter |
4183 |
space indices. Appendix \ref{radialTSF} contains details on the |
459 |
|
|
radial functions. |
460 |
gezelter |
3906 |
|
461 |
gezelter |
3984 |
In the formulas presented in the tables below, the placeholder |
462 |
|
|
function $f(r)$ is used to represent the electrostatic kernel (either |
463 |
|
|
damped or undamped). The main functions that go into the force and |
464 |
gezelter |
3985 |
torque terms, $g_n(r), h_n(r), s_n(r), \mathrm{~and~} t_n(r)$ are |
465 |
|
|
successive derivatives of the shifted electrostatic kernel, $f_n(r)$ |
466 |
|
|
of the same index $n$. The algebra required to evaluate energies, |
467 |
|
|
forces and torques is somewhat tedious, so only the final forms are |
468 |
gezelter |
3986 |
presented in tables \ref{tab:tableenergy} and \ref{tab:tableFORCE}. |
469 |
gezelter |
3996 |
One of the principal findings of our work is that the individual |
470 |
|
|
orientational contributions to the various multipole-multipole |
471 |
|
|
interactions must be treated with distinct radial functions, but each |
472 |
|
|
of these contributions is independently force shifted at the cutoff |
473 |
|
|
radius. |
474 |
gezelter |
3906 |
|
475 |
gezelter |
3982 |
\subsection{Gradient-shifted force (GSF) electrostatics} |
476 |
gezelter |
3985 |
The second, and conceptually simpler approach to force-shifting |
477 |
|
|
maintains only the linear $(r-r_c)$ term in the truncated Taylor |
478 |
|
|
expansion, and has a similar interaction energy for all multipole |
479 |
|
|
orders: |
480 |
gezelter |
3906 |
\begin{equation} |
481 |
gezelter |
3996 |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
482 |
gezelter |
4183 |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) |
483 |
|
|
\hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
484 |
gezelter |
3985 |
\label{generic2} |
485 |
gezelter |
3906 |
\end{equation} |
486 |
gezelter |
4183 |
where $\hat{\mathbf{r}}$ is the unit vector pointing between the two |
487 |
|
|
multipoles, and the sum describes a separate force-shifting that is |
488 |
|
|
applied to each orientational contribution to the energy. Both the |
489 |
|
|
potential and the gradient for force shifting are evaluated for an |
490 |
|
|
image multipole projected onto the surface of the cutoff sphere (see |
491 |
|
|
fig \ref{fig:shiftedMultipoles}). The image multipole retains the |
492 |
gezelter |
3996 |
orientation ($\hat{\mathbf{b}}$) of the interacting multipole. No |
493 |
gezelter |
3990 |
higher order terms $(r-r_c)^n$ appear. The primary difference between |
494 |
|
|
the TSF and GSF methods is the stage at which the Taylor Series is |
495 |
|
|
applied; in the Taylor-shifted approach, it is applied to the kernel |
496 |
|
|
itself. In the Gradient-shifted approach, it is applied to individual |
497 |
gezelter |
4183 |
radial interaction terms in the multipole expansion. Energies from |
498 |
gezelter |
3990 |
this method thus have the general form: |
499 |
gezelter |
3906 |
\begin{equation} |
500 |
gezelter |
3985 |
U= \sum (\text{angular factor}) (\text{radial factor}). |
501 |
|
|
\label{generic3} |
502 |
gezelter |
3906 |
\end{equation} |
503 |
|
|
|
504 |
gezelter |
3986 |
Functional forms for both methods (TSF and GSF) can both be summarized |
505 |
gezelter |
4181 |
using the form of Eq.~\ref{generic3}). The basic forms for the |
506 |
gezelter |
3985 |
energy, force, and torque expressions are tabulated for both shifting |
507 |
gezelter |
3986 |
approaches below -- for each separate orientational contribution, only |
508 |
gezelter |
3985 |
the radial factors differ between the two methods. |
509 |
gezelter |
3906 |
|
510 |
gezelter |
4179 |
\subsection{Generalization of the Wolf shifted potential (SP)} |
511 |
|
|
It is also possible to formulate an extension of the Wolf approach for |
512 |
|
|
multipoles by simply projecting the image multipole onto the surface |
513 |
|
|
of the cutoff sphere, and including the interactions with the central |
514 |
|
|
multipole and the image. This effectively shifts the pair potential |
515 |
|
|
to zero at the cutoff radius, |
516 |
|
|
\begin{equation} |
517 |
|
|
U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
518 |
gezelter |
4183 |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
519 |
gezelter |
4179 |
\label{eq:SP} |
520 |
|
|
\end{equation} |
521 |
|
|
independent of the orientations of the two multipoles. The sum again |
522 |
|
|
describes separate potential shifting that is applied to each |
523 |
|
|
orientational contribution to the energy. |
524 |
|
|
|
525 |
|
|
The shifted potential (SP) method is a simple truncation of the GSF |
526 |
|
|
method for each orientational contribution, leaving out the $(r-r_c)$ |
527 |
|
|
terms that multiply the gradient. Functional forms for the |
528 |
|
|
shifted-potential (SP) method can also be summarized using the form of |
529 |
gezelter |
4181 |
Eq.~\ref{generic3}. The energy, force, and torque expressions are |
530 |
gezelter |
4179 |
tabulated below for all three methods. As in the GSF and TSF methods, |
531 |
|
|
for each separate orientational contribution, only the radial factors |
532 |
|
|
differ between the SP, GSF, and TSF methods. |
533 |
|
|
|
534 |
|
|
|
535 |
gezelter |
3906 |
\subsection{\label{sec:level2}Body and space axes} |
536 |
gezelter |
3989 |
Although objects $\bf a$ and $\bf b$ rotate during a molecular |
537 |
|
|
dynamics (MD) simulation, their multipole tensors remain fixed in |
538 |
|
|
body-frame coordinates. While deriving force and torque expressions, |
539 |
|
|
it is therefore convenient to write the energies, forces, and torques |
540 |
|
|
in intermediate forms involving the vectors of the rotation matrices. |
541 |
|
|
We denote body axes for objects $\bf a$ and $\bf b$ using unit vectors |
542 |
|
|
$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$. |
543 |
gezelter |
4179 |
In a typical simulation, the initial axes are obtained by |
544 |
gezelter |
3989 |
diagonalizing the moment of inertia tensors for the objects. (N.B., |
545 |
|
|
the body axes are generally {\it not} the same as those for which the |
546 |
|
|
quadrupole moment is diagonal.) The rotation matrices are then |
547 |
|
|
propagated during the simulation. |
548 |
gezelter |
3906 |
|
549 |
gezelter |
3989 |
The rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be |
550 |
gezelter |
3985 |
expressed using these unit vectors: |
551 |
gezelter |
3906 |
\begin{eqnarray} |
552 |
|
|
\hat{\mathbf {a}} = |
553 |
|
|
\begin{pmatrix} |
554 |
|
|
\hat{a}_1 \\ |
555 |
|
|
\hat{a}_2 \\ |
556 |
|
|
\hat{a}_3 |
557 |
gezelter |
3989 |
\end{pmatrix}, \qquad |
558 |
gezelter |
3906 |
\hat{\mathbf {b}} = |
559 |
|
|
\begin{pmatrix} |
560 |
|
|
\hat{b}_1 \\ |
561 |
|
|
\hat{b}_2 \\ |
562 |
|
|
\hat{b}_3 |
563 |
|
|
\end{pmatrix} |
564 |
|
|
\end{eqnarray} |
565 |
|
|
% |
566 |
gezelter |
3985 |
These matrices convert from space-fixed $(xyz)$ to body-fixed $(123)$ |
567 |
gezelter |
3989 |
coordinates. |
568 |
|
|
|
569 |
|
|
Allen and Germano,\cite{Allen:2006fk} following earlier work by Price |
570 |
|
|
{\em et al.},\cite{Price:1984fk} showed that if the interaction |
571 |
|
|
energies are written explicitly in terms of $\hat{r}$ and the body |
572 |
|
|
axes ($\hat{a}_m$, $\hat{b}_n$) : |
573 |
gezelter |
3906 |
% |
574 |
gezelter |
3985 |
\begin{equation} |
575 |
gezelter |
3989 |
U(r, \{\hat{a}_m \cdot \hat{r} \}, |
576 |
|
|
\{\hat{b}_n\cdot \hat{r} \}, |
577 |
|
|
\{\hat{a}_m \cdot \hat{b}_n \}) . |
578 |
|
|
\label{ugeneral} |
579 |
|
|
\end{equation} |
580 |
|
|
% |
581 |
|
|
the forces come out relatively cleanly, |
582 |
|
|
% |
583 |
|
|
\begin{equation} |
584 |
|
|
\mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} = \frac{\partial U}{\partial \mathbf{r}} |
585 |
|
|
= \frac{\partial U}{\partial r} \hat{r} |
586 |
|
|
+ \sum_m \left[ |
587 |
|
|
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} |
588 |
|
|
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} |
589 |
|
|
+ \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} |
590 |
|
|
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} |
591 |
|
|
\right] \label{forceequation}. |
592 |
|
|
\end{equation} |
593 |
|
|
|
594 |
|
|
The torques can also be found in a relatively similar |
595 |
|
|
manner, |
596 |
|
|
% |
597 |
|
|
\begin{eqnarray} |
598 |
|
|
\mathbf{\tau}_{\bf a} = |
599 |
|
|
\sum_m |
600 |
|
|
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} |
601 |
|
|
( \hat{r} \times \hat{a}_m ) |
602 |
|
|
-\sum_{mn} |
603 |
|
|
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} |
604 |
|
|
(\hat{a}_m \times \hat{b}_n) \\ |
605 |
|
|
% |
606 |
|
|
\mathbf{\tau}_{\bf b} = |
607 |
|
|
\sum_m |
608 |
|
|
\frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} |
609 |
|
|
( \hat{r} \times \hat{b}_m) |
610 |
|
|
+\sum_{mn} |
611 |
|
|
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} |
612 |
|
|
(\hat{a}_m \times \hat{b}_n) . |
613 |
|
|
\end{eqnarray} |
614 |
|
|
|
615 |
gezelter |
4179 |
Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $ |
616 |
gezelter |
3989 |
is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk} |
617 |
|
|
We also made use of the identities, |
618 |
|
|
% |
619 |
|
|
\begin{align} |
620 |
|
|
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} |
621 |
|
|
=& \frac{1}{r} \left( \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r} |
622 |
|
|
\right) \\ |
623 |
|
|
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} |
624 |
|
|
=& \frac{1}{r} \left( \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r} |
625 |
gezelter |
4179 |
\right). |
626 |
gezelter |
3989 |
\end{align} |
627 |
|
|
|
628 |
|
|
Many of the multipole contractions required can be written in one of |
629 |
|
|
three equivalent forms using the unit vectors $\hat{r}$, $\hat{a}_m$, |
630 |
|
|
and $\hat{b}_n$. In the torque expressions, it is useful to have the |
631 |
|
|
angular-dependent terms available in all three fashions, e.g. for the |
632 |
|
|
dipole-dipole contraction: |
633 |
|
|
% |
634 |
|
|
\begin{equation} |
635 |
gezelter |
3906 |
\mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}} |
636 |
gezelter |
3985 |
= D_{\bf {a}\alpha} D_{\bf {b}\alpha} = |
637 |
gezelter |
4179 |
\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}}. |
638 |
gezelter |
3985 |
\end{equation} |
639 |
gezelter |
3906 |
% |
640 |
gezelter |
3985 |
The first two forms are written using space coordinates. The first |
641 |
|
|
form is standard in the chemistry literature, while the second is |
642 |
|
|
expressed using implied summation notation. The third form shows |
643 |
|
|
explicit sums over body indices and the dot products now indicate |
644 |
|
|
contractions using space indices. |
645 |
gezelter |
3906 |
|
646 |
gezelter |
3989 |
In computing our force and torque expressions, we carried out most of |
647 |
|
|
the work in body coordinates, and have transformed the expressions |
648 |
|
|
back to space-frame coordinates, which are reported below. Interested |
649 |
|
|
readers may consult the supplemental information for this paper for |
650 |
|
|
the intermediate body-frame expressions. |
651 |
gezelter |
3906 |
|
652 |
gezelter |
3982 |
\subsection{The Self-Interaction \label{sec:selfTerm}} |
653 |
|
|
|
654 |
gezelter |
3985 |
In addition to cutoff-sphere neutralization, the Wolf |
655 |
|
|
summation~\cite{Wolf99} and the damped shifted force (DSF) |
656 |
gezelter |
4098 |
extension~\cite{Fennell:2006zl} also include self-interactions that |
657 |
gezelter |
3985 |
are handled separately from the pairwise interactions between |
658 |
|
|
sites. The self-term is normally calculated via a single loop over all |
659 |
|
|
sites in the system, and is relatively cheap to evaluate. The |
660 |
|
|
self-interaction has contributions from two sources. |
661 |
|
|
|
662 |
|
|
First, the neutralization procedure within the cutoff radius requires |
663 |
|
|
a contribution from a charge opposite in sign, but equal in magnitude, |
664 |
|
|
to the central charge, which has been spread out over the surface of |
665 |
|
|
the cutoff sphere. For a system of undamped charges, the total |
666 |
|
|
self-term is |
667 |
gezelter |
3980 |
\begin{equation} |
668 |
gezelter |
4179 |
U_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2}. |
669 |
gezelter |
3980 |
\label{eq:selfTerm} |
670 |
|
|
\end{equation} |
671 |
gezelter |
3985 |
|
672 |
|
|
Second, charge damping with the complementary error function is a |
673 |
|
|
partial analogy to the Ewald procedure which splits the interaction |
674 |
|
|
into real and reciprocal space sums. The real space sum is retained |
675 |
|
|
in the Wolf and DSF methods. The reciprocal space sum is first |
676 |
|
|
minimized by folding the largest contribution (the self-interaction) |
677 |
|
|
into the self-interaction from charge neutralization of the damped |
678 |
|
|
potential. The remainder of the reciprocal space portion is then |
679 |
|
|
discarded (as this contributes the largest computational cost and |
680 |
|
|
complexity to the Ewald sum). For a system containing only damped |
681 |
|
|
charges, the complete self-interaction can be written as |
682 |
gezelter |
3980 |
\begin{equation} |
683 |
gezelter |
4175 |
U_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} + |
684 |
gezelter |
3980 |
\frac{\alpha}{\sqrt{\pi}} \right) \sum_{{\bf a}=1}^N |
685 |
|
|
C_{\bf a}^{2}. |
686 |
|
|
\label{eq:dampSelfTerm} |
687 |
|
|
\end{equation} |
688 |
|
|
|
689 |
|
|
The extension of DSF electrostatics to point multipoles requires |
690 |
gezelter |
4183 |
treatment of the self-neutralization \textit{and} reciprocal |
691 |
gezelter |
3980 |
contributions to the self-interaction for higher order multipoles. In |
692 |
|
|
this section we give formulae for these interactions up to quadrupolar |
693 |
|
|
order. |
694 |
|
|
|
695 |
|
|
The self-neutralization term is computed by taking the {\it |
696 |
|
|
non-shifted} kernel for each interaction, placing a multipole of |
697 |
|
|
equal magnitude (but opposite in polarization) on the surface of the |
698 |
|
|
cutoff sphere, and averaging over the surface of the cutoff sphere. |
699 |
|
|
Because the self term is carried out as a single sum over sites, the |
700 |
|
|
reciprocal-space portion is identical to half of the self-term |
701 |
gezelter |
4179 |
obtained by Smith, and also by Aguado and Madden for the application |
702 |
|
|
of the Ewald sum to multipoles.\cite{Smith82,Smith98,Aguado03} For a |
703 |
|
|
given site which posesses a charge, dipole, and quadrupole, both types |
704 |
|
|
of contribution are given in table \ref{tab:tableSelf}. |
705 |
gezelter |
3980 |
|
706 |
|
|
\begin{table*} |
707 |
|
|
\caption{\label{tab:tableSelf} Self-interaction contributions for |
708 |
|
|
site ({\bf a}) that has a charge $(C_{\bf a})$, dipole |
709 |
|
|
$(\mathbf{D}_{\bf a})$, and quadrupole $(\mathbf{Q}_{\bf a})$} |
710 |
|
|
\begin{ruledtabular} |
711 |
|
|
\begin{tabular}{lccc} |
712 |
|
|
Multipole order & Summed Quantity & Self-neutralization & Reciprocal \\ \hline |
713 |
|
|
Charge & $C_{\bf a}^2$ & $-f(r_c)$ & $-\frac{\alpha}{\sqrt{\pi}}$ \\ |
714 |
|
|
Dipole & $|\mathbf{D}_{\bf a}|^2$ & $\frac{1}{3} \left( h(r_c) + |
715 |
|
|
\frac{2 g(r_c)}{r_c} \right)$ & $-\frac{2 \alpha^3}{3 \sqrt{\pi}}$\\ |
716 |
gezelter |
3989 |
Quadrupole & $2 \mathbf{Q}_{\bf a}:\mathbf{Q}_{\bf a} + \text{Tr}(\mathbf{Q}_{\bf a})^2$ & |
717 |
gezelter |
3980 |
$- \frac{1}{15} \left( t(r_c)+ \frac{4 s(r_c)}{r_c} \right)$ & |
718 |
|
|
$-\frac{4 \alpha^5}{5 \sqrt{\pi}}$ \\ |
719 |
|
|
Charge-Quadrupole & $-2 C_{\bf a} \text{Tr}(\mathbf{Q}_{\bf a})$ & $\frac{1}{3} \left( |
720 |
|
|
h(r_c) + \frac{2 g(r_c)}{r_c} \right)$& $-\frac{2 \alpha^3}{3 \sqrt{\pi}}$ \\ |
721 |
|
|
\end{tabular} |
722 |
|
|
\end{ruledtabular} |
723 |
|
|
\end{table*} |
724 |
|
|
|
725 |
|
|
For sites which simultaneously contain charges and quadrupoles, the |
726 |
|
|
self-interaction includes a cross-interaction between these two |
727 |
|
|
multipole orders. Symmetry prevents the charge-dipole and |
728 |
|
|
dipole-quadrupole interactions from contributing to the |
729 |
|
|
self-interaction. The functions that go into the self-neutralization |
730 |
gezelter |
3985 |
terms, $g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive |
731 |
|
|
derivatives of the electrostatic kernel, $f(r)$ (either the undamped |
732 |
|
|
$1/r$ or the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that |
733 |
|
|
have been evaluated at the cutoff distance. For undamped |
734 |
|
|
interactions, $f(r_c) = 1/r_c$, $g(r_c) = -1/r_c^{2}$, and so on. For |
735 |
|
|
damped interactions, $f(r_c) = B_0(r_c)$, $g(r_c) = B_0'(r_c)$, and so |
736 |
|
|
on. Appendix \ref{SmithFunc} contains recursion relations that allow |
737 |
|
|
rapid evaluation of these derivatives. |
738 |
gezelter |
3980 |
|
739 |
gezelter |
3985 |
\section{Interaction energies, forces, and torques} |
740 |
|
|
The main result of this paper is a set of expressions for the |
741 |
|
|
energies, forces and torques (up to quadrupole-quadrupole order) that |
742 |
gezelter |
4183 |
work for the Taylor-shifted, gradient-shifted, and shifted potential |
743 |
|
|
approximations. These expressions were derived using a set of generic |
744 |
|
|
radial functions. Without using the shifting approximations mentioned |
745 |
|
|
above, some of these radial functions would be identical, and the |
746 |
|
|
expressions coalesce into the familiar forms for unmodified |
747 |
|
|
multipole-multipole interactions. Table \ref{tab:tableenergy} maps |
748 |
|
|
between the generic functions and the radial functions derived for the |
749 |
|
|
three methods. The energy equations are written in terms of lab-frame |
750 |
|
|
representations of the dipoles, quadrupoles, and the unit vector |
751 |
|
|
connecting the two objects, |
752 |
gezelter |
3906 |
|
753 |
|
|
% Energy in space coordinate form ---------------------------------------------------------------------------------------------- |
754 |
|
|
% |
755 |
|
|
% |
756 |
|
|
% u ca cb |
757 |
|
|
% |
758 |
gezelter |
3983 |
\begin{align} |
759 |
|
|
U_{C_{\bf a}C_{\bf b}}(r)=& |
760 |
gezelter |
3985 |
C_{\bf a} C_{\bf b} v_{01}(r) \label{uchch} |
761 |
gezelter |
3983 |
\\ |
762 |
gezelter |
3906 |
% |
763 |
|
|
% u ca db |
764 |
|
|
% |
765 |
gezelter |
3983 |
U_{C_{\bf a}D_{\bf b}}(r)=& |
766 |
gezelter |
3985 |
C_{\bf a} \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right) v_{11}(r) |
767 |
gezelter |
3906 |
\label{uchdip} |
768 |
gezelter |
3983 |
\\ |
769 |
gezelter |
3906 |
% |
770 |
|
|
% u ca qb |
771 |
|
|
% |
772 |
gezelter |
3985 |
U_{C_{\bf a}Q_{\bf b}}(r)=& C_{\bf a } \Bigl[ \text{Tr}Q_{\bf b} |
773 |
|
|
v_{21}(r) + \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot |
774 |
|
|
\hat{r} \right) v_{22}(r) \Bigr] |
775 |
gezelter |
3906 |
\label{uchquad} |
776 |
gezelter |
3983 |
\\ |
777 |
gezelter |
3906 |
% |
778 |
|
|
% u da cb |
779 |
|
|
% |
780 |
gezelter |
3983 |
%U_{D_{\bf a}C_{\bf b}}(r)=& |
781 |
|
|
%-\frac{C_{\bf b}}{4\pi \epsilon_0} |
782 |
|
|
%\left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right) v_{11}(r) \label{udipch} |
783 |
|
|
%\\ |
784 |
gezelter |
3906 |
% |
785 |
|
|
% u da db |
786 |
|
|
% |
787 |
gezelter |
3983 |
U_{D_{\bf a}D_{\bf b}}(r)=& |
788 |
gezelter |
3985 |
-\Bigr[ \left( \mathbf{D}_{\mathbf {a}} \cdot |
789 |
gezelter |
3906 |
\mathbf{D}_{\mathbf{b}} \right) v_{21}(r) |
790 |
|
|
+\left( \mathbf{D}_{\mathbf {a}} \cdot \hat{r} \right) |
791 |
|
|
\left( \mathbf{D}_{\mathbf {b}} \cdot \hat{r} \right) |
792 |
|
|
v_{22}(r) \Bigr] |
793 |
|
|
\label{udipdip} |
794 |
gezelter |
3983 |
\\ |
795 |
gezelter |
3906 |
% |
796 |
|
|
% u da qb |
797 |
|
|
% |
798 |
|
|
\begin{split} |
799 |
|
|
% 1 |
800 |
gezelter |
3983 |
U_{D_{\bf a}Q_{\bf b}}(r) =& |
801 |
gezelter |
3985 |
-\Bigl[ |
802 |
gezelter |
3906 |
\text{Tr}\mathbf{Q}_{\mathbf{b}} |
803 |
|
|
\left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right) |
804 |
|
|
+2 ( \mathbf{D}_{\mathbf{a}} \cdot |
805 |
|
|
\mathbf{Q}_{\mathbf{b}} \cdot \hat{r} ) \Bigr] v_{31}(r) \\ |
806 |
|
|
% 2 |
807 |
gezelter |
3985 |
&- \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right) |
808 |
gezelter |
3906 |
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{32}(r) |
809 |
|
|
\label{udipquad} |
810 |
|
|
\end{split} |
811 |
gezelter |
3983 |
\\ |
812 |
gezelter |
3906 |
% |
813 |
|
|
% u qa cb |
814 |
|
|
% |
815 |
gezelter |
3983 |
%U_{Q_{\bf a}C_{\bf b}}(r)=& |
816 |
|
|
%\frac{C_{\bf b }}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\bf a} v_{21}(r) |
817 |
|
|
%\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) v_{22}(r) \Bigr] |
818 |
|
|
%\label{uquadch} |
819 |
|
|
%\\ |
820 |
gezelter |
3906 |
% |
821 |
|
|
% u qa db |
822 |
|
|
% |
823 |
gezelter |
3983 |
%\begin{split} |
824 |
gezelter |
3906 |
%1 |
825 |
gezelter |
3983 |
%U_{Q_{\bf a}D_{\bf b}}(r)=& |
826 |
|
|
%\frac{1}{4\pi \epsilon_0} \Bigl[ |
827 |
|
|
%\text{Tr}\mathbf{Q}_{\mathbf{a}} |
828 |
|
|
%\left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right) |
829 |
|
|
%+2 ( \mathbf{D}_{\mathbf{b}} \cdot |
830 |
|
|
%\mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r)\\ |
831 |
gezelter |
3906 |
% 2 |
832 |
gezelter |
3983 |
%&+\frac{1}{4\pi \epsilon_0} |
833 |
|
|
%\left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right) |
834 |
|
|
%\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) v_{32}(r) |
835 |
|
|
%\label{uquaddip} |
836 |
|
|
%\end{split} |
837 |
|
|
%\\ |
838 |
gezelter |
3906 |
% |
839 |
|
|
% u qa qb |
840 |
|
|
% |
841 |
|
|
\begin{split} |
842 |
|
|
%1 |
843 |
gezelter |
3983 |
U_{Q_{\bf a}Q_{\bf b}}(r)=& |
844 |
gezelter |
3985 |
\Bigl[ |
845 |
gezelter |
3906 |
\text{Tr} \mathbf{Q}_{\mathbf{a}} \text{Tr} \mathbf{Q}_{\mathbf{b}} |
846 |
gezelter |
3989 |
+2 |
847 |
|
|
\mathbf{Q}_{\mathbf{a}} : \mathbf{Q}_{\mathbf{b}} \Bigr] v_{41}(r) |
848 |
gezelter |
3906 |
\\ |
849 |
|
|
% 2 |
850 |
gezelter |
3985 |
&+\Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}} |
851 |
gezelter |
3906 |
\left( \hat{r} \cdot |
852 |
|
|
\mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) |
853 |
|
|
+\text{Tr}\mathbf{Q}_{\mathbf{b}} |
854 |
|
|
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} |
855 |
|
|
\cdot \hat{r} \right) +4 (\hat{r} \cdot |
856 |
|
|
\mathbf{Q}_{{\mathbf a}}\cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) |
857 |
|
|
\Bigr] v_{42}(r) |
858 |
|
|
\\ |
859 |
|
|
% 4 |
860 |
gezelter |
3985 |
&+ |
861 |
gezelter |
3906 |
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) |
862 |
|
|
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{43}(r). |
863 |
|
|
\label{uquadquad} |
864 |
|
|
\end{split} |
865 |
gezelter |
3983 |
\end{align} |
866 |
gezelter |
3985 |
% |
867 |
gezelter |
3983 |
Note that the energies of multipoles on site $\mathbf{b}$ interacting |
868 |
|
|
with those on site $\mathbf{a}$ can be obtained by swapping indices |
869 |
|
|
along with the sign of the intersite vector, $\hat{r}$. |
870 |
gezelter |
3906 |
|
871 |
|
|
% |
872 |
|
|
% |
873 |
|
|
% TABLE of radial functions ---------------------------------------------------------------------------------------------------------------- |
874 |
|
|
% |
875 |
|
|
|
876 |
gezelter |
3985 |
\begin{sidewaystable} |
877 |
|
|
\caption{\label{tab:tableenergy}Radial functions used in the energy |
878 |
gezelter |
4183 |
and torque equations. The $f, g, h, s, t, \mathrm{and~} u$ |
879 |
|
|
functions used in this table are defined in Appendices |
880 |
|
|
\ref{radialTSF} and \ref{radialGSF}. The gradient shifted (GSF) |
881 |
|
|
functions include the shifted potential (SP) |
882 |
|
|
contributions (\textit{cf.} Eqs. \ref{generic2} and |
883 |
|
|
\ref{eq:SP}).} |
884 |
gezelter |
4177 |
\begin{tabular}{|c|c|l|l|l|} \hline |
885 |
|
|
Generic&Bare Coulomb&Taylor-Shifted (TSF)&Shifted Potential (SP)&Gradient-Shifted (GSF) |
886 |
gezelter |
3906 |
\\ \hline |
887 |
|
|
% |
888 |
|
|
% |
889 |
|
|
% |
890 |
|
|
%Ch-Ch& |
891 |
|
|
$v_{01}(r)$ & |
892 |
|
|
$\frac{1}{r}$ & |
893 |
|
|
$f_0(r)$ & |
894 |
gezelter |
4177 |
$f(r)-f(r_c)$ & |
895 |
|
|
SP $-(r-r_c)g(r_c)$ |
896 |
gezelter |
3906 |
\\ |
897 |
|
|
% |
898 |
|
|
% |
899 |
|
|
% |
900 |
|
|
%Ch-Di& |
901 |
|
|
$v_{11}(r)$ & |
902 |
|
|
$-\frac{1}{r^2}$ & |
903 |
|
|
$g_1(r)$ & |
904 |
gezelter |
4177 |
$g(r)-g(r_c)$ & |
905 |
|
|
SP $-(r-r_c)h(r_c)$ \\ |
906 |
gezelter |
3906 |
% |
907 |
|
|
% |
908 |
|
|
% |
909 |
|
|
%Ch-Qu/Di-Di& |
910 |
|
|
$v_{21}(r)$ & |
911 |
|
|
$-\frac{1}{r^3} $ & |
912 |
|
|
$\frac{g_2(r)}{r} $ & |
913 |
gezelter |
4177 |
$\frac{g(r)}{r}-\frac{g(r_c)}{r_c}$ & |
914 |
gezelter |
4179 |
SP $-(r-r_c) \left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right)$ \\ |
915 |
|
|
% |
916 |
|
|
% |
917 |
|
|
% |
918 |
gezelter |
3906 |
$v_{22}(r)$ & |
919 |
|
|
$\frac{3}{r^3} $ & |
920 |
|
|
$\left(-\frac{g_2(r)}{r} + h_2(r) \right)$ & |
921 |
gezelter |
4179 |
$\left(-\frac{g(r)}{r}+h(r) \right) -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right)$ |
922 |
|
|
& SP $-(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$\\ |
923 |
gezelter |
3906 |
% |
924 |
|
|
% |
925 |
|
|
% |
926 |
|
|
%Di-Qu & |
927 |
|
|
$v_{31}(r)$ & |
928 |
|
|
$\frac{3}{r^4} $ & |
929 |
|
|
$\left(-\frac{g_3(r)}{r^2} + \frac{h_3(r)}{r} \right)$ & |
930 |
gezelter |
4179 |
$\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r}\right)-\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right)$ |
931 |
|
|
& SP $-(r-r_c) \left(\frac{2g(r_c)}{r_c^3}-\frac{2h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ \\ |
932 |
gezelter |
3906 |
% |
933 |
gezelter |
4179 |
% |
934 |
|
|
% |
935 |
gezelter |
3906 |
$v_{32}(r)$ & |
936 |
|
|
$-\frac{15}{r^4} $ & |
937 |
|
|
$\left( \frac{3g_3(r)}{r^2} - \frac{3h_3(r)}{r} + s_3(r) \right)$ & |
938 |
gezelter |
4179 |
$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)$& |
939 |
|
|
SP $-(r-r_c) \left( \frac{-6g(r_c)}{r_c^3}+\frac{6h(r_c)}{r_c^2}\right.$ \\ |
940 |
|
|
&&& $~~~-\left(\frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c)\right)$ & |
941 |
|
|
$\phantom{SP-(r-r_c)}\left.-\frac{3s(r_c)}{r_c}+t(r_c) \right)$\\ |
942 |
gezelter |
3906 |
% |
943 |
|
|
% |
944 |
|
|
% |
945 |
|
|
%Qu-Qu& |
946 |
|
|
$v_{41}(r)$ & |
947 |
|
|
$\frac{3}{r^5} $ & |
948 |
|
|
$\left(-\frac{g_4(r)}{r^3} +\frac{h_4(r)}{r^2} \right) $ & |
949 |
gezelter |
4179 |
$\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right)- \left(-\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2} \right)$ & |
950 |
|
|
SP $-(r-r_c) \left( \frac{3g(r_c)}{r_c^4}-\frac{3h(r_c)}{r_c^3}+\frac{s(r_c)}{r_c^2} \right)$ |
951 |
gezelter |
3906 |
\\ |
952 |
|
|
% 2 |
953 |
|
|
$v_{42}(r)$ & |
954 |
|
|
$- \frac{15}{r^5} $ & |
955 |
|
|
$\left( \frac{3g_4(r)}{r^3} - \frac{3h_4(r)}{r^2}+\frac{s_4(r)}{r} \right)$ & |
956 |
gezelter |
4179 |
$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right)$ & |
957 |
|
|
SP$-(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}\right.$ \\ |
958 |
|
|
&&& $~~~-\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ & |
959 |
|
|
$\phantom{SP-(r-r_c)}\left. -\frac{4s(r_c)}{r_c^2} + \frac{t(r_c)}{r_c}\right)$\\ |
960 |
gezelter |
3906 |
% 3 |
961 |
gezelter |
4179 |
% |
962 |
|
|
% |
963 |
gezelter |
3906 |
$v_{43}(r)$ & |
964 |
|
|
$ \frac{105}{r^5} $ & |
965 |
|
|
$\left(-\frac{15g_4(r)}{r^3}+\frac{15h_4(r)}{r^2}-\frac{6s_4(r)}{r} + t_4(r)\right) $ & |
966 |
gezelter |
4179 |
$ \left(-\frac{15g(r)}{r^3} +\frac{15h(r)}{r^2}-\frac{6s(r)}{r}+t(r)\right) $ & |
967 |
|
|
SP $-(r-r_c)\left(\frac{45g(r_c)}{r_c^4}-\frac{45h(r_c)}{r_c^3}\right.$\\ |
968 |
|
|
&&& $~~~-\left(-\frac{15g(r_c)}{r_c^3}+\frac{15h(r_c)}{r_c^2}-\frac{6s(r_c)}{r_c}+ t(r_c)\right)$ & |
969 |
|
|
$\phantom{SP-(r-r_c)}\left.+\frac{21s(r_c)}{r_c^2}-\frac{6t(r_c)}{r_c}+u(r_c) \right)$\\ |
970 |
gezelter |
4177 |
\hline |
971 |
gezelter |
3906 |
\end{tabular} |
972 |
gezelter |
3985 |
\end{sidewaystable} |
973 |
gezelter |
3906 |
% |
974 |
|
|
% |
975 |
|
|
% FORCE TABLE of radial functions ---------------------------------------------------------------------------------------------------------------- |
976 |
|
|
% |
977 |
|
|
|
978 |
gezelter |
3985 |
\begin{sidewaystable} |
979 |
gezelter |
4179 |
\caption{\label{tab:tableFORCE}Radial functions used in the force |
980 |
|
|
equations. Gradient shifted (GSF) functions are constructed using the shifted |
981 |
|
|
potential (SP) functions. Some of these functions are simple |
982 |
|
|
modifications of the functions found in table \ref{tab:tableenergy}} |
983 |
|
|
\begin{tabular}{|c|c|l|l|l|} \hline |
984 |
|
|
Function&Definition&Taylor-Shifted (TSF)& Shifted Potential (SP) |
985 |
|
|
&Gradient-Shifted (GSF) |
986 |
gezelter |
3906 |
\\ \hline |
987 |
|
|
% |
988 |
|
|
% |
989 |
|
|
% |
990 |
|
|
$w_a(r)$& |
991 |
gezelter |
3985 |
$\frac{d v_{01}}{dr}$& |
992 |
|
|
$g_0(r)$& |
993 |
gezelter |
4179 |
$g(r)$& |
994 |
|
|
SP $-g(r_c)$ \\ |
995 |
gezelter |
3906 |
% |
996 |
|
|
% |
997 |
|
|
$w_b(r)$ & |
998 |
gezelter |
3985 |
$\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $& |
999 |
|
|
$\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ & |
1000 |
gezelter |
4179 |
$h(r) - \frac{v_{11}(r)}{r} $ & |
1001 |
|
|
SP $- h(r_c)$ \\ |
1002 |
gezelter |
3906 |
% |
1003 |
|
|
$w_c(r)$ & |
1004 |
gezelter |
3985 |
$\frac{v_{11}(r)}{r}$ & |
1005 |
|
|
$\frac{g_1(r)}{r} $ & |
1006 |
gezelter |
4179 |
$\frac{v_{11}(r)}{r}$& |
1007 |
gezelter |
3985 |
$\frac{v_{11}(r)}{r}$\\ |
1008 |
gezelter |
3906 |
% |
1009 |
|
|
% |
1010 |
|
|
$w_d(r)$& |
1011 |
gezelter |
3985 |
$\frac{d v_{21}}{dr}$& |
1012 |
|
|
$\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ & |
1013 |
gezelter |
4179 |
$\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)$ & |
1014 |
|
|
SP $-\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $ \\ |
1015 |
gezelter |
3906 |
% |
1016 |
|
|
$w_e(r)$ & |
1017 |
gezelter |
3985 |
$\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ & |
1018 |
|
|
$\frac{v_{22}(r)}{r}$ & |
1019 |
gezelter |
4179 |
$\frac{v_{22}(r)}{r}$ & |
1020 |
gezelter |
3906 |
$\frac{v_{22}(r)}{r}$ \\ |
1021 |
|
|
% |
1022 |
|
|
% |
1023 |
|
|
$w_f(r)$& |
1024 |
gezelter |
3985 |
$\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$& |
1025 |
|
|
$\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ & |
1026 |
gezelter |
4179 |
$ \left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) -\frac{2v_{22}(r)}{r}$& |
1027 |
|
|
SP $- \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$\\ |
1028 |
gezelter |
3906 |
% |
1029 |
|
|
$w_g(r)$& |
1030 |
gezelter |
3985 |
$\frac{v_{31}(r)}{r}$& |
1031 |
|
|
$ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$& |
1032 |
gezelter |
4179 |
$\frac{v_{31}(r)}{r}$& |
1033 |
gezelter |
3906 |
$\frac{v_{31}(r)}{r}$\\ |
1034 |
|
|
% |
1035 |
|
|
$w_h(r)$ & |
1036 |
gezelter |
3985 |
$\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$& |
1037 |
|
|
$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $ & |
1038 |
gezelter |
4179 |
$ \left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) -\frac{v_{31}(r)}{r}$ & |
1039 |
|
|
SP $ - \left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\ |
1040 |
gezelter |
3906 |
% 2 |
1041 |
|
|
$w_i(r)$ & |
1042 |
gezelter |
3985 |
$\frac{v_{32}(r)}{r}$ & |
1043 |
|
|
$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $ & |
1044 |
gezelter |
4179 |
$\frac{v_{32}(r)}{r}$& |
1045 |
gezelter |
3985 |
$\frac{v_{32}(r)}{r}$\\ |
1046 |
gezelter |
3906 |
% |
1047 |
|
|
$w_j(r)$ & |
1048 |
gezelter |
3985 |
$\frac{d v_{32}}{dr} - \frac{3v_{32}}{r}$& |
1049 |
|
|
$\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right) $ & |
1050 |
gezelter |
4179 |
$\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right) -\frac{3v_{32}}{r}$ & |
1051 |
|
|
SP $-\left(\frac{-6g(_cr)}{r_c^3} +\frac{6h(r_c)}{r_c^2} |
1052 |
|
|
-\frac{3s(r_c)}{r_c} +t(r_c) \right)$ \\ |
1053 |
gezelter |
3906 |
% |
1054 |
|
|
$w_k(r)$ & |
1055 |
gezelter |
3985 |
$\frac{d v_{41}}{dr} $ & |
1056 |
|
|
$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2} \right)$ & |
1057 |
gezelter |
4179 |
$\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2} |
1058 |
|
|
\right)$ & |
1059 |
|
|
SP $-\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2} \right)$ \\ |
1060 |
gezelter |
3906 |
% |
1061 |
|
|
$w_l(r)$ & |
1062 |
gezelter |
3985 |
$\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ & |
1063 |
|
|
$\left(-\frac{15g_4(r)}{r^4} +\frac{15h_4(r)}{r^3} -\frac{6s_4(r)}{r^2} +\frac{t_4(r)}{r} \right)$ & |
1064 |
gezelter |
4179 |
$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} |
1065 |
|
|
+\frac{t(r)}{r} \right) -\frac{2v_{42}(r)}{r}$& |
1066 |
|
|
SP$-\left(-\frac{9g(r_c)}{r_c^4} +\frac{9h(r_c)}{r_c^3} -\frac{4s(r_c)}{r_c^2} +\frac{t(r_c)}{r_c} \right)$\\ |
1067 |
gezelter |
3906 |
% |
1068 |
|
|
$w_m(r)$ & |
1069 |
gezelter |
3985 |
$\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$& |
1070 |
gezelter |
4179 |
$\left(\frac{105g_4(r)}{r^4} - \frac{105h_4(r)}{r^3} \right.$ & |
1071 |
|
|
$\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2}\right.$ & |
1072 |
|
|
SP $- \left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}\right.$ \\ |
1073 |
|
|
&& $~~~\left.+ \frac{45s_4(r)}{r^2} - \frac{10t_4(r)}{r} +u_4(r) \right)$ |
1074 |
|
|
& $~~~\left. -\frac{6t(r)}{r} +u(r) \right) -\frac{4v_{43}(r)}{r}$ & |
1075 |
|
|
$\phantom{SP-} \left.+\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $\\ |
1076 |
gezelter |
3906 |
% |
1077 |
|
|
$w_n(r)$ & |
1078 |
gezelter |
3985 |
$\frac{v_{42}(r)}{r}$ & |
1079 |
|
|
$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2} \right)$ & |
1080 |
gezelter |
4179 |
$\frac{v_{42}(r)}{r}$& |
1081 |
gezelter |
3985 |
$\frac{v_{42}(r)}{r}$\\ |
1082 |
gezelter |
3906 |
% |
1083 |
|
|
$w_o(r)$ & |
1084 |
gezelter |
3985 |
$\frac{v_{43}(r)}{r}$& |
1085 |
|
|
$\left(-\frac{15g_4(r)}{r^4} +\frac{15h_4(r)}{r^3} -\frac{6s_4(r)}{r^2} +\frac{t_4(r)}{r} \right)$ & |
1086 |
gezelter |
4179 |
$\frac{v_{43}(r)}{r}$& |
1087 |
gezelter |
3985 |
$\frac{v_{43}(r)}{r}$ \\ \hline |
1088 |
gezelter |
3906 |
% |
1089 |
|
|
|
1090 |
|
|
\end{tabular} |
1091 |
gezelter |
3985 |
\end{sidewaystable} |
1092 |
gezelter |
3906 |
% |
1093 |
|
|
% |
1094 |
|
|
% |
1095 |
|
|
|
1096 |
|
|
\subsection{Forces} |
1097 |
gezelter |
3985 |
The force on object $\bf{a}$, $\mathbf{F}_{\bf a}$, due to object |
1098 |
|
|
$\bf{b}$ is the negative of the force on $\bf{b}$ due to $\bf{a}$. For |
1099 |
|
|
a simple charge-charge interaction, these forces will point along the |
1100 |
|
|
$\pm \hat{r}$ directions, where $\mathbf{r}=\mathbf{r}_b - |
1101 |
|
|
\mathbf{r}_a $. Thus |
1102 |
gezelter |
3906 |
% |
1103 |
|
|
\begin{equation} |
1104 |
|
|
F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}}{\partial r} |
1105 |
|
|
\quad \text{and} \quad F_{\bf b \alpha} |
1106 |
|
|
= - \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}} {\partial r} . |
1107 |
|
|
\end{equation} |
1108 |
|
|
% |
1109 |
gezelter |
3985 |
We list below the force equations written in terms of lab-frame |
1110 |
gezelter |
4183 |
coordinates. The radial functions used in the three methods are listed |
1111 |
gezelter |
3985 |
in Table \ref{tab:tableFORCE} |
1112 |
gezelter |
3906 |
% |
1113 |
gezelter |
3985 |
%SPACE COORDINATES FORCE EQUATIONS |
1114 |
gezelter |
3906 |
% |
1115 |
|
|
% ************************************************************************** |
1116 |
|
|
% f ca cb |
1117 |
|
|
% |
1118 |
gezelter |
3985 |
\begin{align} |
1119 |
|
|
\mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =& |
1120 |
|
|
C_{\bf a} C_{\bf b} w_a(r) \hat{r} \\ |
1121 |
gezelter |
3906 |
% |
1122 |
|
|
% |
1123 |
|
|
% |
1124 |
gezelter |
3985 |
\mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =& |
1125 |
|
|
C_{\bf a} \Bigl[ |
1126 |
gezelter |
3906 |
\left( \hat{r} \cdot \mathbf{D}_{\mathbf{b}} \right) |
1127 |
|
|
w_b(r) \hat{r} |
1128 |
gezelter |
3985 |
+ \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr] \\ |
1129 |
gezelter |
3906 |
% |
1130 |
|
|
% |
1131 |
|
|
% |
1132 |
gezelter |
3985 |
\mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =& |
1133 |
|
|
C_{\bf a } \Bigr[ |
1134 |
gezelter |
3906 |
\text{Tr}\mathbf{Q}_{\bf b} w_d(r) \hat{r} |
1135 |
|
|
+ 2 \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} w_e(r) |
1136 |
gezelter |
3985 |
+ \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} |
1137 |
|
|
\right) w_f(r) \hat{r} \Bigr] \\ |
1138 |
gezelter |
3906 |
% |
1139 |
|
|
% |
1140 |
|
|
% |
1141 |
gezelter |
3985 |
% \begin{equation} |
1142 |
|
|
% \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} = |
1143 |
|
|
% -C_{\bf{b}} \Bigl[ |
1144 |
|
|
% \left( \hat{r} \cdot \mathbf{D}_{\mathbf{a}} \right) w_b(r) \hat{r} |
1145 |
|
|
% + \mathbf{D}_{\mathbf{a}} w_c(r) \Bigr] |
1146 |
|
|
% \end{equation} |
1147 |
gezelter |
3906 |
% |
1148 |
|
|
% |
1149 |
|
|
% |
1150 |
gezelter |
3985 |
\begin{split} |
1151 |
|
|
\mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} =& |
1152 |
gezelter |
3906 |
- \mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}} w_d(r) \hat{r} |
1153 |
|
|
+ \left( \mathbf{D}_{\mathbf {a}} |
1154 |
|
|
\left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right) |
1155 |
gezelter |
3985 |
+ \mathbf{D}_{\mathbf {b}} \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right) \right) w_e(r)\\ |
1156 |
gezelter |
3906 |
% 2 |
1157 |
gezelter |
3985 |
& - \left( \hat{r} \cdot \mathbf{D}_{\mathbf {a}} \right) |
1158 |
|
|
\left( \hat{r} \cdot \mathbf{D}_{\mathbf {b}} \right) w_f(r) \hat{r} |
1159 |
|
|
\end{split}\\ |
1160 |
gezelter |
3906 |
% |
1161 |
|
|
% |
1162 |
|
|
% |
1163 |
|
|
\begin{split} |
1164 |
gezelter |
3985 |
\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =& - \Bigl[ |
1165 |
gezelter |
3906 |
\text{Tr}\mathbf{Q}_{\mathbf{b}} \mathbf{ D}_{\mathbf{a}} |
1166 |
|
|
+2 \mathbf{D}_{\mathbf{a}} \cdot |
1167 |
|
|
\mathbf{Q}_{\mathbf{b}} \Bigr] w_g(r) |
1168 |
gezelter |
3985 |
- \Bigl[ |
1169 |
gezelter |
3906 |
\text{Tr}\mathbf{Q}_{\mathbf{b}} |
1170 |
|
|
\left( \hat{r} \cdot \mathbf{D}_{\mathbf{a}} \right) |
1171 |
|
|
+2 ( \mathbf{D}_{\mathbf{a}} \cdot |
1172 |
|
|
\mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\ |
1173 |
|
|
% 3 |
1174 |
gezelter |
3985 |
& - \Bigl[\mathbf{ D}_{\mathbf{a}} (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) |
1175 |
gezelter |
3906 |
+2 (\hat{r} \cdot \mathbf{D}_{\mathbf{a}} ) (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} ) \Bigr] |
1176 |
|
|
w_i(r) |
1177 |
|
|
% 4 |
1178 |
gezelter |
3985 |
- |
1179 |
gezelter |
3906 |
(\hat{r} \cdot \mathbf{D}_{\mathbf{a}} ) |
1180 |
gezelter |
3985 |
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r} \end{split} \\ |
1181 |
gezelter |
3906 |
% |
1182 |
|
|
% |
1183 |
gezelter |
3985 |
% \begin{equation} |
1184 |
|
|
% \mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} = |
1185 |
|
|
% \frac{C_{\bf b }}{4\pi \epsilon_0} \Bigr[ |
1186 |
|
|
% \text{Tr}\mathbf{Q}_{\bf a} w_d(r) \hat{r} |
1187 |
|
|
% + 2 \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} w_e(r) |
1188 |
|
|
% + \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr] |
1189 |
|
|
% \end{equation} |
1190 |
|
|
% % |
1191 |
|
|
% \begin{equation} |
1192 |
|
|
% \begin{split} |
1193 |
|
|
% \mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} = |
1194 |
|
|
% &\frac{1}{4\pi \epsilon_0} \Bigl[ |
1195 |
|
|
% \text{Tr}\mathbf{Q}_{\mathbf{a}} \mathbf{D}_{\mathbf{b}} |
1196 |
|
|
% +2 \mathbf{D}_{\mathbf{b}} \cdot \mathbf{Q}_{\mathbf{a}} \Bigr] w_g(r) |
1197 |
|
|
% % 2 |
1198 |
|
|
% + \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}} |
1199 |
|
|
% (\hat{r} \cdot \mathbf{D}_{\mathbf{b}}) |
1200 |
|
|
% +2 (\mathbf{D}_{\mathbf{b}} \cdot |
1201 |
|
|
% \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\ |
1202 |
|
|
% % 3 |
1203 |
|
|
% &+ \frac{1}{4\pi \epsilon_0} \Bigl[ \mathbf{D}_{\mathbf{b}} |
1204 |
|
|
% (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) |
1205 |
|
|
% +2 (\hat{r} \cdot \mathbf{D}_{\mathbf{b}}) |
1206 |
|
|
% (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \Bigr] w_i(r) |
1207 |
|
|
% % 4 |
1208 |
|
|
% +\frac{1}{4\pi \epsilon_0} |
1209 |
|
|
% (\hat{r} \cdot \mathbf{D}_{\mathbf{b}}) |
1210 |
|
|
% (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) w_j(r) \hat{r} |
1211 |
|
|
% \end{split} |
1212 |
|
|
% \end{equation} |
1213 |
gezelter |
3906 |
% |
1214 |
|
|
% |
1215 |
|
|
% |
1216 |
|
|
\begin{split} |
1217 |
gezelter |
3985 |
\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =& |
1218 |
|
|
\Bigl[ |
1219 |
gezelter |
3989 |
\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} |
1220 |
|
|
+ 2 \mathbf{Q}_{\mathbf{a}} : \mathbf{Q}_{\mathbf{b}} \Bigr] w_k(r) \hat{r} \\ |
1221 |
gezelter |
3906 |
% 2 |
1222 |
gezelter |
3985 |
&+ \Bigl[ |
1223 |
gezelter |
3906 |
2\text{Tr}\mathbf{Q}_{\mathbf{b}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} ) |
1224 |
|
|
+ 2\text{Tr}\mathbf{Q}_{\mathbf{a}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} ) |
1225 |
|
|
% 3 |
1226 |
|
|
+4 (\mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) |
1227 |
|
|
+ 4(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}}) \Bigr] w_n(r) \\ |
1228 |
|
|
% 4 |
1229 |
gezelter |
3985 |
&+ \Bigl[ |
1230 |
gezelter |
3906 |
\text{Tr}\mathbf{Q}_{\mathbf{a}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) |
1231 |
|
|
+ \text{Tr}\mathbf{Q}_{\mathbf{b}} |
1232 |
|
|
(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) |
1233 |
|
|
% 5 |
1234 |
|
|
+4 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot |
1235 |
|
|
\mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\ |
1236 |
|
|
% |
1237 |
gezelter |
3985 |
&+ \Bigl[ |
1238 |
gezelter |
3906 |
+ 2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} ) |
1239 |
|
|
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) |
1240 |
|
|
%6 |
1241 |
|
|
+2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) |
1242 |
|
|
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_o(r) \\ |
1243 |
|
|
% 7 |
1244 |
gezelter |
3985 |
&+ |
1245 |
gezelter |
3906 |
(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) |
1246 |
gezelter |
3985 |
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) w_m(r) \hat{r} \end{split} |
1247 |
|
|
\end{align} |
1248 |
|
|
Note that the forces for higher multipoles on site $\mathbf{a}$ |
1249 |
|
|
interacting with those of lower order on site $\mathbf{b}$ can be |
1250 |
|
|
obtained by swapping indices in the expressions above. |
1251 |
|
|
|
1252 |
gezelter |
3906 |
% |
1253 |
gezelter |
3985 |
% Torques SECTION ----------------------------------------------------------------------------------------- |
1254 |
gezelter |
3906 |
% |
1255 |
|
|
\subsection{Torques} |
1256 |
gezelter |
3989 |
|
1257 |
gezelter |
3906 |
% |
1258 |
gezelter |
4183 |
The torques for the three methods are given in space-frame |
1259 |
|
|
coordinates: |
1260 |
gezelter |
3906 |
% |
1261 |
|
|
% |
1262 |
gezelter |
3985 |
\begin{align} |
1263 |
|
|
\mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =& |
1264 |
|
|
C_{\bf a} (\hat{r} \times \mathbf{D}_{\mathbf{b}}) v_{11}(r) \\ |
1265 |
gezelter |
3906 |
% |
1266 |
|
|
% |
1267 |
|
|
% |
1268 |
gezelter |
3985 |
\mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =& |
1269 |
|
|
2C_{\bf a} |
1270 |
|
|
\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{22}(r) \\ |
1271 |
gezelter |
3906 |
% |
1272 |
|
|
% |
1273 |
|
|
% |
1274 |
gezelter |
3985 |
% \begin{equation} |
1275 |
|
|
% \mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} = |
1276 |
|
|
% -\frac{C_{\bf b}}{4\pi \epsilon_0} |
1277 |
|
|
% (\hat{r} \times \mathbf{D}_{\mathbf{a}}) v_{11}(r) |
1278 |
|
|
% \end{equation} |
1279 |
gezelter |
3906 |
% |
1280 |
|
|
% |
1281 |
|
|
% |
1282 |
gezelter |
3985 |
\mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} =& |
1283 |
|
|
\mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r) |
1284 |
gezelter |
3906 |
% 2 |
1285 |
gezelter |
3985 |
- |
1286 |
gezelter |
3906 |
(\hat{r} \times \mathbf{D}_{\mathbf {a}} ) |
1287 |
gezelter |
3985 |
(\hat{r} \cdot \mathbf{D}_{\mathbf {b}} ) v_{22}(r)\\ |
1288 |
gezelter |
3906 |
% |
1289 |
|
|
% |
1290 |
|
|
% |
1291 |
gezelter |
3985 |
% \begin{equation} |
1292 |
|
|
% \mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} = |
1293 |
|
|
% -\frac{1}{4\pi \epsilon_0} \mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r) |
1294 |
|
|
% % 2 |
1295 |
|
|
% +\frac{1}{4\pi \epsilon_0} |
1296 |
|
|
% (\hat{r} \cdot \mathbf{D}_{\mathbf {a}} ) |
1297 |
|
|
% (\hat{r} \times \mathbf{D}_{\mathbf {b}} ) v_{22}(r) |
1298 |
|
|
% \end{equation} |
1299 |
gezelter |
3906 |
% |
1300 |
|
|
% |
1301 |
|
|
% |
1302 |
gezelter |
3985 |
\mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} =& |
1303 |
|
|
\Bigl[ |
1304 |
gezelter |
3906 |
-\text{Tr}\mathbf{Q}_{\mathbf{b}} |
1305 |
|
|
(\hat{r} \times \mathbf{D}_{\mathbf{a}} ) |
1306 |
|
|
+2 \mathbf{D}_{\mathbf{a}} \times |
1307 |
|
|
(\mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) |
1308 |
|
|
\Bigr] v_{31}(r) |
1309 |
|
|
% 3 |
1310 |
gezelter |
3985 |
- (\hat{r} \times \mathbf{D}_{\mathbf{a}} ) |
1311 |
|
|
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)\\ |
1312 |
gezelter |
3906 |
% |
1313 |
|
|
% |
1314 |
|
|
% |
1315 |
gezelter |
3985 |
\mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =& |
1316 |
|
|
\Bigl[ |
1317 |
gezelter |
3906 |
+2 ( \mathbf{D}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} ) \times |
1318 |
|
|
\hat{r} |
1319 |
|
|
-2 \mathbf{D}_{\mathbf{a}} \times |
1320 |
|
|
(\mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) |
1321 |
|
|
\Bigr] v_{31}(r) |
1322 |
|
|
% 2 |
1323 |
gezelter |
3985 |
+ |
1324 |
gezelter |
3906 |
(\hat{r} \cdot \mathbf{D}_{\mathbf{a}}) |
1325 |
gezelter |
3985 |
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)\\ |
1326 |
gezelter |
3906 |
% |
1327 |
|
|
% |
1328 |
|
|
% |
1329 |
gezelter |
3985 |
% \begin{equation} |
1330 |
|
|
% \mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} = |
1331 |
|
|
% \frac{1}{4\pi \epsilon_0} \Bigl[ |
1332 |
|
|
% -2 (\mathbf{D}_{\mathbf{b}} \cdot \mathbf{Q}_{\mathbf{a}} ) \times \hat{r} |
1333 |
|
|
% +2 \mathbf{D}_{\mathbf{b}} \times |
1334 |
|
|
% (\mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) |
1335 |
|
|
% \Bigr] v_{31}(r) |
1336 |
|
|
% % 3 |
1337 |
|
|
% - \frac{2}{4\pi \epsilon_0} |
1338 |
|
|
% (\hat{r} \cdot \mathbf{D}_{\mathbf{b}} ) |
1339 |
|
|
% (\hat{r} \cdot \mathbf |
1340 |
|
|
% {Q}_{{\mathbf a}}) \times \hat{r} v_{32}(r) |
1341 |
|
|
% \end{equation} |
1342 |
gezelter |
3906 |
% |
1343 |
|
|
% |
1344 |
|
|
% |
1345 |
gezelter |
3985 |
% \begin{equation} |
1346 |
|
|
% \mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} = |
1347 |
|
|
% \frac{1}{4\pi \epsilon_0} \Bigl[ |
1348 |
|
|
% \text{Tr}\mathbf{Q}_{\mathbf{a}} |
1349 |
|
|
% (\hat{r} \times \mathbf{D}_{\mathbf{b}} ) |
1350 |
|
|
% +2 \mathbf{D}_{\mathbf{b}} \times |
1351 |
|
|
% ( \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r) |
1352 |
|
|
% % 2 |
1353 |
|
|
% +\frac{1}{4\pi \epsilon_0} |
1354 |
|
|
% (\hat{r} \times \mathbf{D}_{\mathbf{b}} ) |
1355 |
|
|
% (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) v_{32}(r) |
1356 |
|
|
% \end{equation} |
1357 |
gezelter |
3906 |
% |
1358 |
|
|
% |
1359 |
|
|
% |
1360 |
|
|
\begin{split} |
1361 |
gezelter |
3985 |
\mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =& |
1362 |
|
|
-4 |
1363 |
gezelter |
3906 |
\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}} |
1364 |
|
|
v_{41}(r) \\ |
1365 |
|
|
% 2 |
1366 |
gezelter |
3985 |
&+ |
1367 |
gezelter |
3906 |
\Bigl[-2\text{Tr}\mathbf{Q}_{\mathbf{b}} |
1368 |
|
|
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times \hat{r} |
1369 |
|
|
+4 \hat{r} \times |
1370 |
|
|
( \mathbf{Q}_{{\mathbf a}} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) |
1371 |
|
|
% 3 |
1372 |
|
|
-4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} )\times |
1373 |
|
|
( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} ) \Bigr] v_{42}(r) \\ |
1374 |
|
|
% 4 |
1375 |
gezelter |
3985 |
&+ 2 |
1376 |
gezelter |
3906 |
\hat{r} \times ( \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) |
1377 |
gezelter |
3985 |
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r) \end{split}\\ |
1378 |
gezelter |
3906 |
% |
1379 |
|
|
% |
1380 |
|
|
% |
1381 |
|
|
\begin{split} |
1382 |
|
|
\mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} = |
1383 |
gezelter |
3985 |
&4 |
1384 |
gezelter |
3906 |
\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}} v_{41}(r) \\ |
1385 |
|
|
% 2 |
1386 |
gezelter |
3985 |
&+ \Bigl[- 2\text{Tr}\mathbf{Q}_{\mathbf{a}} |
1387 |
gezelter |
3906 |
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} ) \times \hat{r} |
1388 |
|
|
-4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot |
1389 |
|
|
\mathbf{Q}_{{\mathbf b}} ) \times |
1390 |
|
|
\hat{r} |
1391 |
|
|
+4 ( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times |
1392 |
|
|
( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) |
1393 |
|
|
\Bigr] v_{42}(r) \\ |
1394 |
|
|
% 4 |
1395 |
gezelter |
3985 |
&+2 |
1396 |
gezelter |
3906 |
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) |
1397 |
gezelter |
3985 |
\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)\end{split} |
1398 |
|
|
\end{align} |
1399 |
|
|
% |
1400 |
|
|
Here, we have defined the matrix cross product in an identical form |
1401 |
|
|
as in Ref. \onlinecite{Smith98}: |
1402 |
|
|
\begin{equation} |
1403 |
|
|
\left[\mathbf{A} \times \mathbf{B}\right]_\alpha = \sum_\beta |
1404 |
|
|
\left[\mathbf{A}_{\alpha+1,\beta} \mathbf{B}_{\alpha+2,\beta} |
1405 |
|
|
-\mathbf{A}_{\alpha+2,\beta} \mathbf{B}_{\alpha+2,\beta} |
1406 |
|
|
\right] |
1407 |
gezelter |
4179 |
\label{eq:matrixCross} |
1408 |
gezelter |
3906 |
\end{equation} |
1409 |
gezelter |
3985 |
where $\alpha+1$ and $\alpha+2$ are regarded as cyclic |
1410 |
|
|
permuations of the matrix indices. |
1411 |
gezelter |
3980 |
|
1412 |
gezelter |
3985 |
All of the radial functions required for torques are identical with |
1413 |
|
|
the radial functions previously computed for the interaction energies. |
1414 |
gezelter |
4183 |
These are tabulated for all three methods in table |
1415 |
gezelter |
3985 |
\ref{tab:tableenergy}. The torques for higher multipoles on site |
1416 |
|
|
$\mathbf{a}$ interacting with those of lower order on site |
1417 |
|
|
$\mathbf{b}$ can be obtained by swapping indices in the expressions |
1418 |
|
|
above. |
1419 |
|
|
|
1420 |
gezelter |
3980 |
\section{Comparison to known multipolar energies} |
1421 |
|
|
|
1422 |
|
|
To understand how these new real-space multipole methods behave in |
1423 |
|
|
computer simulations, it is vital to test against established methods |
1424 |
|
|
for computing electrostatic interactions in periodic systems, and to |
1425 |
|
|
evaluate the size and sources of any errors that arise from the |
1426 |
gezelter |
4183 |
real-space cutoffs. In this paper we test SP, TSF, and GSF |
1427 |
gezelter |
3990 |
electrostatics against analytical methods for computing the energies |
1428 |
|
|
of ordered multipolar arrays. In the following paper, we test the new |
1429 |
|
|
methods against the multipolar Ewald sum for computing the energies, |
1430 |
|
|
forces and torques for a wide range of typical condensed-phase |
1431 |
|
|
(disordered) systems. |
1432 |
gezelter |
3980 |
|
1433 |
|
|
Because long-range electrostatic effects can be significant in |
1434 |
|
|
crystalline materials, ordered multipolar arrays present one of the |
1435 |
|
|
biggest challenges for real-space cutoff methods. The dipolar |
1436 |
|
|
analogues to the Madelung constants were first worked out by Sauer, |
1437 |
|
|
who computed the energies of ordered dipole arrays of zero |
1438 |
|
|
magnetization and obtained a number of these constants.\cite{Sauer} |
1439 |
|
|
This theory was developed more completely by Luttinger and |
1440 |
gezelter |
3986 |
Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays |
1441 |
gezelter |
3990 |
and other periodic structures. |
1442 |
gezelter |
3986 |
|
1443 |
gezelter |
3990 |
To test the new electrostatic methods, we have constructed very large, |
1444 |
|
|
$N=$ 16,000~(bcc) arrays of dipoles in the orientations described in |
1445 |
|
|
Ref. \onlinecite{LT}. These structures include ``A'' lattices with |
1446 |
|
|
nearest neighbor chains of antiparallel dipoles, as well as ``B'' |
1447 |
|
|
lattices with nearest neighbor strings of antiparallel dipoles if the |
1448 |
|
|
dipoles are contained in a plane perpendicular to the dipole direction |
1449 |
|
|
that passes through the dipole. We have also studied the minimum |
1450 |
gezelter |
3980 |
energy structure for the BCC lattice that was found by Luttinger \& |
1451 |
gezelter |
4179 |
Tisza. The total electrostatic energy density for any of the arrays |
1452 |
|
|
is given by: |
1453 |
gezelter |
3980 |
\begin{equation} |
1454 |
|
|
E = C N^2 \mu^2 |
1455 |
|
|
\end{equation} |
1456 |
gezelter |
3990 |
where $C$ is the energy constant (equivalent to the Madelung |
1457 |
|
|
constant), $N$ is the number of dipoles per unit volume, and $\mu$ is |
1458 |
|
|
the strength of the dipole. Energy constants (converged to 1 part in |
1459 |
|
|
$10^9$) are given in the supplemental information. |
1460 |
gezelter |
3980 |
|
1461 |
gezelter |
4172 |
\begin{figure} |
1462 |
|
|
\includegraphics[width=\linewidth]{Dipoles_rCutNew.pdf} |
1463 |
|
|
\caption{Convergence of the lattice energy constants as a function of |
1464 |
|
|
cutoff radius (normalized by the lattice constant, $a$) for the new |
1465 |
|
|
real-space methods. Three dipolar crystal structures were sampled, |
1466 |
|
|
and the analytic energy constants for the three lattices are |
1467 |
|
|
indicated with grey dashed lines. The left panel shows results for |
1468 |
|
|
the undamped kernel ($1/r$), while the damped error function kernel, |
1469 |
|
|
$B_0(r)$ was used in the right panel.} |
1470 |
|
|
\label{fig:Dipoles_rCut} |
1471 |
|
|
\end{figure} |
1472 |
mlamichh |
4163 |
|
1473 |
gezelter |
4172 |
\begin{figure} |
1474 |
|
|
\includegraphics[width=\linewidth]{Dipoles_alphaNew.pdf} |
1475 |
|
|
\caption{Convergence to the lattice energy constants as a function of |
1476 |
|
|
the reduced damping parameter ($\alpha^* = \alpha a$) for the |
1477 |
|
|
different real-space methods in the same three dipolar crystals in |
1478 |
|
|
Figure \ref{fig:Dipoles_rCut}. The left panel shows results for a |
1479 |
|
|
relatively small cutoff radius ($r_c = 4.5 a$) while a larger cutoff |
1480 |
|
|
radius ($r_c = 6 a$) was used in the right panel. } |
1481 |
|
|
\label{fig:Dipoles_alpha} |
1482 |
mlamichh |
4163 |
\end{figure} |
1483 |
gezelter |
4172 |
|
1484 |
gezelter |
3990 |
For the purposes of testing the energy expressions and the |
1485 |
|
|
self-neutralization schemes, the primary quantity of interest is the |
1486 |
|
|
analytic energy constant for the perfect arrays. Convergence to these |
1487 |
|
|
constants are shown as a function of both the cutoff radius, $r_c$, |
1488 |
gezelter |
4174 |
and the damping parameter, $\alpha$ in Figs.\ref{fig:Dipoles_rCut} |
1489 |
gezelter |
4172 |
and \ref{fig:Dipoles_alpha}. We have simultaneously tested a hard |
1490 |
gezelter |
4183 |
cutoff (where the kernel is simply truncated at the cutoff radius) in |
1491 |
|
|
addition to the three new methods. |
1492 |
gezelter |
3986 |
|
1493 |
gezelter |
4174 |
The hard cutoff exhibits oscillations around the analytic energy |
1494 |
gezelter |
3988 |
constants, and converges to incorrect energies when the complementary |
1495 |
gezelter |
4179 |
error function damping kernel is used. The shifted potential (SP) |
1496 |
|
|
converges to the correct energy smoothly by $r_c = 4.5 a$ even for the |
1497 |
|
|
undamped case. This indicates that the shifting and the correction |
1498 |
|
|
provided by the self term are required for obtaining accurate energies. |
1499 |
|
|
The Taylor-shifted force (TSF) approximation appears to perturb the |
1500 |
|
|
potential too much inside the cutoff region to provide accurate |
1501 |
|
|
measures of the energy constants. GSF is a compromise, converging to |
1502 |
|
|
the correct energies within $r_c = 6 a$. |
1503 |
gezelter |
4174 |
|
1504 |
gezelter |
3986 |
{\it Quadrupolar} analogues to the Madelung constants were first |
1505 |
|
|
worked out by Nagai and Nakamura who computed the energies of selected |
1506 |
|
|
quadrupole arrays based on extensions to the Luttinger and Tisza |
1507 |
gezelter |
4174 |
approach.\cite{Nagai01081960,Nagai01091963} |
1508 |
gezelter |
3980 |
|
1509 |
|
|
In analogy to the dipolar arrays, the total electrostatic energy for |
1510 |
|
|
the quadrupolar arrays is: |
1511 |
|
|
\begin{equation} |
1512 |
gezelter |
3996 |
E = C N \frac{3\bar{Q}^2}{4a^5} |
1513 |
gezelter |
3980 |
\end{equation} |
1514 |
gezelter |
3996 |
where $a$ is the lattice parameter, and $\bar{Q}$ is the effective |
1515 |
|
|
quadrupole moment, |
1516 |
|
|
\begin{equation} |
1517 |
gezelter |
4098 |
\bar{Q}^2 = 2 \left(3 Q : Q - (\text{Tr} Q)^2 \right) |
1518 |
gezelter |
3996 |
\end{equation} |
1519 |
|
|
for the primitive quadrupole as defined in Eq. \ref{eq:quadrupole}. |
1520 |
|
|
(For the traceless quadrupole tensor, $\Theta = 3 Q - \text{Tr} Q$, |
1521 |
gezelter |
4098 |
the effective moment, $\bar{Q}^2 = \frac{2}{3} \Theta : \Theta$.) |
1522 |
gezelter |
3980 |
|
1523 |
gezelter |
4174 |
To test the new electrostatic methods for quadrupoles, we have |
1524 |
|
|
constructed very large, $N=$ 8,000~(sc), 16,000~(bcc), and |
1525 |
|
|
32,000~(fcc) arrays of linear quadrupoles in the orientations |
1526 |
|
|
described in Ref. \onlinecite{Nagai01081960}. We have compared the |
1527 |
gezelter |
4179 |
energy constants for these low-energy configurations for linear |
1528 |
|
|
quadrupoles. Convergence to these constants are shown as a function of |
1529 |
|
|
both the cutoff radius, $r_c$, and the damping parameter, $\alpha$ in |
1530 |
|
|
Figs.~\ref{fig:Quadrupoles_rCut} and \ref{fig:Quadrupoles_alpha}. |
1531 |
gezelter |
4174 |
|
1532 |
gezelter |
4172 |
\begin{figure} |
1533 |
gezelter |
4173 |
\includegraphics[width=\linewidth]{Quadrupoles_rcutConvergence.pdf} |
1534 |
|
|
\caption{Convergence of the lattice energy constants as a function of |
1535 |
|
|
cutoff radius (normalized by the lattice constant, $a$) for the new |
1536 |
|
|
real-space methods. Three quadrupolar crystal structures were |
1537 |
|
|
sampled, and the analytic energy constants for the three lattices |
1538 |
|
|
are indicated with grey dashed lines. The left panel shows results |
1539 |
|
|
for the undamped kernel ($1/r$), while the damped error function |
1540 |
|
|
kernel, $B_0(r)$ was used in the right panel.} |
1541 |
gezelter |
4174 |
\label{fig:Quadrupoles_rCut} |
1542 |
gezelter |
4172 |
\end{figure} |
1543 |
|
|
|
1544 |
|
|
|
1545 |
|
|
\begin{figure}[!htbp] |
1546 |
gezelter |
4174 |
\includegraphics[width=\linewidth]{Quadrupoles_newAlpha.pdf} |
1547 |
|
|
\caption{Convergence to the lattice energy constants as a function of |
1548 |
|
|
the reduced damping parameter ($\alpha^* = \alpha a$) for the |
1549 |
|
|
different real-space methods in the same three quadrupolar crystals |
1550 |
|
|
in Figure \ref{fig:Quadrupoles_rCut}. The left panel shows |
1551 |
|
|
results for a relatively small cutoff radius ($r_c = 4.5 a$) while a |
1552 |
|
|
larger cutoff radius ($r_c = 6 a$) was used in the right panel. } |
1553 |
|
|
\label{fig:Quadrupoles_alpha} |
1554 |
gezelter |
4172 |
\end{figure} |
1555 |
|
|
|
1556 |
gezelter |
4174 |
Again, we find that the hard cutoff exhibits oscillations around the |
1557 |
gezelter |
4179 |
analytic energy constants. The shifted potential (SP) approximation |
1558 |
|
|
converges to the correct energy smoothly by $r_c = 3 a$ even for the |
1559 |
|
|
undamped case. The Taylor-shifted force (TSF) approximation again |
1560 |
|
|
appears to perturb the potential too much inside the cutoff region to |
1561 |
|
|
provide accurate measures of the energy constants. GSF again provides |
1562 |
|
|
a compromise between the two methods -- energies are converged by $r_c |
1563 |
|
|
= 4.5 a$, and the approximation is not as perturbative at short range |
1564 |
|
|
as TSF. |
1565 |
gezelter |
4172 |
|
1566 |
gezelter |
4179 |
It is also useful to understand the convergence to the lattice energy |
1567 |
|
|
constants as a function of the reduced damping parameter ($\alpha^* = |
1568 |
|
|
\alpha a$) for the different real-space methods. |
1569 |
|
|
Figures. \ref{fig:Dipoles_alpha} and \ref{fig:Quadrupoles_alpha} show |
1570 |
|
|
this comparison for the dipolar and quadrupolar lattices, |
1571 |
|
|
respectively. All of the methods (except for TSF) have excellent |
1572 |
|
|
behavior for the undamped or weakly-damped cases. All of the methods |
1573 |
|
|
can be forced to converge by increasing the value of $\alpha$, which |
1574 |
gezelter |
4183 |
effectively shortens the overall range of the potential by equalizing |
1575 |
|
|
the truncation effects on the different orientational contributions. |
1576 |
|
|
In the second paper in the series, we discuss how large values of |
1577 |
|
|
$\alpha$ can perturb the force and torque vectors, but both |
1578 |
|
|
weakly-damped or over-damped electrostatics appear to generate |
1579 |
gezelter |
4179 |
reasonable values for the total electrostatic energies under both the |
1580 |
|
|
SP and GSF approximations. |
1581 |
gezelter |
4174 |
|
1582 |
gezelter |
3985 |
\section{Conclusion} |
1583 |
gezelter |
4174 |
We have presented three efficient real-space methods for computing the |
1584 |
gezelter |
4181 |
interactions between point multipoles. One of these (SP) is a |
1585 |
|
|
multipolar generalization of Wolf's method that smoothly shifts |
1586 |
|
|
electrostatic energies to zero at the cutoff radius. Two of these |
1587 |
|
|
methods (GSF and TSF) also smoothly truncate the forces and torques |
1588 |
|
|
(in addition to the energies) at the cutoff radius, making them |
1589 |
|
|
attractive for both molecular dynamics and Monte Carlo simulations. We |
1590 |
|
|
find that the Gradient-Shifted Force (GSF) and the Shifted-Potential |
1591 |
|
|
(SP) methods converge rapidly to the correct lattice energies for |
1592 |
|
|
ordered dipolar and quadrupolar arrays, while the Taylor-Shifted Force |
1593 |
|
|
(TSF) is too severe an approximation to provide accurate convergence |
1594 |
|
|
to lattice energies. |
1595 |
gezelter |
3980 |
|
1596 |
gezelter |
3985 |
In most cases, GSF can obtain nearly quantitative agreement with the |
1597 |
|
|
lattice energy constants with reasonably small cutoff radii. The only |
1598 |
|
|
exception we have observed is for crystals which exhibit a bulk |
1599 |
|
|
macroscopic dipole moment (e.g. Luttinger \& Tisza's $Z_1$ lattice). |
1600 |
|
|
In this particular case, the multipole neutralization scheme can |
1601 |
|
|
interfere with the correct computation of the energies. We note that |
1602 |
|
|
the energies for these arrangements are typically much larger than for |
1603 |
|
|
crystals with net-zero moments, so this is not expected to be an issue |
1604 |
|
|
in most simulations. |
1605 |
gezelter |
3980 |
|
1606 |
gezelter |
4179 |
The techniques used here to derive the force, torque and energy |
1607 |
|
|
expressions can be extended to higher order multipoles, although some |
1608 |
|
|
of the objects (e.g. the matrix cross product in |
1609 |
|
|
Eq. \ref{eq:matrixCross}) will need to be generalized for higher-rank |
1610 |
|
|
tensors. We also note that the definitions of the multipoles used |
1611 |
|
|
here are in a primitive form, and these need some care when comparing |
1612 |
|
|
with experiment or other computational techniques. |
1613 |
|
|
|
1614 |
gezelter |
3985 |
In large systems, these new methods can be made to scale approximately |
1615 |
|
|
linearly with system size, and detailed comparisons with the Ewald sum |
1616 |
|
|
for a wide range of chemical environments follows in the second paper. |
1617 |
gezelter |
3980 |
|
1618 |
gezelter |
3906 |
\begin{acknowledgments} |
1619 |
gezelter |
3985 |
JDG acknowledges helpful discussions with Christopher |
1620 |
|
|
Fennell. Support for this project was provided by the National |
1621 |
gezelter |
4172 |
Science Foundation under grant CHE-1362211. Computational time was |
1622 |
gezelter |
3985 |
provided by the Center for Research Computing (CRC) at the |
1623 |
|
|
University of Notre Dame. |
1624 |
gezelter |
3906 |
\end{acknowledgments} |
1625 |
|
|
|
1626 |
gezelter |
3984 |
\newpage |
1627 |
gezelter |
3906 |
\appendix |
1628 |
|
|
|
1629 |
gezelter |
3984 |
\section{Smith's $B_l(r)$ functions for damped-charge distributions} |
1630 |
gezelter |
3985 |
\label{SmithFunc} |
1631 |
gezelter |
3984 |
The following summarizes Smith's $B_l(r)$ functions and includes |
1632 |
|
|
formulas given in his appendix.\cite{Smith98} The first function |
1633 |
|
|
$B_0(r)$ is defined by |
1634 |
gezelter |
3906 |
% |
1635 |
|
|
\begin{equation} |
1636 |
|
|
B_0(r)=\frac{\text{erfc}(\alpha r)}{r} = \frac{2}{\sqrt{\pi}r}= |
1637 |
|
|
\int_{\alpha r}^{\infty} \text{e}^{-s^2} ds . |
1638 |
|
|
\end{equation} |
1639 |
|
|
% |
1640 |
|
|
The first derivative of this function is |
1641 |
|
|
% |
1642 |
|
|
\begin{equation} |
1643 |
|
|
\frac{dB_0(r)}{dr}=-\frac{1}{r^2}\text{erfc}(\alpha r) |
1644 |
|
|
-\frac{2\alpha}{r\sqrt{\pi}}\text{e}^{-{\alpha}^2r^2} |
1645 |
|
|
\end{equation} |
1646 |
|
|
% |
1647 |
gezelter |
3984 |
which can be used to define a function $B_1(r)$: |
1648 |
gezelter |
3906 |
% |
1649 |
|
|
\begin{equation} |
1650 |
|
|
B_1(r)=-\frac{1}{r}\frac{dB_0(r)}{dr} |
1651 |
|
|
\end{equation} |
1652 |
|
|
% |
1653 |
gezelter |
3984 |
In general, the recurrence relation, |
1654 |
gezelter |
3906 |
\begin{equation} |
1655 |
|
|
B_l(r)=-\frac{1}{r}\frac{dB_{l-1}(r)}{dr} |
1656 |
|
|
= \frac{1}{r^2} \left[ (2l-1)B_{l-1}(r) + \frac {(2\alpha^2)^l}{\alpha \sqrt{\pi}} |
1657 |
|
|
\text{e}^{-{\alpha}^2r^2} |
1658 |
gezelter |
3984 |
\right] , |
1659 |
gezelter |
3906 |
\end{equation} |
1660 |
gezelter |
4183 |
is very useful for building up higher derivatives. As noted by Smith, |
1661 |
|
|
it is possible to approximate the $B_l(r)$ functions, |
1662 |
gezelter |
3906 |
% |
1663 |
|
|
\begin{equation} |
1664 |
|
|
B_l(r)=\frac{(2l)!}{l!2^lr^{2l+1}} - \frac {(2\alpha^2)^{l+1}}{(2l+1)\alpha \sqrt{\pi}} |
1665 |
|
|
+\text{O}(r) . |
1666 |
|
|
\end{equation} |
1667 |
gezelter |
3984 |
\newpage |
1668 |
|
|
\section{The $r$-dependent factors for TSF electrostatics} |
1669 |
gezelter |
4183 |
\label{radialTSF} |
1670 |
gezelter |
3906 |
|
1671 |
|
|
Using the shifted damped functions $f_n(r)$ defined by: |
1672 |
|
|
% |
1673 |
|
|
\begin{equation} |
1674 |
gezelter |
3984 |
f_n(r)= B_0(r) -\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} B_0^{(m)}(r_c) , |
1675 |
gezelter |
3906 |
\end{equation} |
1676 |
|
|
% |
1677 |
gezelter |
3984 |
where the superscript $(m)$ denotes the $m^\mathrm{th}$ derivative. In |
1678 |
|
|
this Appendix, we provide formulas for successive derivatives of this |
1679 |
|
|
function. (If there is no damping, then $B_0(r)$ is replaced by |
1680 |
|
|
$1/r$.) First, we find: |
1681 |
gezelter |
3906 |
% |
1682 |
|
|
\begin{equation} |
1683 |
|
|
\frac{\partial f_n}{\partial r_\alpha}=\hat{r}_\alpha \frac{d f_n}{d r} . |
1684 |
|
|
\end{equation} |
1685 |
|
|
% |
1686 |
gezelter |
3984 |
This formula clearly brings in derivatives of Smith's $B_0(r)$ |
1687 |
|
|
function, and we define higher-order derivatives as follows: |
1688 |
gezelter |
3906 |
% |
1689 |
gezelter |
3984 |
\begin{align} |
1690 |
|
|
g_n(r)=& \frac{d f_n}{d r} = |
1691 |
|
|
B_0^{(1)}(r) -\sum_{m=0}^{n} \frac {(r-r_c)^m}{m!} B_0^{(m+1)}(r_c) \\ |
1692 |
|
|
h_n(r)=& \frac{d^2f_n}{d r^2} = |
1693 |
|
|
B_0^{(2)}(r) -\sum_{m=0}^{n-1} \frac {(r-r_c)^m}{m!} B_0^{(m+2)}(r_c) \\ |
1694 |
|
|
s_n(r)=& \frac{d^3f_n}{d r^3} = |
1695 |
|
|
B_0^{(3)}(r) -\sum_{m=0}^{n-2} \frac {(r-r_c)^m}{m!} B_0^{(m+3)}(r_c) \\ |
1696 |
|
|
t_n(r)=& \frac{d^4f_n}{d r^4} = |
1697 |
|
|
B_0^{(4)}(r) -\sum_{m=0}^{n-3} \frac {(r-r_c)^m}{m!} B_0^{(m+4)}(r_c) \\ |
1698 |
|
|
u_n(r)=& \frac{d^5f_n}{d r^5} = |
1699 |
|
|
B_0^{(5)}(r) -\sum_{m=0}^{n-4} \frac {(r-r_c)^m}{m!} B_0^{(m+5)}(r_c) . |
1700 |
|
|
\end{align} |
1701 |
gezelter |
3906 |
% |
1702 |
gezelter |
3984 |
We note that the last function needed (for quadrupole-quadrupole interactions) is |
1703 |
gezelter |
3906 |
% |
1704 |
|
|
\begin{equation} |
1705 |
gezelter |
3984 |
u_4(r)=B_0^{(5)}(r) - B_0^{(5)}(r_c) . |
1706 |
gezelter |
3906 |
\end{equation} |
1707 |
gezelter |
3989 |
% The functions |
1708 |
|
|
% needed are listed schematically below: |
1709 |
|
|
% % |
1710 |
|
|
% \begin{eqnarray} |
1711 |
|
|
% f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\ |
1712 |
|
|
% g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\ |
1713 |
|
|
% h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\ |
1714 |
|
|
% s_2 \quad s_3 \quad &s_4 \nonumber \\ |
1715 |
|
|
% t_3 \quad &t_4 \nonumber \\ |
1716 |
|
|
% &u_4 \nonumber . |
1717 |
|
|
% \end{eqnarray} |
1718 |
gezelter |
3984 |
The functions $f_n(r)$ to $u_n(r)$ can be computed recursively and |
1719 |
gezelter |
3989 |
stored on a grid for values of $r$ from $0$ to $r_c$. Using these |
1720 |
|
|
functions, we find |
1721 |
gezelter |
3906 |
% |
1722 |
gezelter |
3984 |
\begin{align} |
1723 |
|
|
\frac{\partial f_n}{\partial r_\alpha} =&r_\alpha \frac {g_n}{r} \label{eq:b9}\\ |
1724 |
|
|
\frac{\partial^2 f_n}{\partial r_\alpha \partial r_\beta} =&\delta_{\alpha \beta}\frac {g_n}{r} |
1725 |
|
|
+r_\alpha r_\beta \left( -\frac{g_n}{r^3} +\frac{h_n}{r^2}\right) \\ |
1726 |
gezelter |
3989 |
\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta \partial r_\gamma} =& |
1727 |
gezelter |
3906 |
\left( \delta_{\alpha \beta} r_\gamma + \delta_{\alpha \gamma} r_\beta + |
1728 |
|
|
\delta_{ \beta \gamma} r_\alpha \right) |
1729 |
gezelter |
3989 |
\left( -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right) \nonumber \\ |
1730 |
|
|
& + r_\alpha r_\beta r_\gamma |
1731 |
gezelter |
3984 |
\left( \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \\ |
1732 |
gezelter |
3989 |
\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta \partial |
1733 |
|
|
r_\gamma \partial r_\delta} =& |
1734 |
gezelter |
3906 |
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} |
1735 |
|
|
+ \delta_{\alpha \gamma} \delta_{\beta \delta} |
1736 |
|
|
+\delta_{ \beta \gamma} \delta_{\alpha \delta} \right) |
1737 |
|
|
\left( - \frac{g_n}{r^3} + \frac{h_n}{r^2} \right) \nonumber \\ |
1738 |
gezelter |
3984 |
&+ \left( \delta_{\alpha \beta} r_\gamma r_\delta |
1739 |
|
|
+ \text{5 permutations} |
1740 |
gezelter |
3906 |
\right) \left( \frac{3 g_n}{r^5} - \frac{3h_n}{r^4} + \frac{s_n}{r^3} |
1741 |
|
|
\right) \nonumber \\ |
1742 |
gezelter |
3984 |
&+ r_\alpha r_\beta r_\gamma r_\delta |
1743 |
gezelter |
3906 |
\left( -\frac{15g_n}{r^7} + \frac{15h_n}{r^6} - \frac{6s_n}{r^5} |
1744 |
gezelter |
3984 |
+ \frac{t_n}{r^4} \right)\\ |
1745 |
gezelter |
3906 |
\frac{\partial^5 f_n} |
1746 |
gezelter |
3989 |
{\partial r_\alpha \partial r_\beta \partial r_\gamma \partial |
1747 |
|
|
r_\delta \partial r_\epsilon} =& |
1748 |
gezelter |
3906 |
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} r_\epsilon |
1749 |
gezelter |
3984 |
+ \text{14 permutations} \right) |
1750 |
gezelter |
3906 |
\left( \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \nonumber \\ |
1751 |
gezelter |
3984 |
&+ \left( \delta_{\alpha \beta} r_\gamma r_\delta r_\epsilon |
1752 |
|
|
+ \text{9 permutations} |
1753 |
gezelter |
3906 |
\right) \left(- \frac{15g_n}{r^7}+\frac{15h_n}{r^7} -\frac{6s_n}{r^5} +\frac{t_n}{r^4} |
1754 |
|
|
\right) \nonumber \\ |
1755 |
gezelter |
3984 |
&+ r_\alpha r_\beta r_\gamma r_\delta r_\epsilon |
1756 |
gezelter |
3906 |
\left( \frac{105g_n}{r^9} - \frac{105h_n}{r^8} + \frac{45s_n}{r^7} |
1757 |
gezelter |
3984 |
- \frac{10t_n}{r^6} +\frac{u_n}{r^5} \right) \label{eq:b13} |
1758 |
|
|
\end{align} |
1759 |
gezelter |
3906 |
% |
1760 |
|
|
% |
1761 |
|
|
% |
1762 |
gezelter |
3984 |
\newpage |
1763 |
|
|
\section{The $r$-dependent factors for GSF electrostatics} |
1764 |
gezelter |
4183 |
\label{radialGSF} |
1765 |
gezelter |
3906 |
|
1766 |
gezelter |
3984 |
In Gradient-shifted force electrostatics, the kernel is not expanded, |
1767 |
gezelter |
4183 |
and the expansion is carried out on the individual terms in the |
1768 |
|
|
multipole interaction energies. For damped charges, this still brings |
1769 |
|
|
multiple derivatives of the Smith's $B_0(r)$ function into the |
1770 |
|
|
algebra. To denote these terms, we generalize the notation of the |
1771 |
|
|
previous appendix. For either $f(r)=1/r$ (undamped) or $f(r)=B_0(r)$ |
1772 |
|
|
(damped), |
1773 |
gezelter |
3906 |
% |
1774 |
gezelter |
3984 |
\begin{align} |
1775 |
gezelter |
4179 |
g(r) &= \frac{df}{d r} && &&=-\frac{1}{r^2} |
1776 |
|
|
&&\mathrm{or~~~} -rB_1(r) \\ |
1777 |
|
|
h(r) &= \frac{dg}{d r} &&= \frac{d^2f}{d r^2} &&= \frac{2}{r^3} &&\mathrm{or~~~}-B_1(r) + r^2 B_2(r) \\ |
1778 |
|
|
s(r) &= \frac{dh}{d r} &&= \frac{d^3f}{d r^3} &&=-\frac{6}{r^4}&&\mathrm{or~~~}3rB_2(r) - r^3 B_3(r)\\ |
1779 |
|
|
t(r) &= \frac{ds}{d r} &&= \frac{d^4f}{d r^4} &&= \frac{24}{r^5} &&\mathrm{or~~~} 3 |
1780 |
|
|
B_2(r) - 6r^2 B_3(r) + r^4 B_4(r) \\ |
1781 |
|
|
u(r) &= \frac{dt}{d r} &&= \frac{d^5f}{d r^5} &&=-\frac{120}{r^6} &&\mathrm{or~~~} -15 |
1782 |
|
|
r B_3(r) + 10 r^3B_4(r) -r^5B_5(r). |
1783 |
gezelter |
3984 |
\end{align} |
1784 |
gezelter |
3906 |
% |
1785 |
gezelter |
4179 |
For undamped charges, Table I lists these derivatives under the Bare |
1786 |
|
|
Coulomb column. Equations \ref{eq:b9} to \ref{eq:b13} are still |
1787 |
gezelter |
3989 |
correct for GSF electrostatics if the subscript $n$ is eliminated. |
1788 |
gezelter |
3906 |
|
1789 |
gezelter |
4172 |
\newpage |
1790 |
|
|
|
1791 |
gezelter |
3980 |
\bibliography{multipole} |
1792 |
|
|
|
1793 |
gezelter |
3906 |
\end{document} |
1794 |
|
|
% |
1795 |
|
|
% ****** End of file multipole.tex ****** |