trunk/multipole/multipole1.tex

% ****** Start of file aipsamp.tex ******
%
%   This file is part of the AIP files in the AIP distribution for REVTeX 4.
%   Version 4.1 of REVTeX, October 2009
%
%   Copyright (c) 2009 American Institute of Physics.
%
%   See the AIP README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.1
%
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex  aipsamp
%  2)  bibtex aipsamp
%  3)  latex  aipsamp
%  4)  latex  aipsamp
%
% Use this file as a source of example code for your aip document.
% Use the file aiptemplate.tex as a template for your document.
\documentclass[%
 aip,jcp,
 amsmath,amssymb,
preprint,%
% reprint,%
%author-year,%
%author-numerical,%
jcp]{revtex4-1}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
%\usepackage{bm}% bold math
\usepackage{times}
\usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
\usepackage{url}

%\usepackage[mathlines]{lineno}% Enable numbering of text and display math
%\linenumbers\relax % Commence numbering lines

\begin{document}

\preprint{AIP/123-QED}

\title[Taylor-shifted and Gradient-shifted electrostatics for multipoles]
{Real space alternatives to the Ewald
Sum. I. Taylor-shifted and Gradient-shifted electrostatics for multipoles}

\author{Madan Lamichhane}
\affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}

\author{J. Daniel Gezelter}
 \email{gezelter@nd.edu.}
\affiliation{Department of Chemistry and Biochemistry, University
of Notre Dame, Notre Dame, IN 46556}

\author{Kathie E. Newman}
\affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}


\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
  We have extended the original damped-shifted force (DSF)
  electrostatic kernel and have been able to derive two new
  electrostatic potentials for higher-order multipoles that are based
  on truncated Taylor expansions around the cutoff radius.  For
  multipole-multipole interactions, we find that each of the distinct
  orientational contributions has a separate radial function to ensure
  that the overall forces and torques vanish at the cutoff radius. In
  this paper, we present energy, force, and torque expressions for the
  new models, and compare these real-space interaction models to exact
  results for ordered arrays of multipoles.
\end{abstract}

\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
                             % Classification Scheme.
\keywords{Suggested keywords}%Use showkeys class option if keyword
                              %display desired
\maketitle

\section{Introduction}
There has been increasing interest in real-space methods for
calculating electrostatic interactions in computer simulations of
condensed molecular
systems.\cite{Wolf99,Zahn02,Kast03,BeckD.A.C._bi0486381,Ma05,Fennell:2006zl,Chen:2004du,Chen:2006ii,Rodgers:2006nw,Denesyuk:2008ez,Izvekov:2008wo}
The simplest of these techniques was developed by Wolf {\it et al.}
in their work towards an $\mathcal{O}(N)$ Coulombic sum.\cite{Wolf99}
For systems of point charges, Fennell and Gezelter showed that a
simple damped shifted force (DSF) modification to Wolf's method could
give nearly quantitative agreement with smooth particle mesh Ewald
(SPME)\cite{Essmann95} configurational energy differences as well as
atomic force and molecular torque vectors.\cite{Fennell:2006zl}

The computational efficiency and the accuracy of the DSF method are
surprisingly good, particularly for systems with uniform charge
density. Additionally, dielectric constants obtained using DSF and
similar methods where the force vanishes at $R_\textrm{c}$ are
essentially quantitative.\cite{Izvekov:2008wo} The DSF and other
related methods have now been widely investigated,\cite{Hansen:2012uq}
and DSF is now used routinely in simulations of ionic
liquids,\cite{doi:10.1021/la400226g,McCann:2013fk} flow in carbon
nanotubes,\cite{kannam:094701} gas sorption in metal-organic framework
materials,\cite{Forrest:2012ly} thermal conductivity of methane
hydrates,\cite{English:2008kx} condensation coefficients of
water,\cite{Louden:2013ve} and in tribology at solid-liquid-solid
interfaces.\cite{Tokumasu:2013zr} DSF electrostatics provides a
compromise between the computational speed of real-space cutoffs and
the accuracy of fully-periodic Ewald treatments.

\subsection{Coarse Graining using Point Multipoles}
One common feature of many coarse-graining approaches, which treat
entire molecular subsystems as a single rigid body, is simplification
of the electrostatic interactions between these bodies so that fewer
site-site interactions are required to compute configurational
energies. Notably, the force matching approaches of Voth and coworkers
are an exciting development in their ability to represent realistic
(and {\it reactive}) chemical systems for very large length scales and
long times.  This approach utilized a coarse-graining in interaction
space (CGIS) which fits an effective force for the coarse grained
system using a variational force-matching method to a fine-grained
simulation.\cite{Izvekov:2008wo}

The coarse-graining approaches of Ren \& coworkers,\cite{Golubkov06}
and Essex \&
coworkers,\cite{ISI:000276097500009,ISI:000298664400012}
both utilize Gay-Berne
ellipsoids~\cite{Berne72,Gay81,Luckhurst90,Cleaver96,Berardi98,Ravichandran:1999fk,Berardi99,Pasterny00}
to model dispersive interactions and point multipoles to model
electrostatics for entire molecules or functional groups.

Ichiye and coworkers have recently introduced a number of very fast
water models based on a ``sticky'' multipole model which are
qualitatively better at reproducing the behavior of real water than
the more common point-charge models (SPC/E, TIPnP).  The point charge
models are also substantially more computationally demanding than the
sticky multipole
approach.\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} The
SSDQO model requires the use of an approximate multipole expansion
(AME) as the exact multipole expansion is quite expensive
(particularly when handled via the Ewald sum).\cite{Ichiye:2006qy}
Another particularly important use of point multipoles (and multipole
polarizability) is in the very high-quality AMOEBA water model and
related force fields.\cite{Ponder:2010fk,schnieders:124114,Ren:2011uq}

Higher-order multipoles present a peculiar issue for molecular
dynamics. Multipolar interactions are inherently short-ranged, and
should not need the relatively expensive Ewald treatment.  However,
real-space cutoff methods are normally applied in an orientation-blind
fashion so multipoles which leave and then re-enter a cutoff sphere in
a different orientation can cause energy discontinuities.

This paper outlines an extension of the original DSF electrostatic
kernel to point multipoles.  We have developed two distinct real-space
interaction models for higher-order multipoles based on two truncated
Taylor expansions that are carried out at the cutoff radius.  We are
calling these models {\bf Taylor-shifted} and {\bf Gradient-shifted}
electrostatics.  Because of differences in the initial assumptions,
the two methods yield related, but different expressions for energies,
forces, and torques.  

In this paper we outline the new methodology and give functional forms
for the energies, forces, and torques up to quadrupole-quadrupole
order.  We also compare the new methods to analytic energy constants
for periodic arrays of point multipoles.  In the following paper, we
provide extensive numerical comparisons to Ewald-based electrostatics
in common simulation enviornments.  

\section{Methodology}

\subsection{Self-neutralization, damping, and force-shifting}
The DSF and Wolf methods operate by neutralizing the total charge
contained within the cutoff sphere surrounding each particle.  This is
accomplished by shifting the potential functions to generate image
charges on the surface of the cutoff sphere for each pair interaction
computed within $R_\textrm{c}$. Damping using a complementary error
function is applied to the potential to accelerate convergence. The
potential for the DSF method, where $\alpha$ is the adjustable damping
parameter, is given by
\begin{equation*}
V_\mathrm{DSF}(r) = C_a C_b \Biggr{[}
\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} 
- \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} + \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2} 
+ \frac{2\alpha}{\pi^{1/2}}
\frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}
\right)\left(r_{ij}-R_\mathrm{c}\right)\ \Biggr{]} 
\label{eq:DSFPot}
\end{equation*}

To insure net charge neutrality within each cutoff sphere, an
additional ``self'' term is added to the potential. This term is
constant (as long as the charges and cutoff radius do not change), and
exists outside the normal pair-loop for molecular simulations.  It can
be thought of as a contribution from a charge opposite in sign, but
equal in magnitude, to the central charge, which has been spread out
over the surface of the cutoff sphere.  A portion of the self term is
identical to the self term in the Ewald summation, and comes from the
utilization of the complimentary error function for electrostatic
damping.\cite{deLeeuw80,Wolf99}

There have been recent efforts to extend the Wolf self-neutralization
method to zero out the dipole and higher order multipoles contained
within the cutoff
sphere.\cite{Fukuda:2011jk,Fukuda:2012yu,Fukuda:2013qv}  

In this work, we will be extending the idea of self-neutralization for
the point multipoles in a similar way.  In Figure
\ref{fig:shiftedMultipoles}, the central dipolar site $\mathbf{D}_i$
is interacting with point dipole $\mathbf{D}_j$ and point quadrupole,
$\mathbf{Q}_k$.  The self-neutralization scheme for point multipoles
involves projecting opposing multipoles for sites $j$ and $k$ on the
surface of the cutoff sphere.

\begin{figure}
\includegraphics[width=3in]{SM}
\caption{Reversed multipoles are projected onto the surface of the
  cutoff sphere. The forces, torques, and potential are then smoothly
  shifted to zero as the sites leave the cutoff region.}
\label{fig:shiftedMultipoles}
\end{figure}

As in the point-charge approach, there is a contribution from
self-neutralization of site $i$.  The self term for multipoles is
described in section \ref{sec:selfTerm}.

\subsection{The multipole expansion}

Consider two discrete rigid collections of point charges, denoted as
$\bf a$ and $\bf b$.  In the following, we assume that the two objects
interact via electrostatics only and describe those interactions in
terms of a standard multipole expansion.  Putting the origin of the
coordinate system at the center of mass of $\bf a$, we use vectors
$\mathbf{r}_k$ to denote the positions of all charges $q_k$ in $\bf
a$.  Then the electrostatic potential of object $\bf a$ at
$\mathbf{r}$ is given by
\begin{equation}
V_a(\mathbf r) = \frac{1}{4\pi\epsilon_0}
 \sum_{k \, \text{in \bf a}} \frac{q_k}{\lvert \mathbf{r} - \mathbf{r}_k \rvert}.
