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{rotating}

%\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_{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 a diverse set of chemical
environments.\cite{doi:10.1021/la400226g,McCann:2013fk,kannam:094701,Forrest:2012ly,English:2008kx,Louden:2013ve,Tokumasu:2013zr}
DSF electrostatics provides a compromise between the computational
speed of real-space cutoffs and the accuracy of fully-periodic Ewald
treatments.

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. To do this, the interactions between coarse-grained sites
are typically taken to be point
multipoles.\cite{Golubkov06,ISI:000276097500009,ISI:000298664400012}

Water, in particular, has been modeled recently with point multipoles
up to octupolar
order.\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} For
maximum efficiency, these models require the use of an approximate
multipole expansion as the exact multipole expansion can become quite
expensive (particularly when handled via the Ewald
sum).\cite{Ichiye:2006qy} Point multipoles and multipole
polarizability have also been utilized in the 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 describe 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 somewhat 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 numerical comparisons to Ewald-based electrostatics in common
simulation enviornments.

\section{Methodology}
An efficient real-space electrostatic method involves the use of a
pair-wise functional form,
\begin{equation}
V = \sum_i \sum_{j>i} V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) +
\sum_i V_i^\mathrm{correction}
\end{equation}
that is short-ranged and easily truncated at a cutoff radius,
\begin{equation}
  V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) = \left\{
\begin{array}{ll}
V_\mathrm{approx} (r_{ij}, \Omega_i, \Omega_j) & \quad r \le r_c \\
0 & \quad r > r_c ,
\end{array}
\right.
\end{equation}
along with an easily computed correction term ($\sum_i
V_i^\mathrm{correction}$) which has linear-scaling with the number of
particles.  Here $\Omega_i$ and $\Omega_j$ represent orientational
coordinates of the two sites.  The computational efficiency, energy
conservation, and even some physical properties of a simulation can
depend dramatically on how the $V_\mathrm{approx}$ function behaves at
the cutoff radius. The goal of any approximation method should be to
mimic the real behavior of the electrostatic interactions as closely
as possible without sacrificing the near-linear scaling of a cutoff
method.

\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_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_i C_j \Biggr{[}
\frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}} 
- \frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c} + \left(\frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c^2} 
+ \frac{2\alpha}{\pi^{1/2}}
\frac{\exp\left(-\alpha^2r_c^2\right)}{r_c}
\right)\left(r_{ij}-r_c\right)\ \Biggr{]} 
\label{eq:DSFPot}
\end{equation*}
Note that in this potential and in all electrostatic quantities that
follow, the standard $1/4 \pi \epsilon_{0}$ has been omitted for
clarity.

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 also 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 extend the idea of self-neutralization for the point
multipoles by insuring net charge-neutrality and net-zero moments
within each cutoff sphere.  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.  There are also significant modifications made to make the
forces and torques go smoothly to zero at the cutoff distance.

\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 an additional 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) = 
 \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}) =\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)
= \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)=\hat{M}_{\bf a} \hat{M}_{\bf b} \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 on $\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 to
bare or damped Coulomb kernels (or any other function) 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)= 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}
which clearly vanishes as the $r$ approaches the cutoff radius. 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. For quadrupole-quadrupole interactions, we therefore
require shifted energy expressions that include up to $(r-r_c)^5$ so
that all energies, forces, and torques are zero as $r \rightarrow
r_c$. In each case, we 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)=
C_{\bf a} D_{{\bf b}\alpha}  \frac {\partial f_1(r) }{\partial r_\alpha} 
= C_{\bf a} D_{{\bf b}\alpha}
\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} = C_{\bf a} D_{{\bf b}\alpha}
\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}}{r}
\left[  -\frac{1}{r^2}+\frac{1}{r_c^2}-\frac{2(r-r_c)}{r_c^3} \right] 
+C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta 
\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= (\text{prefactor}) (\text{derivatives}) f_n(r)
\label{generic}
\end{equation}
%
with $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.  For example, in
quadrupole-quadrupole interactions for which the $\text{prefactor}$ is
$Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$, 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, $g_n(r), h_n(r), s_n(r), \mathrm{~and~} t_n(r)$ are
successive derivatives of the shifted electrostatic kernel, $f_n(r)$
of the same index $n$.  The algebra required to evaluate energies,
forces and torques is somewhat tedious, so only the final forms are
presented in tables \ref{tab:tableenergy} and \ref{tab:tableFORCE}.

