--- trunk/multipole/multipole1.tex	2014/01/02 15:46:58	3988
+++ trunk/multipole/multipole1.tex	2014/07/17 18:24:41	4195
@@ -44,7 +44,7 @@ Sum. I. Taylor-shifted and Gradient-shifted electrosta
 %\preprint{AIP/123-QED}
 
 \title{Real space alternatives to the Ewald
-Sum. I. Taylor-shifted and Gradient-shifted electrostatics for multipoles}
+Sum. I. Shifted electrostatics for multipoles}
 
 \author{Madan Lamichhane}
 \affiliation{Department of Physics, University
@@ -65,15 +65,25 @@ of Notre Dame, Notre Dame, IN 46556}
 
 \begin{abstract}
   We have extended the original damped-shifted force (DSF)
-  electrostatic kernel and have been able to derive two new
+  electrostatic kernel and have been able to derive three 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.
+  on truncated Taylor expansions around the cutoff radius.  These
+  include a shifted potential (SP) that generalizes the Wolf method
+  for point multipoles, and Taylor-shifted force (TSF) and
+  gradient-shifted force (GSF) potentials that are both
+  generalizations of DSF electrostatics for multipoles. We find that
+  each of the distinct orientational contributions requires a separate
+  radial function to ensure that pairwise energies, forces and torques
+  all 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.  We find that the GSF and 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.
+  Because real-space methods can be made to scale linearly with system
+  size, SP and GSF are attractive options for large Monte
+  Carlo and molecular dynamics simulations, respectively.
 \end{abstract}
 
 %\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
@@ -133,13 +143,14 @@ kernel to point multipoles.  We describe two distinct 
 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
+kernel to point multipoles.  We describe three distinct real-space
+interaction models for higher-order multipoles based on truncated
 Taylor expansions that are carried out at the cutoff radius.  We are
-calling these models {\bf Taylor-shifted} and {\bf Gradient-shifted}
+calling these models {\bf Taylor-shifted} (TSF), {\bf
+  gradient-shifted} (GSF) and {\bf shifted potential} (SP)
 electrostatics.  Because of differences in the initial assumptions,
-the two methods yield related, but somewhat different expressions for
-energies, forces, and torques.
+the three methods yield related, but distinct 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
@@ -152,24 +163,25 @@ V = \sum_i \sum_{j>i} V_\mathrm{pair}(r_{ij}, \Omega_i
 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}
+U = \sum_i \sum_{j>i} U_\mathrm{pair}(\mathbf{r}_{ij}, \Omega_i, \Omega_j) +
+\sum_i U_i^\mathrm{self}
 \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\{
+  U_\mathrm{pair}(\mathbf{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 ,
+U_\mathrm{approx} (\mathbf{r}_{ij}, \Omega_i, \Omega_j) & \quad \left| \mathbf{r}_{ij} \right| \le r_c \\
+0 & \quad \left| \mathbf{r}_{ij} \right|  > 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
+along with an easily computed self-interaction term ($\sum_i
+U_i^\mathrm{self}$) which scales linearly with the number of
 particles.  Here $\Omega_i$ and $\Omega_j$ represent orientational
-coordinates of the two sites.  The computational efficiency, energy
+coordinates of the two sites, and $\mathbf{r}_{ij}$ is the vector
+between 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
+depend dramatically on how the $U_\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
@@ -180,12 +192,11 @@ computed within $r_c$. Damping using a complementary e
 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
+computed within $r_c$. Damping using a complementary error function is
+applied to the potential to accelerate convergence. The interaction
+for a pair of charges ($C_i$ and $C_j$) in the DSF method,
 \begin{equation*}
-V_\mathrm{DSF}(r) = C_i C_j \Biggr{[}
+U_\mathrm{DSF}(r_{ij}) = 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}}
@@ -193,9 +204,9 @@ Note that in this potential and in all electrostatic q
 \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.
+where $\alpha$ is the adjustable damping parameter.  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
@@ -222,7 +233,7 @@ forces and torques go smoothly to zero at the cutoff d
 forces and torques go smoothly to zero at the cutoff distance.
 
 \begin{figure}
-\includegraphics[width=3in]{SM}
+\includegraphics[width=3in]{SM.eps}
 \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.}
@@ -244,13 +255,13 @@ V_a(\mathbf r) = 
 a$.  Then the electrostatic potential of object $\bf a$ at
 $\mathbf{r}$ is given by
 \begin{equation}
-V_a(\mathbf r) = 
+\phi_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
+labeling specific charges in $\bf a$ and $\bf b$ respectively.  The
 Taylor expansion,
 \begin{equation}
  \frac{1}{\lvert \mathbf{r} - \mathbf{r}_k \rvert} = 
@@ -263,7 +274,7 @@ V_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r} 
 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} 