\end{equation}
The Taylor expansion in $r$ can be written using an implied summation
notation.  Here Greek indices are used to indicate space coordinates
($x$, $y$, $z$) and the subscripts $k$ and $j$ are reserved for
labelling specific charges in $\bf a$ and $\bf b$ respectively.  The
Taylor expansion,
\begin{equation}
 \frac{1}{\lvert \mathbf{r} - \mathbf{r}_k \rvert} = 
\left( 1
-  r_{k\alpha} \frac{\partial}{\partial r_{\alpha}} 
+ \frac{1}{2}  r_{k\alpha} r_{k\beta} \frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} +\dots 
\right)
\frac{1}{r}  ,
\end{equation}
can then be used to express the electrostatic potential on $\bf a$ in
terms of multipole operators,
\begin{equation}
V_{\bf a}(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\hat{M}_{\bf a} \frac{1}{r} 
\end{equation}
where
\begin{equation}
\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}} 
+  Q_{{\bf a}\alpha\beta}
 \frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots
\end{equation}
Here, the point charge, dipole, and quadrupole for object $\bf a$ are
given by $C_{\bf a}$, $D_{{\bf a}\alpha}$, and $Q_{{\bf
    a}\alpha\beta}$, respectively.  These are the primitive multipoles
which can be expressed as a distribution of charges,
\begin{align}
C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \\
D_{{\bf a}\alpha} =&\sum_{k \, \text{in \bf a}} q_k r_{k\alpha} ,\\
Q_{{\bf a}\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k r_{k\alpha}  r_{k\beta} .
\end{align}
Note that the definition of the primitive quadrupole here differs from
the standard traceless form, and contains an additional Taylor-series
based factor of $1/2$. 

It is convenient to locate charges $q_j$ relative to the center of mass of  $\bf b$.  Then with $\bf{r}$ pointing from
$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $), the interaction energy is given by
\begin{equation}
U_{\bf{ab}}(r)
= \frac{1}{4\pi \epsilon_0} \hat{M}_a \sum_{j \, \text{in \bf b}} \frac {q_j}{\vert \bf{r}+\bf{r}_j \vert} .
\end{equation}
This can also be expanded as a Taylor series in $r$.  Using a notation
similar to before to define the multipoles on object {\bf b},
\begin{equation}
\hat{M}_{\bf b} = C_{\bf b} + D_{{\bf b}\alpha} \frac{\partial}{\partial r_{\alpha}} 
+   Q_{{\bf b}\alpha\beta}
 \frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots
\end{equation}
we arrive at the multipole expression for the total interaction energy.
\begin{equation}
U_{\bf{ab}}(r)=\frac{\hat{M}_{\bf a} \hat{M}_{\bf b}}{4\pi \epsilon_0} \frac{1}{r}  \label{kernel}.
\end{equation}
This form has the benefit of separating out the energies of
interaction into contributions from the charge, dipole, and quadrupole
of $\bf a$ interacting with the same multipoles $\bf b$.

\subsection{Damped Coulomb interactions}
In the standard multipole expansion, one typically uses the bare
Coulomb potential, with radial dependence $1/r$, as shown in
Eq.~(\ref{kernel}).  It is also quite common to use a damped Coulomb
interaction, which results from replacing point charges with Gaussian
distributions of charge with width $\alpha$.  In damped multipole
electrostatics, the kernel ($1/r$) of the expansion is replaced with
the function:
\begin{equation}
B_0(r)=\frac{\text{erfc}(\alpha r)}{r} = \frac{2}{\sqrt{\pi}r}
\int_{\alpha r}^{\infty} \text{e}^{-s^2} ds .
\end{equation}
We develop equations below using the function $f(r)$ to represent
either $1/r$ or $B_0(r)$, and all of the techniques can be applied
either to bare or damped Coulomb kernels as long as derivatives of
these functions are known.  Smith's convenient functions $B_l(r)$ are
summarized in Appendix A.


The main goal of this work is to smoothly cut off the interaction
energy as well as forces and torques as $r\rightarrow r_c$.  To
describe how this goal may be met, we use two examples, charge-charge
and charge-dipole, using the bare Coulomb kernel $f(r)=1/r$ to explain
the idea.

\subsection{Shifted-force methods}
In the shifted-force approximation, the interaction energy for two
charges $C_{\bf a}$ and $C_{\bf b}$ separated by a distance $r$ is
written:
\begin{equation}
U_{C_{\bf a}C_{\bf b}}(r)=\frac{1}{4\pi \epsilon_0} C_{\bf a} C_{\bf b}
\left({ \frac{1}{r} - \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2}  }
\right) .
\end{equation}
Two shifting terms appear in this equations, one from the
neutralization procedure ($-1/r_c$), and one that causes the first
derivative to vanish at the cutoff radius.

Since one derivative of the interaction energy is needed for the
force, the minimal perturbation is a term linear in $(r-r_c)$ in the
interaction energy, that is:
\begin{equation}
\frac{d\,}{dr} 
\left( {\frac{1}{r} - \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2}  }
\right) = \left(- \frac{1}{r^2} + \frac{1}{r_c^2} 
\right) .
\end{equation}
There are a number of ways to generalize this derivative shift for
higher-order multipoles.  Below, we present two methods, one based on
higher-order Taylor series for $r$ near $r_c$, and the other based on
linear shift of the kernel gradients at the cutoff itself.

\subsection{Taylor-shifted force (TSF) electrostatics}
In the Taylor-shifted force (TSF) method, the procedure that we follow
is based on a Taylor expansion containing the same number of
derivatives required for each force term to vanish at the cutoff.  For
example, the quadrupole-quadrupole interaction energy requires four
derivatives of the kernel, and the force requires one additional
derivative. We therefore require shifted energy expressions that
include enough terms so that all energies, forces, and torques are
zero as $r \rightarrow r_c$.  In each case, we will subtract off a
function $f_n^{\text{shift}}(r)$ from the kernel $f(r)=1/r$.  The
subscript $n$ indicates the number of derivatives to be taken when
deriving a given multipole energy.  We choose a function with
guaranteed smooth derivatives --- a truncated Taylor series of the
function $f(r)$, e.g.,
%
\begin{equation}
f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)}(r_c)  .
\end{equation}
%
The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$.
Thus, for $f(r)=1/r$, we find
%
\begin{equation}
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} .
\end{equation}
%
Continuing with the example of a charge $\bf a$ interacting with a
dipole $\bf b$, we write
%
\begin{equation}
U_{C_{\bf a}D_{\bf b}}(r)=
\frac{C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}  \frac {\partial f_1(r) }{\partial r_\alpha} 
=\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}   
\frac {r_\alpha}{r} \frac {\partial f_1(r)}{\partial r} .
\end{equation}
%
The force that dipole  $\bf b$ exerts on charge $\bf a$ is
%
\begin{equation}
F_{C_{\bf a}D_{\bf b}\beta} =\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}
\left[ \frac{\delta_{\alpha\beta}}{r} \frac {\partial}{\partial r} + 
\frac{r_\alpha r_\beta}{r^2}
\left( -\frac{1}{r} \frac {\partial} {\partial r} 
+ \frac {\partial ^2} {\partial r^2} \right) \right] f_1(r) .
\end{equation}
%
For undamped coulombic interactions, $f(r)=1/r$, we find
%
\begin{equation}
F_{C_{\bf a}D_{\bf b}\beta} =
\frac{C_{\bf a} D_{{\bf b}\beta} }{4\pi \epsilon_0r}
\left[  -\frac{1}{r^2}+\frac{1}{r_c^2}-\frac{2(r-r_c)}{r_c^3} \right] 
+\frac{C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta }{4\pi \epsilon_0}
\left[ \frac{3}{r^5}-\frac{3}{r^3r_c^2} \right] .
\end{equation}
%
This expansion shows the expected $1/r^3$ dependence of the force.  

In general, we can write
%
\begin{equation}
U=\frac{1}{4\pi \epsilon_0} (\text{prefactor}) (\text{derivatives}) f_n(r)
\label{generic}
\end{equation}
%
where $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole
and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for quadrupole-quadrupole.
An example is the case of quadrupole-quadrupole for which the $\text{prefactor}$ is 
$Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$ and the derivatives are
$\partial^4/\partial r_\alpha \partial r_\beta \partial r_\gamma \partial r_\delta$, with
implied summation combining the space indices.

In the formulas presented in the tables below, the placeholder
function $f(r)$ is used to represent the electrostatic kernel (either
damped or undamped).  The main functions that go into the force and
torque terms, $f_n(r), g_n(r), h_n(r), s_n(r), \mathrm{~and~} t_n(r)$
are successive derivatives of the shifted electrostatic kernel of the
same index $n$.  The algebra required to evaluate energies, forces and
torques is somewhat tedious and are summarized in Appendices A and B.

\subsection{Gradient-shifted force (GSF) electrostatics}
Note the method used in the previous subsection to smoothly shift the
force to zero is a truncated Taylor Series in the radius $r$.  The
second method maintains only the linear $(r-r_c)$ term and has a
similar interaction energy for all multipole orders:
\begin{equation}
U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big \lvert  _{r_c} .
\end{equation}
No higher order terms $(r-r_c)^n$ appear.  The primary difference
between the TSF and GSF methods is the stage at which the Taylor
Series is applied; in the Taylor-shifted approach, it is applied to
the kernel itself.  In the Gradient-shifted approach, it is applied to
individual radial interactions terms in the multipole expansion.
Terms from this method thus have the general form:
\begin{equation}
U=\frac{1}{4\pi \epsilon_0}\sum  (\text{angular factor}) (\text{radial factor}).
\label{generic2}
\end{equation}

Results for both methods can be summarized using the form of
Eq.~(\ref{generic2}) and are listed in Table I below.

\subsection{\label{sec:level2}Body and space axes}

