--- trunk/multipole/multipole1.tex	2014/01/02 15:46:58	3988
+++ trunk/multipole/multipole1.tex	2014/06/09 14:23:48	4177
@@ -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
@@ -152,24 +152,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 +181,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 +193,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
@@ -244,13 +244,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 +263,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,13 +276,20 @@ 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
@@ -306,6 +313,7 @@ In the standard multipole expansion, one typically use
 of $\bf a$ interacting with the same multipoles on $\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 +329,13 @@ 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)$ 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 +432,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}
 %
@@ -439,6 +453,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,19 +465,23 @@ 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(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r}
+\cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert  _{r_c} \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 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 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}
@@ -471,32 +494,20 @@ the radial factors differ between the two methods.
 the radial factors differ between the two methods.
 
 \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.
 
-[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
+The rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be
 expressed using these unit vectors:
 \begin{eqnarray}
 \hat{\mathbf {a}} = 
@@ -504,37 +515,84 @@ 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}_b $
+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}}
@@ -546,12 +604,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 +626,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,7 +641,7 @@ 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}
@@ -598,7 +661,7 @@ site which posesses a charge, dipole, and multipole, b
 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
+site which posesses a charge, dipole, and quadrupole, both types of
 contribution are given in table \ref{tab:tableSelf}.
 
 \begin{table*}
@@ -611,7 +674,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(
@@ -741,8 +804,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}}
@@ -775,8 +838,8 @@ along with the sign of the intersite vector, $\hat{r}$
   \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
+\begin{tabular}{|c|c|l|l|l|} \hline
+Generic&Bare Coulomb&Taylor-Shifted (TSF)&Shifted Potential (SP)&Gradient-Shifted (GSF)
 \\ \hline
 %
 %
@@ -785,7 +848,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 +858,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 +867,14 @@ $\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)
+$\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)$ & SP \\
+&&& $~~~-\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)$\\
 %
 %
 %
@@ -818,18 +882,15 @@ $\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)$ & SP \\
+&&& $-\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)$ 
-\\
+$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)$ &  SP \\
+&&& $~~~- \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)$ \\
 %
 %
 %
@@ -837,27 +898,29 @@ $\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)$ & SP \\
+&&& $~~~- \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)$ 
-\\
+$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r}
+\right)$ & SP \\
+&&& $~~~-\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
+$\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)$ & $~~~~~~~\left.+\frac{21s(r_c)}{r_c^2}
+-\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\
+\hline
 \end{tabular}
 \end{sidewaystable}
 %
@@ -981,49 +1044,8 @@ F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_
 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}
@@ -1134,8 +1156,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,29 +1193,8 @@ 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 both the Taylor-Shifted as well as Gradient-Shifted
 methods are given in space-frame coordinates:
 %
@@ -1355,18 +1356,40 @@ above.
 $\mathbf{b}$ can be obtained by swapping indices in the expressions
 above.
 
+\section{Related real-space methods}
+One can also formulate an extension of the Wolf approach for point
+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 total potential
+to zero at the cutoff radius,
+\begin{equation}
+U^{\text{SP}} = \sum \left[ U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
+U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) \right]
+\label{eq:SP}
+\end{equation}
+where the sum describes separate potential shifting that is applied to
+each orientational contribution to the energy.
+ 
+The energies and torques for the shifted potential (SP) can be easily
+obtained by zeroing out the $(r-r_c)$ terms in the final column of
+table \ref{tab:tableenergy}.  Forces for the SP method retain
+discontinuities at the cutoff sphere, and can be obtained by
+eliminating all functions that depend on $r_c$ in the last column of
+table \ref{tab:tableFORCE}.  The self-energy contributions for the SP
+potential are identical to both the GSF and TSF methods.
+
 \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.
+real-space cutoffs. In this paper we test both 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,117 +1399,134 @@ 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:
 \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.pdf}
+\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.pdf}
+\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
+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), 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.  The SP method is
+essentially an extension of the original Wolf method for multipoles.
+
+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
+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.
 
-
 {\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 the lowest energy configurations for these 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.pdf}
+\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}[!htbp]  
+\includegraphics[width=\linewidth]{Quadrupoles_newAlpha.pdf}
+\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) and
+gradient-shifted force (GSF) approximations converge to the correct
+energy smoothly by $r_c = 4 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.
+
+
 \section{Conclusion}
-We have presented two efficient real-space methods for computing the
+We have presented three 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
@@ -1494,7 +1534,7 @@ provide accurate convergence to lattice energies.  
 (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.  
+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
@@ -1513,7 +1553,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}
@@ -1609,33 +1649,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 +1688,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 \\
@@ -1670,8 +1711,8 @@ we generalize the notation of the previous appendix.  
 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)
+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}\\
@@ -1681,10 +1722,9 @@ For undamped charges, $f(r)=1/r$, Table I lists these 
 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.
+For undamped charges 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
 