+\phi_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r} 
 \end{equation}
 where
 \begin{equation}
@@ -276,16 +287,23 @@ C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \\
     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} .
+C_{\bf a} =&\sum_{k \, \text{in \bf a}} q_k , \label{eq:charge} \\
+D_{{\bf a}\alpha} =&\sum_{k \, \text{in \bf a}} q_k r_{k\alpha}, \label{eq:dipole}\\
+Q_{{\bf a}\alpha\beta} =& \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k
+r_{k\alpha}  r_{k\beta} . \label{eq:quadrupole}
 \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$. 
+based factor of $1/2$.  We are essentially treating the mass
+distribution with higher priority; the moment of inertia tensor,
+$\overleftrightarrow{\mathsf I}$, is diagonalized to obtain body-fixed
+axes, and the charge distribution may result in a quadrupole tensor
+that is not necessarily diagonal in the body frame.  Additional
+reasons for utilizing the primitive quadrupole are discussed in
+section \ref{sec:damped}.
 
 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
+$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $), 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} .
@@ -303,9 +321,10 @@ of $\bf a$ interacting with the same multipoles on $\b
 \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$.
+of $\bf a$ interacting with the same types of multipoles in $\bf b$.
 
 \subsection{Damped Coulomb interactions}
+\label{sec:damped}
 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
@@ -321,7 +340,14 @@ functions $B_l(r)$ are summarized in Appendix A.
 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.
+functions $B_l(r)$, which are used for derivatives of the damped
+kernel, are summarized in Appendix A.  (N.B. there is one important
+distinction between the two kernels, which is the behavior of
+$\nabla^2 \frac{1}{r}$ compared with $\nabla^2 B_0(r)$.  The former is
+zero everywhere except for a delta function evaluated at the origin.
+The latter also has delta function behavior, but is non-zero for $r
+\neq 0$.  Thus the standard justification for using a traceless
+quadrupole tensor fails for the damped case.)
 
 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
@@ -418,7 +444,7 @@ U= (\text{prefactor}) (\text{derivatives}) f_n(r)
 In general, we can write
 %
 \begin{equation}
-U= (\text{prefactor}) (\text{derivatives}) f_n(r)
+U^{\text{TSF}}= (\text{prefactor}) (\text{derivatives}) f_n(r)
 \label{generic}
 \end{equation}
 %
@@ -429,7 +455,8 @@ space indices.
 $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.
+space indices.  Appendix \ref{radialTSF} contains details on the
+radial functions.
 
 In the formulas presented in the tables below, the placeholder
 function $f(r)$ is used to represent the electrostatic kernel (either
@@ -439,6 +466,11 @@ presented in tables \ref{tab:tableenergy} and \ref{tab
 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}.
+One of the principal findings of our work is that the individual
+orientational contributions to the various multipole-multipole
+interactions must be treated with distinct radial functions, but each
+of these contributions is independently force shifted at the cutoff
+radius.  
 
 \subsection{Gradient-shifted force (GSF) electrostatics}
 The second, and conceptually simpler approach to force-shifting
@@ -446,57 +478,75 @@ U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \n
 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} .
+U^{\text{GSF}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
+U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c)
+\hat{\mathbf{r}} \cdot \nabla U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
 \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:
+where $\hat{\mathbf{r}}$ is the unit vector pointing between the two
+multipoles, and the sum describes a separate force-shifting that is
+applied to each orientational contribution to the energy.  Both the
+potential and the gradient for force shifting are evaluated for an
+image multipole projected onto the surface of the cutoff sphere (see
+fig \ref{fig:shiftedMultipoles}).  The image multipole retains the
+orientation ($\hat{\mathbf{b}}$) of the interacting multipole.  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 interaction 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
+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}
+\subsection{Generalization of the Wolf shifted potential (SP)}
+It is also possible to formulate an extension of the Wolf approach for
+multipoles by simply projecting the image multipole onto the surface
+of the cutoff sphere, and including the interactions with the central
+multipole and the image.  This effectively shifts the pair potential
+to zero at the cutoff radius,
+\begin{equation}
+U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
+U(r_c \hat{\mathbf{r}},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
+\label{eq:SP}
+\end{equation}
+independent of the orientations of the two multipoles.  The sum again
+describes separate potential shifting that is applied to each
+orientational contribution to the energy.
+ 
+The shifted potential (SP) method is a simple truncation of the GSF
+method for each orientational contribution, leaving out the $(r-r_c)$
+terms that multiply the gradient. Functional forms for the
+shifted-potential (SP) method can also be summarized using the form of
+Eq.~\ref{generic3}.  The energy, force, and torque expressions are
+tabulated below for all three methods. As in the GSF and TSF methods,
+for each separate orientational contribution, only the radial factors
+differ between the SP, GSF, and TSF methods.
 
-[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}
+\subsection{\label{sec:level2}Body and space axes}
+Although objects $\bf a$ and $\bf b$ rotate during a molecular
+dynamics (MD) simulation, their multipole tensors remain fixed in
+body-frame coordinates. While deriving force and torque expressions,
+it is therefore convenient to write the energies, forces, and torques
+in intermediate forms involving the vectors of the rotation matrices.
+We denote body axes for objects $\bf a$ and $\bf b$ using unit vectors
+$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$.
+In a typical simulation, the initial axes are obtained by
+diagonalizing the moment of inertia tensors for the objects.  (N.B.,
+the body axes are generally {\it not} the same as those for which the
+quadrupole moment is diagonal.)  The rotation matrices are then
+propagated during the simulation.
 
-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
+The rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be
 expressed using these unit vectors:
 \begin{eqnarray}
 \hat{\mathbf {a}} = 
@@ -504,40 +554,87 @@ expressed using these unit vectors:
 \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}\\
+\end{pmatrix}, \qquad
 \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:
+coordinates.
+
+Allen and Germano,\cite{Allen:2006fk} following earlier work by Price
+{\em et al.},\cite{Price:1984fk} showed that 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}
+%
+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}
+
+The torques can also be found in a relatively similar
+manner,
+%
+\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}
+
+Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $
+is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk}
+We also made use of the identities,
+%
+\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}
+
+Many of the multipole contractions required can be written in one of
+three equivalent forms using the 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 all three 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}}
+\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
@@ -546,12 +643,17 @@ contractions using space indices. 
 explicit sums over body indices and the dot products now indicate
 contractions using space indices. 
 
+In computing our force and torque expressions, we carried out most of
+the work in body coordinates, and have transformed the expressions
+back to space-frame coordinates, which are reported below.  Interested
+readers may consult the supplemental information for this paper for
+the intermediate body-frame expressions.
 
 \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
+extension~\cite{Fennell:2006zl} also include 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
@@ -563,7 +665,7 @@ V_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C
 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}
+U_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2}.
 \label{eq:selfTerm}
 \end{equation}
 
@@ -578,14 +680,14 @@ V_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r
 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} + 
+U_\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
+treatment of the self-neutralization \textit{and} reciprocal
 contributions to the self-interaction for higher order multipoles.  In
 this section we give formulae for these interactions up to quadrupolar
 order.
@@ -596,10 +698,10 @@ obtained by Smith and Aguado and Madden for the applic
 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}.