So far, all energies and forces have been written in terms of fixed
space coordinates $x$, $y$, $z$.  Interaction energies are computed
from the generic formulas Eq.~(\ref{generic}) and ~(\ref{generic2})
which combine prefactors with radial functions.  Because objects $\bf
a$ and $\bf b$ both translate and rotate during a molecular dynamics
(MD) simulation, it is desirable to contract all $r$-dependent terms
with dipole and quadrupole moments expressed in terms of their body
axes.  To do so, we follow the methodology of Allen and
Germano,\cite{Allen:2006fk} which was itself based on an earlier paper
by Price {\em et al.}\cite{Price:1984fk}

We denote body axes for objects $\bf a$ and $\bf b$ by unit vectors
$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$
referring to a convenient set of inertial body axes.  (N.B., these
body axes are generally not the same as those for which the quadrupole
moment is diagonal.)  Then,
%
\begin{eqnarray}
\hat{a}_m= a_{mx}\hat{x} + a_{my}\hat{y} + a_{mz}\hat{z}  \\
\hat{b}_m= b_{mx}\hat{x} + b_{my}\hat{y} + b_{mz}\hat{z}  .
\end{eqnarray}
Allen and Germano define matrices $\hat{\mathbf {a}}$ 
and  $\hat{\mathbf {b}}$ using these unit vectors:
\begin{eqnarray}
\hat{\mathbf {a}} = 
\begin{pmatrix}
\hat{a}_1 \\
\hat{a}_2 \\
\hat{a}_3
\end{pmatrix}
=
\begin{pmatrix}
a_{1x} \quad a_{1y} \quad a_{1z} \\
a_{2x} \quad a_{2y} \quad a_{2z} \\
a_{3x} \quad a_{3y} \quad a_{3z}
\end{pmatrix}\\
\hat{\mathbf {b}} = 
\begin{pmatrix}
\hat{b}_1 \\
\hat{b}_2 \\
\hat{b}_3
\end{pmatrix}
=
\begin{pmatrix}
b_{1x}\quad  b_{1y} \quad b_{1z} \\
b_{2x} \quad b_{2y} \quad b_{2z} \\
b_{3x} \quad b_{3y} \quad b_{3z}
\end{pmatrix}  .
\end{eqnarray}
%
These matrices convert from space-fixed $(xyz)$ to object-fixed $(123)$ coordinates.
All contractions of prefactors with derivatives of functions can be written in terms of these matrices.
It proves to be equally convenient to just write any contraction in terms of unit vectors
$\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$.  
We have found it useful to write angular-dependent terms in three different fashions,
illustrated by the following 
three examples from the interaction-energy expressions:
%
\begin{eqnarray}
 \mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}}
=D_{\bf {a}\alpha} D_{\bf {b}\alpha}=
\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}} \\
r^2 \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right)=
r_\alpha Q_{\bf b \alpha \beta} r_\beta = r^2
\sum_{mn}(\hat{r} \cdot \hat{b}_m)  Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \\
r ( \mathbf{D}_{\mathbf{a}} \cdot 
\mathbf{Q}_{\mathbf{b}}  \cdot \hat{r})=
D_{\bf {a}\alpha}  Q_{\bf b \alpha \beta} r_\beta
=r \sum_{lmn} D_{\mathbf{a}l} (\hat{a}_l \cdot \hat{b}_m ) Q_{\mathbf{b}mn}
(\hat{b}_n \cdot \hat{r}) .
\end{eqnarray}
%
[Dan, perhaps there are better examples to show here.]

In each line, the first two terms are written using space coordinates.  The first form is standard
in the chemistry literature, and the second is ``physicist style'' using implied summation notation.  The third 
form shows explicitly sums over body indices and the dot products now indicate contractions using space indices.
We find the first form to be useful in writing equations prior to converting to computer code.  The second form is helpful
in derivations of the interaction energy expressions.  The third one is specifically helpful when deriving forces and torques, as will
be discussed below.  


\subsection{The Self-Interaction \label{sec:selfTerm}}

The Wolf summation~\cite{Wolf99} and the later damped shifted force
(DSF) extension~\cite{Fennell06} included self-interactions that are
handled separately from the pairwise interactions between sites. The
self-term is normally calculated via a single loop over all sites in
the system, and is relatively cheap to evaluate. The self-interaction
has contributions from two sources:
\begin{itemize}
\item The neutralization procedure within the cutoff radius requires a
  contribution from a charge opposite in sign, but equal in magnitude,
  to the central charge, which has been spread out over the surface of
  the cutoff sphere.  For a system of undamped charges, the total
  self-term is
\begin{equation}
V_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2}
\label{eq:selfTerm}
\end{equation}
Note that in this potential and in all electrostatic quantities that
follow, the standard $4 \pi \epsilon_{0}$ has been omitted for
clarity.
\item Charge damping with the complementary error function is a
  partial analogy to the Ewald procedure which splits the interaction
  into real and reciprocal space sums.  The real space sum is retained
  in the Wolf and DSF methods.  The reciprocal space sum is first
  minimized by folding the largest contribution (the self-interaction)
  into the self-interaction from charge neutralization of the damped
  potential.  The remainder of the reciprocal space portion is then
  discarded (as this contributes the largest computational cost and
  complexity to the Ewald sum).  For a system containing only damped
  charges, the complete self-interaction can be written as
\begin{equation}
V_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} + 
  \frac{\alpha}{\sqrt{\pi}} \right) \sum_{{\bf a}=1}^N
      C_{\bf a}^{2}.
\label{eq:dampSelfTerm}
\end{equation}
\end{itemize}

The extension of DSF electrostatics to point multipoles requires
treatment of {\it both} the self-neutralization and reciprocal
contributions to the self-interaction for higher order multipoles.  In
this section we give formulae for these interactions up to quadrupolar
order.

The self-neutralization term is computed by taking the {\it
  non-shifted} kernel for each interaction, placing a multipole of
equal magnitude (but opposite in polarization) on the surface of the
cutoff sphere, and averaging over the surface of the cutoff sphere.
Because the self term is carried out as a single sum over sites, the
reciprocal-space portion is identical to half of the self-term
obtained by Smith and Aguado and Madden for the application of the
Ewald sum to multipoles.\cite{Smith82,Smith98,Aguado03} For a given
site which posesses a charge, dipole, and multipole, both types of
contribution are given in table \ref{tab:tableSelf}.

\begin{table*}
  \caption{\label{tab:tableSelf} Self-interaction contributions for 
    site ({\bf a}) that has a charge $(C_{\bf a})$, dipole
    $(\mathbf{D}_{\bf a})$, and quadrupole $(\mathbf{Q}_{\bf a})$}
\begin{ruledtabular}
\begin{tabular}{lccc}
Multipole order & Summed Quantity & Self-neutralization  & Reciprocal \\ \hline
Charge & $C_{\bf a}^2$ & $-f(r_c)$ & $-\frac{\alpha}{\sqrt{\pi}}$ \\
Dipole & $|\mathbf{D}_{\bf a}|^2$ & $\frac{1}{3} \left( h(r_c) +
  \frac{2 g(r_c)}{r_c} \right)$ & $-\frac{2 \alpha^3}{3 \sqrt{\pi}}$\\
Quadrupole & $2 \text{Tr}(\mathbf{Q}_{\bf a}^2) + \text{Tr}(\mathbf{Q}_{\bf a})^2$ &
$- \frac{1}{15} \left( t(r_c)+ \frac{4 s(r_c)}{r_c} \right)$ &
$-\frac{4 \alpha^5}{5 \sqrt{\pi}}$ \\
Charge-Quadrupole & $-2 C_{\bf a} \text{Tr}(\mathbf{Q}_{\bf a})$ & $\frac{1}{3} \left(
  h(r_c) + \frac{2 g(r_c)}{r_c} \right)$& $-\frac{2 \alpha^3}{3 \sqrt{\pi}}$ \\
\end{tabular}
\end{ruledtabular}
\end{table*}

For sites which simultaneously contain charges and quadrupoles, the
self-interaction includes a cross-interaction between these two
multipole orders.  Symmetry prevents the charge-dipole and
dipole-quadrupole interactions from contributing to the
self-interaction.  The functions that go into the self-neutralization
terms, $f(r), g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive
derivatives of the electrostatic kernel (either the undamped $1/r$ or
the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that are
evaluated at the cutoff distance.  For undamped interactions, $f(r_c)
= 1/r_c$, $g(r_c) = -1/r_c^{2}$, and so on.  For damped interactions,
$f(r_c) = B_0(r_c)$, $g(r_c) = B_0'(r_c)$, and so on.  Appendix XX
contains recursion relations that allow rapid evaluation of these
derivatives.

\section{Energies, forces, and torques}
\subsection{Interaction energies}

We now list multipole interaction energies using a set of generic
radial functions.  Table \ref{tab:tableenergy} maps between the
generic functions and the radial functions derived for both the
Taylor-shifted and Gradient-shifted methods.  This set of equations is
written in terms of space coordinates:

% Energy in space coordinate form ----------------------------------------------------------------------------------------------
%
%
% u ca cb
%
\begin{align}
U_{C_{\bf a}C_{\bf b}}(r)=&
\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r)  \label{uchch}
\\
%
% u ca db
%
U_{C_{\bf a}D_{\bf b}}(r)=&
\frac{C_{\bf a}}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)  v_{11}(r)  
 \label{uchdip}