\subsection{Gradient-shifted force (GSF) electrostatics}
The second, and conceptually simpler approach to force-shifting
maintains only the linear $(r-r_c)$ term in the truncated Taylor
expansion, 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} .
\label{generic2}
\end{equation}
Here the gradient for force shifting is evaluated for an image
multipole projected onto the surface of the cutoff sphere (see fig
\ref{fig:shiftedMultipoles}). 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.  Energies from this method thus have the general
form:
\begin{equation}
U= \sum  (\text{angular factor}) (\text{radial factor}).
\label{generic3}
\end{equation}

Functional forms for both methods (TSF and GSF) can both be summarized
using the form of Eq.~(\ref{generic3}).  The basic forms for the
energy, force, and torque expressions are tabulated for both shifting
approaches below -- for each separate orientational contribution, only
the radial factors differ between the two methods.

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

[XXX Do we need this section in the main paper? or should it go in the
extra materials?]

So far, all energies and forces have been written in terms of fixed
space coordinates.  Interaction energies are computed from the generic
formulas Eq.~(\ref{generic}) and ~(\ref{generic2}) which combine
orientational 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 have followed the methodology of Allen and
Germano,\cite{Allen:2006fk} which was itself based on earlier work 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}
Rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be
expressed 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 body-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$. In the torque
expressions, it is useful to have the angular-dependent terms
available in three different fashions, e.g. for the dipole-dipole
contraction:
%
\begin{equation}
 \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}}
\end{equation}
%
The first two forms are written using space coordinates.  The first
form is standard in the chemistry literature, while the second is
expressed using implied summation notation.  The third form shows
explicit sums over body indices and the dot products now indicate
contractions using space indices. 


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

In addition to cutoff-sphere neutralization, the Wolf
summation~\cite{Wolf99} and the damped shifted force (DSF)
extension~\cite{Fennell:2006zl} also 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.

First, 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}

Second, 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}

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, $g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive
derivatives of the electrostatic kernel, $f(r)$ (either the undamped
$1/r$ or the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that
have been 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 \ref{SmithFunc} contains recursion relations that allow
rapid evaluation of these derivatives.

\section{Interaction energies, forces, and torques}
The main result of this paper is a set of expressions for the
energies, forces and torques (up to quadrupole-quadrupole order) that
work for both the Taylor-shifted and Gradient-shifted approximations.
These expressions were derived using a set of generic radial
functions.  Without using the shifting approximations mentioned above,
some of these radial functions would be identical, and the expressions
coalesce into the familiar forms for unmodified multipole-multipole
interactions.  Table \ref{tab:tableenergy} maps between the generic
functions and the radial functions derived for both the Taylor-shifted
and Gradient-shifted methods.  The energy equations are written in
terms of lab-frame representations of the dipoles, quadrupoles, and
the unit vector connecting the two objects,

% Energy in space coordinate form ----------------------------------------------------------------------------------------------
%
%
% u ca cb
%
\begin{align}
U_{C_{\bf a}C_{\bf b}}(r)=&
C_{\bf a} C_{\bf b}  v_{01}(r)  \label{uchch}
\\
%
% u ca db
%
U_{C_{\bf a}D_{\bf b}}(r)=&
C_{\bf a} \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)=& C_{\bf a } \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)=&
-\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) =&
-\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
&- \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)=&
\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
&+\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
&+ 
\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{sidewaystable}
  \caption{\label{tab:tableenergy}Radial functions used in the energy
    and torque equations.  The $f, g, h, s, t, \mathrm{and} u$
    functions used in this table are defined in Appendices B and C.}
\begin{tabular}{|c|c|l|l|} \hline
Generic&Bare 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)$ \\ \hline
\end{tabular}
\end{sidewaystable}
%
%
% FORCE  TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
%

