| 44 |
|
%\preprint{AIP/123-QED} |
| 45 |
|
|
| 46 |
|
\title{Real space alternatives to the Ewald |
| 47 |
< |
Sum. I. Taylor-shifted and Gradient-shifted electrostatics for multipoles} |
| 47 |
> |
Sum. I. Shifted electrostatics for multipoles} |
| 48 |
|
|
| 49 |
|
\author{Madan Lamichhane} |
| 50 |
|
\affiliation{Department of Physics, University |
| 65 |
|
|
| 66 |
|
\begin{abstract} |
| 67 |
|
We have extended the original damped-shifted force (DSF) |
| 68 |
< |
electrostatic kernel and have been able to derive two new |
| 68 |
> |
electrostatic kernel and have been able to derive three new |
| 69 |
|
electrostatic potentials for higher-order multipoles that are based |
| 70 |
< |
on truncated Taylor expansions around the cutoff radius. For |
| 71 |
< |
multipole-multipole interactions, we find that each of the distinct |
| 72 |
< |
orientational contributions has a separate radial function to ensure |
| 73 |
< |
that the overall forces and torques vanish at the cutoff radius. In |
| 74 |
< |
this paper, we present energy, force, and torque expressions for the |
| 75 |
< |
new models, and compare these real-space interaction models to exact |
| 76 |
< |
results for ordered arrays of multipoles. |
| 70 |
> |
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 |
> |
size, SP and GSF are attractive options for large Monte |
| 86 |
> |
Carlo and molecular dynamics simulations, respectively. |
| 87 |
|
\end{abstract} |
| 88 |
|
|
| 89 |
|
%\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy |
| 143 |
|
a different orientation can cause energy discontinuities. |
| 144 |
|
|
| 145 |
|
This paper outlines an extension of the original DSF electrostatic |
| 146 |
< |
kernel to point multipoles. We describe two distinct real-space |
| 147 |
< |
interaction models for higher-order multipoles based on two truncated |
| 146 |
> |
kernel to point multipoles. We describe three distinct real-space |
| 147 |
> |
interaction models for higher-order multipoles based on truncated |
| 148 |
|
Taylor expansions that are carried out at the cutoff radius. We are |
| 149 |
< |
calling these models {\bf Taylor-shifted} and {\bf Gradient-shifted} |
| 149 |
> |
calling these models {\bf Taylor-shifted} (TSF), {\bf |
| 150 |
> |
gradient-shifted} (GSF) and {\bf shifted potential} (SP) |
| 151 |
|
electrostatics. Because of differences in the initial assumptions, |
| 152 |
< |
the two methods yield related, but somewhat different expressions for |
| 153 |
< |
energies, forces, and torques. |
| 152 |
> |
the two methods yield related, but distinct expressions for energies, |
| 153 |
> |
forces, and torques. |
| 154 |
|
|
| 155 |
|
In this paper we outline the new methodology and give functional forms |
| 156 |
|
for the energies, forces, and torques up to quadrupole-quadrupole |
| 163 |
|
An efficient real-space electrostatic method involves the use of a |
| 164 |
|
pair-wise functional form, |
| 165 |
|
\begin{equation} |
| 166 |
< |
V = \sum_i \sum_{j>i} V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) + |
| 167 |
< |
\sum_i V_i^\mathrm{correction} |
| 166 |
> |
U = \sum_i \sum_{j>i} U_\mathrm{pair}(\mathbf{r}_{ij}, \Omega_i, \Omega_j) + |
| 167 |
> |
\sum_i U_i^\mathrm{self} |
| 168 |
|
\end{equation} |
| 169 |
|
that is short-ranged and easily truncated at a cutoff radius, |
| 170 |
|
\begin{equation} |
| 171 |
< |
V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) = \left\{ |
| 171 |
> |
U_\mathrm{pair}(\mathbf{r}_{ij},\Omega_i, \Omega_j) = \left\{ |
| 172 |
|
\begin{array}{ll} |
| 173 |
< |
V_\mathrm{approx} (r_{ij}, \Omega_i, \Omega_j) & \quad r \le r_c \\ |
| 174 |
< |
0 & \quad r > r_c , |
| 173 |
> |
U_\mathrm{approx} (\mathbf{r}_{ij}, \Omega_i, \Omega_j) & \quad \left| \mathbf{r}_{ij} \right| \le r_c \\ |
| 174 |
> |
0 & \quad \left| \mathbf{r}_{ij} \right| > r_c , |
| 175 |
|
\end{array} |
| 176 |
|
\right. |
| 177 |
|
\end{equation} |
| 178 |
< |
along with an easily computed correction term ($\sum_i |
| 179 |
< |
V_i^\mathrm{correction}$) which has linear-scaling with the number of |
| 178 |
> |
along with an easily computed self-interaction term ($\sum_i |
| 179 |
> |
U_i^\mathrm{self}$) which scales linearly with the number of |
| 180 |
|
particles. Here $\Omega_i$ and $\Omega_j$ represent orientational |
| 181 |
< |
coordinates of the two sites. The computational efficiency, energy |
| 181 |
> |
coordinates of the two sites, and $\mathbf{r}_{ij}$ is the vector |
| 182 |
> |
between the two sites. The computational efficiency, energy |
| 183 |
|
conservation, and even some physical properties of a simulation can |
| 184 |
< |
depend dramatically on how the $V_\mathrm{approx}$ function behaves at |
| 184 |
> |
depend dramatically on how the $U_\mathrm{approx}$ function behaves at |
| 185 |
|
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 |
| 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 |
< |
computed within $r_c$. Damping using a complementary error |
| 196 |
< |
function is applied to the potential to accelerate convergence. The |
| 197 |
< |
potential for the DSF method, where $\alpha$ is the adjustable damping |
| 186 |
< |
parameter, is given by |
| 195 |
> |
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 |
|
\begin{equation*} |
| 199 |
< |
V_\mathrm{DSF}(r) = C_i C_j \Biggr{[} |
| 199 |
> |
U_\mathrm{DSF}(r_{ij}) = C_i C_j \Biggr{[} |
| 200 |
|
\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} |
| 201 |
|
- \frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c} + \left(\frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c^2} |
| 202 |
|
+ \frac{2\alpha}{\pi^{1/2}} |
| 204 |
|
\right)\left(r_{ij}-r_c\right)\ \Biggr{]} |
| 205 |
|
\label{eq:DSFPot} |
| 206 |
|
\end{equation*} |
| 207 |
< |
Note that in this potential and in all electrostatic quantities that |
| 208 |
< |
follow, the standard $1/4 \pi \epsilon_{0}$ has been omitted for |
| 209 |
< |
clarity. |
| 207 |
> |
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 |
|
|
| 211 |
|
To insure net charge neutrality within each cutoff sphere, an |
| 212 |
|
additional ``self'' term is added to the potential. This term is |
| 255 |
|
a$. Then the electrostatic potential of object $\bf a$ at |
| 256 |
|
$\mathbf{r}$ is given by |
| 257 |
|
\begin{equation} |
| 258 |
< |
V_a(\mathbf r) = |
| 258 |
> |
\phi_a(\mathbf r) = |
| 259 |
|
\sum_{k \, \text{in \bf a}} \frac{q_k}{\lvert \mathbf{r} - \mathbf{r}_k \rvert}. |
| 260 |
|
\end{equation} |
| 261 |
|
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 |
< |
labelling specific charges in $\bf a$ and $\bf b$ respectively. The |
| 264 |
> |
labeling specific charges in $\bf a$ and $\bf b$ respectively. The |
| 265 |
|
Taylor expansion, |
| 266 |
|
\begin{equation} |
| 267 |
|
\frac{1}{\lvert \mathbf{r} - \mathbf{r}_k \rvert} = |
| 274 |
|
can then be used to express the electrostatic potential on $\bf a$ in |
| 275 |
|
terms of multipole operators, |
| 276 |
|
\begin{equation} |
| 277 |
< |
V_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r} |
| 277 |
> |
\phi_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r} |
| 278 |
|
\end{equation} |
| 279 |
|
where |
| 280 |
|
\begin{equation} |
| 287 |
|
a}\alpha\beta}$, respectively. These are the primitive multipoles |
| 288 |
|
which can be expressed as a distribution of charges, |
| 289 |
|
\begin{align} |
| 290 |
< |
C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \\ |
| 291 |
< |
D_{{\bf a}\alpha} =&\sum_{k \, \text{in \bf a}} q_k r_{k\alpha} ,\\ |
| 292 |
< |
Q_{{\bf a}\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k r_{k\alpha} r_{k\beta} . |
| 290 |
> |
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 |
|
\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 |
< |
based factor of $1/2$. |
| 297 |
> |
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 |
|
|
| 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 |
< |
$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $), the interaction energy is given by |
| 306 |
> |
$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $), the interaction energy is given by |
| 307 |
|
\begin{equation} |
| 308 |
|
U_{\bf{ab}}(r) |
| 309 |
|
= \hat{M}_a \sum_{j \, \text{in \bf b}} \frac {q_j}{\vert \bf{r}+\bf{r}_j \vert} . |
| 321 |
|
\end{equation} |
| 322 |
|
This form has the benefit of separating out the energies of |
| 323 |
|
interaction into contributions from the charge, dipole, and quadrupole |
| 324 |
< |
of $\bf a$ interacting with the same multipoles on $\bf b$. |
| 324 |
> |
of $\bf a$ interacting with the same multipoles in $\bf b$. |
| 325 |
|
|
| 326 |
|
\subsection{Damped Coulomb interactions} |
| 327 |
+ |
\label{sec:damped} |
| 328 |
|
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 |
| 340 |
|
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 |
< |
functions $B_l(r)$ are summarized in Appendix A. |
| 343 |
> |
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 |
> |
quadrupole tensor fails for the damped case.) |
| 351 |
|
|
| 352 |
|
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 |
| 444 |
|
In general, we can write |
| 445 |
|
% |
| 446 |
|
\begin{equation} |
| 447 |
< |
U= (\text{prefactor}) (\text{derivatives}) f_n(r) |
| 447 |
> |
U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r) |
| 448 |
|
\label{generic} |
| 449 |
|
\end{equation} |
| 450 |
|
% |
| 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 |
< |
space indices. |
| 458 |
> |
space indices. Appendix \ref{radialTSF} contains details on the |
| 459 |
> |
radial functions. |
| 460 |
|
|
| 461 |
|
In the formulas presented in the tables below, the placeholder |
| 462 |
|
function $f(r)$ is used to represent the electrostatic kernel (either |
| 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 |
|
presented in tables \ref{tab:tableenergy} and \ref{tab:tableFORCE}. |
| 469 |
+ |
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 |
|
|
| 475 |
|
\subsection{Gradient-shifted force (GSF) electrostatics} |
| 476 |
|
The second, and conceptually simpler approach to force-shifting |
| 478 |
|
expansion, and has a similar interaction energy for all multipole |
| 479 |
|
orders: |
| 480 |
|
\begin{equation} |
| 481 |
< |
U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big |
| 482 |
< |
\lvert _{r_c} . |
| 481 |
> |
U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) - |
| 482 |
> |
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 |
|
\label{generic2} |
| 485 |
|
\end{equation} |
| 486 |
< |
Here the gradient for force shifting is evaluated for an image |
| 487 |
< |
multipole projected onto the surface of the cutoff sphere (see fig |
| 488 |
< |
\ref{fig:shiftedMultipoles}). No higher order terms $(r-r_c)^n$ |
| 489 |
< |
appear. The primary difference between the TSF and GSF methods is the |
| 490 |
< |
stage at which the Taylor Series is applied; in the Taylor-shifted |
| 491 |
< |
approach, it is applied to the kernel itself. In the Gradient-shifted |
| 492 |
< |
approach, it is applied to individual radial interactions terms in the |
| 493 |
< |
multipole expansion. Energies from this method thus have the general |
| 494 |
< |
form: |
| 486 |
> |
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 |
> |
orientation ($\hat{\mathbf{b}}$) of the interacting multipole. No |
| 493 |
> |
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 |
> |
radial interaction terms in the multipole expansion. Energies from |
| 498 |
> |
this method thus have the general form: |
| 499 |
|
\begin{equation} |
| 500 |
|
U= \sum (\text{angular factor}) (\text{radial factor}). |
| 501 |
|
\label{generic3} |
| 502 |
|
\end{equation} |
| 503 |
|
|
| 504 |
|
Functional forms for both methods (TSF and GSF) can both be summarized |
| 505 |
< |
using the form of Eq.~(\ref{generic3}). The basic forms for the |
| 505 |
> |
using the form of Eq.~\ref{generic3}). The basic forms for the |
| 506 |
|
energy, force, and torque expressions are tabulated for both shifting |
| 507 |
|
approaches below -- for each separate orientational contribution, only |
| 508 |
|
the radial factors differ between the two methods. |
| 509 |
|
|
| 510 |
< |
\subsection{\label{sec:level2}Body and space axes} |
| 510 |
> |
\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 |
> |
U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right] |
| 519 |
> |
\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 |
> |
Eq.~\ref{generic3}. The energy, force, and torque expressions are |
| 530 |
> |
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 |
|
|
| 475 |
– |
[XXX Do we need this section in the main paper? or should it go in the |
| 476 |
– |
extra materials?] |
| 534 |
|
|
| 535 |
< |
So far, all energies and forces have been written in terms of fixed |
| 536 |
< |
space coordinates. Interaction energies are computed from the generic |
| 537 |
< |
formulas Eq.~(\ref{generic}) and ~(\ref{generic2}) which combine |
| 538 |
< |
orientational prefactors with radial functions. Because objects $\bf |
| 539 |
< |
a$ and $\bf b$ both translate and rotate during a molecular dynamics |
| 540 |
< |
(MD) simulation, it is desirable to contract all $r$-dependent terms |
| 541 |
< |
with dipole and quadrupole moments expressed in terms of their body |
| 542 |
< |
axes. To do so, we have followed the methodology of Allen and |
| 543 |
< |
Germano,\cite{Allen:2006fk} which was itself based on earlier work by |
| 544 |
< |
Price {\em et al.}\cite{Price:1984fk} |
| 535 |
> |
\subsection{\label{sec:level2}Body and space axes} |
| 536 |
> |
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 |
> |
In a typical simulation, the initial axes are obtained by |
| 544 |
> |
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 |
|
|
| 549 |
< |
We denote body axes for objects $\bf a$ and $\bf b$ by unit vectors |
| 490 |
< |
$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$ |
| 491 |
< |
referring to a convenient set of inertial body axes. (N.B., these |
| 492 |
< |
body axes are generally not the same as those for which the quadrupole |
| 493 |
< |
moment is diagonal.) Then, |
| 494 |
< |
% |
| 495 |
< |
\begin{eqnarray} |
| 496 |
< |
\hat{a}_m= a_{mx}\hat{x} + a_{my}\hat{y} + a_{mz}\hat{z} \\ |
| 497 |
< |
\hat{b}_m= b_{mx}\hat{x} + b_{my}\hat{y} + b_{mz}\hat{z} . |
| 498 |
< |
\end{eqnarray} |
| 499 |
< |
Rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be |
| 549 |
> |
The rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be |
| 550 |
|
expressed using these unit vectors: |
| 551 |
|
\begin{eqnarray} |
| 552 |
|
\hat{\mathbf {a}} = |
| 554 |
|
\hat{a}_1 \\ |
| 555 |
|
\hat{a}_2 \\ |
| 556 |
|
\hat{a}_3 |
| 557 |
< |
\end{pmatrix} |
| 508 |
< |
= |
| 509 |
< |
\begin{pmatrix} |
| 510 |
< |
a_{1x} \quad a_{1y} \quad a_{1z} \\ |
| 511 |
< |
a_{2x} \quad a_{2y} \quad a_{2z} \\ |
| 512 |
< |
a_{3x} \quad a_{3y} \quad a_{3z} |
| 513 |
< |
\end{pmatrix}\\ |
| 557 |
> |
\end{pmatrix}, \qquad |
| 558 |
|
\hat{\mathbf {b}} = |
| 559 |
|
\begin{pmatrix} |
| 560 |
|
\hat{b}_1 \\ |
| 561 |
|
\hat{b}_2 \\ |
| 562 |
|
\hat{b}_3 |
| 563 |
|
\end{pmatrix} |
| 520 |
– |
= |
| 521 |
– |
\begin{pmatrix} |
| 522 |
– |
b_{1x} \quad b_{1y} \quad b_{1z} \\ |
| 523 |
– |
b_{2x} \quad b_{2y} \quad b_{2z} \\ |
| 524 |
– |
b_{3x} \quad b_{3y} \quad b_{3z} |
| 525 |
– |
\end{pmatrix} . |
| 564 |
|
\end{eqnarray} |
| 565 |
|
% |
| 566 |
|
These matrices convert from space-fixed $(xyz)$ to body-fixed $(123)$ |
| 567 |
< |
coordinates. All contractions of prefactors with derivatives of |
| 568 |
< |
functions can be written in terms of these matrices. It proves to be |
| 569 |
< |
equally convenient to just write any contraction in terms of unit |
| 570 |
< |
vectors $\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$. In the torque |
| 571 |
< |
expressions, it is useful to have the angular-dependent terms |
| 572 |
< |
available in three different fashions, e.g. for the dipole-dipole |
| 573 |
< |
contraction: |
| 567 |
> |
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 |
> |
% |
| 574 |
> |
\begin{equation} |
| 575 |
> |
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 |
> |
Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $ |
| 616 |
> |
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 |
> |
\right). |
| 626 |
> |
\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 |
|
\mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}} |
| 636 |
|
= D_{\bf {a}\alpha} D_{\bf {b}\alpha} = |
| 637 |
< |
\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}} |
| 637 |
> |
\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}}. |
| 638 |
|
\end{equation} |
| 639 |
|
% |
| 640 |
|
The first two forms are written using space coordinates. The first |
| 643 |
|
explicit sums over body indices and the dot products now indicate |
| 644 |
|
contractions using space indices. |
| 645 |
|
|
| 646 |
+ |
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 |
|
|
| 652 |
|
\subsection{The Self-Interaction \label{sec:selfTerm}} |
| 653 |
|
|
| 654 |
|
In addition to cutoff-sphere neutralization, the Wolf |
| 655 |
|
summation~\cite{Wolf99} and the damped shifted force (DSF) |
| 656 |
< |
extension~\cite{Fennell:2006zl} also included self-interactions that |
| 656 |
> |
extension~\cite{Fennell:2006zl} also include self-interactions that |
| 657 |
|
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 |
| 665 |
|
the cutoff sphere. For a system of undamped charges, the total |
| 666 |
|
self-term is |
| 667 |
|
\begin{equation} |
| 668 |
< |
V_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2} |
| 668 |
> |
U_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2}. |
| 669 |
|
\label{eq:selfTerm} |
| 670 |
|
\end{equation} |
| 671 |
|
|
| 680 |
|
complexity to the Ewald sum). For a system containing only damped |
| 681 |
|
charges, the complete self-interaction can be written as |
| 682 |
|
\begin{equation} |
| 683 |
< |
V_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} + |
| 683 |
> |
U_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} + |
| 684 |
|
\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 |
< |
treatment of {\it both} the self-neutralization and reciprocal |
| 690 |
> |
treatment of the self-neutralization \textit{and} reciprocal |
| 691 |
|
contributions to the self-interaction for higher order multipoles. In |
| 692 |
|
this section we give formulae for these interactions up to quadrupolar |
| 693 |
|
order. |
| 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 |
< |
obtained by Smith and Aguado and Madden for the application of the |
| 702 |
< |
Ewald sum to multipoles.\cite{Smith82,Smith98,Aguado03} For a given |
| 703 |
< |
site which posesses a charge, dipole, and multipole, both types of |
| 704 |
< |
contribution are given in table \ref{tab:tableSelf}. |
| 701 |
> |
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 |
|
|
| 706 |
|
\begin{table*} |
| 707 |
|
\caption{\label{tab:tableSelf} Self-interaction contributions for |
| 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 |
< |
Quadrupole & $2 \text{Tr}(\mathbf{Q}_{\bf a}^2) + \text{Tr}(\mathbf{Q}_{\bf a})^2$ & |
| 716 |
> |
Quadrupole & $2 \mathbf{Q}_{\bf a}:\mathbf{Q}_{\bf a} + \text{Tr}(\mathbf{Q}_{\bf a})^2$ & |
| 717 |
|
$- \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( |
| 739 |
|
\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 |
< |
work for both the Taylor-shifted and Gradient-shifted approximations. |
| 743 |
< |
These expressions were derived using a set of generic radial |
| 744 |
< |
functions. Without using the shifting approximations mentioned above, |
| 745 |
< |
some of these radial functions would be identical, and the expressions |
| 746 |
< |
coalesce into the familiar forms for unmodified multipole-multipole |
| 747 |
< |
interactions. Table \ref{tab:tableenergy} maps between the generic |
| 748 |
< |
functions and the radial functions derived for both the Taylor-shifted |
| 749 |
< |
and Gradient-shifted methods. The energy equations are written in |
| 750 |
< |
terms of lab-frame representations of the dipoles, quadrupoles, and |
| 751 |
< |
the unit vector connecting the two objects, |
| 742 |
> |
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 |
|
|
| 753 |
|
% Energy in space coordinate form ---------------------------------------------------------------------------------------------- |
| 754 |
|
% |
| 843 |
|
U_{Q_{\bf a}Q_{\bf b}}(r)=& |
| 844 |
|
\Bigl[ |
| 845 |
|
\text{Tr} \mathbf{Q}_{\mathbf{a}} \text{Tr} \mathbf{Q}_{\mathbf{b}} |
| 846 |
< |
+2 \text{Tr} \left( |
| 847 |
< |
\mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} \right) \Bigr] v_{41}(r) |
| 846 |
> |
+2 |
| 847 |
> |
\mathbf{Q}_{\mathbf{a}} : \mathbf{Q}_{\mathbf{b}} \Bigr] v_{41}(r) |
| 848 |
|
\\ |
| 849 |
|
% 2 |
| 850 |
|
&+\Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}} |
| 875 |
|
|
| 876 |
|
\begin{sidewaystable} |
| 877 |
|
\caption{\label{tab:tableenergy}Radial functions used in the energy |
| 878 |
< |
and torque equations. The $f, g, h, s, t, \mathrm{and} u$ |
| 879 |
< |
functions used in this table are defined in Appendices B and C.} |
| 880 |
< |
\begin{tabular}{|c|c|l|l|} \hline |
| 881 |
< |
Generic&Bare Coulomb&Taylor-Shifted&Gradient-Shifted |
| 878 |
> |
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 |
> |
\begin{tabular}{|c|c|l|l|l|} \hline |
| 885 |
> |
Generic&Bare Coulomb&Taylor-Shifted (TSF)&Shifted Potential (SP)&Gradient-Shifted (GSF) |
| 886 |
|
\\ \hline |
| 887 |
|
% |
| 888 |
|
% |
| 891 |
|
$v_{01}(r)$ & |
| 892 |
|
$\frac{1}{r}$ & |
| 893 |
|
$f_0(r)$ & |
| 894 |
< |
$f(r)-f(r_c)-(r-r_c)g(r_c)$ |
| 894 |
> |
$f(r)-f(r_c)$ & |
| 895 |
> |
SP $-(r-r_c)g(r_c)$ |
| 896 |
|
\\ |
| 897 |
|
% |
| 898 |
|
% |
| 901 |
|
$v_{11}(r)$ & |
| 902 |
|
$-\frac{1}{r^2}$ & |
| 903 |
|
$g_1(r)$ & |
| 904 |
< |
$g(r)-g(r_c)-(r-r_c)h(r_c)$ \\ |
| 904 |
> |
$g(r)-g(r_c)$ & |
| 905 |
> |
SP $-(r-r_c)h(r_c)$ \\ |
| 906 |
|
% |
| 907 |
|
% |
| 908 |
|
% |
| 910 |
|
$v_{21}(r)$ & |
| 911 |
|
$-\frac{1}{r^3} $ & |
| 912 |
|
$\frac{g_2(r)}{r} $ & |
| 913 |
< |
$\frac{g(r)}{r}-\frac{g(r_c)}{r_c} -(r-r_c) |
| 914 |
< |
\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right)$ \\ |
| 913 |
> |
$\frac{g(r)}{r}-\frac{g(r_c)}{r_c}$ & |
| 914 |
> |
SP $-(r-r_c) \left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right)$ \\ |
| 915 |
> |
% |
| 916 |
> |
% |
| 917 |
> |
% |
| 918 |
|
$v_{22}(r)$ & |
| 919 |
|
$\frac{3}{r^3} $ & |
| 920 |
|
$\left(-\frac{g_2(r)}{r} + h_2(r) \right)$ & |
| 921 |
< |
$\left(-\frac{g(r)}{r}+h(r) \right) |
| 922 |
< |
-\left(-\frac{g(r_c)}{r_c}+h(r_c) \right)$ \\ |
| 812 |
< |
&&& $ ~~~-(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$ |
| 813 |
< |
\\ |
| 921 |
> |
$\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 |
|
% |
| 924 |
|
% |
| 925 |
|
% |
| 927 |
|
$v_{31}(r)$ & |
| 928 |
|
$\frac{3}{r^4} $ & |
| 929 |
|
$\left(-\frac{g_3(r)}{r^2} + \frac{h_3(r)}{r} \right)$ & |
| 930 |
< |
$\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r} \right) |
| 931 |
< |
-\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right) $\\ |
| 823 |
< |
&&&$ ~~~ -(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)$ |
| 824 |
< |
\\ |
| 930 |
> |
$\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 |
|
% |
| 933 |
+ |
% |
| 934 |
+ |
% |
| 935 |
|
$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 |
< |
$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right) |
| 939 |
< |
- \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c) \right)$ \\ |
| 940 |
< |
&&&$ ~~~ -(r-r_c) \left( \frac{-6g(r_c)}{r_c^3}+\frac{6h(r_c)}{r_c^2}-\frac{3s(r_c)}{r_c}+t(r_c) \right)$ |
| 941 |
< |
\\ |
| 938 |
> |
$\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 |
|
% |
| 943 |
|
% |
| 944 |
|
% |
| 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 |
< |
$\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right) |
| 950 |
< |
- \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2} \right)$ \\ |
| 842 |
< |
&&&$ ~~~ -(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)$ |
| 949 |
> |
$\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 |
|
\\ |
| 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 |
< |
$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right) |
| 957 |
< |
-\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ \\ |
| 958 |
< |
&&&$ ~~~ -(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3} |
| 959 |
< |
-\frac{4s(r_c)}{r_c^2} + \frac{t(r_c)}{r_c}\right)$ |
| 852 |
< |
\\ |
| 956 |
> |
$\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 |
|
% 3 |
| 961 |
+ |
% |
| 962 |
+ |
% |
| 963 |
|
$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 |
< |
$\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s(r)}{r} + t(r)\right)$ \\ |
| 967 |
< |
&&&$~~~ -\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)$ \\ |
| 968 |
< |
&&&$~~~ -(r-r_c)\left(\frac{45g(r_c)}{r_c^4}-\frac{45h(r_c)}{r_c^3}+\frac{21s(r_c)}{r_c^2} |
| 969 |
< |
-\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\ \hline |
| 966 |
> |
$ \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 |
> |
\hline |
| 971 |
|
\end{tabular} |
| 972 |
|
\end{sidewaystable} |
| 973 |
|
% |
| 976 |
|
% |
| 977 |
|
|
| 978 |
|
\begin{sidewaystable} |
| 979 |
< |
\caption{\label{tab:tableFORCE}Radial functions used in the force equations.} |
| 980 |
< |
\begin{tabular}{|c|c|l|l|} \hline |
| 981 |
< |
Function&Definition&Taylor-Shifted&Gradient-Shifted |
| 979 |
> |
\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 |
|
\\ \hline |
| 987 |
|
% |
| 988 |
|
% |
| 990 |
|
$w_a(r)$& |
| 991 |
|
$\frac{d v_{01}}{dr}$& |
| 992 |
|
$g_0(r)$& |
| 993 |
< |
$g(r)-g(r_c)$ \\ |
| 993 |
> |
$g(r)$& |
| 994 |
> |
SP $-g(r_c)$ \\ |
| 995 |
|
% |
| 996 |
|
% |
| 997 |
|
$w_b(r)$ & |
| 998 |
|
$\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $& |
| 999 |
|
$\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ & |
| 1000 |
< |
$h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\ |
| 1000 |
> |
$h(r) - \frac{v_{11}(r)}{r} $ & |
| 1001 |
> |
SP $- h(r_c)$ \\ |
| 1002 |
|
% |
| 1003 |
|
$w_c(r)$ & |
| 1004 |
|
$\frac{v_{11}(r)}{r}$ & |
| 1005 |
|
$\frac{g_1(r)}{r} $ & |
| 1006 |
+ |
$\frac{v_{11}(r)}{r}$& |
| 1007 |
|
$\frac{v_{11}(r)}{r}$\\ |
| 1008 |
|
% |
| 1009 |
|
% |
| 1010 |
|
$w_d(r)$& |
| 1011 |
|
$\frac{d v_{21}}{dr}$& |
| 1012 |
|
$\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ & |
| 1013 |
< |
$\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right) |
| 1014 |
< |
-\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $ \\ |
| 1013 |
> |
$\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 |
|
% |
| 1016 |
|
$w_e(r)$ & |
| 1017 |
|
$\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ & |
| 1018 |
|
$\frac{v_{22}(r)}{r}$ & |
| 1019 |
+ |
$\frac{v_{22}(r)}{r}$ & |
| 1020 |
|
$\frac{v_{22}(r)}{r}$ \\ |
| 1021 |
|
% |
| 1022 |
|
% |
| 1023 |
|
$w_f(r)$& |
| 1024 |
|
$\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 |
< |
$ \left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) $ \\ |
| 1027 |
< |
&&& $ ~~~- \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) |
| 910 |
< |
\right)-\frac{2v_{22}(r)}{r}$\\ |
| 1026 |
> |
$ \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 |
|
% |
| 1029 |
|
$w_g(r)$& |
| 1030 |
|
$\frac{v_{31}(r)}{r}$& |
| 1031 |
|
$ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$& |
| 1032 |
+ |
$\frac{v_{31}(r)}{r}$& |
| 1033 |
|
$\frac{v_{31}(r)}{r}$\\ |
| 1034 |
|
% |
| 1035 |
|
$w_h(r)$ & |
| 1036 |
|
$\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 |
< |
$ \left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) - \left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\ |
| 1039 |
< |
&&& $ ~~~ -\frac{v_{31}(r)}{r}$ \\ |
| 1038 |
> |
$ \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 |
|
% 2 |
| 1041 |
|
$w_i(r)$ & |
| 1042 |
|
$\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 |
+ |
$\frac{v_{32}(r)}{r}$& |
| 1045 |
|
$\frac{v_{32}(r)}{r}$\\ |
| 1046 |
|
% |
| 1047 |
|
$w_j(r)$ & |
| 1048 |
|
$\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 |
< |
$\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right)$ \\ |
| 1051 |
< |
&&& $~~~-\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) -\frac{3v_{32}}{r}$ \\ |
| 1050 |
> |
$\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 |
|
% |
| 1054 |
|
$w_k(r)$ & |
| 1055 |
|
$\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 |
< |
$\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2} \right) |
| 1058 |
< |
-\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2} \right)$ \\ |
| 1057 |
> |
$\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 |
|
% |
| 1061 |
|
$w_l(r)$ & |
| 1062 |
|
$\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 |
< |
$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\ |
| 1065 |
< |
&&& $~~~ -\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) |
| 1066 |
< |
-\frac{2v_{42}(r)}{r}$\\ |
| 1064 |
> |
$\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 |
|
% |
| 1068 |
|
$w_m(r)$ & |
| 1069 |
|
$\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$& |
| 1070 |
< |
$\left(\frac{105g_4(r)}{r^4} - \frac{105h_4(r)}{r^3} + \frac{45s_4(r)}{r^2} - \frac{10t_4(r)}{r} +u_4(r) \right)$ & |
| 1071 |
< |
$\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$\\ |
| 1072 |
< |
&&& $~~~- \left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3} |
| 1073 |
< |
+\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $\\ |
| 1074 |
< |
&&& $~~~-\frac{4v_{43}(r)}{r}$ \\ |
| 1070 |
> |
$\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 |
|
% |
| 1077 |
|
$w_n(r)$ & |
| 1078 |
|
$\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 |
+ |
$\frac{v_{42}(r)}{r}$& |
| 1081 |
|
$\frac{v_{42}(r)}{r}$\\ |
| 1082 |
|
% |
| 1083 |
|
$w_o(r)$ & |
| 1084 |
|
$\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 |
+ |
$\frac{v_{43}(r)}{r}$& |
| 1087 |
|
$\frac{v_{43}(r)}{r}$ \\ \hline |
| 1088 |
|
% |
| 1089 |
|
|
| 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} |
| 985 |
– |
% |
| 986 |
– |
Obtaining the force from the interaction energy expressions is the |
| 987 |
– |
same for higher-order multipole interactions -- the trick is to make |
| 988 |
– |
sure that all $r$-dependent derivatives are considered. This is |
| 989 |
– |
straighforward if the interaction energies are written explicitly in |
| 990 |
– |
terms of $\hat{r}$ and the body axes ($\hat{a}_m$, |
| 991 |
– |
$\hat{b}_n$) : |
| 992 |
– |
% |
| 993 |
– |
\begin{equation} |
| 994 |
– |
U(r,\{\hat{a}_m \cdot \hat{r} \}, |
| 995 |
– |
\{\hat{b}_n\cdot \hat{r} \}, |
| 996 |
– |
\{\hat{a}_m \cdot \hat{b}_n \}) . |
| 997 |
– |
\label{ugeneral} |
| 998 |
– |
\end{equation} |
| 999 |
– |
% |
| 1000 |
– |
Allen and Germano,\cite{Allen:2006fk} showed that if the energy is |
| 1001 |
– |
written in this form, the forces come out relatively cleanly, |
| 1002 |
– |
% |
| 1003 |
– |
\begin{equation} |
| 1004 |
– |
\mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} = \frac{\partial U}{\partial \mathbf{r}} |
| 1005 |
– |
= \frac{\partial U}{\partial r} \hat{r} |
| 1006 |
– |
+ \sum_m \left[ |
| 1007 |
– |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} |
| 1008 |
– |
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} |
| 1009 |
– |
+ \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} |
| 1010 |
– |
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} |
| 1011 |
– |
\right] \label{forceequation}. |
| 1012 |
– |
\end{equation} |
| 1013 |
– |
% |
| 1014 |
– |
Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ |
| 1015 |
– |
is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk} |
| 1016 |
– |
In simplifying the algebra, we have also used: |
| 1017 |
– |
% |
| 1018 |
– |
\begin{align} |
| 1019 |
– |
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} |
| 1020 |
– |
=& \frac{1}{r} \left( \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r} |
| 1021 |
– |
\right) \\ |
| 1022 |
– |
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} |
| 1023 |
– |
=& \frac{1}{r} \left( \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r} |
| 1024 |
– |
\right) . |
| 1025 |
– |
\end{align} |
| 1108 |
|
% |
| 1109 |
|
We list below the force equations written in terms of lab-frame |
| 1110 |
< |
coordinates. The radial functions used in the two methods are listed |
| 1110 |
> |
coordinates. The radial functions used in the three methods are listed |
| 1111 |
|
in Table \ref{tab:tableFORCE} |
| 1112 |
|
% |
| 1113 |
|
%SPACE COORDINATES FORCE EQUATIONS |
| 1216 |
|
\begin{split} |
| 1217 |
|
\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =& |
| 1218 |
|
\Bigl[ |
| 1219 |
< |
\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} \hat{r} |
| 1220 |
< |
+ 2 \text{Tr} ( \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_k(r) \hat{r} \\ |
| 1219 |
> |
\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 |
|
% 2 |
| 1222 |
|
&+ \Bigl[ |
| 1223 |
|
2\text{Tr}\mathbf{Q}_{\mathbf{b}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} ) |
| 1253 |
|
% Torques SECTION ----------------------------------------------------------------------------------------- |
| 1254 |
|
% |
| 1255 |
|
\subsection{Torques} |
| 1256 |
< |
When energies are written in the form of Eq.~({\ref{ugeneral}), then |
| 1175 |
< |
torques can be found in a relatively straightforward |
| 1176 |
< |
manner,\cite{Allen:2006fk} |
| 1256 |
> |
|
| 1257 |
|
% |
| 1258 |
< |
\begin{eqnarray} |
| 1259 |
< |
\mathbf{\tau}_{\bf a} = |
| 1180 |
< |
\sum_m |
| 1181 |
< |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} |
| 1182 |
< |
( \hat{r} \times \hat{a}_m ) |
| 1183 |
< |
-\sum_{mn} |
| 1184 |
< |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} |
| 1185 |
< |
(\hat{a}_m \times \hat{b}_n) \\ |
| 1186 |
< |
% |
| 1187 |
< |
\mathbf{\tau}_{\bf b} = |
| 1188 |
< |
\sum_m |
| 1189 |
< |
\frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} |
| 1190 |
< |
( \hat{r} \times \hat{b}_m) |
| 1191 |
< |
+\sum_{mn} |
| 1192 |
< |
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} |
| 1193 |
< |
(\hat{a}_m \times \hat{b}_n) . |
| 1194 |
< |
\end{eqnarray} |
| 1258 |
> |
The torques for the three methods are given in space-frame |
| 1259 |
> |
coordinates: |
| 1260 |
|
% |
| 1261 |
|
% |
| 1197 |
– |
The torques for both the Taylor-Shifted as well as Gradient-Shifted |
| 1198 |
– |
methods are given in space-frame coordinates: |
| 1199 |
– |
% |
| 1200 |
– |
% |
| 1262 |
|
\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) \\ |
| 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 |
+ |
\label{eq:matrixCross} |
| 1408 |
|
\end{equation} |
| 1409 |
|
where $\alpha+1$ and $\alpha+2$ are regarded as cyclic |
| 1410 |
|
permuations of the matrix indices. |
| 1411 |
|
|
| 1412 |
|
All of the radial functions required for torques are identical with |
| 1413 |
|
the radial functions previously computed for the interaction energies. |
| 1414 |
< |
These are tabulated for both shifted force methods in table |
| 1414 |
> |
These are tabulated for all three methods in table |
| 1415 |
|
\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 |
| 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 |
< |
real-space cutoffs. In this paper we test Taylor-shifted and |
| 1427 |
< |
Gradient-shifted electrostatics against analytical methods for |
| 1428 |
< |
computing the energies of ordered multipolar arrays. In the following |
| 1429 |
< |
paper, we test the new methods against the multipolar Ewald sum for |
| 1430 |
< |
computing the energies, forces and torques for a wide range of typical |
| 1431 |
< |
condensed-phase (disordered) systems. |
| 1426 |
> |
real-space cutoffs. In this paper we test SP, TSF, and GSF |
| 1427 |
> |
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 |
|
|
| 1433 |
|
Because long-range electrostatic effects can be significant in |
| 1434 |
|
crystalline materials, ordered multipolar arrays present one of the |
| 1438 |
|
magnetization and obtained a number of these constants.\cite{Sauer} |
| 1439 |
|
This theory was developed more completely by Luttinger and |
| 1440 |
|
Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays |
| 1441 |
< |
and other periodic structures. We have repeated the Luttinger \& |
| 1380 |
< |
Tisza series summations to much higher order and obtained the energy |
| 1381 |
< |
constants (converged to one part in $10^9$) in table \ref{tab:LT}. |
| 1441 |
> |
and other periodic structures. |
| 1442 |
|
|
| 1443 |
< |
\begin{table*}[h] |
| 1444 |
< |
\centering{ |
| 1445 |
< |
\caption{Luttinger \& Tisza arrays and their associated |
| 1446 |
< |
energy constants. Type "A" arrays have nearest neighbor strings of |
| 1447 |
< |
antiparallel dipoles. Type "B" arrays have nearest neighbor |
| 1448 |
< |
strings of antiparallel dipoles if the dipoles are contained in a |
| 1449 |
< |
plane perpendicular to the dipole direction that passes through |
| 1390 |
< |
the dipole.} |
| 1391 |
< |
} |
| 1392 |
< |
\label{tab:LT} |
| 1393 |
< |
\begin{ruledtabular} |
| 1394 |
< |
\begin{tabular}{cccc} |
| 1395 |
< |
Array Type & Lattice & Dipole Direction & Energy constants \\ \hline |
| 1396 |
< |
A & SC & 001 & -2.676788684 \\ |
| 1397 |
< |
A & BCC & 001 & 0 \\ |
| 1398 |
< |
A & BCC & 111 & -1.770078733 \\ |
| 1399 |
< |
A & FCC & 001 & 2.166932835 \\ |
| 1400 |
< |
A & FCC & 011 & -1.083466417 \\ |
| 1401 |
< |
B & SC & 001 & -2.676788684 \\ |
| 1402 |
< |
B & BCC & 001 & -1.338394342 \\ |
| 1403 |
< |
B & BCC & 111 & -1.770078733 \\ |
| 1404 |
< |
B & FCC & 001 & -1.083466417 \\ |
| 1405 |
< |
B & FCC & 011 & -1.807573634 \\ |
| 1406 |
< |
-- & BCC & minimum & -1.985920929 \\ |
| 1407 |
< |
\end{tabular} |
| 1408 |
< |
\end{ruledtabular} |
| 1409 |
< |
\end{table*} |
| 1410 |
< |
|
| 1411 |
< |
In addition to the A and B arrays, there is an additional minimum |
| 1443 |
> |
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 |
|
energy structure for the BCC lattice that was found by Luttinger \& |
| 1451 |
< |
Tisza. The total electrostatic energy for any of the arrays is given |
| 1452 |
< |
by: |
| 1451 |
> |
Tisza. The total electrostatic energy density for any of the arrays |
| 1452 |
> |
is given by: |
| 1453 |
|
\begin{equation} |
| 1454 |
|
E = C N^2 \mu^2 |
| 1455 |
|
\end{equation} |
| 1456 |
< |
where $C$ is the energy constant given in table \ref{tab:LT}, $N$ is |
| 1457 |
< |
the number of dipoles per unit volume, and $\mu$ is the strength of |
| 1458 |
< |
the dipole. |
| 1456 |
> |
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 |
|
|
| 1461 |
< |
To test the new electrostatic methods, we have constructed very large, |
| 1462 |
< |
$N$ = 8,000~(sc), 16,000~(bcc), or 32,000~(fcc) arrays of dipoles in |
| 1463 |
< |
the orientations described in table \ref{tab:LT}. For the purposes of |
| 1464 |
< |
testing the energy expressions and the self-neutralization schemes, |
| 1465 |
< |
the primary quantity of interest is the analytic energy constant for |
| 1466 |
< |
the perfect arrays. Convergence to these constants are shown as a |
| 1467 |
< |
function of both the cutoff radius, $r_c$, and the damping parameter, |
| 1468 |
< |
$\alpha$ in Figs. \ref{fig:energyConstVsCutoff} and XXX. We have |
| 1469 |
< |
simultaneously tested a hard cutoff (where the kernel is simply |
| 1470 |
< |
truncated at the cutoff radius), as well as a shifted potential (SP) |
| 1471 |
< |
form which includes a potential-shifting and self-interaction term, |
| 1433 |
< |
but does not shift the forces and torques smoothly at the cutoff |
| 1434 |
< |
radius. |
| 1461 |
> |
\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 |
|
|
| 1473 |
|
\begin{figure} |
| 1474 |
< |
\includegraphics[width=4.5in]{energyConstVsCutoff} |
| 1475 |
< |
\caption{Convergence to the analytic energy constants as a function of |
| 1476 |
< |
cutoff radius (normalized by the lattice constant) for the different |
| 1477 |
< |
real-space methods. The two crystals shown here are the ``B'' array |
| 1478 |
< |
for bcc crystals with the dipoles along the 001 direction (upper), |
| 1479 |
< |
as well as the minimum energy bcc lattice (lower). The analytic |
| 1480 |
< |
energy constants are shown as a grey dashed line. The left panel |
| 1481 |
< |
shows results for the undamped kernel ($1/r$), while the damped |
| 1445 |
< |
error function kernel, $B_0(r)$ was used in the right panel. } |
| 1446 |
< |
\label{fig:energyConstVsCutoff} |
| 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 |
|
\end{figure} |
| 1483 |
|
|
| 1484 |
< |
The Hard cutoff exhibits oscillations around the analytic energy |
| 1485 |
< |
constants, and converges to incorrect energies when the complementary |
| 1486 |
< |
error function damping kernel is used. The shifted potential (SP) and |
| 1487 |
< |
gradient-shifted force (GSF) approximations converge to the correct |
| 1488 |
< |
energy smoothly by $r_c / 6 a$ even for the undamped case. This |
| 1489 |
< |
indicates that the correction provided by the self term is required |
| 1490 |
< |
for obtaining accurate energies. The Taylor-shifted force (TSF) |
| 1491 |
< |
approximation appears to perturb the potential too much inside the |
| 1457 |
< |
cutoff region to provide accurate measures of the energy constants. |
| 1484 |
> |
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 |
> |
and the damping parameter, $\alpha$ in Figs.\ref{fig:Dipoles_rCut} |
| 1489 |
> |
and \ref{fig:Dipoles_alpha}. We have simultaneously tested a hard |
| 1490 |
> |
cutoff (where the kernel is simply truncated at the cutoff radius) in |
| 1491 |
> |
addition to the three new methods. |
| 1492 |
|
|
| 1493 |
+ |
The hard cutoff exhibits oscillations around the analytic energy |
| 1494 |
+ |
constants, and converges to incorrect energies when the complementary |
| 1495 |
+ |
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 |
|
|
| 1504 |
|
{\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 |
< |
approach.\cite{Nagai01081960,Nagai01091963} We have compared the |
| 1464 |
< |
energy constants for the lowest energy configurations for linear |
| 1465 |
< |
quadrupoles shown in table \ref{tab:NNQ} |
| 1507 |
> |
approach.\cite{Nagai01081960,Nagai01091963} |
| 1508 |
|
|
| 1467 |
– |
\begin{table*} |
| 1468 |
– |
\centering{ |
| 1469 |
– |
\caption{Nagai and Nakamura Quadurpolar arrays}} |
| 1470 |
– |
\label{tab:NNQ} |
| 1471 |
– |
\begin{ruledtabular} |
| 1472 |
– |
\begin{tabular}{ccc} |
| 1473 |
– |
Lattice & Quadrupole Direction & Energy constants \\ \hline |
| 1474 |
– |
SC & 111 & -8.3 \\ |
| 1475 |
– |
BCC & 011 & -21.7 \\ |
| 1476 |
– |
FCC & 111 & -80.5 |
| 1477 |
– |
\end{tabular} |
| 1478 |
– |
\end{ruledtabular} |
| 1479 |
– |
\end{table*} |
| 1480 |
– |
|
| 1509 |
|
In analogy to the dipolar arrays, the total electrostatic energy for |
| 1510 |
|
the quadrupolar arrays is: |
| 1511 |
|
\begin{equation} |
| 1512 |
< |
E = C \frac{3}{4} N^2 Q^2 |
| 1512 |
> |
E = C N \frac{3\bar{Q}^2}{4a^5} |
| 1513 |
|
\end{equation} |
| 1514 |
< |
where $Q$ is the quadrupole moment. |
| 1514 |
> |
where $a$ is the lattice parameter, and $\bar{Q}$ is the effective |
| 1515 |
> |
quadrupole moment, |
| 1516 |
> |
\begin{equation} |
| 1517 |
> |
\bar{Q}^2 = 2 \left(3 Q : Q - (\text{Tr} Q)^2 \right) |
| 1518 |
> |
\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 |
> |
the effective moment, $\bar{Q}^2 = \frac{2}{3} \Theta : \Theta$.) |
| 1522 |
|
|
| 1523 |
+ |
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 |
+ |
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 |
+ |
|
| 1532 |
+ |
\begin{figure} |
| 1533 |
+ |
\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 |
+ |
\label{fig:Quadrupoles_rCut} |
| 1542 |
+ |
\end{figure} |
| 1543 |
+ |
|
| 1544 |
+ |
|
| 1545 |
+ |
\begin{figure}[!htbp] |
| 1546 |
+ |
\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 |
+ |
\end{figure} |
| 1555 |
+ |
|
| 1556 |
+ |
Again, we find that the hard cutoff exhibits oscillations around the |
| 1557 |
+ |
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 |
+ |
|
| 1566 |
+ |
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 |
+ |
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 |
+ |
reasonable values for the total electrostatic energies under both the |
| 1580 |
+ |
SP and GSF approximations. |
| 1581 |
+ |
|
| 1582 |
|
\section{Conclusion} |
| 1583 |
< |
We have presented two efficient real-space methods for computing the |
| 1584 |
< |
interactions between point multipoles. These methods have the benefit |
| 1585 |
< |
of smoothly truncating the energies, forces, and torques at the cutoff |
| 1586 |
< |
radius, making them attractive for both molecular dynamics (MD) and |
| 1587 |
< |
Monte Carlo (MC) simulations. We find that the Gradient-Shifted Force |
| 1588 |
< |
(GSF) and the Shifted-Potential (SP) methods converge rapidly to the |
| 1589 |
< |
correct lattice energies for ordered dipolar and quadrupolar arrays, |
| 1590 |
< |
while the Taylor-Shifted Force (TSF) is too severe an approximation to |
| 1591 |
< |
provide accurate convergence to lattice energies. |
| 1583 |
> |
We have presented three efficient real-space methods for computing the |
| 1584 |
> |
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 |
|
|
| 1596 |
|
In most cases, GSF can obtain nearly quantitative agreement with the |
| 1597 |
|
lattice energy constants with reasonably small cutoff radii. The only |
| 1603 |
|
crystals with net-zero moments, so this is not expected to be an issue |
| 1604 |
|
in most simulations. |
| 1605 |
|
|
| 1606 |
+ |
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 |
|
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. |
| 1618 |
|
\begin{acknowledgments} |
| 1619 |
|
JDG acknowledges helpful discussions with Christopher |
| 1620 |
|
Fennell. Support for this project was provided by the National |
| 1621 |
< |
Science Foundation under grant CHE-0848243. Computational time was |
| 1621 |
> |
Science Foundation under grant CHE-1362211. Computational time was |
| 1622 |
|
provided by the Center for Research Computing (CRC) at the |
| 1623 |
|
University of Notre Dame. |
| 1624 |
|
\end{acknowledgments} |
| 1657 |
|
\text{e}^{-{\alpha}^2r^2} |
| 1658 |
|
\right] , |
| 1659 |
|
\end{equation} |
| 1660 |
< |
is very useful for building up higher derivatives. Using these |
| 1661 |
< |
formulas, we find: |
| 1660 |
> |
is very useful for building up higher derivatives. As noted by Smith, |
| 1661 |
> |
it is possible to approximate the $B_l(r)$ functions, |
| 1662 |
|
% |
| 1558 |
– |
\begin{align} |
| 1559 |
– |
\frac{dB_0}{dr}=&-rB_1(r) \\ |
| 1560 |
– |
\frac{d^2B_0}{dr^2}=& - B_1(r) + r^2 B_2(r) \\ |
| 1561 |
– |
\frac{d^3B_0}{dr^3}=& 3 r B_2(r) - r^3 B_3(r) \\ |
| 1562 |
– |
\frac{d^4B_0}{dr^4}=& 3 B_2(r) - 6 r^2 B_3(r) + r^4 B_4(r) \\ |
| 1563 |
– |
\frac{d^5B_0}{dr^5}=& - 15 r B_3(r) + 10 r^3 B_4(r) - r^5 B_5(r) . |
| 1564 |
– |
\end{align} |
| 1565 |
– |
% |
| 1566 |
– |
As noted by Smith, it is possible to approximate the $B_l(r)$ |
| 1567 |
– |
functions, |
| 1568 |
– |
% |
| 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 |
|
\newpage |
| 1668 |
|
\section{The $r$-dependent factors for TSF electrostatics} |
| 1669 |
+ |
\label{radialTSF} |
| 1670 |
|
|
| 1671 |
|
Using the shifted damped functions $f_n(r)$ defined by: |
| 1672 |
|
% |
| 1704 |
|
\begin{equation} |
| 1705 |
|
u_4(r)=B_0^{(5)}(r) - B_0^{(5)}(r_c) . |
| 1706 |
|
\end{equation} |
| 1707 |
< |
|
| 1707 |
> |
% 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 |
|
The functions $f_n(r)$ to $u_n(r)$ can be computed recursively and |
| 1719 |
< |
stored on a grid for values of $r$ from $0$ to $r_c$. The functions |
| 1720 |
< |
needed are listed schematically below: |
| 1719 |
> |
stored on a grid for values of $r$ from $0$ to $r_c$. Using these |
| 1720 |
> |
functions, we find |
| 1721 |
|
% |
| 1617 |
– |
\begin{eqnarray} |
| 1618 |
– |
f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\ |
| 1619 |
– |
g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\ |
| 1620 |
– |
h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\ |
| 1621 |
– |
s_2 \quad s_3 \quad &s_4 \nonumber \\ |
| 1622 |
– |
t_3 \quad &t_4 \nonumber \\ |
| 1623 |
– |
&u_4 \nonumber . |
| 1624 |
– |
\end{eqnarray} |
| 1625 |
– |
|
| 1626 |
– |
Using these functions, we find |
| 1627 |
– |
% |
| 1722 |
|
\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 |
< |
\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta r_\gamma} =& |
| 1726 |
> |
\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta \partial r_\gamma} =& |
| 1727 |
|
\left( \delta_{\alpha \beta} r_\gamma + \delta_{\alpha \gamma} r_\beta + |
| 1728 |
|
\delta_{ \beta \gamma} r_\alpha \right) |
| 1729 |
< |
\left( -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right) |
| 1730 |
< |
+ r_\alpha r_\beta r_\gamma |
| 1729 |
> |
\left( -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right) \nonumber \\ |
| 1730 |
> |
& + r_\alpha r_\beta r_\gamma |
| 1731 |
|
\left( \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \\ |
| 1732 |
< |
\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta r_\gamma r_\delta} =& |
| 1732 |
> |
\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta \partial |
| 1733 |
> |
r_\gamma \partial r_\delta} =& |
| 1734 |
|
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} |
| 1735 |
|
+ \delta_{\alpha \gamma} \delta_{\beta \delta} |
| 1736 |
|
+\delta_{ \beta \gamma} \delta_{\alpha \delta} \right) |
| 1743 |
|
\left( -\frac{15g_n}{r^7} + \frac{15h_n}{r^6} - \frac{6s_n}{r^5} |
| 1744 |
|
+ \frac{t_n}{r^4} \right)\\ |
| 1745 |
|
\frac{\partial^5 f_n} |
| 1746 |
< |
{\partial r_\alpha \partial r_\beta r_\gamma r_\delta r_\epsilon} =& |
| 1746 |
> |
{\partial r_\alpha \partial r_\beta \partial r_\gamma \partial |
| 1747 |
> |
r_\delta \partial r_\epsilon} =& |
| 1748 |
|
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} r_\epsilon |
| 1749 |
|
+ \text{14 permutations} \right) |
| 1750 |
|
\left( \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \nonumber \\ |
| 1761 |
|
% |
| 1762 |
|
\newpage |
| 1763 |
|
\section{The $r$-dependent factors for GSF electrostatics} |
| 1764 |
+ |
\label{radialGSF} |
| 1765 |
|
|
| 1766 |
|
In Gradient-shifted force electrostatics, the kernel is not expanded, |
| 1767 |
< |
rather the individual terms in the multipole interaction energies. |
| 1768 |
< |
For damped charges , this still brings into the algebra multiple |
| 1769 |
< |
derivatives of the Smith's $B_0(r)$ function. To denote these terms, |
| 1770 |
< |
we generalize the notation of the previous appendix. For $f(r)=1/r$ |
| 1771 |
< |
(bare Coulomb) or $f(r)=B_0(r)$ (smeared charge) |
| 1767 |
> |
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 |
|
% |
| 1774 |
|
\begin{align} |
| 1775 |
< |
g(r)=& \frac{df}{d r}\\ |
| 1776 |
< |
h(r)=& \frac{dg}{d r} = \frac{d^2f}{d r^2} \\ |
| 1777 |
< |
s(r)=& \frac{dh}{d r} = \frac{d^3f}{d r^3} \\ |
| 1778 |
< |
t(r)=& \frac{ds}{d r} = \frac{d^4f}{d r^4} \\ |
| 1779 |
< |
u(r)=& \frac{dt}{d r} = \frac{d^5f}{d r^5} . |
| 1775 |
> |
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 |
|
\end{align} |
| 1784 |
|
% |
| 1785 |
< |
For undamped charges, $f(r)=1/r$, Table I lists these derivatives |
| 1786 |
< |
under the column ``Bare Coulomb.'' Equations \ref{eq:b9} to |
| 1787 |
< |
\ref{eq:b13} are still correct for GSF electrostatics if the subscript |
| 1687 |
< |
$n$ is eliminated. |
| 1785 |
> |
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 |
> |
correct for GSF electrostatics if the subscript $n$ is eliminated. |
| 1788 |
|
|
| 1789 |
|
\newpage |
| 1790 |
|
|