+obtained by Smith, and also by 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 quadrupole, both types
+of contribution are given in table \ref{tab:tableSelf}.
 
 \begin{table*}
   \caption{\label{tab:tableSelf} Self-interaction contributions for 
@@ -611,7 +713,7 @@ Quadrupole & $2 \text{Tr}(\mathbf{Q}_{\bf a}^2) + \tex
 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$ &
+Quadrupole & $2 \mathbf{Q}_{\bf a}:\mathbf{Q}_{\bf a} + \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(
@@ -637,16 +739,16 @@ work for both the Taylor-shifted and Gradient-shifted 
 \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,
+work for the Taylor-shifted, gradient-shifted, and shifted potential
+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 the
+three 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 ----------------------------------------------------------------------------------------------
 %
@@ -741,8 +843,8 @@ +2 \text{Tr} \left(
 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
+\mathbf{Q}_{\mathbf{a}} : \mathbf{Q}_{\mathbf{b}} \Bigr] v_{41}(r)
 \\
 % 2
 &+\Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
@@ -773,10 +875,14 @@ along with the sign of the intersite vector, $\hat{r}$
 
 \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
+    and torque equations.  The $f, g, h, s, t, \mathrm{and~} u$
+    functions used in this table are defined in Appendices
+    \ref{radialTSF} and \ref{radialGSF}.  The gradient shifted (GSF)
+    functions include the shifted potential (SP)
+    contributions (\textit{cf.} Eqs. \ref{generic2} and
+    \ref{eq:SP}).}
+\begin{tabular}{|c|c|l|l|l|} \hline
+Generic&Bare Coulomb&Taylor-Shifted (TSF)&Shifted Potential (SP)&Gradient-Shifted (GSF)
 \\ \hline
 %
 %
@@ -785,7 +891,8 @@ $f(r)-f(r_c)-(r-r_c)g(r_c)$
 $v_{01}(r)$ &
 $\frac{1}{r}$ &
 $f_0(r)$ &
-$f(r)-f(r_c)-(r-r_c)g(r_c)$
+$f(r)-f(r_c)$ &
+SP $-(r-r_c)g(r_c)$
 \\
 %
 %
@@ -794,7 +901,8 @@ $g(r)-g(r_c)-(r-r_c)h(r_c)$ \\
 $v_{11}(r)$ &
 $-\frac{1}{r^2}$ &
 $g_1(r)$ &
-$g(r)-g(r_c)-(r-r_c)h(r_c)$ \\
+$g(r)-g(r_c)$ &
+SP $-(r-r_c)h(r_c)$ \\
 %
 %
 %
@@ -802,15 +910,16 @@ $\frac{g(r)}{r}-\frac{g(r_c)}{r_c} -(r-r_c)
 $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)$ \\
+$\frac{g(r)}{r}-\frac{g(r_c)}{r_c}$ &
+SP $-(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)$ 
-\\
+$\left(-\frac{g(r)}{r}+h(r) \right) -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right)$ 
+& SP $-(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$\\
 %
 %
 %
@@ -818,18 +927,18 @@ $\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r} \right)
 $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)$ 
-\\
+$\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)$ 
+& SP $-(r-r_c) \left(\frac{2g(r_c)}{r_c^3}-\frac{2h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$  \\
 % 
+%
+%
 $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)$ 
-\\
+$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)$&
+SP $-(r-r_c) \left( \frac{-6g(r_c)}{r_c^3}+\frac{6h(r_c)}{r_c^2}\right.$ \\
+&&& $~~~-\left(\frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c)\right)$ & 
+$\phantom{SP-(r-r_c)}\left.-\frac{3s(r_c)}{r_c}+t(r_c) \right)$\\
 %
 %
 %
@@ -837,27 +946,28 @@ $\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right)
 $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)$ 
+$\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)$ & 
+SP $-(r-r_c) \left( \frac{3g(r_c)}{r_c^4}-\frac{3h(r_c)}{r_c^3}+\frac{s(r_c)}{r_c^2} \right)$ 
 \\
 % 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)$ 
-\\
+$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right)$ & 
+SP$-(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}\right.$  \\
+&&& $~~~-\left( \frac{3g(r_c)}{r_c^3} -  \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ & 
+$\phantom{SP-(r-r_c)}\left. -\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
+$ \left(-\frac{15g(r)}{r^3} +\frac{15h(r)}{r^2}-\frac{6s(r)}{r}+t(r)\right) $ &
+SP $-(r-r_c)\left(\frac{45g(r_c)}{r_c^4}-\frac{45h(r_c)}{r_c^3}\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)$ & 
+$\phantom{SP-(r-r_c)}\left.+\frac{21s(r_c)}{r_c^2}-\frac{6t(r_c)}{r_c}+u(r_c) \right)$\\
+\hline
 \end{tabular}
 \end{sidewaystable}
 %
@@ -866,9 +976,13 @@ -\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\ \hline
 %
 
 \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