\begin{sidewaystable}
\caption{\label{tab:tableFORCE}Radial functions used in the force equations.}
\begin{tabular}{|c|c|l|l|} \hline
Function&Definition&Taylor-Shifted&Gradient-Shifted
\\ \hline
%
%
%
$w_a(r)$&
$\frac{d v_{01}}{dr}$&
$g_0(r)$&
$g(r)-g(r_c)$ \\
%
%
$w_b(r)$ &
$\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{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{v_{11}(r)}{r}$ &
$\frac{g_1(r)}{r} $ &
$\frac{v_{11}(r)}{r}$\\
%
%
$w_d(r)$&
$\frac{d v_{21}}{dr}$&
$\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}$ &
$\frac{v_{22}(r)}{r}$ \\
% 
%
$w_f(r)$&
$\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{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)$& 
$\frac{v_{31}(r)}{r}$& 
$ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$& 
$\frac{v_{31}(r)}{r}$\\
%
$w_h(r)$ &
$\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{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)$ &
$\frac{v_{32}(r)}{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)$ &
$\frac{d v_{32}}{dr}  - \frac{3v_{32}}{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)$ &
$\frac{d v_{41}}{dr} $ &
$\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)$ &
$\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{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_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)
-\frac{2v_{42}(r)}{r}$\\
%
$w_m(r)$ &
$\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{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)$ &
$\frac{v_{42}(r)}{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)$ &
$\frac{v_{43}(r)}{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}$  \\ \hline
%

\end{tabular}
\end{sidewaystable}
%
%
%

\subsection{Forces}
The force on object $\bf{a}$, $\mathbf{F}_{\bf a}$, due to object
$\bf{b}$ is the negative of the force 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}
%
Obtaining the force from the interaction energy expressions is the
same for higher-order multipole interactions -- the trick is to make
sure that all $r$-dependent derivatives are considered.  This is
straighforward if the interaction energies are written explicitly in
terms of $\hat{r}$ and the body axes ($\hat{a}_m$,
$\hat{b}_n$) :
%
\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}
%
Allen and Germano,\cite{Allen:2006fk} showed that if the energy is
written in this form, the forces come out relatively cleanly,
%
\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 that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $
is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk}
In simplifying the algebra, we have also used:
%
\begin{align}
\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{align}
%
We list below the force equations written in terms of lab-frame
coordinates.  The radial functions used in the two methods are listed
in Table \ref{tab:tableFORCE}
%
%SPACE COORDINATES FORCE EQUATIONS
%
% **************************************************************************
% f ca cb
%
\begin{align}
\mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =&
C_{\bf a} C_{\bf b}  w_a(r) \hat{r} \\
%
%
%
\mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =&
C_{\bf a} \Bigl[
 \left( \hat{r} \cdot \mathbf{D}_{\mathbf{b}} \right) 
w_b(r) \hat{r}  
+ \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr] \\
%
%
%
\mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =&
C_{\bf a } \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] \\
%
%
%
% \begin{equation}
% \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} = 
% -C_{\bf{b}} \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{split}
\mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} =&
- \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}
\end{split}\\
%
%
%
\begin{split}
\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =& - \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)
 - \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
& - \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
 -
(\hat{r} \cdot \mathbf{D}_{\mathbf{a}} )
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r} \end{split} \\
%
%
% \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{split}
\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =& 
 \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
&+ \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
&+  \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} \\
%
&+ \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
&+ 
(\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{align}
Note that the forces for higher multipoles on site $\mathbf{a}$
interacting with those of lower order on site $\mathbf{b}$ can be
obtained by swapping indices in the expressions above.

%
% Torques SECTION -----------------------------------------------------------------------------------------
%
\subsection{Torques}
When energies are written in the form of Eq.~({\ref{ugeneral}), then
  torques can be found in a relatively straightforward
  manner,\cite{Allen:2006fk}
%
\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}
%
%
The torques for both the Taylor-Shifted as well as Gradient-Shifted
methods are given in space-frame coordinates:
%
%
\begin{align}
\mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =& 
C_{\bf a}  (\hat{r} \times  \mathbf{D}_{\mathbf{b}}) v_{11}(r) \\
%
%
%
\mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =&
2C_{\bf a}
\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{22}(r) \\
%
%
%
% \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}
%
%
%
\mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} =&
 \mathbf{D}_{\mathbf {a}}  \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
% 2
-
(\hat{r} \times \mathbf{D}_{\mathbf {a}} )
(\hat{r} \cdot \mathbf{D}_{\mathbf {b}} )  v_{22}(r)\\
%
%
%
% \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}
%
%
%
\mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} =& 
 \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