\\
%
% u ca qb
%
U_{C_{\bf a}Q_{\bf b}}(r)=&
\frac{C_{\bf a }}{4\pi \epsilon_0} \Bigl[ \text{Tr}Q_{\bf b} v_{21}(r) 
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{22}(r) \Bigr]
\label{uchquad}
\\
%
% u da cb
%
%U_{D_{\bf a}C_{\bf b}}(r)=&
%-\frac{C_{\bf b}}{4\pi \epsilon_0}  
%\left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)   v_{11}(r) \label{udipch}
%\\
%
% u da db
%
U_{D_{\bf a}D_{\bf b}}(r)=&
-\frac{1}{4\pi \epsilon_0} \Bigr[ \left( \mathbf{D}_{\mathbf {a}} \cdot
\mathbf{D}_{\mathbf{b}} \right)  v_{21}(r)
+\left( \mathbf{D}_{\mathbf {a}} \cdot \hat{r} \right)
\left( \mathbf{D}_{\mathbf {b}} \cdot \hat{r} \right)  
v_{22}(r) \Bigr]
\label{udipdip}
\\
%
% u da qb
%
\begin{split}
% 1
U_{D_{\bf a}Q_{\bf b}}(r) =&
-\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{b}}
\left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
+2 ( \mathbf{D}_{\mathbf{a}} \cdot
\mathbf{Q}_{\mathbf{b}} \cdot \hat{r} ) \Bigr] v_{31}(r) \\
% 2
&-\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
 \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{32}(r)
\label{udipquad}
\end{split}
\\
%
% u qa cb
%
%U_{Q_{\bf a}C_{\bf b}}(r)=&
%\frac{C_{\bf b }}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\bf a}  v_{21}(r)
%\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right)  v_{22}(r)  \Bigr]
%\label{uquadch}
%\\
%
% u qa db
%
%\begin{split}
%1
%U_{Q_{\bf a}D_{\bf b}}(r)=&
%\frac{1}{4\pi \epsilon_0} \Bigl[
%\text{Tr}\mathbf{Q}_{\mathbf{a}}
%\left(  \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)
%+2 ( \mathbf{D}_{\mathbf{b}} \cdot 
%\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r}) \Bigr] v_{31}(r)\\
% 2
%&+\frac{1}{4\pi \epsilon_0}
%\left(  \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)
%\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) v_{32}(r)
%\label{uquaddip}
%\end{split}
%\\
%
% u qa qb
%
\begin{split}
%1
U_{Q_{\bf a}Q_{\bf b}}(r)=&
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr} \mathbf{Q}_{\mathbf{a}} \text{Tr} \mathbf{Q}_{\mathbf{b}}
+2 \text{Tr} \left(
\mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} \right) \Bigr] v_{41}(r)
\\
% 2
&+ \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
 \left( \hat{r} \cdot 
\mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right)
+\text{Tr}\mathbf{Q}_{\mathbf{b}}
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}}
 \cdot \hat{r} \right)  +4 (\hat{r}  \cdot
\mathbf{Q}_{{\mathbf a}}\cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
\Bigr] v_{42}(r)
 \\
% 4
&+ \frac{1}{4\pi \epsilon_0}
\left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right)
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}}  \cdot \hat{r} \right) v_{43}(r).
\label{uquadquad}
\end{split}
\end{align}

Note that the energies of multipoles on site $\mathbf{b}$ interacting
with those on site $\mathbf{a}$ can be obtained by swapping indices
along with the sign of the intersite vector, $\hat{r}$.

%
%
% TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
%

\begin{table*}
\caption{\label{tab:tableenergy}Radial functions used in the energy and torque equations.  Functions
used in this table are defined in Appendices B and C.}
\begin{ruledtabular}
\begin{tabular}{|l|c|l|l}
Generic&Coulomb&Taylor-Shifted&Gradient-Shifted
\\ \hline
%
%
%
%Ch-Ch&
$v_{01}(r)$ &
$\frac{1}{r}$ &
$f_0(r)$ &
$f(r)-f(r_c)-(r-r_c)g(r_c)$
\\
%
%
%
%Ch-Di&
$v_{11}(r)$ &
$-\frac{1}{r^2}$ &
$g_1(r)$ &
$g(r)-g(r_c)-(r-r_c)h(r_c)$ \\
%
%
%
%Ch-Qu/Di-Di&
$v_{21}(r)$ &
$-\frac{1}{r^3}  $ & 
$\frac{g_2(r)}{r} $ &
$\frac{g(r)}{r}-\frac{g(r_c)}{r_c} -(r-r_c)
\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right)$ \\
$v_{22}(r)$ &
$\frac{3}{r^3}  $ &
$\left(-\frac{g_2(r)}{r} + h_2(r) \right)$ &
$\left(-\frac{g(r)}{r}+h(r) \right)
-\left(-\frac{g(r_c)}{r_c}+h(r_c) \right) $ \\
&&&$  -(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$ 
\\
%
%
%
%Di-Qu &
$v_{31}(r)$ &
$\frac{3}{r^4}  $ &
$\left(-\frac{g_3(r)}{r^2} + \frac{h_3(r)}{r} \right)$ &
$\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) $\\
&&&$  -(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)$ 
\\
% 
$v_{32}(r)$ &
$-\frac{15}{r^4}  $ &
$\left( \frac{3g_3(r)}{r^2} - \frac{3h_3(r)}{r} + s_3(r) \right)$ &
$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)
- \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c) \right)$ \\
&&&$  -(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)$ 
\\
%
%
%
%Qu-Qu&
$v_{41}(r)$ &
$\frac{3}{r^5} $ &
$\left(-\frac{g_4(r)}{r^3} +\frac{h_4(r)}{r^2} \right) $  &
$\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)$ \\
&&&$  -(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)$ 
\\
% 2
$v_{42}(r)$ &
$- \frac{15}{r^5}   $ &
$\left( \frac{3g_4(r)}{r^3} - \frac{3h_4(r)}{r^2}+\frac{s_4(r)}{r} \right)$ &
$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right)
-\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ \\
&&&$  -(r-r_c) \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)$ 
\\
% 3
$v_{43}(r)$ &
$ \frac{105}{r^5}  $ &
$\left(-\frac{15g_4(r)}{r^3}+\frac{15h_4(r)}{r^2}-\frac{6s_4(r)}{r} + t_4(r)\right) $ &
$\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s(r)}{r} + t(r)\right)$ \\
&&&$ -\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)$ \\
&&&$ -(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}
-\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\
\end{tabular}
\end{ruledtabular}
\end{table*}
%
%
% FORCE  TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
%

\begin{table}
\caption{\label{tab:tableFORCE}Radial functions used in the force equations.}
\begin{ruledtabular}
\begin{tabular}{cc}
Generic&Method 1 or Method 2
\\ \hline
%
%
%
$w_a(r)$&
$\frac{d v_{01}}{dr}$ \\
%
%
$w_b(r)$ &
$\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $ \\
%
$w_c(r)$ &
$\frac{v_{11}(r)}{r}$ \\
%
%
$w_d(r)$&
$\frac{d v_{21}}{dr}$ \\
%
$w_e(r)$ &
$\frac{v_{22}(r)}{r}$ \\
% 
%
$w_f(r)$&
$\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$\\
%
$w_g(r)$& 
$\frac{v_{31}(r)}{r}$\\
%
$w_h(r)$ &
$\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$\\
% 2
$w_i(r)$ &
$\frac{v_{32}(r)}{r}$ \\
% 
$w_j(r)$ &
$\frac{d v_{32}}{dr}  - \frac{3v_{32}}{r}$ \\
%
$w_k(r)$ &
$\frac{d v_{41}}{dr} $  \\
%
$w_l(r)$ &
$\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ \\
%
$w_m(r)$ &
$\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$ \\
%
$w_n(r)$ &
$\frac{v_{42}(r)}{r}$ \\
%
$w_o(r)$ &
$\frac{v_{43}(r)}{r}$ \\
%

\end{tabular}
\end{ruledtabular}
\end{table}
%
%
%

\subsection{Forces}

The force $\mathbf{F}_{\bf a}$ on $\bf{a}$ due to $\bf{b}$ is the negative of
the force  $\mathbf{F}_{\bf b}$ on $\bf{b}$ due to $\bf{a}$.  For a simple charge-charge
interaction, these forces will point along the $\pm \hat{r}$ directions, where 
$\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $.  Thus
%
\begin{equation}
F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}}{\partial r} 
\quad \text{and} \quad  F_{\bf b \alpha} 
= - \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}} {\partial r}  .
\end{equation}
%
The concept of obtaining a force from an energy by taking a gradient is the same for
higher-order multipole interactions, the trick is to make sure that all
$r$-dependent derivatives are considered.
As is pointed out by Allen and Germano,\cite{Allen:2006fk} this is straightforward if the
interaction energies are written recognizing explicit 
$\hat{r}$ and body axes ($\hat{a}_m$, $\hat{b}_n$) dependences:
%
\begin{equation}
U(r,\{\hat{a}_m \cdot \hat{r} \}, 
\{\hat{b}_n\cdot \hat{r} \} 
\{\hat{a}_m \cdot \hat{b}_n \}) .
\label{ugeneral}
\end{equation}
%
Then,
%
\begin{equation}
\mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} =  \frac{\partial U}{\partial \mathbf{r}} 
= \frac{\partial U}{\partial r} \hat{r} 
 + \sum_m \left[ 
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} 
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} 
+ \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} 
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} 
\right] \label{forceequation}.
\end{equation}
%
Note our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ is opposite
that of Allen and Germano.\cite{Allen:2006fk}  In simplifying the algebra, we also use:
%
\begin{eqnarray}
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} 
= \frac{1}{r} \left(  \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r}
\right) \\
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} 
= \frac{1}{r} \left(  \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r} 
\right) .
\end{eqnarray}
%
We list below the force equations written in terms  of space coordinates.  The
radial functions used in the two methods are listed in Table II.
%
%SPACE COORDINATES FORCE EQUTIONS
%
% **************************************************************************
% f ca cb
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} = 
\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r} 
\end{equation}
%
%
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} = 
\frac{C_{\bf a}}{4\pi \epsilon_0} \Bigl[
 \left( \hat{r} \cdot \mathbf{D}_{\mathbf{b}} \right) 
w_b(r) \hat{r}  
+ \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr]
\end{equation}
%
%
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} = 
\frac{C_{\bf a }}{4\pi \epsilon_0} \Bigr[
\text{Tr}\mathbf{Q}_{\bf b} w_d(r) \hat{r}
+ 2  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} w_e(r)
 + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