+\caption{\label{tab:tableFORCE}Radial functions used in the force
+  equations. Gradient shifted (GSF) functions are constructed using the shifted
+    potential (SP) functions.  Some of these functions are simple
+    modifications of the functions found in table \ref{tab:tableenergy}}
+\begin{tabular}{|c|c|l|l|l|} \hline
+Function&Definition&Taylor-Shifted (TSF)& Shifted Potential (SP)
+&Gradient-Shifted (GSF)
 \\ \hline
 %
 %
@@ -876,91 +990,100 @@ $g(r)-g(r_c)$ \\
 $w_a(r)$&
 $\frac{d v_{01}}{dr}$&
 $g_0(r)$&
-$g(r)-g(r_c)$ \\
+$g(r)$&
+SP $-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} $ \\
+$h(r) - \frac{v_{11}(r)}{r} $ &
+SP $- h(r_c)$ \\
 %
 $w_c(r)$ &
 $\frac{v_{11}(r)}{r}$ &
 $\frac{g_1(r)}{r} $ &
+$\frac{v_{11}(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) $ \\
+$\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)$ &
+SP $-\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}$ &
 $\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}$\\
+  $ \left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) -\frac{2v_{22}(r)}{r}$&
+SP $- \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$\\
 %
 $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}$&
 $\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}$ \\
+$ \left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) -\frac{v_{31}(r)}{r}$ &
+SP $ - \left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\
 % 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}$&
 $\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}$ \\
+$\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right) -\frac{3v_{32}}{r}$ &
+SP $-\left(\frac{-6g(_cr)}{r_c^3} +\frac{6h(r_c)}{r_c^2}
+  -\frac{3s(r_c)}{r_c} +t(r_c) \right)$ \\
 %
 $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)$ \\
+$\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}
+\right)$ &
+SP $-\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}$\\
+$\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}$&
+SP$-\left(-\frac{9g(r_c)}{r_c^4} +\frac{9h(r_c)}{r_c^3} -\frac{4s(r_c)}{r_c^2} +\frac{t(r_c)}{r_c} \right)$\\
 %
 $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}$  \\
+$\left(\frac{105g_4(r)}{r^4} - \frac{105h_4(r)}{r^3} \right.$ &
+$\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2}\right.$ &
+SP $- \left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}\right.$ \\
+&& $~~~\left.+ \frac{45s_4(r)}{r^2} - \frac{10t_4(r)}{r} +u_4(r) \right)$
+& $~~~\left. -\frac{6t(r)}{r} +u(r) \right) -\frac{4v_{43}(r)}{r}$ &
+$\phantom{SP-} \left.+\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $\\ 
 %
 $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}$&
 $\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}$&
 $\frac{v_{43}(r)}{r}$  \\ \hline
 %
 
@@ -982,50 +1105,9 @@ F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_
 \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
+coordinates.  The radial functions used in the three methods are listed
 in Table \ref{tab:tableFORCE}
 %
 %SPACE COORDINATES FORCE EQUATIONS
@@ -1134,8 +1216,8 @@ w_i(r)
 \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} \\