- (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)\\
%
%
%
\mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =&
 \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
+
(\hat{r} \cdot \mathbf{D}_{\mathbf{a}})
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)\\
%
%
%
% \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{split}
\mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =&
-4 
\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}}
v_{41}(r) \\
% 2
&+ 
\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
&+ 2
\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}\\
%
%
%
\begin{split}
\mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} =  
&4
\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}} v_{41}(r) \\
% 2
&+  \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
&+2
(\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{align}
%
Here, we have defined the matrix cross product in an identical form
as in Ref. \onlinecite{Smith98}:
\begin{equation}
\left[\mathbf{A} \times \mathbf{B}\right]_\alpha = \sum_\beta
\left[\mathbf{A}_{\alpha+1,\beta} \mathbf{B}_{\alpha+2,\beta}
  -\mathbf{A}_{\alpha+2,\beta} \mathbf{B}_{\alpha+2,\beta} 
\right]
\end{equation}
where $\alpha+1$ and $\alpha+2$ are regarded as cyclic
permuations of the matrix indices.

All of the radial functions required for torques are identical with
the radial functions previously computed for the interaction energies.
These are tabulated for both shifted force methods in table
\ref{tab:tableenergy}.  The torques for higher multipoles on site
$\mathbf{a}$ interacting with those of lower order on site
$\mathbf{b}$ can be obtained by swapping indices in the expressions
above.

\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 energy
constants (converged to one part in $10^9$) in table \ref{tab:LT}.

\begin{table*}[h]
\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 \\
   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 \\
  --     &      BCC      &    minimum       &      -1.985920929 \\
\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 any of the arrays is given
by:
\begin{equation}
  E = C N^2 \mu^2
\end{equation}
where $C$ is the energy constant given in table \ref{tab:LT}, $N$ is
the number of dipoles per unit volume, and $\mu$ is the strength of
the dipole.

To test the new electrostatic methods, we have constructed very large
($N = 8,000, 16,000, 32,000$) arrays of dipoles in the orientations
described in \ref{tab:LT}. For the purposes of this paper, the primary
quantity of interest is the analytic energy constant for the perfect
arrays. Convergence to these constants are shown as a function of both
the cutoff radius, $r_c$, and the damping parameter, $\alpha$ in Fig.
XXX. We have simultaneously tested a hard cutoff (where the kernel is
simply truncated at the cutoff radius), as well as a shifted potential
(SP) form which includes a potential-shifting and self-interaction
term, but does not shift the forces and torques smoothly at the cutoff
radius.

{\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.

\section{Conclusion}
We have presented two efficient real-space methods for computing the
interactions between point multipoles.  These methods have the benefit
of smoothly truncating the energies, forces, and torques at the cutoff
radius, making them attractive for both molecular dynamics (MD) and
Monte Carlo (MC) simulations.  We find that the Gradient-Shifted Force
(GSF) and the Shifted-Potential (SP) methods converge rapidly to the
correct lattice energies for ordered dipolar and quadrupolar arrays,
while the Taylor-Shifted Force (TSF) is too severe an approximation to
provide accurate convergence to lattice energies.  

In most cases, GSF can obtain nearly quantitative agreement with the
lattice energy constants with reasonably small cutoff radii.  The only
exception we have observed is for crystals which exhibit a bulk
macroscopic dipole moment (e.g. Luttinger \& Tisza's $Z_1$ lattice).
In this particular case, the multipole neutralization scheme can
interfere with the correct computation of the energies.  We note that
the energies for these arrangements are typically much larger than for
crystals with net-zero moments, so this is not expected to be an issue
in most simulations.

In large systems, these new methods can be made to scale approximately
linearly with system size, and detailed comparisons with the Ewald sum
for a wide range of chemical environments follows in the second paper.

\begin{acknowledgments}
  JDG acknowledges helpful discussions with Christopher
  Fennell. 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}
\label{SmithFunc}
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.

\newpage

\bibliography{multipole}

\end{document}
%
% ****** End of file multipole.tex ******
Revision:	3986
Committed:	Tue Dec 31 17:28:43 2013 UTC (10 years, 5 months ago) by gezelter
Content type:	application/x-tex
File size:	64549 byte(s)
Log Message:	Migrated much to Supplemental, nearing completion on manuscript 1.