\end{equation}
%
%
%
\begin{equation}
\mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} = 
-\frac{C_{\bf{b}}}{4\pi \epsilon_0} \Bigl[
\left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) w_b(r) \hat{r}
+ \mathbf{D}_{\mathbf{a}} w_c(r) \Bigr]
\end{equation}
%
%
%
\begin{equation}
\mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
- \mathbf{D}_{\mathbf {a}} \cdot  \mathbf{D}_{\mathbf{b}} w_d(r) \hat{r}
+ \left( \mathbf{D}_{\mathbf {a}} 
\left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)
+ \mathbf{D}_{\mathbf {b}} \left( \mathbf{D}_{\mathbf{a}}  \cdot \hat{r} \right) \right) w_e(r)
% 2
 - \left( \hat{r} \cdot \mathbf{D}_{\mathbf {a}} \right) 
\left( \hat{r} \cdot \mathbf{D}_{\mathbf {b}} \right) w_f(r) \hat{r} \Bigr]
\end{equation}
%
%
%
\begin{equation}
\begin{split}
\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} = 
& - \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{b}} \mathbf{ D}_{\mathbf{a}}
+2 \mathbf{D}_{\mathbf{a}} \cdot 
\mathbf{Q}_{\mathbf{b}} \Bigr] w_g(r)
 - \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{b}}
\left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) 
+2 ( \mathbf{D}_{\mathbf{a}} \cdot 
\mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
% 3
& - \frac{1}{4\pi \epsilon_0} \Bigl[\mathbf{ D}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) 
+2  (\hat{r} \cdot \mathbf{D}_{\mathbf{a}} ) (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} )  \Bigr]
w_i(r)
% 4
 -\frac{1}{4\pi \epsilon_0} 
(\hat{r} \cdot \mathbf{D}_{\mathbf{a}} )
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r}  
\end{split}
\end{equation}
%
%
\begin{equation}
\mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} = 
\frac{C_{\bf b }}{4\pi \epsilon_0} \Bigr[
\text{Tr}\mathbf{Q}_{\bf a} w_d(r) \hat{r}
+ 2  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} w_e(r)
 + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
\end{equation}
%
\begin{equation}
\begin{split}
\mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} = 
&\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{a}} \mathbf{D}_{\mathbf{b}} 
+2 \mathbf{D}_{\mathbf{b}} \cdot \mathbf{Q}_{\mathbf{a}}  \Bigr] w_g(r)
% 2
+ \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
(\hat{r} \cdot  \mathbf{D}_{\mathbf{b}})
+2 (\mathbf{D}_{\mathbf{b}} \cdot 
\mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
% 3
&+ \frac{1}{4\pi \epsilon_0} \Bigl[ \mathbf{D}_{\mathbf{b}} 
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) 
+2  (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
(\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} ) \Bigr]   w_i(r)
% 4
+\frac{1}{4\pi \epsilon_0} 
(\hat{r} \cdot \mathbf{D}_{\mathbf{b}}) 
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}}  \cdot \hat{r}) w_j(r) \hat{r} 
\end{split}
\end{equation}
%
%
%
\begin{equation}
\begin{split}
\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} = 
+\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} \hat{r}
+ 2 \text{Tr} ( \mathbf{Q}_{\mathbf{a}} \cdot  \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_k(r) \hat{r} \\
% 2
+\frac{1}{4\pi \epsilon_0} \Bigl[
2\text{Tr}\mathbf{Q}_{\mathbf{b}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )  
+ 2\text{Tr}\mathbf{Q}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} ) 
% 3
+4 (\mathbf{Q}_{\mathbf{a}}  \cdot  \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})  
+  4(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}}) \Bigr] w_n(r) \\
% 4
+ \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{a}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) 
+ \text{Tr}\mathbf{Q}_{\mathbf{b}}
(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})  
% 5
+4 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot  
\mathbf{Q}_{\mathbf{b}}   \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
%
+ \frac{1}{4\pi \epsilon_0} \Bigl[
+ 2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
%6
+2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \hat{r})
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_o(r) \\
%  7
+ \frac{1}{4\pi \epsilon_0} 
(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r}) 
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) w_m(r) \hat{r} 
\end{split}
\end{equation}
%
%
% TORQUES SECTION -----------------------------------------------------------------------------------------
%
\subsection{Torques}

Following again Allen and Germano,\cite{Allen:2006fk} when energies are written in the form
of Eq.~({\ref{ugeneral}), then torques can be expressed as:
%
\begin{eqnarray}
\mathbf{\tau}_{\bf a} =
 \sum_m 
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} 
( \hat{r} \times \hat{a}_m )
-\sum_{mn}
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} 
(\hat{a}_m \times \hat{b}_n) \\
%
\mathbf{\tau}_{\bf b} =
 \sum_m 
\frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} 
( \hat{r} \times \hat{b}_m)
+\sum_{mn}
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} 
(\hat{a}_m \times \hat{b}_n) .
\end{eqnarray}
%
%
Here we list the torque equations written in terms of space coordinates.
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} = 
\frac{C_{\bf a}}{4\pi \epsilon_0}  (\hat{r} \times  \mathbf{D}_{\mathbf{b}}) v_{11}(r) 
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} = 
\frac{2C_{\bf a}}{4\pi \epsilon_0}
\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{22}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =  
-\frac{C_{\bf b}}{4\pi \epsilon_0}  
(\hat{r} \times \mathbf{D}_{\mathbf{a}})  v_{11}(r) 
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} = 
\frac{1}{4\pi \epsilon_0}  \mathbf{D}_{\mathbf {a}}  \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
% 2
-\frac{1}{4\pi \epsilon_0} 
(\hat{r} \times \mathbf{D}_{\mathbf {a}} )
(\hat{r} \cdot \mathbf{D}_{\mathbf {b}} )  v_{22}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} = 
-\frac{1}{4\pi \epsilon_0} \mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
% 2
+\frac{1}{4\pi \epsilon_0} 
(\hat{r} \cdot \mathbf{D}_{\mathbf {a}} )
(\hat{r} \times \mathbf{D}_{\mathbf {b}} ) v_{22}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
-\text{Tr}\mathbf{Q}_{\mathbf{b}}
(\hat{r} \times \mathbf{D}_{\mathbf{a}} )
+2 \mathbf{D}_{\mathbf{a}}  \times 
(\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
\Bigr] v_{31}(r)
% 3
-\frac{1}{4\pi \epsilon_0}
\ (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =
 \frac{1}{4\pi \epsilon_0} \Bigl[
+2 ( \mathbf{D}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} ) \times
\hat{r} 
-2 \mathbf{D}_{\mathbf{a}}  \times 
(\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
\Bigr] v_{31}(r)
% 2
+\frac{2}{4\pi \epsilon_0}
(\hat{r} \cdot \mathbf{D}_{\mathbf{a}})
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} =
\frac{1}{4\pi \epsilon_0} \Bigl[
-2 (\mathbf{D}_{\mathbf{b}}  \cdot \mathbf{Q}_{\mathbf{a}} ) \times \hat{r} 
+2 \mathbf{D}_{\mathbf{b}}  \times 
(\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
\Bigr] v_{31}(r)
% 3
- \frac{2}{4\pi \epsilon_0}
(\hat{r} \cdot \mathbf{D}_{\mathbf{b}} )
(\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}}) \times \hat{r} v_{32}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{a}}
(\hat{r} \times \mathbf{D}_{\mathbf{b}} )
+2 \mathbf{D}_{\mathbf{b}}  \times
( \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r)
% 2
+\frac{1}{4\pi \epsilon_0}
(\hat{r} \times \mathbf{D}_{\mathbf{b}} )
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) v_{32}(r)
\end{equation}
%
%
%
\begin{equation}
\begin{split}
\mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
&-\frac{4}{4\pi \epsilon_0} 
\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}}
v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[-2\text{Tr}\mathbf{Q}_{\mathbf{b}}
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times \hat{r}
+4 \hat{r} \times 
( \mathbf{Q}_{{\mathbf a}} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) 
% 3
-4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} )\times 
( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} ) \Bigr] v_{42}(r) \\
% 4
&+ \frac{2}{4\pi \epsilon_0}
\hat{r} \times ( \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
\end{split}
\end{equation}
%
%
%
\begin{equation}
\begin{split}
\mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} =  
&\frac{4}{4\pi \epsilon_0} 
\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}} v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} \Bigl[- 2\text{Tr}\mathbf{Q}_{\mathbf{a}}
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} ) \times \hat{r}
-4  (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot 
\mathbf{Q}_{{\mathbf b}} ) \times
\hat{r} 
+4 ( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times 
( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
\Bigr] v_{42}(r) \\
% 4
&+ \frac{2}{4\pi \epsilon_0}
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
\end{split}
\end{equation}
%

\section{Comparison to known multipolar energies}

To understand how these new real-space multipole methods behave in
computer simulations, it is vital to test against established methods
for computing electrostatic interactions in periodic systems, and to
evaluate the size and sources of any errors that arise from the
real-space cutoffs. In this paper we test Taylor-shifted and
Gradient-shifted electrostatics against analytical methods for
computing the energies of ordered multipolar arrays.  In the following
paper, we test the new methods against the multipolar Ewald sum for
computing the energies, forces and torques for a wide range of typical
condensed-phase (disordered) systems.

Because long-range electrostatic effects can be significant in
crystalline materials, ordered multipolar arrays present one of the
biggest challenges for real-space cutoff methods.  The dipolar
analogues to the Madelung constants were first worked out by Sauer,
who computed the energies of ordered dipole arrays of zero
magnetization and obtained a number of these constants.\cite{Sauer}
This theory was developed more completely by Luttinger and
Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays and
other periodic structures.  We have repeated the Luttinger \& Tisza
series summations to much higher order and obtained the following
energy constants (converged to one part in $10^9$):
\begin{table*}
\centering{
  \caption{Luttinger \& Tisza arrays and their associated
    energy constants. Type "A" arrays have nearest neighbor strings of
    antiparallel dipoles.  Type "B" arrays have nearest neighbor
    strings of antiparallel dipoles if the dipoles are contained in a
    plane perpendicular to the dipole direction that passes through
    the dipole.}
}
\label{tab:LT}
\begin{ruledtabular}
\begin{tabular}{cccc}
Array Type &  Lattice &   Dipole Direction &    Energy constants \\ \hline
   A     &      SC       &      001         &      -2.676788684 \\
   A     &      BCC      &      001         &       0 \\
   A     &      BCC      &      111         &      -1.770078733 \\
   A     &      FCC      &      001         &       2.166932835 \\
   A     &      FCC      &      011         &      -1.083466417 \\

   *     &      BCC      &    minimum       &      -1.985920929 \\

   B     &      SC        &     001         &      -2.676788684 \\
   B     &      BCC       &     001         &      -1.338394342 \\
   B     &      BCC       &     111         &     -1.770078733 \\
   B     &      FCC       &     001         &      -1.083466417 \\
   B     &      FCC       &     011         &      -1.807573634 
\end{tabular}
\end{ruledtabular}
\end{table*}

In addition to the A and B arrays, there is an additional minimum
energy structure for the BCC lattice that was found by Luttinger \&
Tisza.  The total electrostatic energy for an array is given by:
\begin{equation}
  E = C N^2 \mu^2
\end{equation}
where $C$ is the energy constant given above, $N$ is the number of
dipoles per unit volume, and $\mu$ is the strength of the dipole.

{\it Quadrupolar} analogues to the Madelung constants were first worked out by Nagai and Nakamura who
computed the energies of selected quadrupole arrays based on
extensions to the Luttinger and Tisza
approach.\cite{Nagai01081960,Nagai01091963}  We have compared the
energy constants for the lowest energy configurations for linear
quadrupoles shown in table \ref{tab:NNQ}

\begin{table*}
\centering{
  \caption{Nagai and Nakamura Quadurpolar arrays}}
\label{tab:NNQ}
\begin{ruledtabular}
\begin{tabular}{ccc}
 Lattice &   Quadrupole Direction &    Energy constants \\ \hline
  SC       &      111         &      -8.3 \\
  BCC      &      011         &      -21.7 \\
  FCC      &      111         &      -80.5 
\end{tabular}
\end{ruledtabular}
\end{table*}

In analogy to the dipolar arrays, the total electrostatic energy for
the quadrupolar arrays is:
\begin{equation}
 E = C \frac{3}{4} N^2 Q^2
\end{equation} 
where $Q$ is the quadrupole moment.


\begin{acknowledgments}
Support for this project was provided by the National Science
Foundation under grant CHE-0848243. Computational time was provided by
the Center for Research Computing (CRC) at the University of Notre
Dame.
\end{acknowledgments}

\newpage
\appendix

\section{Smith's $B_l(r)$ functions for damped-charge distributions}

The following summarizes Smith's $B_l(r)$ functions and includes
formulas given in his appendix.\cite{Smith98} The first function
$B_0(r)$ is defined by
%
\begin{equation}
B_0(r)=\frac{\text{erfc}(\alpha r)}{r} = \frac{2}{\sqrt{\pi}r}=
\int_{\alpha r}^{\infty} \text{e}^{-s^2} ds .
\end{equation}
%
The first derivative of this function is
%
\begin{equation}
\frac{dB_0(r)}{dr}=-\frac{1}{r^2}\text{erfc}(\alpha r)
-\frac{2\alpha}{r\sqrt{\pi}}\text{e}^{-{\alpha}^2r^2}
\end{equation}
%
which can be used to define a function $B_1(r)$:
%
\begin{equation}
B_1(r)=-\frac{1}{r}\frac{dB_0(r)}{dr}
\end{equation}
%
In general, the recurrence relation,
\begin{equation}
B_l(r)=-\frac{1}{r}\frac{dB_{l-1}(r)}{dr} 
= \frac{1}{r^2} \left[ (2l-1)B_{l-1}(r) + \frac {(2\alpha^2)^l}{\alpha \sqrt{\pi}}
\text{e}^{-{\alpha}^2r^2}
\right] ,
\end{equation}
is very useful for building up higher derivatives.  Using these
formulas, we find:
%
\begin{align}
\frac{dB_0}{dr}=&-rB_1(r) \\
\frac{d^2B_0}{dr^2}=&    - B_1(r) + r^2 B_2(r) \\
\frac{d^3B_0}{dr^3}=&   3 r B_2(r) - r^3 B_3(r) \\
\frac{d^4B_0}{dr^4}=&   3 B_2(r) - 6 r^2 B_3(r) + r^4 B_4(r) \\
\frac{d^5B_0}{dr^5}=& - 15 r B_3(r) + 10 r^3 B_4(r) - r^5 B_5(r) .
\end{align}
%
As noted by Smith, it is possible to approximate the $B_l(r)$
functions, 
%
\begin{equation}
B_l(r)=\frac{(2l)!}{l!2^lr^{2l+1}} - \frac {(2\alpha^2)^{l+1}}{(2l+1)\alpha \sqrt{\pi}}
+\text{O}(r) .
\end{equation}
\newpage
\section{The $r$-dependent factors for TSF electrostatics}

Using the shifted damped functions $f_n(r)$ defined by:
%
\begin{equation}
f_n(r)= B_0(r) -\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} B_0^{(m)}(r_c)  ,
\end{equation}
%
where the superscript $(m)$ denotes the $m^\mathrm{th}$ derivative. In
this Appendix, we provide formulas for successive derivatives of this
function.  (If there is no damping, then $B_0(r)$ is replaced by
$1/r$.)  First, we find:
%
\begin{equation}
\frac{\partial f_n}{\partial r_\alpha}=\hat{r}_\alpha \frac{d f_n}{d r} .
\end{equation}
%
This formula clearly brings in derivatives of Smith's $B_0(r)$
function, and we define higher-order derivatives as follows:
%
\begin{align}
g_n(r)=& \frac{d f_n}{d r} =
B_0^{(1)}(r) -\sum_{m=0}^{n} \frac {(r-r_c)^m}{m!} B_0^{(m+1)}(r_c) \\
h_n(r)=& \frac{d^2f_n}{d r^2} =
B_0^{(2)}(r) -\sum_{m=0}^{n-1} \frac {(r-r_c)^m}{m!} B_0^{(m+2)}(r_c) \\
s_n(r)=& \frac{d^3f_n}{d r^3} =
B_0^{(3)}(r) -\sum_{m=0}^{n-2} \frac {(r-r_c)^m}{m!} B_0^{(m+3)}(r_c) \\
t_n(r)=& \frac{d^4f_n}{d r^4} =
B_0^{(4)}(r) -\sum_{m=0}^{n-3} \frac {(r-r_c)^m}{m!} B_0^{(m+4)}(r_c) \\
u_n(r)=& \frac{d^5f_n}{d r^5} =
B_0^{(5)}(r) -\sum_{m=0}^{n-4} \frac {(r-r_c)^m}{m!} B_0^{(m+5)}(r_c)  .
\end{align}
%
We note that the last function needed (for quadrupole-quadrupole interactions) is
%
\begin{equation}
u_4(r)=B_0^{(5)}(r) - B_0^{(5)}(r_c) .
\end{equation}

The functions $f_n(r)$ to $u_n(r)$ can be computed recursively and
stored on a grid for values of $r$ from $0$ to $r_c$.  The functions
needed are listed schematically below:
%
\begin{eqnarray}
f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\
g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\
h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\
s_2 \quad s_3 \quad &s_4 \nonumber \\
t_3 \quad &t_4 \nonumber \\
&u_4 \nonumber .
\end{eqnarray}

Using these functions, we find
%
\begin{align}
\frac{\partial f_n}{\partial r_\alpha} =&r_\alpha \frac {g_n}{r} \label{eq:b9}\\
\frac{\partial^2 f_n}{\partial r_\alpha \partial r_\beta} =&\delta_{\alpha \beta}\frac {g_n}{r} 
+r_\alpha r_\beta \left( -\frac{g_n}{r^3} +\frac{h_n}{r^2}\right) \\
\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta r_\gamma} =&
\left( \delta_{\alpha \beta} r_\gamma + \delta_{\alpha \gamma} r_\beta + 
\delta_{ \beta \gamma} r_\alpha \right)  
\left(  -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right)
+ r_\alpha r_\beta r_\gamma 
\left(  \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \\
\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta r_\gamma r_\delta} =& 
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} 
+ \delta_{\alpha \gamma} \delta_{\beta \delta}
 +\delta_{ \beta \gamma} \delta_{\alpha \delta} \right)