+\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} 
++ 2  \mathbf{Q}_{\mathbf{a}} :  \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}} )  
@@ -1171,33 +1253,12 @@ When energies are written in the form of Eq.~({\ref{ug
 % 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 the three methods are given in space-frame
+coordinates:
 %
 %
-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) \\
@@ -1343,13 +1404,14 @@ as in Ref. \onlinecite{Smith98}:
 \left[\mathbf{A}_{\alpha+1,\beta} \mathbf{B}_{\alpha+2,\beta}
   -\mathbf{A}_{\alpha+2,\beta} \mathbf{B}_{\alpha+2,\beta} 
 \right]
+\label{eq:matrixCross}
 \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
+These are tabulated for all three 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
@@ -1361,12 +1423,12 @@ real-space cutoffs. In this paper we test Taylor-shift
 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.
+real-space cutoffs. In this paper we test SP, TSF, and GSF
+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
@@ -1376,125 +1438,160 @@ and other periodic structures.  We have repeated the L
 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}.
+and other periodic structures.  
 
-\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
+To test the new electrostatic methods, we have constructed very large,
+$N=$ 16,000~(bcc) arrays of dipoles in the orientations described in
+Ref. \onlinecite{LT}.  These structures include ``A'' lattices with
+nearest neighbor chains of antiparallel dipoles, as well as ``B''
+lattices with nearest neighbor strings of antiparallel dipoles if the
+dipoles are contained in a plane perpendicular to the dipole direction
+that passes through the dipole.  We have also studied the 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:
+Tisza.  The total electrostatic energy density 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.
+where $C$ is the energy constant (equivalent to the Madelung
+constant), $N$ is the number of dipoles per unit volume, and $\mu$ is
+the strength of the dipole. Energy constants (converged to 1 part in
+$10^9$) are given in the supplemental information.
 
-To test the new electrostatic methods, we have constructed very large,
-$N$ = 8,000~(sc), 16,000~(bcc), or 32,000~(fcc) arrays of dipoles in
-the orientations described in table \ref{tab:LT}.  For the purposes of
-testing the energy expressions and the self-neutralization schemes,
-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 Figs.  \ref{fig:energyConstVsCutoff} and 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.
+\begin{figure}
+\includegraphics[width=\linewidth]{Dipoles_rCutNew.eps}
+\caption{Convergence of the lattice energy constants as a function of
+  cutoff radius (normalized by the lattice constant, $a$) for the new
+  real-space methods.  Three dipolar crystal structures were sampled,
+  and the analytic energy constants for the three lattices are
+  indicated with grey dashed lines.  The left panel shows results for
+  the undamped kernel ($1/r$), while the damped error function kernel,
+  $B_0(r)$ was used in the right panel.}
+\label{fig:Dipoles_rCut}
+\end{figure}
 
 \begin{figure}
-\includegraphics[width=4.5in]{energyConstVsCutoff}
-\caption{Convergence to the analytic energy constants as a function of
-  cutoff radius (normalized by the lattice constant) for the different
-  real-space methods. The two crystals shown here are the ``B'' array
-  for bcc crystals with the dipoles along the 001 direction (upper),
-  as well as the minimum energy bcc lattice (lower).  The analytic
-  energy constants are shown as a grey dashed line.  The left panel
-  shows results for the undamped kernel ($1/r$), while the damped
-  error function kernel, $B_0(r)$ was used in the right panel. }
-\label{fig:energyConstVsCutoff}
+\includegraphics[width=\linewidth]{Dipoles_alphaNew.eps}
+\caption{Convergence to the lattice energy constants as a function of
+  the reduced damping parameter ($\alpha^* = \alpha a$) for the
+  different real-space methods in the same three dipolar crystals in
+  Figure \ref{fig:Dipoles_rCut}.  The left panel shows results for a
+  relatively small cutoff radius ($r_c = 4.5 a$) while a larger cutoff
+  radius ($r_c = 6 a$) was used in the right panel. }
+\label{fig:Dipoles_alpha}
 \end{figure}
 
-The Hard cutoff exhibits oscillations around the analytic energy
-constants, and converges to incorrect energies when the complementary
-error function damping kernel is used.  The shifted potential (SP) and
-gradient-shifted force (GSF) approximations converge to the correct
-energy smoothly by $r_c / 6 a$ even for the undamped case.  This
-indicates that the correction provided by the self term is required
-for obtaining accurate energies.  The Taylor-shifted force (TSF)
-approximation appears to perturb the potential too much inside the
-cutoff region to provide accurate measures of the energy constants.
+For the purposes of testing the energy expressions and the
+self-neutralization schemes, 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 Figs.\ref{fig:Dipoles_rCut}
+and \ref{fig:Dipoles_alpha}. We have simultaneously tested a hard
+cutoff (where the kernel is simply truncated at the cutoff radius) in
+addition to the three new methods.
 
+The hard cutoff exhibits oscillations around the analytic energy
+constants, and converges to incorrect energies when the complementary
+error function damping kernel is used.  The shifted potential (SP)
+converges to the correct energy smoothly by $r_c = 4.5 a$ even for the
+undamped case.  This indicates that the shifting and the correction
+provided by the self term are  required for obtaining accurate energies.
+The Taylor-shifted force (TSF) approximation appears to perturb the
+potential too much inside the cutoff region to provide accurate
+measures of the energy constants.  GSF is a compromise, converging to
+the correct energies within $r_c = 6 a$. 
 
 {\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}
+approach.\cite{Nagai01081960,Nagai01091963} 
 
-\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
+ E = C N \frac{3\bar{Q}^2}{4a^5} 
 \end{equation} 
-where $Q$ is the quadrupole moment.
+where $a$ is the lattice parameter, and $\bar{Q}$ is the effective
+quadrupole moment,
+\begin{equation}
+\bar{Q}^2 = 2 \left(3 Q : Q - (\text{Tr} Q)^2 \right)
+\end{equation}
+for the primitive quadrupole as defined in Eq. \ref{eq:quadrupole}.
+(For the traceless quadrupole tensor, $\Theta = 3 Q - \text{Tr} Q$,
+the effective moment, $\bar{Q}^2 = \frac{2}{3} \Theta : \Theta$.)
 
+To test the new electrostatic methods for quadrupoles, we have
+constructed very large, $N=$ 8,000~(sc), 16,000~(bcc), and
+32,000~(fcc) arrays of linear quadrupoles in the orientations
+described in Ref. \onlinecite{Nagai01081960}.  We have compared the
+energy constants for these low-energy configurations for linear
+quadrupoles. Convergence to these constants are shown as a function of
+both the cutoff radius, $r_c$, and the damping parameter, $\alpha$ in
+Figs.~\ref{fig:Quadrupoles_rCut} and \ref{fig:Quadrupoles_alpha}.
+
+\begin{figure}
+\includegraphics[width=\linewidth]{Quadrupoles_rcutConvergence.eps}
+\caption{Convergence of the lattice energy constants as a function of
+  cutoff radius (normalized by the lattice constant, $a$) for the new
+  real-space methods.  Three quadrupolar crystal structures were
+  sampled, and the analytic energy constants for the three lattices
+  are indicated with grey dashed lines.  The left panel shows results
+  for the undamped kernel ($1/r$), while the damped error function
+  kernel, $B_0(r)$ was used in the right panel.}
+\label{fig:Quadrupoles_rCut}
+\end{figure}
+
+
+\begin{figure}
+\includegraphics[width=\linewidth]{Quadrupoles_newAlpha.eps}
+\caption{Convergence to the lattice energy constants as a function of
+  the reduced damping parameter ($\alpha^* = \alpha a$) for the
+  different real-space methods in the same three quadrupolar crystals
+  in Figure \ref{fig:Quadrupoles_rCut}.  The left panel shows
+  results for a relatively small cutoff radius ($r_c = 4.5 a$) while a
+  larger cutoff radius ($r_c = 6 a$) was used in the right panel.  }
+\label{fig:Quadrupoles_alpha}
+\end{figure}
+
+Again, we find that the hard cutoff exhibits oscillations around the
+analytic energy constants.  The shifted potential (SP) approximation
+converges to the correct energy smoothly by $r_c = 3 a$ even for the
+undamped case.  The Taylor-shifted force (TSF) approximation again
+appears to perturb the potential too much inside the cutoff region to
+provide accurate measures of the energy constants.  GSF again provides
+a compromise between the two methods -- energies are converged by $r_c
+= 4.5 a$, and the approximation is not as perturbative at short range
+as TSF.
+
+It is also useful to understand the convergence to the lattice energy
+constants as a function of the reduced damping parameter ($\alpha^* =
+\alpha a$) for the different real-space methods.
+Figures \ref{fig:Dipoles_alpha} and \ref{fig:Quadrupoles_alpha} show
+this comparison for the dipolar and quadrupolar lattices,
+respectively.  All of the methods (except for TSF) have excellent
+behavior for the undamped or weakly-damped cases.  All of the methods
+can be forced to converge by increasing the value of $\alpha$, which
+effectively shortens the overall range of the potential by equalizing
+the truncation effects on the different orientational contributions.
+In the second paper in the series, we discuss how large values of
+$\alpha$ can perturb the force and torque vectors, but both
+weakly-damped or over-damped electrostatics appear to generate
+reasonable values for the total electrostatic energies under both the
+SP and GSF approximations.
+
 \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.  
+We have presented three efficient real-space methods for computing the
+interactions between point multipoles.  One of these (SP) is a
+multipolar generalization of Wolf's method that smoothly shifts
+electrostatic energies to zero at the cutoff radius.  Two of these
+methods (GSF and TSF) also smoothly truncate the forces and torques
+(in addition to the energies) at the cutoff radius, making them
+attractive for both molecular dynamics and Monte Carlo 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
@@ -1506,6 +1603,14 @@ In large systems, these new methods can be made to sca
 crystals with net-zero moments, so this is not expected to be an issue
 in most simulations.
 
+The techniques used here to derive the force, torque and energy
+expressions can be extended to higher order multipoles, although some
+of the objects (e.g. the matrix cross product in
+Eq. \ref{eq:matrixCross}) will need to be generalized for higher-rank
+tensors.  We also note that the definitions of the multipoles used
+here are in a primitive form, and these need some care when comparing
+with experiment or other computational techniques.
+
 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.
@@ -1513,7 +1618,7 @@ for a wide range of chemical environments follows in t
 \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
+  Science Foundation under grant CHE-1362211. Computational time was
   provided by the Center for Research Computing (CRC) at the
   University of Notre Dame.
 \end{acknowledgments}
@@ -1552,26 +1657,16 @@ is very useful for building up higher derivatives.  Us
 \text{e}^{-{\alpha}^2r^2}
 \right] ,
 \end{equation}