\left( - \frac{g_n}{r^3} + \frac{h_n}{r^2} \right)  \nonumber \\ 
&+ \left( \delta_{\alpha \beta} r_\gamma r_\delta
+ \text{5 permutations}
\right) \left( \frac{3 g_n}{r^5} - \frac{3h_n}{r^4} + \frac{s_n}{r^3} 
\right) \nonumber \\
&+ r_\alpha r_\beta r_\gamma r_\delta
\left(  -\frac{15g_n}{r^7} + \frac{15h_n}{r^6} - \frac{6s_n}{r^5}
+ \frac{t_n}{r^4} \right)\\
\frac{\partial^5 f_n}
{\partial r_\alpha \partial r_\beta r_\gamma r_\delta r_\epsilon} =& 
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} r_\epsilon
+ \text{14 permutations} \right) 
\left(  \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \nonumber \\
&+ \left( \delta_{\alpha \beta} r_\gamma r_\delta r_\epsilon
+ \text{9 permutations}
\right) \left(- \frac{15g_n}{r^7}+\frac{15h_n}{r^7} -\frac{6s_n}{r^5} +\frac{t_n}{r^4} 
\right) \nonumber \\
&+ r_\alpha r_\beta r_\gamma r_\delta r_\epsilon
\left(  \frac{105g_n}{r^9} - \frac{105h_n}{r^8} + \frac{45s_n}{r^7}
- \frac{10t_n}{r^6} +\frac{u_n}{r^5} \right) \label{eq:b13}
\end{align}
%
%
%
\newpage
\section{The $r$-dependent factors for GSF electrostatics}

In Gradient-shifted force electrostatics, the kernel is not expanded,
rather the individual terms in the multipole interaction energies.
For damped charges , this still brings into the algebra multiple
derivatives of the Smith's $B_0(r)$ function.  To denote these terms,
we generalize the notation of the previous appendix.  For $f(r)=1/r$
(bare Coulomb) or $f(r)=B_0(r)$ (smeared charge)
%
\begin{align}
g(r)=& \frac{df}{d r}\\
h(r)=& \frac{dg}{d r} = \frac{d^2f}{d r^2} \\
s(r)=& \frac{dh}{d r} = \frac{d^3f}{d r^3} \\
t(r)=& \frac{ds}{d r} = \frac{d^4f}{d r^4} \\
u(r)=& \frac{dt}{d r} = \frac{d^5f}{d r^5} .
\end{align}
%
For undamped charges, $f(r)=1/r$, Table I lists these derivatives
under the column ``Bare Coulomb.''  Equations \ref{eq:b9} to
\ref{eq:b13} are still correct for GSF electrostatics if the subscript
$n$ is eliminated.

\section{Extra Material}
%
%
%Energy in body coordinate form ---------------------------------------------------------------
%
Here are the interaction energies written in terms of the body coordinates:

%
% u ca cb
%
\begin{equation}
U_{C_{\bf a}C_{\bf b}}(r)=
\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r) 
\end{equation}
%
% u ca db
%
\begin{equation}
U_{C_{\bf a}D_{\bf b}}(r)=
\frac{C_{\bf a}}{4\pi \epsilon_0} 
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} \,  v_{11}(r) 
\end{equation}
%
% u ca qb
%
\begin{equation}
U_{C_{\bf a}Q_{\bf b}}(r)=
\frac{C_{\bf a }\text{Tr}Q_{\bf b}}{4\pi \epsilon_0}  
v_{21}(r) \nonumber \\
+\frac{C_{\bf a}}{4\pi \epsilon_0}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r})
v_{22}(r) 
\end{equation}
%
% u da cb
%
\begin{equation}
U_{D_{\bf a}C_{\bf b}}(r)=
-\frac{C_{\bf b}}{4\pi \epsilon_0}  
\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} \,  v_{11}(r) 
\end{equation}
%
% u da db
%
\begin{equation}
\begin{split}
% 1
U_{D_{\bf a}D_{\bf b}}(r)&=
-\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m} 
(\hat{a}_m \cdot   \hat{b}_n)
D_{\mathbf{b}n} v_{21}(r) \\
% 2
&-\frac{1}{4\pi \epsilon_0} 
\sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m} 
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} 
v_{22}(r)
\end{split}
\end{equation}
%
% u da qb
%
\begin{equation}
\begin{split}
% 1
U_{D_{\bf a}Q_{\bf b}}(r)&=
-\frac{1}{4\pi \epsilon_0} \left(
\text{Tr}Q_{\mathbf{b}}
\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
+2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \cdot \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
\right)  v_{31}(r) \\
% 2
&-\frac{1}{4\pi \epsilon_0}
\sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{b}_n \cdot \hat{r}) v_{32}(r)
\end{split}
\end{equation}
%
% u qa cb
%
\begin{equation}
U_{Q_{\bf a}C_{\bf b}}(r)=
\frac{C_{\bf b }\text{Tr}Q_{\bf a}}{4\pi \epsilon_0}  v_{21}(r)
+\frac{C_{\bf b}}{4\pi \epsilon_0}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) v_{22}(r)
\end{equation}
%
% u qa db
%
\begin{equation}
\begin{split}
%1
U_{Q_{\bf a}D_{\bf b}}(r)&=
\frac{1}{4\pi \epsilon_0} \left(
\text{Tr}Q_{\mathbf{a}}
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n}
+2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{b}_l \cdot \hat{a}_m)
Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
\right) v_{31}(r)  \\
% 2
&+\frac{1}{4\pi \epsilon_0}
\sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{a}_n \cdot \hat{r}) v_{32}(r)
\end{split}
\end{equation}
%
% u qa qb
%
\begin{equation}
\begin{split}
%1
U_{Q_{\bf a}Q_{\bf b}}(r)&=
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}}
+2\sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p) 
Q_{\mathbf{a}lm}  Q_{\mathbf{b}np} 
(\hat{a}_m \cdot \hat{b}_n) \Bigr] 
v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[ \text{Tr}Q_{\mathbf{a}}
\sum_{lm} (\hat{r} \cdot \hat{b}_l )
Q_{{\mathbf b}lm}
(\hat{b}_m \cdot \hat{r})
+\text{Tr}Q_{\mathbf{b}}
\sum_{lm} (\hat{r} \cdot \hat{a}_l )
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{r}) \\
% 3
&+4 \sum_{lmnp}
(\hat{r} \cdot \hat{a}_l ) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{b}_n)
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r})
\Bigr] v_{42}(r)  \\
% 4
&+ \frac{1}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \cdot \hat{a}_l) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{r})
\sum_{np}  (\hat{r} \cdot \hat{b}_n) 
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r})  v_{43}(r).
\end{split}
\end{equation}
%


% BODY coordinates force equations --------------------------------------------
%
%
Here are the force equations written in terms of body coordinates.
%
% f ca cb
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} = 
\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r}
\end{equation}
%
% f ca db
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} = 
\frac{C_{\bf a}}{4\pi \epsilon_0}  
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} w_b(r) \hat{r}
+\frac{C_{\bf a}}{4\pi \epsilon_0}  
\sum_n  D_{\mathbf{b}n} \hat{b}_n w_c(r)
\end{equation}
%
% f ca qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} = 
\frac{1}{4\pi \epsilon_0}  
C_{\bf a }\text{Tr}Q_{\bf b} w_d(r) \hat{r}
+ 2C_{\bf a } \sum_l  \hat{b}_l Q_{{\mathbf b}ln} (\hat{b}_n \cdot \hat{r}) w_e(r) \\
% 2
+\frac{C_{\bf a}}{4\pi \epsilon_0}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r}) w_f(r) \hat{r} 
\end{split}
\end{equation}
%
% f da cb
%
\begin{equation}
\mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} = 
-\frac{C_{\bf{b}}}{4\pi \epsilon_0} 
\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} w_b(r)  \hat{r}
-\frac{C_{\bf{b}}}{4\pi \epsilon_0} 
\sum_n  D_{\mathbf{a}n} \hat{a}_n w_c(r)
\end{equation}
% 
% f da db
%
\begin{equation}
\begin{split}
% 1
\mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} &= 
-\frac{1}{4\pi \epsilon_0}
 \sum_{mn} D_{\mathbf {a}m} 
(\hat{a}_m \cdot   \hat{b}_n)
D_{\mathbf{b}n}  w_d(r) \hat{r}
-\frac{1}{4\pi \epsilon_0} 
\sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m} 
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} w_f(r) \hat{r} \\
% 2
& \quad + \frac{1}{4\pi \epsilon_0} 
\Bigl[ \sum_m D_{\mathbf {a}m} 
\hat{a}_m \sum_n D_{\mathbf{b}n} 
(\hat{b}_n \cdot \hat{r}) 
+ \sum_m D_{\mathbf {b}m} 
\hat{b}_m \sum_n D_{\mathbf{a}n} 
(\hat{a}_n \cdot \hat{r}) \Bigr] w_e(r)  \\
\end{split}
\end{equation}
%
% f da qb
%
\begin{equation}
\begin{split}
% 1
&\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} = 
 - \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{b}}
\sum_l  D_{\mathbf{a}l} \hat{a}_l
+2\sum_{lmn} D_{\mathbf{a}l} 
(\hat{a}_l \cdot \hat{b}_m)
Q_{\mathbf{b}mn} \hat{b}_n  \Bigr] w_g(r) \\
% 3
&  - \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{b}}
\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
+2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \cdot \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
% 4
&+ \frac{1}{4\pi \epsilon_0}
\Bigl[\sum_l  D_{\mathbf{a}l} \hat{a}_l
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{b}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{a}_l)
 D_{\mathbf{a}l}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn} \hat{b}_n \Bigr]   w_i(r)\\
% 6
&  -\frac{1}{4\pi \epsilon_0} 
\sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l} 
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{b}_n \cdot \hat{r})  w_j(r)  \hat{r}
\end{split}
\end{equation}
%
% force qa cb
%
\begin{equation}
\begin{split}
% 1
\mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0}  
C_{\bf b }\text{Tr}Q_{\bf a} \hat{r} w_d(r)
+ \frac{2C_{\bf b }}{4\pi \epsilon_0}  \sum_l  \hat{a}_l Q_{{\mathbf a}ln} (\hat{a}_n \cdot \hat{r}) w_e(r) \\
% 2
&  +\frac{C_{\bf b}}{4\pi \epsilon_0}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) w_f(r) \hat{r} 
\end{split}
\end{equation}
%
% f qa db 
%
\begin{equation}
\begin{split}
% 1
&\mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{a}}
\sum_l  D_{\mathbf{b}l} \hat{b}_l
+2\sum_{lmn} D_{\mathbf{b}l} 
(\hat{b}_l \cdot \hat{a}_m)
Q_{\mathbf{a}mn} \hat{a}_n  \Bigr] 
w_g(r)\\
% 3
&  + \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{a}}
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} 
+2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{b}_l \cdot \hat{a}_m)
Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
% 4
&  + \frac{1}{4\pi \epsilon_0} \Bigl[ \sum_l  D_{\mathbf{b}l} \hat{b}_l
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{a}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{b}_l)
 D_{\mathbf{b}l}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn} \hat{a}_n \Bigr]   w_i(r) \\
% 6
&  +\frac{1}{4\pi \epsilon_0} 
\sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l} 
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{a}_n \cdot \hat{r})  w_j(r)  \hat{r}
\end{split}
\end{equation}
%
% f qa qb
%
\begin{equation}
\begin{split}
&\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}} 
+ 2 \sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p) 
Q_{\mathbf{a}lm}  Q_{\mathbf{b}np} 
(\hat{a}_m \cdot \hat{b}_n) \Bigr] w_k(r) \hat{r}\\
&+\frac{1}{4\pi \epsilon_0} \Bigl[
2\text{Tr}Q_{\mathbf{b}} \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  \hat{a}_m  
+ 2\text{Tr}Q_{\mathbf{a}} \sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm}  \hat{b}_m \\
&+ 4\sum_{lmnp} \hat{a}_l Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_p \cdot \hat{r})   
+ 4\sum_{lmnp} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{b}_n) Q_{\mathbf{b}np} \hat{b}_p
\Bigr] w_n(r) \\
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[ \text{Tr}Q_{\mathbf{a}}
\sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm} (\hat{b}_m \cdot \hat{r}) 
+ \text{Tr}Q_{\mathbf{b}}
\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  (\hat{a}_m \cdot \hat{r}) \\
&+4\sum_{lmnp} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{b}_n) 
Q_{\mathbf{b}np}  (\hat{b}_p \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
%
&+\frac{1}{4\pi \epsilon_0} \Bigl[
2\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} \hat{a}_m
\sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_n \cdot \hat{r}) \\
&+2 \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
\sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} \hat{b}_n \Bigr] w_o(r) \hat{r} \\
&  + \frac{1}{4\pi \epsilon_0} 
\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
\sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_p \cdot \hat{r}) w_m(r) \hat{r} 
\end{split}
\end{equation}
%
Here we list the form of the non-zero damped shifted multipole torques showing
explicitly dependences on body axes:
%
%  t ca db
%
\begin{equation}
\mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} = 
\frac{C_{\bf a}}{4\pi \epsilon_0}  
\sum_n  (\hat{r} \times \hat{b}_n)  D_{\mathbf{b}n} \,  v_{11}(r) 
\end{equation}
%
% t ca qb
%
\begin{equation}
\mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} = 
\frac{2C_{\bf a}}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \times \hat{b}_l) Q_{{\mathbf b}lm} (\hat{b}_m \cdot \hat{r}) v_{22}(r)
\end{equation}
%
%  t da cb
%
\begin{equation}
\mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} = 
-\frac{C_{\bf b}}{4\pi \epsilon_0}  
\sum_n  (\hat{r} \times \hat{a}_n)  D_{\mathbf{a}n} \,  v_{11}(r) 
\end{equation}%
%
%
%  ta da db
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m} 
(\hat{a}_m \times  \hat{b}_n)
D_{\mathbf{b}n} v_{21}(r) \\
% 2
&-\frac{1}{4\pi \epsilon_0} 
\sum_m (\hat{r} \times \hat{a}_m) D_{\mathbf {a}m} 
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
\end{split}
\end{equation}
%
%  tb da db
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} &= 
-\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m} 
(\hat{a}_m \times  \hat{b}_n)
D_{\mathbf{b}n} v_{21}(r) \\
% 2
&+\frac{1}{4\pi \epsilon_0} 
\sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m} 
\sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf {b}n} v_{22}(r) 
\end{split}
\end{equation}
%
% ta da qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0} \left(
-\text{Tr}Q_{\mathbf{b}}
\sum_n (\hat{r} \times \hat{a}_n) D_{\mathbf{a}n}
+2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \times \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
\right) v_{31}(r)\\
% 2
&-\frac{1}{4\pi \epsilon_0}
\sum_l (\hat{r} \times \hat{a}_l) D_{\mathbf{a}l}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{b}_n \cdot \hat{r}) v_{32}(r)
\end{split}
\end{equation}
%
% tb da qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0} \left(
-2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \cdot \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{r} \times \hat{b}_n)
-2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \times \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
\right) v_{31}(r) \\
% 2
&-\frac{2}{4\pi \epsilon_0}
\sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{r}\times \hat{b}_n) v_{32}(r)
\end{split}
\end{equation}
%
% ta qa cb
%
\begin{equation}
\mathbf{\tau}_{{\bf a}Q_{\bf a}C_{\bf b}} = 
\frac{2C_{\bf a}}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{{\mathbf a}lm} (\hat{r} \times \hat{a}_m) v_{22}(r)
\end{equation}
%
% ta qa db
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0} \left(
2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{b}_l \cdot \hat{a}_m)
Q_{\mathbf{a}mn} (\hat{r} \times \hat{a}_n)
+2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{a}_l \times \hat{b}_m)
Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
\right) v_{31}(r) \\
% 2
&+\frac{2}{4\pi \epsilon_0}
\sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{r}\times \hat{a}_n) v_{32}(r)
\end{split}
\end{equation}
%
% tb qa db
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0} \left(
\text{Tr}Q_{\mathbf{a}}
\sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf{b}n}
+2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{a}_l \times \hat{b}_m)
Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
\right) v_{31}(r)\\
% 2
&\frac{1}{4\pi \epsilon_0}
\sum_l (\hat{r} \times \hat{b}_l) D_{\mathbf{b}l}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{a}_n \cdot \hat{r}) v_{32}(r)
\end{split}
\end{equation}
%
% ta qa qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} &= 
-\frac{4}{4\pi \epsilon_0} 
\sum_{lmnp} (\hat{a}_l \times \hat{b}_p)
Q_{\mathbf{a}lm}  Q_{\mathbf{b}np} 
(\hat{a}_m \cdot \hat{b}_n) v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[
2\text{Tr}Q_{\mathbf{b}}
\sum_{lm} (\hat{r} \cdot \hat{a}_l )
Q_{{\mathbf a}lm}
(\hat{r} \times \hat{a}_m)
+4 \sum_{lmnp}
(\hat{r} \times \hat{a}_l ) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{b}_n)
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r}) \\
% 3
&-4 \sum_{lmnp}
(\hat{r} \cdot \hat{a}_l ) 
Q_{{\mathbf a}lm}
(\hat{a}_m \times \hat{b}_n)
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r})
\Bigr] v_{42}(r) \\
% 4
&+ \frac{2}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \times \hat{a}_l) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{r})
\sum_{np}  (\hat{r} \cdot \hat{b}_n) 
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r})  v_{43}(r)\\
\end{split}
\end{equation}
%
% tb qa qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} &= 
\frac{4}{4\pi \epsilon_0} 
\sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
Q_{\mathbf{a}lm}  Q_{\mathbf{b}np} 
(\hat{a}_m \times \hat{b}_n) v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[
2\text{Tr}Q_{\mathbf{a}}
\sum_{lm} (\hat{r} \cdot \hat{b}_l )
Q_{{\mathbf b}lm}
(\hat{r} \times \hat{b}_m)
+4 \sum_{lmnp}
(\hat{r} \cdot \hat{a}_l ) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{b}_n)
Q_{{\mathbf b}np}
(\hat{r} \times \hat{b}_p) \\
% 3
&+4 \sum_{lmnp}
(\hat{r} \cdot \hat{a}_l ) 
Q_{{\mathbf a}lm}
(\hat{a}_m \times \hat{b}_n)
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r})
\Bigr] v_{42}(r)  \\
% 4
&+ \frac{2}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \cdot \hat{a}_l) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{r})
\sum_{np}  (\hat{r} \times \hat{b}_n) 
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r}) v_{43}(r). \\
\end{split}
\end{equation}
%
\begin{table*}
\caption{\label{tab:tableFORCE2}Radial functions used in the force equations.}
\begin{ruledtabular}
\begin{tabular}{|l|l|l|}
Generic&Taylor-shifted Force&Gradient-shifted Force
\\ \hline
%
%
%
$w_a(r)$&
$g_0(r)$&
$g(r)-g(r_c)$ \\
%
%
$w_b(r)$ &
$\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ &
$h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\
%
$w_c(r)$ &
$\frac{g_1(r)}{r} $ &
$\frac{v_{11}(r)}{r}$ \\
%
%
$w_d(r)$&
$\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ &
$\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) $\\
%
$w_e(r)$ &
$\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ &
$\frac{v_{22}(r)}{r}$ \\
% 
%
$w_f(r)$&
$\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ &
$\left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) - $ \\
&&$\left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)-\frac{2v_{22}(r)}{r}$\\
%
$w_g(r)$& $ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$& 
$\frac{v_{31}(r)}{r}$\\
%
$w_h(r)$ &
$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
$\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) $ \\
&&$-\frac{v_{31}(r)}{r}$\\
% 2
$w_i(r)$ &
$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
$\frac{v_{32}(r)}{r}$ \\
% 
$w_j(r)$ &
$\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right)  $ &
$\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right) $  \\
&&$\left(\frac{-6g(_cr)}{r_c^3} +\frac{6h(r_c)}{r_c^2} -\frac{3s(r_c)}{r_c} +t(r_c) \right) -\frac{3v_{32}}{r}$ \\
%
$w_k(r)$ &
$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
$\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}  \right)$  \\
&&$\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2}  \right)$ \\
%
$w_l(r)$ &
$\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)$ &
$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\
&&$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)
-\frac{2v_{42}(r)}{r}$ \\
%
$w_m(r)$ &
$\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)$ &
$\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$ \\
&&$\left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3} 
+\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $ \\ 
&&$-\frac{4v_{43}(r)}{r}$ \\
%
$w_n(r)$ &
$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
$\frac{v_{42}(r)}{r}$ \\
%
$w_o(r)$ &
$\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)$ &
$\frac{v_{43}(r)}{r}$ \\
%
\end{tabular}
\end{ruledtabular}
\end{table*}

\newpage

\bibliography{multipole}

\end{document}
%
% ****** End of file multipole.tex ******
Revision:	3984
Committed:	Sat Dec 28 18:41:19 2013 UTC (11 years, 10 months ago) by gezelter
Content type:	application/x-tex
File size:	77999 byte(s)
Log Message:	more edits