-is very useful for building up higher derivatives.  Using these
-formulas, we find:
+is very useful for building up higher derivatives. As noted by Smith,
+it is possible to approximate the $B_l(r)$ functions,
 %
-\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}
+\label{radialTSF}
 
 Using the shifted damped functions $f_n(r)$ defined by:
 %
@@ -1609,33 +1704,33 @@ u_4(r)=B_0^{(5)}(r) - B_0^{(5)}(r_c) .
 \begin{equation}
 u_4(r)=B_0^{(5)}(r) - B_0^{(5)}(r_c) .
 \end{equation}
-
+% 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}
 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:
+stored on a grid for values of $r$ from $0$ to $r_c$.  Using these
+functions, we find
 %
-\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} =&
+\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta \partial 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{g_n}{r^3} +\frac{h_n}{r^2} \right) \nonumber \\
+& + 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} =& 
+\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta \partial
+  r_\gamma \partial r_\delta} =& 
 \left( \delta_{\alpha \beta} \delta_{\gamma \delta} 
 + \delta_{\alpha \gamma} \delta_{\beta \delta}
  +\delta_{ \beta \gamma} \delta_{\alpha \delta} \right)
@@ -1648,7 +1743,8 @@ + \frac{t_n}{r^4} \right)\\
 \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} =& 
+{\partial r_\alpha \partial r_\beta \partial r_\gamma \partial
+  r_\delta \partial 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 \\
@@ -1665,26 +1761,30 @@ - \frac{10t_n}{r^6} +\frac{u_n}{r^5} \right) \label{eq
 %
 \newpage
 \section{The $r$-dependent factors for GSF electrostatics}
+\label{radialGSF}
 
 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)
+and the expansion is carried out on the individual terms in the
+multipole interaction energies. For damped charges, this still brings
+multiple derivatives of the Smith's $B_0(r)$ function into the
+algebra. To denote these terms, we generalize the notation of the
+previous appendix. For either $f(r)=1/r$ (undamped) or $f(r)=B_0(r)$
+(damped),
 %
 \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} .
+g(r) &= \frac{df}{d r} &&                      &&=-\frac{1}{r^2}
+&&\mathrm{or~~~} -rB_1(r) \\
+h(r) &= \frac{dg}{d r} &&= \frac{d^2f}{d r^2} &&= \frac{2}{r^3} &&\mathrm{or~~~}-B_1(r) + r^2 B_2(r) \\
+s(r) &= \frac{dh}{d r} &&= \frac{d^3f}{d r^3} &&=-\frac{6}{r^4}&&\mathrm{or~~~}3rB_2(r) - r^3 B_3(r)\\
+t(r) &= \frac{ds}{d r} &&= \frac{d^4f}{d r^4} &&= \frac{24}{r^5} &&\mathrm{or~~~} 3
+B_2(r) - 6r^2 B_3(r) + r^4 B_4(r) \\
+u(r) &= \frac{dt}{d r} &&= \frac{d^5f}{d r^5} &&=-\frac{120}{r^6} &&\mathrm{or~~~} -15
+r B_3(r) + 10 r^3B_4(r) -r^5B_5(r).
 \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.
+For undamped charges, Table I lists these derivatives under the Bare
+Coulomb column. Equations \ref{eq:b9} to \ref{eq:b13} are still
+correct for GSF electrostatics if the subscript $n$ is eliminated.
 
 \newpage