ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/multipole/multipole1.tex
(Generate patch)

Comparing trunk/multipole/multipole1.tex (file contents):
Revision 3983 by gezelter, Fri Dec 27 18:09:59 2013 UTC vs.
Revision 3990 by gezelter, Fri Jan 3 18:28:01 2014 UTC

# Line 26 | Line 26 | preprint,%
26   % reprint,%
27   %author-year,%
28   %author-numerical,%
29 < ]{revtex4-1}
29 > jcp]{revtex4-1}
30  
31   \usepackage{graphicx}% Include figure files
32   \usepackage{dcolumn}% Align table columns on decimal point
# Line 34 | Line 34 | preprint,%
34   \usepackage{times}
35   \usepackage[version=3]{mhchem}  % this is a great package for formatting chemical reactions
36   \usepackage{url}
37 + \usepackage{rotating}
38  
39   %\usepackage[mathlines]{lineno}% Enable numbering of text and display math
40   %\linenumbers\relax % Commence numbering lines
41  
42   \begin{document}
43  
44 < \preprint{AIP/123-QED}
44 > %\preprint{AIP/123-QED}
45  
46 < \title[Taylor-shifted and Gradient-shifted electrostatics for multipoles]
46 < {Real space alternatives to the Ewald
46 > \title{Real space alternatives to the Ewald
47   Sum. I. Taylor-shifted and Gradient-shifted electrostatics for multipoles}
48  
49   \author{Madan Lamichhane}
# Line 76 | Line 76 | of Notre Dame, Notre Dame, IN 46556}
76    results for ordered arrays of multipoles.
77   \end{abstract}
78  
79 < \pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
79 > %\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
80                               % Classification Scheme.
81 < \keywords{Suggested keywords}%Use showkeys class option if keyword
81 > %\keywords{Suggested keywords}%Use showkeys class option if keyword
82                                %display desired
83   \maketitle
84  
# Line 98 | Line 98 | similar methods where the force vanishes at $R_\textrm
98   The computational efficiency and the accuracy of the DSF method are
99   surprisingly good, particularly for systems with uniform charge
100   density. Additionally, dielectric constants obtained using DSF and
101 < similar methods where the force vanishes at $R_\textrm{c}$ are
101 > similar methods where the force vanishes at $r_{c}$ are
102   essentially quantitative.\cite{Izvekov:2008wo} The DSF and other
103   related methods have now been widely investigated,\cite{Hansen:2012uq}
104 < and DSF is now used routinely in simulations of ionic
105 < liquids,\cite{doi:10.1021/la400226g,McCann:2013fk} flow in carbon
106 < nanotubes,\cite{kannam:094701} gas sorption in metal-organic framework
107 < materials,\cite{Forrest:2012ly} thermal conductivity of methane
108 < hydrates,\cite{English:2008kx} condensation coefficients of
109 < water,\cite{Louden:2013ve} and in tribology at solid-liquid-solid
110 < interfaces.\cite{Tokumasu:2013zr} DSF electrostatics provides a
111 < compromise between the computational speed of real-space cutoffs and
112 < the accuracy of fully-periodic Ewald treatments.
104 > and DSF is now used routinely in a diverse set of chemical
105 > environments.\cite{doi:10.1021/la400226g,McCann:2013fk,kannam:094701,Forrest:2012ly,English:2008kx,Louden:2013ve,Tokumasu:2013zr}
106 > DSF electrostatics provides a compromise between the computational
107 > speed of real-space cutoffs and the accuracy of fully-periodic Ewald
108 > treatments.
109  
114 \subsection{Coarse Graining using Point Multipoles}
110   One common feature of many coarse-graining approaches, which treat
111   entire molecular subsystems as a single rigid body, is simplification
112   of the electrostatic interactions between these bodies so that fewer
113   site-site interactions are required to compute configurational
114 < energies. Notably, the force matching approaches of Voth and coworkers
115 < are an exciting development in their ability to represent realistic
116 < (and {\it reactive}) chemical systems for very large length scales and
122 < long times.  This approach utilized a coarse-graining in interaction
123 < space (CGIS) which fits an effective force for the coarse grained
124 < system using a variational force-matching method to a fine-grained
125 < simulation.\cite{Izvekov:2008wo}
114 > energies. To do this, the interactions between coarse-grained sites
115 > are typically taken to be point
116 > multipoles.\cite{Golubkov06,ISI:000276097500009,ISI:000298664400012}
117  
118 < The coarse-graining approaches of Ren \& coworkers,\cite{Golubkov06}
119 < and Essex \&
120 < coworkers,\cite{ISI:000276097500009,ISI:000298664400012}
121 < both utilize Gay-Berne
122 < ellipsoids~\cite{Berne72,Gay81,Luckhurst90,Cleaver96,Berardi98,Ravichandran:1999fk,Berardi99,Pasterny00}
123 < to model dispersive interactions and point multipoles to model
124 < electrostatics for entire molecules or functional groups.
125 <
135 < Ichiye and coworkers have recently introduced a number of very fast
136 < water models based on a ``sticky'' multipole model which are
137 < qualitatively better at reproducing the behavior of real water than
138 < the more common point-charge models (SPC/E, TIPnP).  The point charge
139 < models are also substantially more computationally demanding than the
140 < sticky multipole
141 < approach.\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} The
142 < SSDQO model requires the use of an approximate multipole expansion
143 < (AME) as the exact multipole expansion is quite expensive
144 < (particularly when handled via the Ewald sum).\cite{Ichiye:2006qy}
145 < Another particularly important use of point multipoles (and multipole
146 < polarizability) is in the very high-quality AMOEBA water model and
118 > Water, in particular, has been modeled recently with point multipoles
119 > up to octupolar
120 > order.\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} For
121 > maximum efficiency, these models require the use of an approximate
122 > multipole expansion as the exact multipole expansion can become quite
123 > expensive (particularly when handled via the Ewald
124 > sum).\cite{Ichiye:2006qy} Point multipoles and multipole
125 > polarizability have also been utilized in the AMOEBA water model and
126   related force fields.\cite{Ponder:2010fk,schnieders:124114,Ren:2011uq}
127  
128   Higher-order multipoles present a peculiar issue for molecular
# Line 154 | Line 133 | kernel to point multipoles.  We have developed two dis
133   a different orientation can cause energy discontinuities.
134  
135   This paper outlines an extension of the original DSF electrostatic
136 < kernel to point multipoles.  We have developed two distinct real-space
136 > kernel to point multipoles.  We describe two distinct real-space
137   interaction models for higher-order multipoles based on two truncated
138   Taylor expansions that are carried out at the cutoff radius.  We are
139   calling these models {\bf Taylor-shifted} and {\bf Gradient-shifted}
140   electrostatics.  Because of differences in the initial assumptions,
141 < the two methods yield related, but different expressions for energies,
142 < forces, and torques.  
141 > the two methods yield related, but somewhat different expressions for
142 > energies, forces, and torques.
143  
144   In this paper we outline the new methodology and give functional forms
145   for the energies, forces, and torques up to quadrupole-quadrupole
146   order.  We also compare the new methods to analytic energy constants
147 < for periodic arrays of point multipoles.  In the following paper, we
148 < provide extensive numerical comparisons to Ewald-based electrostatics
149 < in common simulation enviornments.  
147 > for periodic arrays of point multipoles. In the following paper, we
148 > provide numerical comparisons to Ewald-based electrostatics in common
149 > simulation enviornments.
150  
151   \section{Methodology}
152 + An efficient real-space electrostatic method involves the use of a
153 + pair-wise functional form,
154 + \begin{equation}
155 + V = \sum_i \sum_{j>i} V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) +
156 + \sum_i V_i^\mathrm{self}
157 + \end{equation}
158 + that is short-ranged and easily truncated at a cutoff radius,
159 + \begin{equation}
160 +  V_\mathrm{pair}(r_{ij},\Omega_i, \Omega_j) = \left\{
161 + \begin{array}{ll}
162 + V_\mathrm{approx} (r_{ij}, \Omega_i, \Omega_j) & \quad r \le r_c \\
163 + 0 & \quad r > r_c ,
164 + \end{array}
165 + \right.
166 + \end{equation}
167 + along with an easily computed self-interaction term ($\sum_i
168 + V_i^\mathrm{self}$) which has linear-scaling with the number of
169 + particles.  Here $\Omega_i$ and $\Omega_j$ represent orientational
170 + coordinates of the two sites.  The computational efficiency, energy
171 + conservation, and even some physical properties of a simulation can
172 + depend dramatically on how the $V_\mathrm{approx}$ function behaves at
173 + the cutoff radius. The goal of any approximation method should be to
174 + mimic the real behavior of the electrostatic interactions as closely
175 + as possible without sacrificing the near-linear scaling of a cutoff
176 + method.
177  
178   \subsection{Self-neutralization, damping, and force-shifting}
179   The DSF and Wolf methods operate by neutralizing the total charge
180   contained within the cutoff sphere surrounding each particle.  This is
181   accomplished by shifting the potential functions to generate image
182   charges on the surface of the cutoff sphere for each pair interaction
183 < computed within $R_\textrm{c}$. Damping using a complementary error
183 > computed within $r_c$. Damping using a complementary error
184   function is applied to the potential to accelerate convergence. The
185   potential for the DSF method, where $\alpha$ is the adjustable damping
186   parameter, is given by
187   \begin{equation*}
188 < V_\mathrm{DSF}(r) = C_a C_b \Biggr{[}
188 > V_\mathrm{DSF}(r) = C_i C_j \Biggr{[}
189   \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}
190 < - \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} + \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}
190 > - \frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c} + \left(\frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c^2}
191   + \frac{2\alpha}{\pi^{1/2}}
192 < \frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}
193 < \right)\left(r_{ij}-R_\mathrm{c}\right)\ \Biggr{]}
192 > \frac{\exp\left(-\alpha^2r_c^2\right)}{r_c}
193 > \right)\left(r_{ij}-r_c\right)\ \Biggr{]}
194   \label{eq:DSFPot}
195   \end{equation*}
196 + Note that in this potential and in all electrostatic quantities that
197 + follow, the standard $1/4 \pi \epsilon_{0}$ has been omitted for
198 + clarity.
199  
200   To insure net charge neutrality within each cutoff sphere, an
201   additional ``self'' term is added to the potential. This term is
# Line 199 | Line 206 | damping.\cite{deLeeuw80,Wolf99}
206   over the surface of the cutoff sphere.  A portion of the self term is
207   identical to the self term in the Ewald summation, and comes from the
208   utilization of the complimentary error function for electrostatic
209 < damping.\cite{deLeeuw80,Wolf99}
209 > damping.\cite{deLeeuw80,Wolf99} There have also been recent efforts to
210 > extend the Wolf self-neutralization method to zero out the dipole and
211 > higher order multipoles contained within the cutoff
212 > sphere.\cite{Fukuda:2011jk,Fukuda:2012yu,Fukuda:2013qv}
213  
214 < There have been recent efforts to extend the Wolf self-neutralization
215 < method to zero out the dipole and higher order multipoles contained
216 < within the cutoff
217 < sphere.\cite{Fukuda:2011jk,Fukuda:2012yu,Fukuda:2013qv}  
218 <
219 < In this work, we will be extending the idea of self-neutralization for
220 < the point multipoles in a similar way.  In Figure
221 < \ref{fig:shiftedMultipoles}, the central dipolar site $\mathbf{D}_i$
222 < is interacting with point dipole $\mathbf{D}_j$ and point quadrupole,
213 < $\mathbf{Q}_k$.  The self-neutralization scheme for point multipoles
214 < involves projecting opposing multipoles for sites $j$ and $k$ on the
215 < surface of the cutoff sphere.
214 > In this work, we extend the idea of self-neutralization for the point
215 > multipoles by insuring net charge-neutrality and net-zero moments
216 > within each cutoff sphere.  In Figure \ref{fig:shiftedMultipoles}, the
217 > central dipolar site $\mathbf{D}_i$ is interacting with point dipole
218 > $\mathbf{D}_j$ and point quadrupole, $\mathbf{Q}_k$.  The
219 > self-neutralization scheme for point multipoles involves projecting
220 > opposing multipoles for sites $j$ and $k$ on the surface of the cutoff
221 > sphere.  There are also significant modifications made to make the
222 > forces and torques go smoothly to zero at the cutoff distance.
223  
224   \begin{figure}
225   \includegraphics[width=3in]{SM}
# Line 222 | Line 229 | As in the point-charge approach, there is a contributi
229   \label{fig:shiftedMultipoles}
230   \end{figure}
231  
232 < As in the point-charge approach, there is a contribution from
233 < self-neutralization of site $i$.  The self term for multipoles is
232 > As in the point-charge approach, there is an additional contribution
233 > from self-neutralization of site $i$.  The self term for multipoles is
234   described in section \ref{sec:selfTerm}.
235  
236   \subsection{The multipole expansion}
# Line 237 | Line 244 | V_a(\mathbf r) = \frac{1}{4\pi\epsilon_0}
244   a$.  Then the electrostatic potential of object $\bf a$ at
245   $\mathbf{r}$ is given by
246   \begin{equation}
247 < V_a(\mathbf r) = \frac{1}{4\pi\epsilon_0}
247 > V_a(\mathbf r) =
248   \sum_{k \, \text{in \bf a}} \frac{q_k}{\lvert \mathbf{r} - \mathbf{r}_k \rvert}.
249   \end{equation}
250   The Taylor expansion in $r$ can be written using an implied summation
# Line 256 | Line 263 | V_{\bf a}(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\hat{M
263   can then be used to express the electrostatic potential on $\bf a$ in
264   terms of multipole operators,
265   \begin{equation}
266 < V_{\bf a}(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\hat{M}_{\bf a} \frac{1}{r}
266 > V_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r}
267   \end{equation}
268   where
269   \begin{equation}
# Line 281 | Line 288 | U_{\bf{ab}}(r)
288   $\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $), the interaction energy is given by
289   \begin{equation}
290   U_{\bf{ab}}(r)
291 < = \frac{1}{4\pi \epsilon_0} \hat{M}_a \sum_{j \, \text{in \bf b}} \frac {q_j}{\vert \bf{r}+\bf{r}_j \vert} .
291 > = \hat{M}_a \sum_{j \, \text{in \bf b}} \frac {q_j}{\vert \bf{r}+\bf{r}_j \vert} .
292   \end{equation}
293   This can also be expanded as a Taylor series in $r$.  Using a notation
294   similar to before to define the multipoles on object {\bf b},
# Line 292 | Line 299 | U_{\bf{ab}}(r)=\frac{\hat{M}_{\bf a} \hat{M}_{\bf b}}{
299   \end{equation}
300   we arrive at the multipole expression for the total interaction energy.
301   \begin{equation}
302 < U_{\bf{ab}}(r)=\frac{\hat{M}_{\bf a} \hat{M}_{\bf b}}{4\pi \epsilon_0} \frac{1}{r}  \label{kernel}.
302 > U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r}  \label{kernel}.
303   \end{equation}
304   This form has the benefit of separating out the energies of
305   interaction into contributions from the charge, dipole, and quadrupole
306 < of $\bf a$ interacting with the same multipoles $\bf b$.
306 > of $\bf a$ interacting with the same multipoles on $\bf b$.
307  
308   \subsection{Damped Coulomb interactions}
309   In the standard multipole expansion, one typically uses the bare
# Line 311 | Line 318 | either $1/r$ or $B_0(r)$, and all of the techniques ca
318   \int_{\alpha r}^{\infty} \text{e}^{-s^2} ds .
319   \end{equation}
320   We develop equations below using the function $f(r)$ to represent
321 < either $1/r$ or $B_0(r)$, and all of the techniques can be applied
322 < either to bare or damped Coulomb kernels as long as derivatives of
323 < these functions are known.  Smith's convenient functions $B_l(r)$ are
324 < summarized in Appendix A.
321 > either $1/r$ or $B_0(r)$, and all of the techniques can be applied to
322 > bare or damped Coulomb kernels (or any other function) as long as
323 > derivatives of these functions are known.  Smith's convenient
324 > functions $B_l(r)$ are summarized in Appendix A.
325  
319 \subsection{Taylor-shifted force (TSF) electrostatics}
320
326   The main goal of this work is to smoothly cut off the interaction
327   energy as well as forces and torques as $r\rightarrow r_c$.  To
328   describe how this goal may be met, we use two examples, charge-charge
329 < and charge-dipole, using the bare Coulomb kernel $f(r)=1/r$ to explain
330 < the idea.
329 > and charge-dipole, using the bare Coulomb kernel, $f(r)=1/r$, to
330 > explain the idea.
331  
332 + \subsection{Shifted-force methods}
333   In the shifted-force approximation, the interaction energy for two
334   charges $C_{\bf a}$ and $C_{\bf b}$ separated by a distance $r$ is
335   written:
336   \begin{equation}
337 < U_{C_{\bf a}C_{\bf b}}(r)=\frac{1}{4\pi \epsilon_0} C_{\bf a} C_{\bf b}
337 > U_{C_{\bf a}C_{\bf b}}(r)= C_{\bf a} C_{\bf b}
338   \left({ \frac{1}{r} - \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2}  }
339   \right) .
340   \end{equation}
341   Two shifting terms appear in this equations, one from the
342 < neutralization procedure ($-1/r_c$), and one that will make the first
343 < derivative also vanish at the cutoff radius.  
342 > neutralization procedure ($-1/r_c$), and one that causes the first
343 > derivative to vanish at the cutoff radius.
344  
345   Since one derivative of the interaction energy is needed for the
346   force, the minimal perturbation is a term linear in $(r-r_c)$ in the
# Line 345 | Line 351 | There are a number of ways to generalize this derivati
351   \right) = \left(- \frac{1}{r^2} + \frac{1}{r_c^2}
352   \right) .
353   \end{equation}
354 < There are a number of ways to generalize this derivative shift for
355 < higher-order multipoles.
354 > which clearly vanishes as the $r$ approaches the cutoff radius. There
355 > are a number of ways to generalize this derivative shift for
356 > higher-order multipoles.  Below, we present two methods, one based on
357 > higher-order Taylor series for $r$ near $r_c$, and the other based on
358 > linear shift of the kernel gradients at the cutoff itself.
359  
360 + \subsection{Taylor-shifted force (TSF) electrostatics}
361   In the Taylor-shifted force (TSF) method, the procedure that we follow
362   is based on a Taylor expansion containing the same number of
363   derivatives required for each force term to vanish at the cutoff.  For
364   example, the quadrupole-quadrupole interaction energy requires four
365   derivatives of the kernel, and the force requires one additional
366 < derivative. We therefore require shifted energy expressions that
367 < include enough terms so that all energies, forces, and torques are
368 < zero as $r \rightarrow r_c$.  In each case, we will subtract off a
369 < function $f_n^{\text{shift}}(r)$ from the kernel $f(r)=1/r$.  The
370 < index $n$ indicates the number of derivatives to be taken when
371 < deriving a given multipole energy.  We choose a function with
372 < guaranteed smooth derivatives --- a truncated Taylor series of the
373 < function $f(r)$, e.g.,
366 > derivative. For quadrupole-quadrupole interactions, we therefore
367 > require shifted energy expressions that include up to $(r-r_c)^5$ so
368 > that all energies, forces, and torques are zero as $r \rightarrow
369 > r_c$. In each case, we subtract off a function $f_n^{\text{shift}}(r)$
370 > from the kernel $f(r)=1/r$.  The subscript $n$ indicates the number of
371 > derivatives to be taken when deriving a given multipole energy.  We
372 > choose a function with guaranteed smooth derivatives -- a truncated
373 > Taylor series of the function $f(r)$, e.g.,
374   %
375   \begin{equation}
376 < f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert  _{r_c}  .
376 > f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)}(r_c)  .
377   \end{equation}
378   %
379   The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$.
# Line 378 | Line 388 | U_{C_{\bf a}D_{\bf b}}(r)=
388   %
389   \begin{equation}
390   U_{C_{\bf a}D_{\bf b}}(r)=
391 < \frac{C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}  \frac {\partial f_1(r) }{\partial r_\alpha}
392 < =\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}  
391 > C_{\bf a} D_{{\bf b}\alpha}  \frac {\partial f_1(r) }{\partial r_\alpha}
392 > = C_{\bf a} D_{{\bf b}\alpha}
393   \frac {r_\alpha}{r} \frac {\partial f_1(r)}{\partial r} .
394   \end{equation}
395   %
396 < The force that dipole  $\bf b$ puts on charge $\bf a$ is
396 > The force that dipole  $\bf b$ exerts on charge $\bf a$ is
397   %
398   \begin{equation}
399 < F_{C_{\bf a}D_{\bf b}\beta} =\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}
399 > F_{C_{\bf a}D_{\bf b}\beta} = C_{\bf a} D_{{\bf b}\alpha}
400   \left[ \frac{\delta_{\alpha\beta}}{r} \frac {\partial}{\partial r} +
401   \frac{r_\alpha r_\beta}{r^2}
402   \left( -\frac{1}{r} \frac {\partial} {\partial r}
403   + \frac {\partial ^2} {\partial r^2} \right) \right] f_1(r) .
404   \end{equation}
405   %
406 < For $f(r)=1/r$, we find
406 > For undamped coulombic interactions, $f(r)=1/r$, we find
407   %
408   \begin{equation}
409   F_{C_{\bf a}D_{\bf b}\beta} =
410 < \frac{C_{\bf a} D_{{\bf b}\beta} }{4\pi \epsilon_0r}
410 > \frac{C_{\bf a} D_{{\bf b}\beta}}{r}
411   \left[  -\frac{1}{r^2}+\frac{1}{r_c^2}-\frac{2(r-r_c)}{r_c^3} \right]
412 < +\frac{C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta }{4\pi \epsilon_0}
412 > +C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta
413   \left[ \frac{3}{r^5}-\frac{3}{r^3r_c^2} \right] .
414   \end{equation}
415   %
416   This expansion shows the expected $1/r^3$ dependence of the force.  
417  
418 < In general, we write
418 > In general, we can write
419   %
420   \begin{equation}
421 < U=\frac{1}{4\pi \epsilon_0} (\text{prefactor}) (\text{derivatives}) f_n(r)
421 > U= (\text{prefactor}) (\text{derivatives}) f_n(r)
422   \label{generic}
423   \end{equation}
424   %
425 < where $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole
426 < and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for quadrupole-quadrupole.
427 < An example is the case of quadrupole-quadrupole for which the $\text{prefactor}$ is
428 < $Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$ and the derivatives are
429 < $\partial^4/\partial r_\alpha \partial r_\beta \partial r_\gamma \partial r_\delta$, with
430 < implied summation combining the space indices.
425 > with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for
426 > charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and
427 > $n=4$ for quadrupole-quadrupole.  For example, in
428 > quadrupole-quadrupole interactions for which the $\text{prefactor}$ is
429 > $Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$, the derivatives are
430 > $\partial^4/\partial r_\alpha \partial r_\beta \partial
431 > r_\gamma \partial r_\delta$, with implied summation combining the
432 > space indices.
433  
434 < To apply this method  to the smeared-charge approach,  
435 < we write $f(r)=\text{erfc}(\alpha r)/r$.  By using one function $f(r)$ for both
436 < approaches, we simplify the tabulation of equations used.  Because
437 < of the many derivatives that are taken, the algebra is tedious and are summarized
438 < in Appendices A and B.  
434 > In the formulas presented in the tables below, the placeholder
435 > function $f(r)$ is used to represent the electrostatic kernel (either
436 > damped or undamped).  The main functions that go into the force and
437 > torque terms, $g_n(r), h_n(r), s_n(r), \mathrm{~and~} t_n(r)$ are
438 > successive derivatives of the shifted electrostatic kernel, $f_n(r)$
439 > of the same index $n$.  The algebra required to evaluate energies,
440 > forces and torques is somewhat tedious, so only the final forms are
441 > presented in tables \ref{tab:tableenergy} and \ref{tab:tableFORCE}.
442  
443   \subsection{Gradient-shifted force (GSF) electrostatics}
444 <
445 < Note the method used in the previous subsection to shift a force is basically that of using
446 < a truncated Taylor Series in the radius $r$.  An alternate method exists, best explained by
447 < writing one shifted formula for all interaction energies $U(r)$:
444 > The second, and conceptually simpler approach to force-shifting
445 > maintains only the linear $(r-r_c)$ term in the truncated Taylor
446 > expansion, and has a similar interaction energy for all multipole
447 > orders:
448   \begin{equation}
449 < U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big \lvert  _{r_c} .
449 > U^{\text{GSF}} =
450 > U(\mathbf{r}, \hat{\mathbf{a}}, \hat{\mathbf{b}}) -
451 > U(\mathbf{r}_c,\hat{\mathbf{a}}, \hat{\mathbf{b}}) - (r-r_c) \hat{r}
452 > \cdot \nabla U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) \Big \lvert  _{r_c} .
453 > \label{generic2}
454   \end{equation}
455 < Note that this method uses only the linear term, $(r-r_c)$ in the Taylor series, no higher order terms
456 < $(r-r_c)^n$ appear.   The primary difference between methods 1 and 2 originates
457 < with the stage in the derivation where the Taylor Series is applied; in method 1, it is applied to the
458 < kernel.  In method 2, it is applied to individual interaction energies of the multipole expansion.
459 < Terms from this method thus have the general form:
455 > Both the potential and the gradient for force shifting are evaluated
456 > for an image multipole projected onto the surface of the cutoff sphere
457 > (see fig \ref{fig:shiftedMultipoles}).  The image multipole retains
458 > the orientation ($\hat{\mathbf{b}}$) of the interacting multipole.  No
459 > higher order terms $(r-r_c)^n$ appear.  The primary difference between
460 > the TSF and GSF methods is the stage at which the Taylor Series is
461 > applied; in the Taylor-shifted approach, it is applied to the kernel
462 > itself.  In the Gradient-shifted approach, it is applied to individual
463 > radial interactions terms in the multipole expansion.  Energies from
464 > this method thus have the general form:
465   \begin{equation}
466 < U=\frac{1}{4\pi \epsilon_0}\sum  (\text{angular factor}) (\text{radial factor}).
467 < \label{generic2}
466 > U= \sum  (\text{angular factor}) (\text{radial factor}).
467 > \label{generic3}
468   \end{equation}
469  
470 < Results for both methods can be summarized using the form of Eq.~(\ref{generic2})
471 < and are listed in Table I below.
470 > Functional forms for both methods (TSF and GSF) can both be summarized
471 > using the form of Eq.~(\ref{generic3}).  The basic forms for the
472 > energy, force, and torque expressions are tabulated for both shifting
473 > approaches below -- for each separate orientational contribution, only
474 > the radial factors differ between the two methods.
475  
476   \subsection{\label{sec:level2}Body and space axes}
477 + Although objects $\bf a$ and $\bf b$ rotate during a molecular
478 + dynamics (MD) simulation, their multipole tensors remain fixed in
479 + body-frame coordinates. While deriving force and torque expressions,
480 + it is therefore convenient to write the energies, forces, and torques
481 + in intermediate forms involving the vectors of the rotation matrices.
482 + We denote body axes for objects $\bf a$ and $\bf b$ using unit vectors
483 + $\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$.
484 + In a typical simulation , the initial axes are obtained by
485 + diagonalizing the moment of inertia tensors for the objects.  (N.B.,
486 + the body axes are generally {\it not} the same as those for which the
487 + quadrupole moment is diagonal.)  The rotation matrices are then
488 + propagated during the simulation.
489  
490 < Up to this point, all energies and forces have been written in terms of fixed space
491 < coordinates $x$, $y$, $z$.  Interaction energies are computed from the generic formulas Eq.~(\ref{generic}) and ~(\ref{generic2}) which
453 < combine prefactors with radial functions.  But because objects
454 < $\bf a$ and  $\bf b$ both translate and rotate as part of a MD simulation,
455 < it is desirable to contract all $r$-dependent terms with dipole and quadrupole
456 < moments expressed in terms of their body axes.  
457 < Since the interaction energy expressions will be used to derive both forces and torques,
458 < we follow here the development of Allen and Germano, which was itself based on an
459 < earlier paper by Price \em et al.\em
460 <
461 < Denote body axes for objects  $\bf a$ and  $\bf b$ by unit vectors
462 < $\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$ referring to a convenient
463 < set of inertial body axes.  (Note, these body axes are generally not the same as those for which the
464 < quadrupole moment is diagonal.)  Then,
465 < %
490 > The rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be
491 > expressed using these unit vectors:
492   \begin{eqnarray}
467 \hat{a}_m= a_{mx}\hat{x} + a_{my}\hat{y} + a_{mz}\hat{z}  \\
468 \hat{b}_m= b_{mx}\hat{x} + b_{my}\hat{y} + b_{mz}\hat{z}  .
469 \end{eqnarray}
470 Allen and Germano define matrices $\hat{\mathbf {a}}$
471 and  $\hat{\mathbf {b}}$ using these unit vectors:
472 \begin{eqnarray}
493   \hat{\mathbf {a}} =
494   \begin{pmatrix}
495   \hat{a}_1 \\
496   \hat{a}_2 \\
497   \hat{a}_3
498 < \end{pmatrix}
479 < =
480 < \begin{pmatrix}
481 < a_{1x} \quad a_{1y} \quad a_{1z} \\
482 < a_{2x} \quad a_{2y} \quad a_{2z} \\
483 < a_{3x} \quad a_{3y} \quad a_{3z}
484 < \end{pmatrix}\\
498 > \end{pmatrix}, \qquad
499   \hat{\mathbf {b}} =
500   \begin{pmatrix}
501   \hat{b}_1 \\
502   \hat{b}_2 \\
503   \hat{b}_3
504   \end{pmatrix}
491 =
492 \begin{pmatrix}
493 b_{1x}\quad  b_{1y} \quad b_{1z} \\
494 b_{2x} \quad b_{2y} \quad b_{2z} \\
495 b_{3x} \quad b_{3y} \quad b_{3z}
496 \end{pmatrix}  .
505   \end{eqnarray}
506   %
507 < These matrices convert from space-fixed $(xyz)$ to object-fixed $(123)$ coordinates.
508 < All contractions of prefactors with derivatives of functions can be written in terms of these matrices.
509 < It proves to be equally convenient to just write any contraction in terms of unit vectors
510 < $\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$.  
511 < We have found it useful to write angular-dependent terms in three different fashions,
512 < illustrated by the following
513 < three examples from the interaction-energy expressions:
507 > These matrices convert from space-fixed $(xyz)$ to body-fixed $(123)$
508 > coordinates.
509 >
510 > Allen and Germano,\cite{Allen:2006fk} following earlier work by Price
511 > {\em et al.},\cite{Price:1984fk} showed that if the interaction
512 > energies are written explicitly in terms of $\hat{r}$ and the body
513 > axes ($\hat{a}_m$, $\hat{b}_n$) :
514   %
515 + \begin{equation}
516 + U(r, \{\hat{a}_m \cdot \hat{r} \},
517 + \{\hat{b}_n\cdot \hat{r} \},
518 + \{\hat{a}_m \cdot \hat{b}_n \}) .
519 + \label{ugeneral}
520 + \end{equation}
521 + %
522 + the forces come out relatively cleanly,
523 + %
524 + \begin{equation}
525 + \mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} =  \frac{\partial U}{\partial \mathbf{r}}
526 + = \frac{\partial U}{\partial r} \hat{r}
527 + + \sum_m \left[
528 + \frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})}
529 + \frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}}
530 + + \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})}
531 + \frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}}
532 + \right] \label{forceequation}.
533 + \end{equation}
534 +
535 + The torques can also be found in a relatively similar
536 + manner,
537 + %
538   \begin{eqnarray}
539 < \mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}}
540 < =D_{\bf {a}\alpha} D_{\bf {b}\alpha}=
541 < \sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}} \\
542 < r^2 \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right)=
543 < r_\alpha Q_{\bf b \alpha \beta} r_\beta = r^2
544 < \sum_{mn}(\hat{r} \cdot \hat{b}_m)  Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \\
545 < r ( \mathbf{D}_{\mathbf{a}} \cdot
546 < \mathbf{Q}_{\mathbf{b}}  \cdot \hat{r})=
547 < D_{\bf {a}\alpha}  Q_{\bf b \alpha \beta} r_\beta
548 < =r \sum_{lmn} D_{\mathbf{a}l} (\hat{a}_l \cdot \hat{b}_m ) Q_{\mathbf{b}mn}
549 < (\hat{b}_n \cdot \hat{r}) .
539 > \mathbf{\tau}_{\bf a} =
540 > \sum_m
541 > \frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})}
542 > ( \hat{r} \times \hat{a}_m )
543 > -\sum_{mn}
544 > \frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)}
545 > (\hat{a}_m \times \hat{b}_n) \\
546 > %
547 > \mathbf{\tau}_{\bf b} =
548 > \sum_m
549 > \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})}
550 > ( \hat{r} \times \hat{b}_m)
551 > +\sum_{mn}
552 > \frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)}
553 > (\hat{a}_m \times \hat{b}_n) .
554   \end{eqnarray}
555 +
556 + Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $
557 + is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk}
558 + We also made use of the identities,
559   %
560 < [Dan, perhaps there are better examples to show here.]
560 > \begin{align}
561 > \frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}}
562 > =& \frac{1}{r} \left(  \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r}
563 > \right) \\
564 > \frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}}
565 > =& \frac{1}{r} \left(  \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r}
566 > \right) .
567 > \end{align}
568  
569 < In each line, the first two terms are written using space coordinates.  The first form is standard
570 < in the chemistry literature, and the second is ``physicist style'' using implied summation notation.  The third
571 < form shows explicitly sums over body indices and the dot products now indicate contractions using space indices.
572 < We find the first form to be useful in writing equations prior to converting to computer code.  The second form is helpful
573 < in derivations of the interaction energy expressions.  The third one is specifically helpful when deriving forces and torques, as will
574 < be discussed below.  
569 > Many of the multipole contractions required can be written in one of
570 > three equivalent forms using the unit vectors $\hat{r}$, $\hat{a}_m$,
571 > and $\hat{b}_n$. In the torque expressions, it is useful to have the
572 > angular-dependent terms available in all three fashions, e.g. for the
573 > dipole-dipole contraction:
574 > %
575 > \begin{equation}
576 > \mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}}
577 > = D_{\bf {a}\alpha} D_{\bf {b}\alpha} =
578 > \sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}}
579 > \end{equation}
580 > %
581 > The first two forms are written using space coordinates.  The first
582 > form is standard in the chemistry literature, while the second is
583 > expressed using implied summation notation.  The third form shows
584 > explicit sums over body indices and the dot products now indicate
585 > contractions using space indices.
586  
587 + In computing our force and torque expressions, we carried out most of
588 + the work in body coordinates, and have transformed the expressions
589 + back to space-frame coordinates, which are reported below.  Interested
590 + readers may consult the supplemental information for this paper for
591 + the intermediate body-frame expressions.
592  
593   \subsection{The Self-Interaction \label{sec:selfTerm}}
594  
595 < The Wolf summation~\cite{Wolf99} and the later damped shifted force
596 < (DSF) extension~\cite{Fennell06} included self-interactions that are
597 < handled separately from the pairwise interactions between sites. The
598 < self-term is normally calculated via a single loop over all sites in
599 < the system, and is relatively cheap to evaluate. The self-interaction
600 < has contributions from two sources:
601 < \begin{itemize}
602 < \item The neutralization procedure within the cutoff radius requires a
603 <  contribution from a charge opposite in sign, but equal in magnitude,
604 <  to the central charge, which has been spread out over the surface of
605 <  the cutoff sphere.  For a system of undamped charges, the total
606 <  self-term is
595 > In addition to cutoff-sphere neutralization, the Wolf
596 > summation~\cite{Wolf99} and the damped shifted force (DSF)
597 > extension~\cite{Fennell:2006zl} also included self-interactions that
598 > are handled separately from the pairwise interactions between
599 > sites. The self-term is normally calculated via a single loop over all
600 > sites in the system, and is relatively cheap to evaluate. The
601 > self-interaction has contributions from two sources.
602 >
603 > First, the neutralization procedure within the cutoff radius requires
604 > a contribution from a charge opposite in sign, but equal in magnitude,
605 > to the central charge, which has been spread out over the surface of
606 > the cutoff sphere.  For a system of undamped charges, the total
607 > self-term is
608   \begin{equation}
609   V_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2}
610   \label{eq:selfTerm}
611   \end{equation}
612 < Note that in this potential and in all electrostatic quantities that
613 < follow, the standard $4 \pi \epsilon_{0}$ has been omitted for
614 < clarity.
615 < \item Charge damping with the complementary error function is a
616 <  partial analogy to the Ewald procedure which splits the interaction
617 <  into real and reciprocal space sums.  The real space sum is retained
618 <  in the Wolf and DSF methods.  The reciprocal space sum is first
619 <  minimized by folding the largest contribution (the self-interaction)
620 <  into the self-interaction from charge neutralization of the damped
621 <  potential.  The remainder of the reciprocal space portion is then
622 <  discarded (as this contributes the largest computational cost and
560 <  complexity to the Ewald sum).  For a system containing only damped
561 <  charges, the complete self-interaction can be written as
612 >
613 > Second, charge damping with the complementary error function is a
614 > partial analogy to the Ewald procedure which splits the interaction
615 > into real and reciprocal space sums.  The real space sum is retained
616 > in the Wolf and DSF methods.  The reciprocal space sum is first
617 > minimized by folding the largest contribution (the self-interaction)
618 > into the self-interaction from charge neutralization of the damped
619 > potential.  The remainder of the reciprocal space portion is then
620 > discarded (as this contributes the largest computational cost and
621 > complexity to the Ewald sum).  For a system containing only damped
622 > charges, the complete self-interaction can be written as
623   \begin{equation}
624   V_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} +
625    \frac{\alpha}{\sqrt{\pi}} \right) \sum_{{\bf a}=1}^N
626        C_{\bf a}^{2}.
627   \label{eq:dampSelfTerm}
628   \end{equation}
568 \end{itemize}
629  
630   The extension of DSF electrostatics to point multipoles requires
631   treatment of {\it both} the self-neutralization and reciprocal
# Line 594 | Line 654 | Quadrupole & $2 \text{Tr}(\mathbf{Q}_{\bf a}^2) + \tex
654   Charge & $C_{\bf a}^2$ & $-f(r_c)$ & $-\frac{\alpha}{\sqrt{\pi}}$ \\
655   Dipole & $|\mathbf{D}_{\bf a}|^2$ & $\frac{1}{3} \left( h(r_c) +
656    \frac{2 g(r_c)}{r_c} \right)$ & $-\frac{2 \alpha^3}{3 \sqrt{\pi}}$\\
657 < Quadrupole & $2 \text{Tr}(\mathbf{Q}_{\bf a}^2) + \text{Tr}(\mathbf{Q}_{\bf a})^2$ &
657 > Quadrupole & $2 \mathbf{Q}_{\bf a}:\mathbf{Q}_{\bf a} + \text{Tr}(\mathbf{Q}_{\bf a})^2$ &
658   $- \frac{1}{15} \left( t(r_c)+ \frac{4 s(r_c)}{r_c} \right)$ &
659   $-\frac{4 \alpha^5}{5 \sqrt{\pi}}$ \\
660   Charge-Quadrupole & $-2 C_{\bf a} \text{Tr}(\mathbf{Q}_{\bf a})$ & $\frac{1}{3} \left(
# Line 608 | Line 668 | terms, $f(r), g(r), h(r), s(r), \mathrm{~and~} t(r)$ a
668   multipole orders.  Symmetry prevents the charge-dipole and
669   dipole-quadrupole interactions from contributing to the
670   self-interaction.  The functions that go into the self-neutralization
671 < terms, $f(r), g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive
672 < derivatives of the electrostatic kernel (either the undamped $1/r$ or
673 < the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that are
674 < evaluated at the cutoff distance.  For undamped interactions, $f(r_c)
675 < = 1/r_c$, $g(r_c) = -1/r_c^{2}$, and so on.  For damped interactions,
676 < $f(r_c) = B_0(r_c)$, $g(r_c) = B_0'(r_c)$, and so on.  Appendix XX
677 < contains recursion relations that allow rapid evaluation of these
678 < derivatives.
671 > terms, $g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive
672 > derivatives of the electrostatic kernel, $f(r)$ (either the undamped
673 > $1/r$ or the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that
674 > have been evaluated at the cutoff distance.  For undamped
675 > interactions, $f(r_c) = 1/r_c$, $g(r_c) = -1/r_c^{2}$, and so on.  For
676 > damped interactions, $f(r_c) = B_0(r_c)$, $g(r_c) = B_0'(r_c)$, and so
677 > on.  Appendix \ref{SmithFunc} contains recursion relations that allow
678 > rapid evaluation of these derivatives.
679  
680 < \section{Energies, forces, and torques}
681 < \subsection{Interaction energies}
682 <
683 < We now list multipole interaction energies using a set of generic
684 < radial functions.  Table \ref{tab:tableenergy} maps between the
685 < generic functions and the radial functions derived for both the
686 < Taylor-shifted and Gradient-shifted methods.  This set of equations is
687 < written in terms of space coordinates:
680 > \section{Interaction energies, forces, and torques}
681 > The main result of this paper is a set of expressions for the
682 > energies, forces and torques (up to quadrupole-quadrupole order) that
683 > work for both the Taylor-shifted and Gradient-shifted approximations.
684 > These expressions were derived using a set of generic radial
685 > functions.  Without using the shifting approximations mentioned above,
686 > some of these radial functions would be identical, and the expressions
687 > coalesce into the familiar forms for unmodified multipole-multipole
688 > interactions.  Table \ref{tab:tableenergy} maps between the generic
689 > functions and the radial functions derived for both the Taylor-shifted
690 > and Gradient-shifted methods.  The energy equations are written in
691 > terms of lab-frame representations of the dipoles, quadrupoles, and
692 > the unit vector connecting the two objects,
693  
694   % Energy in space coordinate form ----------------------------------------------------------------------------------------------
695   %
# Line 633 | Line 698 | U_{C_{\bf a}C_{\bf b}}(r)=&
698   %
699   \begin{align}
700   U_{C_{\bf a}C_{\bf b}}(r)=&
701 < \frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r)  \label{uchch}
701 > C_{\bf a} C_{\bf b}  v_{01}(r)  \label{uchch}
702   \\
703   %
704   % u ca db
705   %
706   U_{C_{\bf a}D_{\bf b}}(r)=&
707 < \frac{C_{\bf a}}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)  v_{11}(r)  
707 > C_{\bf a} \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)  v_{11}(r)  
708   \label{uchdip}
709   \\
710   %
711   % u ca qb
712   %
713 < U_{C_{\bf a}Q_{\bf b}}(r)=&
714 < \frac{C_{\bf a }}{4\pi \epsilon_0} \Bigl[ \text{Tr}Q_{\bf b} v_{21}(r)
715 < \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{22}(r) \Bigr]
713 > U_{C_{\bf a}Q_{\bf b}}(r)=& C_{\bf a } \Bigl[ \text{Tr}Q_{\bf b}
714 > v_{21}(r) + \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot
715 >  \hat{r} \right) v_{22}(r) \Bigr]
716   \label{uchquad}
717   \\
718   %
# Line 661 | Line 726 | -\frac{1}{4\pi \epsilon_0} \Bigr[ \left( \mathbf{D}_{\
726   % u da db
727   %
728   U_{D_{\bf a}D_{\bf b}}(r)=&
729 < -\frac{1}{4\pi \epsilon_0} \Bigr[ \left( \mathbf{D}_{\mathbf {a}} \cdot
729 > -\Bigr[ \left( \mathbf{D}_{\mathbf {a}} \cdot
730   \mathbf{D}_{\mathbf{b}} \right)  v_{21}(r)
731   +\left( \mathbf{D}_{\mathbf {a}} \cdot \hat{r} \right)
732   \left( \mathbf{D}_{\mathbf {b}} \cdot \hat{r} \right)  
# Line 674 | Line 739 | -\frac{1}{4\pi \epsilon_0} \Bigl[
739   \begin{split}
740   % 1
741   U_{D_{\bf a}Q_{\bf b}}(r) =&
742 < -\frac{1}{4\pi \epsilon_0} \Bigl[
742 > -\Bigl[
743   \text{Tr}\mathbf{Q}_{\mathbf{b}}
744   \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
745   +2 ( \mathbf{D}_{\mathbf{a}} \cdot
746   \mathbf{Q}_{\mathbf{b}} \cdot \hat{r} ) \Bigr] v_{31}(r) \\
747   % 2
748 < &-\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
748 > &- \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
749   \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{32}(r)
750   \label{udipquad}
751   \end{split}
# Line 717 | Line 782 | U_{Q_{\bf a}Q_{\bf b}}(r)=&
782   \begin{split}
783   %1
784   U_{Q_{\bf a}Q_{\bf b}}(r)=&
785 < \frac{1}{4\pi \epsilon_0} \Bigl[
785 > \Bigl[
786   \text{Tr} \mathbf{Q}_{\mathbf{a}} \text{Tr} \mathbf{Q}_{\mathbf{b}}
787 < +2 \text{Tr} \left(
788 < \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} \right) \Bigr] v_{41}(r)
787 > +2
788 > \mathbf{Q}_{\mathbf{a}} : \mathbf{Q}_{\mathbf{b}} \Bigr] v_{41}(r)
789   \\
790   % 2
791 < &+ \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
791 > &+\Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
792   \left( \hat{r} \cdot
793   \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right)
794   +\text{Tr}\mathbf{Q}_{\mathbf{b}}
# Line 733 | Line 798 | +\text{Tr}\mathbf{Q}_{\mathbf{b}}
798   \Bigr] v_{42}(r)
799   \\
800   % 4
801 < &+ \frac{1}{4\pi \epsilon_0}
801 > &+
802   \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right)
803   \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}}  \cdot \hat{r} \right) v_{43}(r).
804   \label{uquadquad}
805   \end{split}
806   \end{align}
807 <
807 > %
808   Note that the energies of multipoles on site $\mathbf{b}$ interacting
809   with those on site $\mathbf{a}$ can be obtained by swapping indices
810   along with the sign of the intersite vector, $\hat{r}$.
# Line 749 | Line 814 | along with the sign of the intersite vector, $\hat{r}$
814   % TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
815   %
816  
817 < \begin{table*}
818 < \caption{\label{tab:tableenergy}Radial functions used in the energy and torque equations.  Functions
819 < used in this table are defined in Appendices B and C.}
820 < \begin{ruledtabular}
821 < \begin{tabular}{|l|c|l|l}
822 < Generic&Coulomb&Taylor-Shifted&Gradient-Shifted
817 > \begin{sidewaystable}
818 >  \caption{\label{tab:tableenergy}Radial functions used in the energy
819 >    and torque equations.  The $f, g, h, s, t, \mathrm{and} u$
820 >    functions used in this table are defined in Appendices B and C.}
821 > \begin{tabular}{|c|c|l|l|} \hline
822 > Generic&Bare Coulomb&Taylor-Shifted&Gradient-Shifted
823   \\ \hline
824   %
825   %
# Line 786 | Line 851 | -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right) $ \\
851   $\frac{3}{r^3}  $ &
852   $\left(-\frac{g_2(r)}{r} + h_2(r) \right)$ &
853   $\left(-\frac{g(r)}{r}+h(r) \right)
854 < -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right) $ \\
855 < &&&$  -(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$
854 > -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right)$  \\
855 > &&& $ ~~~-(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$
856   \\
857   %
858   %
# Line 798 | Line 863 | -\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right
863   $\left(-\frac{g_3(r)}{r^2} + \frac{h_3(r)}{r} \right)$ &
864   $\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r} \right)
865   -\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right) $\\
866 < &&&$  -(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)$
866 > &&&$ ~~~ -(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)$
867   \\
868   %
869   $v_{32}(r)$ &
# Line 806 | Line 871 | - \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} +
871   $\left( \frac{3g_3(r)}{r^2} - \frac{3h_3(r)}{r} + s_3(r) \right)$ &
872   $\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)
873   - \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c) \right)$ \\
874 < &&&$  -(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)$
874 > &&&$ ~~~ -(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)$
875   \\
876   %
877   %
# Line 817 | Line 882 | - \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2}
882   $\left(-\frac{g_4(r)}{r^3} +\frac{h_4(r)}{r^2} \right) $  &
883   $\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right)
884   - \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2} \right)$ \\
885 < &&&$  -(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)$
885 > &&&$ ~~~ -(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)$
886   \\
887   % 2
888   $v_{42}(r)$ &
# Line 825 | Line 890 | -\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+
890   $\left( \frac{3g_4(r)}{r^3} - \frac{3h_4(r)}{r^2}+\frac{s_4(r)}{r} \right)$ &
891   $\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right)
892   -\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ \\
893 < &&&$  -(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}
893 > &&&$ ~~~ -(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}
894   -\frac{4s(r_c)}{r_c^2} + \frac{t(r_c)}{r_c}\right)$
895   \\
896   % 3
# Line 833 | Line 898 | $\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s
898   $ \frac{105}{r^5}  $ &
899   $\left(-\frac{15g_4(r)}{r^3}+\frac{15h_4(r)}{r^2}-\frac{6s_4(r)}{r} + t_4(r)\right) $ &
900   $\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s(r)}{r} + t(r)\right)$ \\
901 < &&&$ -\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)$ \\
902 < &&&$ -(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}
903 < -\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\
901 > &&&$~~~ -\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)$ \\
902 > &&&$~~~ -(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}
903 > -\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\ \hline
904   \end{tabular}
905 < \end{ruledtabular}
841 < \end{table*}
905 > \end{sidewaystable}
906   %
907   %
908   % FORCE  TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
909   %
910  
911 < \begin{table}
911 > \begin{sidewaystable}
912   \caption{\label{tab:tableFORCE}Radial functions used in the force equations.}
913 < \begin{ruledtabular}
914 < \begin{tabular}{cc}
851 < Generic&Method 1 or Method 2
913 > \begin{tabular}{|c|c|l|l|} \hline
914 > Function&Definition&Taylor-Shifted&Gradient-Shifted
915   \\ \hline
916   %
917   %
918   %
919   $w_a(r)$&
920 < $\frac{d v_{01}}{dr}$ \\
920 > $\frac{d v_{01}}{dr}$&
921 > $g_0(r)$&
922 > $g(r)-g(r_c)$ \\
923   %
924   %
925   $w_b(r)$ &
926 < $\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $ \\
926 > $\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $&
927 > $\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ &
928 > $h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\
929   %
930   $w_c(r)$ &
931 < $\frac{v_{11}(r)}{r}$ \\
931 > $\frac{v_{11}(r)}{r}$ &
932 > $\frac{g_1(r)}{r} $ &
933 > $\frac{v_{11}(r)}{r}$\\
934   %
935   %
936   $w_d(r)$&
937 < $\frac{d v_{21}}{dr}$ \\
937 > $\frac{d v_{21}}{dr}$&
938 > $\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ &
939 > $\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)
940 > -\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $ \\
941   %
942   $w_e(r)$ &
943 + $\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ &
944 + $\frac{v_{22}(r)}{r}$ &
945   $\frac{v_{22}(r)}{r}$ \\
946   %
947   %
948   $w_f(r)$&
949 < $\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$\\
949 > $\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$&
950 > $\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ &
951 >  $ \left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) $  \\
952 >  &&& $ ~~~- \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c)
953 >     \right)-\frac{2v_{22}(r)}{r}$\\
954   %
955   $w_g(r)$&
956 + $\frac{v_{31}(r)}{r}$&
957 + $ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$&
958   $\frac{v_{31}(r)}{r}$\\
959   %
960   $w_h(r)$ &
961 < $\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$\\
961 > $\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$&
962 > $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
963 > $ \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) $ \\
964 > &&& $ ~~~ -\frac{v_{31}(r)}{r}$ \\
965   % 2
966   $w_i(r)$ &
967 < $\frac{v_{32}(r)}{r}$ \\
967 > $\frac{v_{32}(r)}{r}$ &
968 > $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
969 > $\frac{v_{32}(r)}{r}$\\
970   %
971   $w_j(r)$ &
972 < $\frac{d v_{32}}{dr}  - \frac{3v_{32}}{r}$ \\
972 > $\frac{d v_{32}}{dr}  - \frac{3v_{32}}{r}$&
973 > $\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right)  $ &
974 > $\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right)$ \\
975 > &&& $~~~-\left(\frac{-6g(_cr)}{r_c^3} +\frac{6h(r_c)}{r_c^2}
976 >  -\frac{3s(r_c)}{r_c} +t(r_c) \right) -\frac{3v_{32}}{r}$ \\
977   %
978   $w_k(r)$ &
979 < $\frac{d v_{41}}{dr} $  \\
979 > $\frac{d v_{41}}{dr} $ &
980 > $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
981 > $\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}  \right)
982 > -\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2}  \right)$ \\
983   %
984   $w_l(r)$ &
985 < $\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ \\
985 > $\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ &
986 > $\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)$ &
987 > $\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\
988 > &&& $~~~ -\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)
989 > -\frac{2v_{42}(r)}{r}$\\
990   %
991   $w_m(r)$ &
992 < $\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$ \\
992 > $\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$&
993 > $\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)$ &
994 > $\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$\\
995 > &&& $~~~- \left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}
996 > +\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $\\
997 > &&& $~~~-\frac{4v_{43}(r)}{r}$  \\
998   %
999   $w_n(r)$ &
1000 < $\frac{v_{42}(r)}{r}$ \\
1000 > $\frac{v_{42}(r)}{r}$ &
1001 > $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
1002 > $\frac{v_{42}(r)}{r}$\\
1003   %
1004   $w_o(r)$ &
1005 < $\frac{v_{43}(r)}{r}$ \\
1005 > $\frac{v_{43}(r)}{r}$&
1006 > $\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)$ &
1007 > $\frac{v_{43}(r)}{r}$  \\ \hline
1008   %
1009  
1010   \end{tabular}
1011 < \end{ruledtabular}
907 < \end{table}
1011 > \end{sidewaystable}
1012   %
1013   %
1014   %
1015  
1016   \subsection{Forces}
1017 <
1018 < The force $\mathbf{F}_{\bf a}$ on $\bf{a}$ due to $\bf{b}$ is the negative of
1019 < the force  $\mathbf{F}_{\bf b}$ on $\bf{b}$ due to $\bf{a}$.  For a simple charge-charge
1020 < interaction, these forces will point along the $\pm \hat{r}$ directions, where
1021 < $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $.  Thus
1017 > The force on object $\bf{a}$, $\mathbf{F}_{\bf a}$, due to object
1018 > $\bf{b}$ is the negative of the force on $\bf{b}$ due to $\bf{a}$. For
1019 > a simple charge-charge interaction, these forces will point along the
1020 > $\pm \hat{r}$ directions, where $\mathbf{r}=\mathbf{r}_b -
1021 > \mathbf{r}_a $.  Thus
1022   %
1023   \begin{equation}
1024   F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}}{\partial r}
# Line 922 | Line 1026 | The concept of obtaining a force from an energy by tak
1026   = - \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}} {\partial r}  .
1027   \end{equation}
1028   %
1029 < The concept of obtaining a force from an energy by taking a gradient is the same for
1030 < higher-order multipole interactions, the trick is to make sure that all
1031 < $r$-dependent derivatives are considered.
928 < As is pointed out by Allen and Germano, this is straightforward if the
929 < interaction energies are written recognizing explicit
930 < $\hat{r}$ and body axes ($\hat{a}_m$, $\hat{b}_n$) dependences:
1029 > We list below the force equations written in terms of lab-frame
1030 > coordinates.  The radial functions used in the two methods are listed
1031 > in Table \ref{tab:tableFORCE}
1032   %
1033 < \begin{equation}
933 < U(r,\{\hat{a}_m \cdot \hat{r} \},
934 < \{\hat{b}_n\cdot \hat{r} \}
935 < \{\hat{a}_m \cdot \hat{b}_n \}) .
936 < \label{ugeneral}
937 < \end{equation}
1033 > %SPACE COORDINATES FORCE EQUATIONS
1034   %
939 Then,
940 %
941 \begin{equation}
942 \mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} =  \frac{\partial U}{\partial \mathbf{r}}
943 = \frac{\partial U}{\partial r} \hat{r}
944 + \sum_m \left[
945 \frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})}
946 \frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}}
947 + \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})}
948 \frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}}
949 \right] \label{forceequation}.
950 \end{equation}
951 %
952 Note our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ is opposite
953 that of Allen and Germano.  In simplifying the algebra, we also use:
954 %
955 \begin{eqnarray}
956 \frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}}
957 = \frac{1}{r} \left(  \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r}
958 \right) \\
959 \frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}}
960 = \frac{1}{r} \left(  \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r}
961 \right) .
962 \end{eqnarray}
963 %
964 We list below the force equations written in terms  of space coordinates.  The
965 radial functions used in the two methods are listed in Table II.
966 %
967 %SPACE COORDINATES FORCE EQUTIONS
968 %
1035   % **************************************************************************
1036   % f ca cb
1037   %
1038 < \begin{equation}
1039 < \mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =
1040 < \frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r}
975 < \end{equation}
1038 > \begin{align}
1039 > \mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =&
1040 > C_{\bf a} C_{\bf b}  w_a(r) \hat{r} \\
1041   %
1042   %
1043   %
1044 < \begin{equation}
1045 < \mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =
981 < \frac{C_{\bf a}}{4\pi \epsilon_0} \Bigl[
1044 > \mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =&
1045 > C_{\bf a} \Bigl[
1046   \left( \hat{r} \cdot \mathbf{D}_{\mathbf{b}} \right)
1047   w_b(r) \hat{r}  
1048 < + \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr]
985 < \end{equation}
1048 > + \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr] \\
1049   %
1050   %
1051   %
1052 < \begin{equation}
1053 < \mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =
991 < \frac{C_{\bf a }}{4\pi \epsilon_0} \Bigr[
1052 > \mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =&
1053 > C_{\bf a } \Bigr[
1054   \text{Tr}\mathbf{Q}_{\bf b} w_d(r) \hat{r}
1055   + 2  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} w_e(r)
1056 < + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
1057 < \end{equation}
1056 > + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}
1057 > \right) w_f(r) \hat{r} \Bigr] \\
1058   %
1059   %
1060   %
1061 < \begin{equation}
1062 < \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1063 < -\frac{C_{\bf{b}}}{4\pi \epsilon_0} \Bigl[
1064 < \left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) w_b(r) \hat{r}
1065 < + \mathbf{D}_{\mathbf{a}} w_c(r) \Bigr]
1066 < \end{equation}
1061 > % \begin{equation}
1062 > % \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1063 > % -C_{\bf{b}} \Bigl[
1064 > % \left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) w_b(r) \hat{r}
1065 > % + \mathbf{D}_{\mathbf{a}} w_c(r) \Bigr]
1066 > % \end{equation}
1067   %
1068   %
1069   %
1070 < \begin{equation}
1071 < \mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} =
1010 < \frac{1}{4\pi \epsilon_0} \Bigl[
1070 > \begin{split}
1071 > \mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} =&
1072   - \mathbf{D}_{\mathbf {a}} \cdot  \mathbf{D}_{\mathbf{b}} w_d(r) \hat{r}
1073   + \left( \mathbf{D}_{\mathbf {a}}
1074   \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)
1075 < + \mathbf{D}_{\mathbf {b}} \left( \mathbf{D}_{\mathbf{a}}  \cdot \hat{r} \right) \right) w_e(r)
1075 > + \mathbf{D}_{\mathbf {b}} \left( \mathbf{D}_{\mathbf{a}}  \cdot \hat{r} \right) \right) w_e(r)\\
1076   % 2
1077 < - \left( \hat{r} \cdot \mathbf{D}_{\mathbf {a}} \right)
1078 < \left( \hat{r} \cdot \mathbf{D}_{\mathbf {b}} \right) w_f(r) \hat{r} \Bigr]
1079 < \end{equation}
1077 > & - \left( \hat{r} \cdot \mathbf{D}_{\mathbf {a}} \right)
1078 > \left( \hat{r} \cdot \mathbf{D}_{\mathbf {b}} \right) w_f(r) \hat{r}
1079 > \end{split}\\
1080   %
1081   %
1082   %
1022 \begin{equation}
1083   \begin{split}
1084 < \mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1025 < & - \frac{1}{4\pi \epsilon_0} \Bigl[
1084 > \mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =& - \Bigl[
1085   \text{Tr}\mathbf{Q}_{\mathbf{b}} \mathbf{ D}_{\mathbf{a}}
1086   +2 \mathbf{D}_{\mathbf{a}} \cdot
1087   \mathbf{Q}_{\mathbf{b}} \Bigr] w_g(r)
1088 < - \frac{1}{4\pi \epsilon_0} \Bigl[
1088 > - \Bigl[
1089   \text{Tr}\mathbf{Q}_{\mathbf{b}}
1090   \left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right)
1091   +2 ( \mathbf{D}_{\mathbf{a}} \cdot
1092   \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
1093   % 3
1094 < & - \frac{1}{4\pi \epsilon_0} \Bigl[\mathbf{ D}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
1094 > & - \Bigl[\mathbf{ D}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
1095   +2  (\hat{r} \cdot \mathbf{D}_{\mathbf{a}} ) (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} )  \Bigr]
1096   w_i(r)
1097   % 4
1098 < -\frac{1}{4\pi \epsilon_0}
1098 > -
1099   (\hat{r} \cdot \mathbf{D}_{\mathbf{a}} )
1100 < (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r}  
1042 < \end{split}
1043 < \end{equation}
1100 > (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r} \end{split} \\
1101   %
1102   %
1103 < \begin{equation}
1104 < \mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} =
1105 < \frac{C_{\bf b }}{4\pi \epsilon_0} \Bigr[
1106 < \text{Tr}\mathbf{Q}_{\bf a} w_d(r) \hat{r}
1107 < + 2  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} w_e(r)
1108 < + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
1109 < \end{equation}
1103 > % \begin{equation}
1104 > % \mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} =
1105 > % \frac{C_{\bf b }}{4\pi \epsilon_0} \Bigr[
1106 > % \text{Tr}\mathbf{Q}_{\bf a} w_d(r) \hat{r}
1107 > % + 2  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} w_e(r)
1108 > %  + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
1109 > % \end{equation}
1110 > % %
1111 > % \begin{equation}
1112 > % \begin{split}
1113 > % \mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1114 > % &\frac{1}{4\pi \epsilon_0} \Bigl[
1115 > % \text{Tr}\mathbf{Q}_{\mathbf{a}} \mathbf{D}_{\mathbf{b}}
1116 > % +2 \mathbf{D}_{\mathbf{b}} \cdot \mathbf{Q}_{\mathbf{a}}  \Bigr] w_g(r)
1117 > % % 2
1118 > % + \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
1119 > % (\hat{r} \cdot  \mathbf{D}_{\mathbf{b}})
1120 > % +2 (\mathbf{D}_{\mathbf{b}} \cdot
1121 > % \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
1122 > % % 3
1123 > % &+ \frac{1}{4\pi \epsilon_0} \Bigl[ \mathbf{D}_{\mathbf{b}}
1124 > % (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1125 > % +2  (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1126 > % (\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} ) \Bigr]   w_i(r)
1127 > % % 4
1128 > % +\frac{1}{4\pi \epsilon_0}
1129 > % (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1130 > % (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}}  \cdot \hat{r}) w_j(r) \hat{r}
1131 > % \end{split}
1132 > % \end{equation}
1133   %
1054 \begin{equation}
1055 \begin{split}
1056 \mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1057 &\frac{1}{4\pi \epsilon_0} \Bigl[
1058 \text{Tr}\mathbf{Q}_{\mathbf{a}} \mathbf{D}_{\mathbf{b}}
1059 +2 \mathbf{D}_{\mathbf{b}} \cdot \mathbf{Q}_{\mathbf{a}}  \Bigr] w_g(r)
1060 % 2
1061 + \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
1062 (\hat{r} \cdot  \mathbf{D}_{\mathbf{b}})
1063 +2 (\mathbf{D}_{\mathbf{b}} \cdot
1064 \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
1065 % 3
1066 &+ \frac{1}{4\pi \epsilon_0} \Bigl[ \mathbf{D}_{\mathbf{b}}
1067 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1068 +2  (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1069 (\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} ) \Bigr]   w_i(r)
1070 % 4
1071 +\frac{1}{4\pi \epsilon_0}
1072 (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1073 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}}  \cdot \hat{r}) w_j(r) \hat{r}
1074 \end{split}
1075 \end{equation}
1134   %
1135   %
1078 %
1079 \begin{equation}
1136   \begin{split}
1137 < \mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1138 < +\frac{1}{4\pi \epsilon_0} \Bigl[
1139 < \text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} \hat{r}
1140 < + 2 \text{Tr} ( \mathbf{Q}_{\mathbf{a}} \cdot  \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_k(r) \hat{r} \\
1137 > \mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =&
1138 > \Bigl[
1139 > \text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}}
1140 > + 2  \mathbf{Q}_{\mathbf{a}} :  \mathbf{Q}_{\mathbf{b}} \Bigr] w_k(r) \hat{r} \\
1141   % 2
1142 < +\frac{1}{4\pi \epsilon_0} \Bigl[
1142 > &+ \Bigl[
1143   2\text{Tr}\mathbf{Q}_{\mathbf{b}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )  
1144   + 2\text{Tr}\mathbf{Q}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} )
1145   % 3
1146   +4 (\mathbf{Q}_{\mathbf{a}}  \cdot  \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})  
1147   +  4(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}}) \Bigr] w_n(r) \\
1148   % 4
1149 < + \frac{1}{4\pi \epsilon_0} \Bigl[
1149 > &+  \Bigl[
1150   \text{Tr}\mathbf{Q}_{\mathbf{a}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1151   + \text{Tr}\mathbf{Q}_{\mathbf{b}}
1152   (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})  
# Line 1098 | Line 1154 | + \frac{1}{4\pi \epsilon_0} \Bigl[
1154   +4 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot  
1155   \mathbf{Q}_{\mathbf{b}}   \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
1156   %
1157 < + \frac{1}{4\pi \epsilon_0} \Bigl[
1157 > &+ \Bigl[
1158   + 2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )
1159   (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1160   %6
1161   +2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \hat{r})
1162   (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_o(r) \\
1163   %  7
1164 < + \frac{1}{4\pi \epsilon_0}
1164 > &+
1165   (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
1166 < (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) w_m(r) \hat{r}
1167 < \end{split}
1168 < \end{equation}
1166 > (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) w_m(r) \hat{r} \end{split}
1167 > \end{align}
1168 > Note that the forces for higher multipoles on site $\mathbf{a}$
1169 > interacting with those of lower order on site $\mathbf{b}$ can be
1170 > obtained by swapping indices in the expressions above.
1171 >
1172   %
1173 + % Torques SECTION -----------------------------------------------------------------------------------------
1174   %
1115 % TORQUES SECTION -----------------------------------------------------------------------------------------
1116 %
1175   \subsection{Torques}
1176  
1119 Following again Allen and Germano, when energies are written in the form
1120 of Eq.~({\ref{ugeneral}), then torques can be expressed as:
1177   %
1178 < \begin{eqnarray}
1179 < \mathbf{\tau}_{\bf a} =
1124 < \sum_m
1125 < \frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})}
1126 < ( \hat{r} \times \hat{a}_m )
1127 < -\sum_{mn}
1128 < \frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)}
1129 < (\hat{a}_m \times \hat{b}_n) \\
1178 > The torques for both the Taylor-Shifted as well as Gradient-Shifted
1179 > methods are given in space-frame coordinates:
1180   %
1131 \mathbf{\tau}_{\bf b} =
1132 \sum_m
1133 \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})}
1134 ( \hat{r} \times \hat{b}_m)
1135 +\sum_{mn}
1136 \frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)}
1137 (\hat{a}_m \times \hat{b}_n) .
1138 \end{eqnarray}
1181   %
1182 + \begin{align}
1183 + \mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =&
1184 + C_{\bf a}  (\hat{r} \times  \mathbf{D}_{\mathbf{b}}) v_{11}(r) \\
1185   %
1141 Here we list the torque equations written in terms of space coordinates.
1186   %
1187   %
1188 + \mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =&
1189 + 2C_{\bf a}
1190 + \hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{22}(r) \\
1191   %
1145 \begin{equation}
1146 \mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =
1147 \frac{C_{\bf a}}{4\pi \epsilon_0}  (\hat{r} \times  \mathbf{D}_{\mathbf{b}}) v_{11}(r)
1148 \end{equation}
1192   %
1193   %
1194 + % \begin{equation}
1195 + % \mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =  
1196 + % -\frac{C_{\bf b}}{4\pi \epsilon_0}  
1197 + % (\hat{r} \times \mathbf{D}_{\mathbf{a}})  v_{11}(r)
1198 + % \end{equation}
1199   %
1152 \begin{equation}
1153 \mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =
1154 \frac{2C_{\bf a}}{4\pi \epsilon_0}
1155 \hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{22}(r)
1156 \end{equation}
1200   %
1201   %
1202 < %
1203 < \begin{equation}
1161 < \mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =  
1162 < -\frac{C_{\bf b}}{4\pi \epsilon_0}  
1163 < (\hat{r} \times \mathbf{D}_{\mathbf{a}})  v_{11}(r)
1164 < \end{equation}
1165 < %
1166 < %
1167 < %
1168 < \begin{equation}
1169 < \mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} =
1170 < \frac{1}{4\pi \epsilon_0}  \mathbf{D}_{\mathbf {a}}  \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1202 > \mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} =&
1203 > \mathbf{D}_{\mathbf {a}}  \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1204   % 2
1205 < -\frac{1}{4\pi \epsilon_0}
1205 > -
1206   (\hat{r} \times \mathbf{D}_{\mathbf {a}} )
1207 < (\hat{r} \cdot \mathbf{D}_{\mathbf {b}} )  v_{22}(r)
1175 < \end{equation}
1207 > (\hat{r} \cdot \mathbf{D}_{\mathbf {b}} )  v_{22}(r)\\
1208   %
1209   %
1210   %
1211 < \begin{equation}
1212 < \mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} =
1213 < -\frac{1}{4\pi \epsilon_0} \mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1214 < % 2
1215 < +\frac{1}{4\pi \epsilon_0}
1216 < (\hat{r} \cdot \mathbf{D}_{\mathbf {a}} )
1217 < (\hat{r} \times \mathbf{D}_{\mathbf {b}} ) v_{22}(r)
1218 < \end{equation}
1211 > % \begin{equation}
1212 > % \mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} =
1213 > % -\frac{1}{4\pi \epsilon_0} \mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1214 > % % 2
1215 > % +\frac{1}{4\pi \epsilon_0}
1216 > % (\hat{r} \cdot \mathbf{D}_{\mathbf {a}} )
1217 > % (\hat{r} \times \mathbf{D}_{\mathbf {b}} ) v_{22}(r)
1218 > % \end{equation}
1219   %
1220   %
1221   %
1222 < \begin{equation}
1223 < \mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1192 < \frac{1}{4\pi \epsilon_0} \Bigl[
1222 > \mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} =&
1223 > \Bigl[
1224   -\text{Tr}\mathbf{Q}_{\mathbf{b}}
1225   (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
1226   +2 \mathbf{D}_{\mathbf{a}}  \times
1227   (\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1228   \Bigr] v_{31}(r)
1229   % 3
1230 < -\frac{1}{4\pi \epsilon_0}
1231 < \ (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
1201 < (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)
1202 < \end{equation}
1230 > - (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
1231 > (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)\\
1232   %
1233   %
1234   %
1235 < \begin{equation}
1236 < \mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =
1208 < \frac{1}{4\pi \epsilon_0} \Bigl[
1235 > \mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =&
1236 > \Bigl[
1237   +2 ( \mathbf{D}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} ) \times
1238   \hat{r}
1239   -2 \mathbf{D}_{\mathbf{a}}  \times
1240   (\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1241   \Bigr] v_{31}(r)
1242   % 2
1243 < +\frac{2}{4\pi \epsilon_0}
1243 > +
1244   (\hat{r} \cdot \mathbf{D}_{\mathbf{a}})
1245 < (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)
1218 < \end{equation}
1245 > (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)\\
1246   %
1247   %
1248   %
1249 < \begin{equation}
1250 < \mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1251 < \frac{1}{4\pi \epsilon_0} \Bigl[
1252 < -2 (\mathbf{D}_{\mathbf{b}}  \cdot \mathbf{Q}_{\mathbf{a}} ) \times \hat{r}
1253 < +2 \mathbf{D}_{\mathbf{b}}  \times
1254 < (\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
1255 < \Bigr] v_{31}(r)
1256 < % 3
1257 < - \frac{2}{4\pi \epsilon_0}
1258 < (\hat{r} \cdot \mathbf{D}_{\mathbf{b}} )
1259 < (\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}}) \times \hat{r} v_{32}(r)
1260 < \end{equation}
1249 > % \begin{equation}
1250 > % \mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1251 > % \frac{1}{4\pi \epsilon_0} \Bigl[
1252 > % -2 (\mathbf{D}_{\mathbf{b}}  \cdot \mathbf{Q}_{\mathbf{a}} ) \times \hat{r}
1253 > % +2 \mathbf{D}_{\mathbf{b}}  \times
1254 > % (\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
1255 > % \Bigr] v_{31}(r)
1256 > % % 3
1257 > % - \frac{2}{4\pi \epsilon_0}
1258 > % (\hat{r} \cdot \mathbf{D}_{\mathbf{b}} )
1259 > % (\hat{r} \cdot  \mathbf
1260 > % {Q}_{{\mathbf a}}) \times \hat{r} v_{32}(r)
1261 > % \end{equation}
1262   %
1263   %
1264   %
1265 < \begin{equation}
1266 < \mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} =
1267 < \frac{1}{4\pi \epsilon_0} \Bigl[
1268 < \text{Tr}\mathbf{Q}_{\mathbf{a}}
1269 < (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1270 < +2 \mathbf{D}_{\mathbf{b}}  \times
1271 < ( \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r)
1272 < % 2
1273 < +\frac{1}{4\pi \epsilon_0}
1274 < (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1275 < (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) v_{32}(r)
1276 < \end{equation}
1265 > % \begin{equation}
1266 > % \mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} =
1267 > % \frac{1}{4\pi \epsilon_0} \Bigl[
1268 > % \text{Tr}\mathbf{Q}_{\mathbf{a}}
1269 > % (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1270 > % +2 \mathbf{D}_{\mathbf{b}}  \times
1271 > % ( \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r)
1272 > % % 2
1273 > % +\frac{1}{4\pi \epsilon_0}
1274 > % (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1275 > % (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) v_{32}(r)
1276 > % \end{equation}
1277   %
1278   %
1279   %
1252 \begin{equation}
1280   \begin{split}
1281 < \mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1282 < &-\frac{4}{4\pi \epsilon_0}
1281 > \mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =&
1282 > -4
1283   \mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}}
1284   v_{41}(r) \\
1285   % 2
1286 < &+ \frac{1}{4\pi \epsilon_0}
1286 > &+
1287   \Bigl[-2\text{Tr}\mathbf{Q}_{\mathbf{b}}
1288   (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times \hat{r}
1289   +4 \hat{r} \times
# Line 1265 | Line 1292 | -4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} )\times
1292   -4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} )\times
1293   ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} ) \Bigr] v_{42}(r) \\
1294   % 4
1295 < &+ \frac{2}{4\pi \epsilon_0}
1295 > &+ 2
1296   \hat{r} \times ( \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1297 < (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
1271 < \end{split}
1272 < \end{equation}
1297 > (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r) \end{split}\\
1298   %
1299   %
1300   %
1276 \begin{equation}
1301   \begin{split}
1302   \mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} =  
1303 < &\frac{4}{4\pi \epsilon_0}
1303 > &4
1304   \mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}} v_{41}(r) \\
1305   % 2
1306 < &+ \frac{1}{4\pi \epsilon_0} \Bigl[- 2\text{Tr}\mathbf{Q}_{\mathbf{a}}
1306 > &+  \Bigl[- 2\text{Tr}\mathbf{Q}_{\mathbf{a}}
1307   (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} ) \times \hat{r}
1308   -4  (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot
1309   \mathbf{Q}_{{\mathbf b}} ) \times
# Line 1288 | Line 1312 | +4 ( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times
1312   ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
1313   \Bigr] v_{42}(r) \\
1314   % 4
1315 < &+ \frac{2}{4\pi \epsilon_0}
1315 > &+2
1316   (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1317 < \hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
1318 < \end{split}
1295 < \end{equation}
1317 > \hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)\end{split}
1318 > \end{align}
1319   %
1320 + Here, we have defined the matrix cross product in an identical form
1321 + as in Ref. \onlinecite{Smith98}:
1322 + \begin{equation}
1323 + \left[\mathbf{A} \times \mathbf{B}\right]_\alpha = \sum_\beta
1324 + \left[\mathbf{A}_{\alpha+1,\beta} \mathbf{B}_{\alpha+2,\beta}
1325 +  -\mathbf{A}_{\alpha+2,\beta} \mathbf{B}_{\alpha+2,\beta}
1326 + \right]
1327 + \end{equation}
1328 + where $\alpha+1$ and $\alpha+2$ are regarded as cyclic
1329 + permuations of the matrix indices.
1330  
1331 + All of the radial functions required for torques are identical with
1332 + the radial functions previously computed for the interaction energies.
1333 + These are tabulated for both shifted force methods in table
1334 + \ref{tab:tableenergy}.  The torques for higher multipoles on site
1335 + $\mathbf{a}$ interacting with those of lower order on site
1336 + $\mathbf{b}$ can be obtained by swapping indices in the expressions
1337 + above.
1338 +
1339 + \section{Related real-space methods}
1340 + One can also formulate a shifted potential,
1341 + \begin{equation}
1342 + U^{\text{SP}} = U(\mathbf{r},\hat{\mathbf{a}}, \hat{\mathbf{b}}) -
1343 + U(\mathbf{r}_c, \hat{\mathbf{a}}, \hat{\mathbf{b}}),
1344 + \label{eq:SP}
1345 + \end{equation}
1346 + obtained by projecting the image multipole onto the surface of the
1347 + cutoff sphere.  The shifted potential (SP) can be thought of as a
1348 + simple extension to the original Wolf method.  The energies and
1349 + torques for the SP can be easily obtained by zeroing out the $(r-r_c)$
1350 + terms in the final column of table \ref{tab:tableenergy}.  SP forces
1351 + (which retain discontinuities at the cutoff sphere) can be obtained by
1352 + eliminating all functions that depend on $r_c$ in the last column of
1353 + table \ref{tab:tableFORCE}.  The self-energy contributions to the SP
1354 + potential are identical to both the GSF and TSF methods.
1355 +
1356   \section{Comparison to known multipolar energies}
1357  
1358   To understand how these new real-space multipole methods behave in
1359   computer simulations, it is vital to test against established methods
1360   for computing electrostatic interactions in periodic systems, and to
1361   evaluate the size and sources of any errors that arise from the
1362 < real-space cutoffs. In this paper we test Taylor-shifted and
1363 < Gradient-shifted electrostatics against analytical methods for
1364 < computing the energies of ordered multipolar arrays.  In the following
1365 < paper, we test the new methods against the multipolar Ewald sum for
1366 < computing the energies, forces and torques for a wide range of typical
1367 < condensed-phase (disordered) systems.
1362 > real-space cutoffs. In this paper we test both TSF and GSF
1363 > electrostatics against analytical methods for computing the energies
1364 > of ordered multipolar arrays.  In the following paper, we test the new
1365 > methods against the multipolar Ewald sum for computing the energies,
1366 > forces and torques for a wide range of typical condensed-phase
1367 > (disordered) systems.
1368  
1369   Because long-range electrostatic effects can be significant in
1370   crystalline materials, ordered multipolar arrays present one of the
# Line 1315 | Line 1373 | Tisza\cite{LT,LT2} who tabulated energy constants for
1373   who computed the energies of ordered dipole arrays of zero
1374   magnetization and obtained a number of these constants.\cite{Sauer}
1375   This theory was developed more completely by Luttinger and
1376 < Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays and
1377 < other periodic structures.  We have repeated the Luttinger \& Tisza
1320 < series summations to much higher order and obtained the following
1321 < energy constants (converged to one part in $10^9$):
1322 < \begin{table*}
1323 < \centering{
1324 <  \caption{Luttinger \& Tisza arrays and their associated
1325 <    energy constants. Type "A" arrays have nearest neighbor strings of
1326 <    antiparallel dipoles.  Type "B" arrays have nearest neighbor
1327 <    strings of antiparallel dipoles if the dipoles are contained in a
1328 <    plane perpendicular to the dipole direction that passes through
1329 <    the dipole.}
1330 < }
1331 < \label{tab:LT}
1332 < \begin{ruledtabular}
1333 < \begin{tabular}{cccc}
1334 < Array Type &  Lattice &   Dipole Direction &    Energy constants \\ \hline
1335 <   A     &      SC       &      001         &      -2.676788684 \\
1336 <   A     &      BCC      &      001         &       0 \\
1337 <   A     &      BCC      &      111         &      -1.770078733 \\
1338 <   A     &      FCC      &      001         &       2.166932835 \\
1339 <   A     &      FCC      &      011         &      -1.083466417 \\
1376 > Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays
1377 > and other periodic structures.  
1378  
1379 <   *     &      BCC      &    minimum       &      -1.985920929 \\
1380 <
1381 <   B     &      SC        &     001         &      -2.676788684 \\
1382 <   B     &      BCC       &     001         &      -1.338394342 \\
1383 <   B     &      BCC       &     111         &     -1.770078733 \\
1384 <   B     &      FCC       &     001         &      -1.083466417 \\
1385 <   B     &      FCC       &     011         &      -1.807573634
1348 < \end{tabular}
1349 < \end{ruledtabular}
1350 < \end{table*}
1351 <
1352 < In addition to the A and B arrays, there is an additional minimum
1379 > To test the new electrostatic methods, we have constructed very large,
1380 > $N=$ 16,000~(bcc) arrays of dipoles in the orientations described in
1381 > Ref. \onlinecite{LT}.  These structures include ``A'' lattices with
1382 > nearest neighbor chains of antiparallel dipoles, as well as ``B''
1383 > lattices with nearest neighbor strings of antiparallel dipoles if the
1384 > dipoles are contained in a plane perpendicular to the dipole direction
1385 > that passes through the dipole.  We have also studied the minimum
1386   energy structure for the BCC lattice that was found by Luttinger \&
1387 < Tisza.  The total electrostatic energy for an array is given by:
1387 > Tisza.  The total electrostatic energy for any of the arrays is given
1388 > by:
1389   \begin{equation}
1390    E = C N^2 \mu^2
1391   \end{equation}
1392 < where $C$ is the energy constant given above, $N$ is the number of
1393 < dipoles per unit volume, and $\mu$ is the strength of the dipole.
1392 > where $C$ is the energy constant (equivalent to the Madelung
1393 > constant), $N$ is the number of dipoles per unit volume, and $\mu$ is
1394 > the strength of the dipole. Energy constants (converged to 1 part in
1395 > $10^9$) are given in the supplemental information.
1396  
1397 < {\it Quadrupolar} analogues to the Madelung constants were first worked out by Nagai and Nakamura who
1398 < computed the energies of selected quadrupole arrays based on
1399 < extensions to the Luttinger and Tisza
1400 < approach.\cite{Nagai01081960,Nagai01091963}  We have compared the
1401 < energy constants for the lowest energy configurations for linear
1402 < quadrupoles shown in table \ref{tab:NNQ}
1397 > For the purposes of testing the energy expressions and the
1398 > self-neutralization schemes, the primary quantity of interest is the
1399 > analytic energy constant for the perfect arrays.  Convergence to these
1400 > constants are shown as a function of both the cutoff radius, $r_c$,
1401 > and the damping parameter, $\alpha$ in Figs.
1402 > \ref{fig:energyConstVsCutoff} and XXX. We have simultaneously tested a
1403 > hard cutoff (where the kernel is simply truncated at the cutoff
1404 > radius), as well as a shifted potential (SP) form which includes a
1405 > potential-shifting and self-interaction term, but does not shift the
1406 > forces and torques smoothly at the cutoff radius.  The SP method is
1407 > essentially an extension of the original Wolf method for multipoles.
1408  
1409 < \begin{table*}
1410 < \centering{
1411 <  \caption{Nagai and Nakamura Quadurpolar arrays}}
1412 < \label{tab:NNQ}
1413 < \begin{ruledtabular}
1414 < \begin{tabular}{ccc}
1415 < Lattice &   Quadrupole Direction &    Energy constants \\ \hline
1416 <  SC       &      111         &      -8.3 \\
1417 <  BCC      &      011         &      -21.7 \\
1418 <  FCC      &      111         &      -80.5
1419 < \end{tabular}
1420 < \end{ruledtabular}
1380 < \end{table*}
1409 > \begin{figure}[!htbp]
1410 > \includegraphics[width=4.5in]{energyConstVsCutoff}
1411 > \caption{Convergence to the analytic energy constants as a function of
1412 >  cutoff radius (normalized by the lattice constant) for the different
1413 >  real-space methods. The two crystals shown here are the ``B'' array
1414 >  for bcc crystals with the dipoles along the 001 direction (upper),
1415 >  as well as the minimum energy bcc lattice (lower).  The analytic
1416 >  energy constants are shown as a grey dashed line.  The left panel
1417 >  shows results for the undamped kernel ($1/r$), while the damped
1418 >  error function kernel, $B_0(r)$ was used in the right panel. }
1419 > \label{fig:energyConstVsCutoff}
1420 > \end{figure}
1421  
1422 + The Hard cutoff exhibits oscillations around the analytic energy
1423 + constants, and converges to incorrect energies when the complementary
1424 + error function damping kernel is used.  The shifted potential (SP) and
1425 + gradient-shifted force (GSF) approximations converge to the correct
1426 + energy smoothly by $r_c / 6 a$ even for the undamped case.  This
1427 + indicates that the correction provided by the self term is required
1428 + for obtaining accurate energies.  The Taylor-shifted force (TSF)
1429 + approximation appears to perturb the potential too much inside the
1430 + cutoff region to provide accurate measures of the energy constants.
1431 +
1432 + {\it Quadrupolar} analogues to the Madelung constants were first
1433 + worked out by Nagai and Nakamura who computed the energies of selected
1434 + quadrupole arrays based on extensions to the Luttinger and Tisza
1435 + approach.\cite{Nagai01081960,Nagai01091963} We have compared the
1436 + energy constants for the lowest energy configurations for linear
1437 + quadrupoles.  
1438 +
1439   In analogy to the dipolar arrays, the total electrostatic energy for
1440   the quadrupolar arrays is:
1441   \begin{equation}
1442   E = C \frac{3}{4} N^2 Q^2
1443   \end{equation}
1444 < where $Q$ is the quadrupole moment.
1444 > where $Q$ is the quadrupole moment.  The lowest energy
1445  
1446 + \section{Conclusion}
1447 + We have presented two efficient real-space methods for computing the
1448 + interactions between point multipoles.  These methods have the benefit
1449 + of smoothly truncating the energies, forces, and torques at the cutoff
1450 + radius, making them attractive for both molecular dynamics (MD) and
1451 + Monte Carlo (MC) simulations.  We find that the Gradient-Shifted Force
1452 + (GSF) and the Shifted-Potential (SP) methods converge rapidly to the
1453 + correct lattice energies for ordered dipolar and quadrupolar arrays,
1454 + while the Taylor-Shifted Force (TSF) is too severe an approximation to
1455 + provide accurate convergence to lattice energies.  
1456  
1457 + In most cases, GSF can obtain nearly quantitative agreement with the
1458 + lattice energy constants with reasonably small cutoff radii.  The only
1459 + exception we have observed is for crystals which exhibit a bulk
1460 + macroscopic dipole moment (e.g. Luttinger \& Tisza's $Z_1$ lattice).
1461 + In this particular case, the multipole neutralization scheme can
1462 + interfere with the correct computation of the energies.  We note that
1463 + the energies for these arrangements are typically much larger than for
1464 + crystals with net-zero moments, so this is not expected to be an issue
1465 + in most simulations.
1466  
1467 + In large systems, these new methods can be made to scale approximately
1468 + linearly with system size, and detailed comparisons with the Ewald sum
1469 + for a wide range of chemical environments follows in the second paper.
1470  
1392
1393
1394
1395
1471   \begin{acknowledgments}
1472 < Support for this project was provided by the National Science
1473 < Foundation under grant CHE-0848243. Computational time was provided by
1474 < the Center for Research Computing (CRC) at the University of Notre
1475 < Dame.
1472 >  JDG acknowledges helpful discussions with Christopher
1473 >  Fennell. Support for this project was provided by the National
1474 >  Science Foundation under grant CHE-0848243. Computational time was
1475 >  provided by the Center for Research Computing (CRC) at the
1476 >  University of Notre Dame.
1477   \end{acknowledgments}
1478  
1479 + \newpage
1480   \appendix
1481  
1482 < \section{Smith's $B_l(r)$ functions for smeared-charge distributions}
1483 <
1484 < The following summarizes Smith's $B_l(r)$ functions and
1485 < includes formulas given in his appendix.  
1486 <
1410 < The first function $B_0(r)$ is defined by
1482 > \section{Smith's $B_l(r)$ functions for damped-charge distributions}
1483 > \label{SmithFunc}
1484 > The following summarizes Smith's $B_l(r)$ functions and includes
1485 > formulas given in his appendix.\cite{Smith98} The first function
1486 > $B_0(r)$ is defined by
1487   %
1488   \begin{equation}
1489   B_0(r)=\frac{\text{erfc}(\alpha r)}{r} = \frac{2}{\sqrt{\pi}r}=
# Line 1421 | Line 1497 | and can be rewritten in terms of a function $B_1(r)$:
1497   -\frac{2\alpha}{r\sqrt{\pi}}\text{e}^{-{\alpha}^2r^2}
1498   \end{equation}
1499   %
1500 < and can be rewritten in terms of a function $B_1(r)$:
1500 > which can be used to define a function $B_1(r)$:
1501   %
1502   \begin{equation}
1503   B_1(r)=-\frac{1}{r}\frac{dB_0(r)}{dr}
1504   \end{equation}
1505   %
1506 < In general,
1506 > In general, the recurrence relation,
1507   \begin{equation}
1508   B_l(r)=-\frac{1}{r}\frac{dB_{l-1}(r)}{dr}
1509   = \frac{1}{r^2} \left[ (2l-1)B_{l-1}(r) + \frac {(2\alpha^2)^l}{\alpha \sqrt{\pi}}
1510   \text{e}^{-{\alpha}^2r^2}
1511 < \right] .
1511 > \right] ,
1512   \end{equation}
1513 + is very useful for building up higher derivatives.  Using these
1514 + formulas, we find:
1515   %
1516 < Using these formulas, we find
1516 > \begin{align}
1517 > \frac{dB_0}{dr}=&-rB_1(r) \\
1518 > \frac{d^2B_0}{dr^2}=&    - B_1(r) + r^2 B_2(r) \\
1519 > \frac{d^3B_0}{dr^3}=&   3 r B_2(r) - r^3 B_3(r) \\
1520 > \frac{d^4B_0}{dr^4}=&   3 B_2(r) - 6 r^2 B_3(r) + r^4 B_4(r) \\
1521 > \frac{d^5B_0}{dr^5}=& - 15 r B_3(r) + 10 r^3 B_4(r) - r^5 B_5(r) .
1522 > \end{align}
1523   %
1524 < \begin{eqnarray}
1525 < \frac{dB_0}{dr}=-rB_1(r) \\
1442 < \frac{d^2B_0}{dr^2}=-B_1(r) + r^2B_2(r) \\
1443 < \frac{d^3B_0}{dr^3}=3rB_2(r) - r^3B_3(r) \\
1444 < \frac{d^4B_0}{dr^4}=3B_2(r) - 6r^2B_3(r)+r^4B_4(r) \\
1445 < \frac{d^5B_0}{dr^5}=-15rB_3(r) + 10r^3B_4(r) -r^5B_5(r) .
1446 < \end{eqnarray}
1524 > As noted by Smith, it is possible to approximate the $B_l(r)$
1525 > functions,
1526   %
1448 As noted by Smith,
1449 %
1527   \begin{equation}
1528   B_l(r)=\frac{(2l)!}{l!2^lr^{2l+1}} - \frac {(2\alpha^2)^{l+1}}{(2l+1)\alpha \sqrt{\pi}}
1529   +\text{O}(r) .
1530   \end{equation}
1531 + \newpage
1532 + \section{The $r$-dependent factors for TSF electrostatics}
1533  
1455 \section{Method 1, the $r$-dependent factors}
1456
1534   Using the shifted damped functions $f_n(r)$ defined by:
1535   %
1536   \begin{equation}
1537 < f_n(r)= B_0 \Big \lvert  _r -\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} B_0^{(m)} \Big \lvert  _{r_c}  ,
1537 > f_n(r)= B_0(r) -\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} B_0^{(m)}(r_c)  ,
1538   \end{equation}
1539   %
1540 < we first provide formulas for successive derivatives of this function.  (If there is
1541 < no damping, then $B_0(r)$ is replaced by $1/r$, as discussed in Section~\ref{damped???}.)  First, we find:
1540 > where the superscript $(m)$ denotes the $m^\mathrm{th}$ derivative. In
1541 > this Appendix, we provide formulas for successive derivatives of this
1542 > function.  (If there is no damping, then $B_0(r)$ is replaced by
1543 > $1/r$.)  First, we find:
1544   %
1545   \begin{equation}
1546   \frac{\partial f_n}{\partial r_\alpha}=\hat{r}_\alpha \frac{d f_n}{d r} .
1547   \end{equation}
1548   %
1549 < This formula clearly brings in derivatives of Smith's $B_0(r)$ function, motivating us to
1550 < define higher-order derivatives as  follows:
1549 > This formula clearly brings in derivatives of Smith's $B_0(r)$
1550 > function, and we define higher-order derivatives as follows:
1551   %
1552 < \begin{eqnarray}
1553 < g_n(r)= \frac{d f_n}{d r} =
1554 < B_0^{(1)} \Big \lvert  _r -\sum_{m=0}^{n} \frac {(r-r_c)^m}{m!} B_0^{(m+1)} \Big \lvert  _{r_c} \\
1555 < h_n(r)= \frac{d^2f_n}{d r^2} =
1556 < B_0^{(2)} \Big \lvert  _r -\sum_{m=0}^{n-1} \frac {(r-r_c)^m}{m!} B_0^{(m+2)} \Big \lvert  _{r_c} \\
1557 < s_n(r)= \frac{d^3f_n}{d r^3} =
1558 < B_0^{(3)} \Big \lvert  _r -\sum_{m=0}^{n-2} \frac {(r-r_c)^m}{m!} B_0^{(m+3)} \Big \lvert  _{r_c} \\
1559 < t_n(r)= \frac{d^4f_n}{d r^4} =
1560 < B_0^{(4)} \Big \lvert  _r -\sum_{m=0}^{n-3} \frac {(r-r_c)^m}{m!} B_0^{(m+4)} \Big \lvert  _{r_c} \\
1561 < u_n(r)= \frac{d^5f_n}{d r^5} =
1562 < B_0^{(5)} \Big \lvert  _r -\sum_{m=0}^{n-4} \frac {(r-r_c)^m}{m!} B_0^{(m+5)} \Big \lvert  _{r_c}  .
1563 < \end{eqnarray}
1552 > \begin{align}
1553 > g_n(r)=& \frac{d f_n}{d r} =
1554 > B_0^{(1)}(r) -\sum_{m=0}^{n} \frac {(r-r_c)^m}{m!} B_0^{(m+1)}(r_c) \\
1555 > h_n(r)=& \frac{d^2f_n}{d r^2} =
1556 > B_0^{(2)}(r) -\sum_{m=0}^{n-1} \frac {(r-r_c)^m}{m!} B_0^{(m+2)}(r_c) \\
1557 > s_n(r)=& \frac{d^3f_n}{d r^3} =
1558 > B_0^{(3)}(r) -\sum_{m=0}^{n-2} \frac {(r-r_c)^m}{m!} B_0^{(m+3)}(r_c) \\
1559 > t_n(r)=& \frac{d^4f_n}{d r^4} =
1560 > B_0^{(4)}(r) -\sum_{m=0}^{n-3} \frac {(r-r_c)^m}{m!} B_0^{(m+4)}(r_c) \\
1561 > u_n(r)=& \frac{d^5f_n}{d r^5} =
1562 > B_0^{(5)}(r) -\sum_{m=0}^{n-4} \frac {(r-r_c)^m}{m!} B_0^{(m+5)}(r_c)  .
1563 > \end{align}
1564   %
1565 < We note that the last function needed (for quadrupole-quadrupole) is
1565 > We note that the last function needed (for quadrupole-quadrupole interactions) is
1566   %
1567   \begin{equation}
1568 < u_4(r)=B_0^{(5)} \Big \lvert  _r - B_0^{(5)} \Big \lvert  _{r_c} .
1568 > u_4(r)=B_0^{(5)}(r) - B_0^{(5)}(r_c) .
1569   \end{equation}
1570 <
1571 < The functions $f_n(r)$ to $u_n(r)$ are recursively computed and stored for values of $r$
1572 < from $0$ to $r_c$.  The functions needed are listed schematically below:
1570 > % The functions
1571 > % needed are listed schematically below:
1572 > % %
1573 > % \begin{eqnarray}
1574 > % f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\
1575 > % g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\
1576 > % h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\
1577 > % s_2 \quad s_3 \quad &s_4 \nonumber \\
1578 > % t_3 \quad &t_4 \nonumber \\
1579 > % &u_4 \nonumber .
1580 > % \end{eqnarray}
1581 > The functions $f_n(r)$ to $u_n(r)$ can be computed recursively and
1582 > stored on a grid for values of $r$ from $0$ to $r_c$.  Using these
1583 > functions, we find
1584   %
1585 < \begin{eqnarray}
1586 < f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\
1587 < g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\
1588 < h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\
1589 < s_2 \quad s_3 \quad &s_4 \nonumber \\
1500 < t_3 \quad &t_4 \nonumber \\
1501 < &u_4 \nonumber .
1502 < \end{eqnarray}
1503 <
1504 < Using these functions, we find
1505 < %
1506 < \begin{equation}
1507 < \frac{\partial f_n}{\partial r_\alpha} =r_\alpha \frac {g_n}{r}
1508 < \end{equation}
1509 < %
1510 < \begin{equation}
1511 < \frac{\partial^2 f_n}{\partial r_\alpha \partial r_\beta} =\delta_{\alpha \beta}\frac {g_n}{r}
1512 < +r_\alpha r_\beta \left( -\frac{g_n}{r^3} +\frac{h_n}{r^2}\right)
1513 < \end{equation}
1514 < %
1515 < \begin{equation}
1516 < \frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta r_\gamma} =
1585 > \begin{align}
1586 > \frac{\partial f_n}{\partial r_\alpha} =&r_\alpha \frac {g_n}{r} \label{eq:b9}\\
1587 > \frac{\partial^2 f_n}{\partial r_\alpha \partial r_\beta} =&\delta_{\alpha \beta}\frac {g_n}{r}
1588 > +r_\alpha r_\beta \left( -\frac{g_n}{r^3} +\frac{h_n}{r^2}\right) \\
1589 > \frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta \partial r_\gamma} =&
1590   \left( \delta_{\alpha \beta} r_\gamma + \delta_{\alpha \gamma} r_\beta +
1591   \delta_{ \beta \gamma} r_\alpha \right)  
1592 < \left(  -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right)
1593 < + r_\alpha r_\beta r_\gamma
1594 < \left(  \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right)
1595 < \end{equation}
1596 < %
1524 < \begin{eqnarray}
1525 < \frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta r_\gamma r_\delta} =
1592 > \left(  -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right) \nonumber \\
1593 > & + r_\alpha r_\beta r_\gamma
1594 > \left(  \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \\
1595 > \frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta \partial
1596 >  r_\gamma \partial r_\delta} =&
1597   \left( \delta_{\alpha \beta} \delta_{\gamma \delta}
1598   + \delta_{\alpha \gamma} \delta_{\beta \delta}
1599   +\delta_{ \beta \gamma} \delta_{\alpha \delta} \right)
1600   \left( - \frac{g_n}{r^3} + \frac{h_n}{r^2} \right)  \nonumber \\
1601 < + \left( \delta_{\alpha \beta} r_\gamma r_\delta
1602 < + \text{5 perm}
1601 > &+ \left( \delta_{\alpha \beta} r_\gamma r_\delta
1602 > + \text{5 permutations}
1603   \right) \left( \frac{3 g_n}{r^5} - \frac{3h_n}{r^4} + \frac{s_n}{r^3}
1604   \right) \nonumber \\
1605 < + r_\alpha r_\beta r_\gamma r_\delta
1605 > &+ r_\alpha r_\beta r_\gamma r_\delta
1606   \left(  -\frac{15g_n}{r^7} + \frac{15h_n}{r^6} - \frac{6s_n}{r^5}
1607 < + \frac{t_n}{r^4} \right)
1537 < \end{eqnarray}
1538 < %
1539 < \begin{eqnarray}
1607 > + \frac{t_n}{r^4} \right)\\
1608   \frac{\partial^5 f_n}
1609 < {\partial r_\alpha \partial r_\beta r_\gamma r_\delta r_\epsilon} =
1609 > {\partial r_\alpha \partial r_\beta \partial r_\gamma \partial
1610 >  r_\delta \partial r_\epsilon} =&
1611   \left( \delta_{\alpha \beta} \delta_{\gamma \delta} r_\epsilon
1612 < + \text{14 perm} \right)
1612 > + \text{14 permutations} \right)
1613   \left(  \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \nonumber \\
1614 < + \left( \delta_{\alpha \beta} r_\gamma r_\delta r_\epsilon
1615 < + \text{9 perm}
1614 > &+ \left( \delta_{\alpha \beta} r_\gamma r_\delta r_\epsilon
1615 > + \text{9 permutations}
1616   \right) \left(- \frac{15g_n}{r^7}+\frac{15h_n}{r^7} -\frac{6s_n}{r^5} +\frac{t_n}{r^4}
1617   \right) \nonumber \\
1618 < + r_\alpha r_\beta r_\gamma r_\delta r_\epsilon
1618 > &+ r_\alpha r_\beta r_\gamma r_\delta r_\epsilon
1619   \left(  \frac{105g_n}{r^9} - \frac{105h_n}{r^8} + \frac{45s_n}{r^7}
1620 < - \frac{10t_n}{r^6} +\frac{u_n}{r^5} \right)
1621 < \end{eqnarray}
1620 > - \frac{10t_n}{r^6} +\frac{u_n}{r^5} \right) \label{eq:b13}
1621 > \end{align}
1622   %
1623   %
1624   %
1625 < \section{Method 2, the $r$-dependent factors}
1625 > \newpage
1626 > \section{The $r$-dependent factors for GSF electrostatics}
1627  
1628 < In method 2, the kernel is not expanded, rather the individual terms in the multipole interaction energies,
1629 < see Eq. (20?).  For a smeared-charge distribution, this still brings into the algebra multiple derivatives
1630 < of the kernel $B_0(r)$.  To denote these terms, we generalize the notation of the previous appendix.
1631 < For $f(r)=1/r$ (bare Coulomb) or $f(r)=B_0(r)$ (smeared charge)
1628 > In Gradient-shifted force electrostatics, the kernel is not expanded,
1629 > rather the individual terms in the multipole interaction energies.
1630 > For damped charges , this still brings into the algebra multiple
1631 > derivatives of the Smith's $B_0(r)$ function.  To denote these terms,
1632 > we generalize the notation of the previous appendix.  For either
1633 > $f(r)=1/r$ (undamped) or $f(r)=B_0(r)$ (damped),
1634   %
1635 < \begin{eqnarray}
1636 < g(r)= \frac{df}{d r}\\
1637 < h(r)= \frac{dg}{d r} = \frac{d^2f}{d r^2} \\
1638 < s(r)= \frac{dh}{d r} = \frac{d^3f}{d r^3} \\
1639 < t(r)= \frac{ds}{d r} = \frac{d^4f}{d r^4} \\
1640 < u(r)= \frac{dt}{d r} =\frac{d^5f}{d r^5} .
1641 < \end{eqnarray}
1635 > \begin{align}
1636 > g(r)=& \frac{df}{d r}\\
1637 > h(r)=& \frac{dg}{d r} = \frac{d^2f}{d r^2} \\
1638 > s(r)=& \frac{dh}{d r} = \frac{d^3f}{d r^3} \\
1639 > t(r)=& \frac{ds}{d r} = \frac{d^4f}{d r^4} \\
1640 > u(r)=& \frac{dt}{d r} = \frac{d^5f}{d r^5} .
1641 > \end{align}
1642   %
1643 < For $f(r)=1/r$, Table I lists these derivatives under the column ``Bare Coulomb.''  Checks of algebra can be made by using limiting forms
1644 < of equations, e.g., the leading term in the function $g_n(r)$ has $r$ dependence given by $g(r)$.  Equations (B9) to B(13)
1645 < are correct for method 2 if one just eliminates the subscript $n$.
1574 <
1575 < \section{Extra Material}
1576 < %
1577 < %
1578 < %Energy in body coordinate form ---------------------------------------------------------------
1579 < %
1580 < Here are the interaction energies written in terms of the body coordinates:
1581 <
1582 < %
1583 < % u ca cb
1584 < %
1585 < \begin{equation}
1586 < U_{C_{\bf a}C_{\bf b}}(r)=
1587 < \frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r)
1588 < \end{equation}
1589 < %
1590 < % u ca db
1591 < %
1592 < \begin{equation}
1593 < U_{C_{\bf a}D_{\bf b}}(r)=
1594 < \frac{C_{\bf a}}{4\pi \epsilon_0}
1595 < \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} \,  v_{11}(r)
1596 < \end{equation}
1597 < %
1598 < % u ca qb
1599 < %
1600 < \begin{equation}
1601 < U_{C_{\bf a}Q_{\bf b}}(r)=
1602 < \frac{C_{\bf a }\text{Tr}Q_{\bf b}}{4\pi \epsilon_0}  
1603 < v_{21}(r) \nonumber \\
1604 < +\frac{C_{\bf a}}{4\pi \epsilon_0}
1605 < \sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r})
1606 < v_{22}(r)
1607 < \end{equation}
1608 < %
1609 < % u da cb
1610 < %
1611 < \begin{equation}
1612 < U_{D_{\bf a}C_{\bf b}}(r)=
1613 < -\frac{C_{\bf b}}{4\pi \epsilon_0}  
1614 < \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} \,  v_{11}(r)
1615 < \end{equation}
1616 < %
1617 < % u da db
1618 < %
1619 < \begin{equation}
1620 < \begin{split}
1621 < % 1
1622 < U_{D_{\bf a}D_{\bf b}}(r)&=
1623 < -\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
1624 < (\hat{a}_m \cdot   \hat{b}_n)
1625 < D_{\mathbf{b}n} v_{21}(r) \\
1626 < % 2
1627 < &-\frac{1}{4\pi \epsilon_0}
1628 < \sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1629 < \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n}
1630 < v_{22}(r)
1631 < \end{split}
1632 < \end{equation}
1633 < %
1634 < % u da qb
1635 < %
1636 < \begin{equation}
1637 < \begin{split}
1638 < % 1
1639 < U_{D_{\bf a}Q_{\bf b}}(r)&=
1640 < -\frac{1}{4\pi \epsilon_0} \left(
1641 < \text{Tr}Q_{\mathbf{b}}
1642 < \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
1643 < +2\sum_{lmn}D_{\mathbf{a}l}
1644 < (\hat{a}_l \cdot \hat{b}_m)
1645 < Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
1646 < \right)  v_{31}(r) \\
1647 < % 2
1648 < &-\frac{1}{4\pi \epsilon_0}
1649 < \sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
1650 < \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1651 < Q_{{\mathbf b}mn}
1652 < (\hat{b}_n \cdot \hat{r}) v_{32}(r)
1653 < \end{split}
1654 < \end{equation}
1655 < %
1656 < % u qa cb
1657 < %
1658 < \begin{equation}
1659 < U_{Q_{\bf a}C_{\bf b}}(r)=
1660 < \frac{C_{\bf b }\text{Tr}Q_{\bf a}}{4\pi \epsilon_0}  v_{21}(r)
1661 < +\frac{C_{\bf b}}{4\pi \epsilon_0}
1662 < \sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) v_{22}(r)
1663 < \end{equation}
1664 < %
1665 < % u qa db
1666 < %
1667 < \begin{equation}
1668 < \begin{split}
1669 < %1
1670 < U_{Q_{\bf a}D_{\bf b}}(r)&=
1671 < \frac{1}{4\pi \epsilon_0} \left(
1672 < \text{Tr}Q_{\mathbf{a}}
1673 < \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n}
1674 < +2\sum_{lmn}D_{\mathbf{b}l}
1675 < (\hat{b}_l \cdot \hat{a}_m)
1676 < Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
1677 < \right) v_{31}(r)  \\
1678 < % 2
1679 < &+\frac{1}{4\pi \epsilon_0}
1680 < \sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
1681 < \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1682 < Q_{{\mathbf a}mn}
1683 < (\hat{a}_n \cdot \hat{r}) v_{32}(r)
1684 < \end{split}
1685 < \end{equation}
1686 < %
1687 < % u qa qb
1688 < %
1689 < \begin{equation}
1690 < \begin{split}
1691 < %1
1692 < U_{Q_{\bf a}Q_{\bf b}}(r)&=
1693 < \frac{1}{4\pi \epsilon_0} \Bigl[
1694 < \text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}}
1695 < +2\sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
1696 < Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
1697 < (\hat{a}_m \cdot \hat{b}_n) \Bigr]
1698 < v_{41}(r) \\
1699 < % 2
1700 < &+ \frac{1}{4\pi \epsilon_0}
1701 < \Bigl[ \text{Tr}Q_{\mathbf{a}}
1702 < \sum_{lm} (\hat{r} \cdot \hat{b}_l )
1703 < Q_{{\mathbf b}lm}
1704 < (\hat{b}_m \cdot \hat{r})
1705 < +\text{Tr}Q_{\mathbf{b}}
1706 < \sum_{lm} (\hat{r} \cdot \hat{a}_l )
1707 < Q_{{\mathbf a}lm}
1708 < (\hat{a}_m \cdot \hat{r}) \\
1709 < % 3
1710 < &+4 \sum_{lmnp}
1711 < (\hat{r} \cdot \hat{a}_l )
1712 < Q_{{\mathbf a}lm}
1713 < (\hat{a}_m \cdot \hat{b}_n)
1714 < Q_{{\mathbf b}np}
1715 < (\hat{b}_p \cdot \hat{r})
1716 < \Bigr] v_{42}(r)  \\
1717 < % 4
1718 < &+ \frac{1}{4\pi \epsilon_0}
1719 < \sum_{lm} (\hat{r} \cdot \hat{a}_l)
1720 < Q_{{\mathbf a}lm}
1721 < (\hat{a}_m \cdot \hat{r})
1722 < \sum_{np}  (\hat{r} \cdot \hat{b}_n)
1723 < Q_{{\mathbf b}np}
1724 < (\hat{b}_p \cdot \hat{r})  v_{43}(r).
1725 < \end{split}
1726 < \end{equation}
1727 < %
1728 <
1729 <
1730 < % BODY coordinates force equations --------------------------------------------
1731 < %
1732 < %
1733 < Here are the force equations written in terms of body coordinates.
1734 < %
1735 < % f ca cb
1736 < %
1737 < \begin{equation}
1738 < \mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =
1739 < \frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r}
1740 < \end{equation}
1741 < %
1742 < % f ca db
1743 < %
1744 < \begin{equation}
1745 < \mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =
1746 < \frac{C_{\bf a}}{4\pi \epsilon_0}  
1747 < \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} w_b(r) \hat{r}
1748 < +\frac{C_{\bf a}}{4\pi \epsilon_0}  
1749 < \sum_n  D_{\mathbf{b}n} \hat{b}_n w_c(r)
1750 < \end{equation}
1751 < %
1752 < % f ca qb
1753 < %
1754 < \begin{equation}
1755 < \begin{split}
1756 < % 1
1757 < \mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =
1758 < \frac{1}{4\pi \epsilon_0}  
1759 < C_{\bf a }\text{Tr}Q_{\bf b} w_d(r) \hat{r}
1760 < + 2C_{\bf a } \sum_l  \hat{b}_l Q_{{\mathbf b}ln} (\hat{b}_n \cdot \hat{r}) w_e(r) \\
1761 < % 2
1762 < +\frac{C_{\bf a}}{4\pi \epsilon_0}
1763 < \sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r}) w_f(r) \hat{r}
1764 < \end{split}
1765 < \end{equation}
1766 < %
1767 < % f da cb
1768 < %
1769 < \begin{equation}
1770 < \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1771 < -\frac{C_{\bf{b}}}{4\pi \epsilon_0}
1772 < \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} w_b(r)  \hat{r}
1773 < -\frac{C_{\bf{b}}}{4\pi \epsilon_0}
1774 < \sum_n  D_{\mathbf{a}n} \hat{a}_n w_c(r)
1775 < \end{equation}
1776 < %
1777 < % f da db
1778 < %
1779 < \begin{equation}
1780 < \begin{split}
1781 < % 1
1782 < \mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} &=
1783 < -\frac{1}{4\pi \epsilon_0}
1784 < \sum_{mn} D_{\mathbf {a}m}
1785 < (\hat{a}_m \cdot   \hat{b}_n)
1786 < D_{\mathbf{b}n}  w_d(r) \hat{r}
1787 < -\frac{1}{4\pi \epsilon_0}
1788 < \sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1789 < \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} w_f(r) \hat{r} \\
1790 < % 2
1791 < & \quad + \frac{1}{4\pi \epsilon_0}
1792 < \Bigl[ \sum_m D_{\mathbf {a}m}
1793 < \hat{a}_m \sum_n D_{\mathbf{b}n}
1794 < (\hat{b}_n \cdot \hat{r})
1795 < + \sum_m D_{\mathbf {b}m}
1796 < \hat{b}_m \sum_n D_{\mathbf{a}n}
1797 < (\hat{a}_n \cdot \hat{r}) \Bigr] w_e(r)  \\
1798 < \end{split}
1799 < \end{equation}
1800 < %
1801 < % f da qb
1802 < %
1803 < \begin{equation}
1804 < \begin{split}
1805 < % 1
1806 < &\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1807 < - \frac{1}{4\pi \epsilon_0} \Bigl[
1808 < \text{Tr}Q_{\mathbf{b}}
1809 < \sum_l  D_{\mathbf{a}l} \hat{a}_l
1810 < +2\sum_{lmn} D_{\mathbf{a}l}
1811 < (\hat{a}_l \cdot \hat{b}_m)
1812 < Q_{\mathbf{b}mn} \hat{b}_n  \Bigr] w_g(r) \\
1813 < % 3
1814 < &  - \frac{1}{4\pi \epsilon_0} \Bigl[
1815 < \text{Tr}Q_{\mathbf{b}}
1816 < \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
1817 < +2\sum_{lmn}D_{\mathbf{a}l}
1818 < (\hat{a}_l \cdot \hat{b}_m)
1819 < Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
1820 < % 4
1821 < &+ \frac{1}{4\pi \epsilon_0}
1822 < \Bigl[\sum_l  D_{\mathbf{a}l} \hat{a}_l
1823 < \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1824 < Q_{{\mathbf b}mn}
1825 < (\hat{b}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{a}_l)
1826 < D_{\mathbf{a}l}
1827 < \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1828 < Q_{{\mathbf b}mn} \hat{b}_n \Bigr]   w_i(r)\\
1829 < % 6
1830 < &  -\frac{1}{4\pi \epsilon_0}
1831 < \sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
1832 < \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1833 < Q_{{\mathbf b}mn}
1834 < (\hat{b}_n \cdot \hat{r})  w_j(r)  \hat{r}
1835 < \end{split}
1836 < \end{equation}
1837 < %
1838 < % force qa cb
1839 < %
1840 < \begin{equation}
1841 < \begin{split}
1842 < % 1
1843 < \mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} &=
1844 < \frac{1}{4\pi \epsilon_0}  
1845 < C_{\bf b }\text{Tr}Q_{\bf a} \hat{r} w_d(r)
1846 < + \frac{2C_{\bf b }}{4\pi \epsilon_0}  \sum_l  \hat{a}_l Q_{{\mathbf a}ln} (\hat{a}_n \cdot \hat{r}) w_e(r) \\
1847 < % 2
1848 < &  +\frac{C_{\bf b}}{4\pi \epsilon_0}
1849 < \sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) w_f(r) \hat{r}
1850 < \end{split}
1851 < \end{equation}
1852 < %
1853 < % f qa db
1854 < %
1855 < \begin{equation}
1856 < \begin{split}
1857 < % 1
1858 < &\mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1859 < \frac{1}{4\pi \epsilon_0} \Bigl[
1860 < \text{Tr}Q_{\mathbf{a}}
1861 < \sum_l  D_{\mathbf{b}l} \hat{b}_l
1862 < +2\sum_{lmn} D_{\mathbf{b}l}
1863 < (\hat{b}_l \cdot \hat{a}_m)
1864 < Q_{\mathbf{a}mn} \hat{a}_n  \Bigr]
1865 < w_g(r)\\
1866 < % 3
1867 < &  + \frac{1}{4\pi \epsilon_0} \Bigl[
1868 < \text{Tr}Q_{\mathbf{a}}
1869 < \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n}
1870 < +2\sum_{lmn}D_{\mathbf{b}l}
1871 < (\hat{b}_l \cdot \hat{a}_m)
1872 < Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
1873 < % 4
1874 < &  + \frac{1}{4\pi \epsilon_0} \Bigl[ \sum_l  D_{\mathbf{b}l} \hat{b}_l
1875 < \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1876 < Q_{{\mathbf a}mn}
1877 < (\hat{a}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{b}_l)
1878 < D_{\mathbf{b}l}
1879 < \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1880 < Q_{{\mathbf a}mn} \hat{a}_n \Bigr]   w_i(r) \\
1881 < % 6
1882 < &  +\frac{1}{4\pi \epsilon_0}
1883 < \sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
1884 < \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1885 < Q_{{\mathbf a}mn}
1886 < (\hat{a}_n \cdot \hat{r})  w_j(r)  \hat{r}
1887 < \end{split}
1888 < \end{equation}
1889 < %
1890 < % f qa qb
1891 < %
1892 < \begin{equation}
1893 < \begin{split}
1894 < &\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1895 < \frac{1}{4\pi \epsilon_0} \Bigl[
1896 < \text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}}
1897 < + 2 \sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
1898 < Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
1899 < (\hat{a}_m \cdot \hat{b}_n) \Bigr] w_k(r) \hat{r}\\
1900 < &+\frac{1}{4\pi \epsilon_0} \Bigl[
1901 < 2\text{Tr}Q_{\mathbf{b}} \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  \hat{a}_m  
1902 < + 2\text{Tr}Q_{\mathbf{a}} \sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm}  \hat{b}_m \\
1903 < &+ 4\sum_{lmnp} \hat{a}_l Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_p \cdot \hat{r})  
1904 < + 4\sum_{lmnp} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{b}_n) Q_{\mathbf{b}np} \hat{b}_p
1905 < \Bigr] w_n(r) \\
1906 < &+ \frac{1}{4\pi \epsilon_0}
1907 < \Bigl[ \text{Tr}Q_{\mathbf{a}}
1908 < \sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm} (\hat{b}_m \cdot \hat{r})
1909 < + \text{Tr}Q_{\mathbf{b}}
1910 < \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  (\hat{a}_m \cdot \hat{r}) \\
1911 < &+4\sum_{lmnp} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{b}_n)
1912 < Q_{\mathbf{b}np}  (\hat{b}_p \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
1913 < %
1914 < &+\frac{1}{4\pi \epsilon_0} \Bigl[
1915 < 2\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} \hat{a}_m
1916 < \sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_n \cdot \hat{r}) \\
1917 < &+2 \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
1918 < \sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} \hat{b}_n \Bigr] w_o(r) \hat{r} \\
1919 < &  + \frac{1}{4\pi \epsilon_0}
1920 < \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
1921 < \sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_p \cdot \hat{r}) w_m(r) \hat{r}
1922 < \end{split}
1923 < \end{equation}
1924 < %
1925 < Here we list the form of the non-zero damped shifted multipole torques showing
1926 < explicitly dependences on body axes:
1927 < %
1928 < %  t ca db
1929 < %
1930 < \begin{equation}
1931 < \mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =
1932 < \frac{C_{\bf a}}{4\pi \epsilon_0}  
1933 < \sum_n  (\hat{r} \times \hat{b}_n)  D_{\mathbf{b}n} \,  v_{11}(r)
1934 < \end{equation}
1935 < %
1936 < % t ca qb
1937 < %
1938 < \begin{equation}
1939 < \mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =
1940 < \frac{2C_{\bf a}}{4\pi \epsilon_0}
1941 < \sum_{lm} (\hat{r} \times \hat{b}_l) Q_{{\mathbf b}lm} (\hat{b}_m \cdot \hat{r}) v_{22}(r)
1942 < \end{equation}
1943 < %
1944 < %  t da cb
1945 < %
1946 < \begin{equation}
1947 < \mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =
1948 < -\frac{C_{\bf b}}{4\pi \epsilon_0}  
1949 < \sum_n  (\hat{r} \times \hat{a}_n)  D_{\mathbf{a}n} \,  v_{11}(r)
1950 < \end{equation}%
1951 < %
1952 < %
1953 < %  ta da db
1954 < %
1955 < \begin{equation}
1956 < \begin{split}
1957 < % 1
1958 < \mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} &=
1959 < \frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
1960 < (\hat{a}_m \times  \hat{b}_n)
1961 < D_{\mathbf{b}n} v_{21}(r) \\
1962 < % 2
1963 < &-\frac{1}{4\pi \epsilon_0}
1964 < \sum_m (\hat{r} \times \hat{a}_m) D_{\mathbf {a}m}
1965 < \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
1966 < \end{split}
1967 < \end{equation}
1968 < %
1969 < %  tb da db
1970 < %
1971 < \begin{equation}
1972 < \begin{split}
1973 < % 1
1974 < \mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} &=
1975 < -\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
1976 < (\hat{a}_m \times  \hat{b}_n)
1977 < D_{\mathbf{b}n} v_{21}(r) \\
1978 < % 2
1979 < &+\frac{1}{4\pi \epsilon_0}
1980 < \sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1981 < \sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
1982 < \end{split}
1983 < \end{equation}
1984 < %
1985 < % ta da qb
1986 < %
1987 < \begin{equation}
1988 < \begin{split}
1989 < % 1
1990 < \mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} &=
1991 < \frac{1}{4\pi \epsilon_0} \left(
1992 < -\text{Tr}Q_{\mathbf{b}}
1993 < \sum_n (\hat{r} \times \hat{a}_n) D_{\mathbf{a}n}
1994 < +2\sum_{lmn}D_{\mathbf{a}l}
1995 < (\hat{a}_l \times \hat{b}_m)
1996 < Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
1997 < \right) v_{31}(r)\\
1998 < % 2
1999 < &-\frac{1}{4\pi \epsilon_0}
2000 < \sum_l (\hat{r} \times \hat{a}_l) D_{\mathbf{a}l}
2001 < \sum_{mn} (\hat{r} \cdot \hat{b}_m)
2002 < Q_{{\mathbf b}mn}
2003 < (\hat{b}_n \cdot \hat{r}) v_{32}(r)
2004 < \end{split}
2005 < \end{equation}
2006 < %
2007 < % tb da qb
2008 < %
2009 < \begin{equation}
2010 < \begin{split}
2011 < % 1
2012 < \mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} &=
2013 < \frac{1}{4\pi \epsilon_0} \left(
2014 < -2\sum_{lmn}D_{\mathbf{a}l}
2015 < (\hat{a}_l \cdot \hat{b}_m)
2016 < Q_{\mathbf{b}mn} (\hat{r} \times \hat{b}_n)
2017 < -2\sum_{lmn}D_{\mathbf{a}l}
2018 < (\hat{a}_l \times \hat{b}_m)
2019 < Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
2020 < \right) v_{31}(r) \\
2021 < % 2
2022 < &-\frac{2}{4\pi \epsilon_0}
2023 < \sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
2024 < \sum_{mn} (\hat{r} \cdot \hat{b}_m)
2025 < Q_{{\mathbf b}mn}
2026 < (\hat{r}\times \hat{b}_n) v_{32}(r)
2027 < \end{split}
2028 < \end{equation}
2029 < %
2030 < % ta qa cb
2031 < %
2032 < \begin{equation}
2033 < \mathbf{\tau}_{{\bf a}Q_{\bf a}C_{\bf b}} =
2034 < \frac{2C_{\bf a}}{4\pi \epsilon_0}
2035 < \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{{\mathbf a}lm} (\hat{r} \times \hat{a}_m) v_{22}(r)
2036 < \end{equation}
2037 < %
2038 < % ta qa db
2039 < %
2040 < \begin{equation}
2041 < \begin{split}
2042 < % 1
2043 < \mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} &=
2044 < \frac{1}{4\pi \epsilon_0} \left(
2045 < 2\sum_{lmn}D_{\mathbf{b}l}
2046 < (\hat{b}_l \cdot \hat{a}_m)
2047 < Q_{\mathbf{a}mn} (\hat{r} \times \hat{a}_n)
2048 < +2\sum_{lmn}D_{\mathbf{b}l}
2049 < (\hat{a}_l \times \hat{b}_m)
2050 < Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
2051 < \right) v_{31}(r) \\
2052 < % 2
2053 < &+\frac{2}{4\pi \epsilon_0}
2054 < \sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
2055 < \sum_{mn} (\hat{r} \cdot \hat{a}_m)
2056 < Q_{{\mathbf a}mn}
2057 < (\hat{r}\times \hat{a}_n) v_{32}(r)
2058 < \end{split}
2059 < \end{equation}
2060 < %
2061 < % tb qa db
2062 < %
2063 < \begin{equation}
2064 < \begin{split}
2065 < % 1
2066 < \mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} &=
2067 < \frac{1}{4\pi \epsilon_0} \left(
2068 < \text{Tr}Q_{\mathbf{a}}
2069 < \sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf{b}n}
2070 < +2\sum_{lmn}D_{\mathbf{b}l}
2071 < (\hat{a}_l \times \hat{b}_m)
2072 < Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
2073 < \right) v_{31}(r)\\
2074 < % 2
2075 < &\frac{1}{4\pi \epsilon_0}
2076 < \sum_l (\hat{r} \times \hat{b}_l) D_{\mathbf{b}l}
2077 < \sum_{mn} (\hat{r} \cdot \hat{a}_m)
2078 < Q_{{\mathbf a}mn}
2079 < (\hat{a}_n \cdot \hat{r}) v_{32}(r)
2080 < \end{split}
2081 < \end{equation}
2082 < %
2083 < % ta qa qb
2084 < %
2085 < \begin{equation}
2086 < \begin{split}
2087 < % 1
2088 < \mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} &=
2089 < -\frac{4}{4\pi \epsilon_0}
2090 < \sum_{lmnp} (\hat{a}_l \times \hat{b}_p)
2091 < Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
2092 < (\hat{a}_m \cdot \hat{b}_n) v_{41}(r) \\
2093 < % 2
2094 < &+ \frac{1}{4\pi \epsilon_0}
2095 < \Bigl[
2096 < 2\text{Tr}Q_{\mathbf{b}}
2097 < \sum_{lm} (\hat{r} \cdot \hat{a}_l )
2098 < Q_{{\mathbf a}lm}
2099 < (\hat{r} \times \hat{a}_m)
2100 < +4 \sum_{lmnp}
2101 < (\hat{r} \times \hat{a}_l )
2102 < Q_{{\mathbf a}lm}
2103 < (\hat{a}_m \cdot \hat{b}_n)
2104 < Q_{{\mathbf b}np}
2105 < (\hat{b}_p \cdot \hat{r}) \\
2106 < % 3
2107 < &-4 \sum_{lmnp}
2108 < (\hat{r} \cdot \hat{a}_l )
2109 < Q_{{\mathbf a}lm}
2110 < (\hat{a}_m \times \hat{b}_n)
2111 < Q_{{\mathbf b}np}
2112 < (\hat{b}_p \cdot \hat{r})
2113 < \Bigr] v_{42}(r) \\
2114 < % 4
2115 < &+ \frac{2}{4\pi \epsilon_0}
2116 < \sum_{lm} (\hat{r} \times \hat{a}_l)
2117 < Q_{{\mathbf a}lm}
2118 < (\hat{a}_m \cdot \hat{r})
2119 < \sum_{np}  (\hat{r} \cdot \hat{b}_n)
2120 < Q_{{\mathbf b}np}
2121 < (\hat{b}_p \cdot \hat{r})  v_{43}(r)\\
2122 < \end{split}
2123 < \end{equation}
2124 < %
2125 < % tb qa qb
2126 < %
2127 < \begin{equation}
2128 < \begin{split}
2129 < % 1
2130 < \mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} &=
2131 < \frac{4}{4\pi \epsilon_0}
2132 < \sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
2133 < Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
2134 < (\hat{a}_m \times \hat{b}_n) v_{41}(r) \\
2135 < % 2
2136 < &+ \frac{1}{4\pi \epsilon_0}
2137 < \Bigl[
2138 < 2\text{Tr}Q_{\mathbf{a}}
2139 < \sum_{lm} (\hat{r} \cdot \hat{b}_l )
2140 < Q_{{\mathbf b}lm}
2141 < (\hat{r} \times \hat{b}_m)
2142 < +4 \sum_{lmnp}
2143 < (\hat{r} \cdot \hat{a}_l )
2144 < Q_{{\mathbf a}lm}
2145 < (\hat{a}_m \cdot \hat{b}_n)
2146 < Q_{{\mathbf b}np}
2147 < (\hat{r} \times \hat{b}_p) \\
2148 < % 3
2149 < &+4 \sum_{lmnp}
2150 < (\hat{r} \cdot \hat{a}_l )
2151 < Q_{{\mathbf a}lm}
2152 < (\hat{a}_m \times \hat{b}_n)
2153 < Q_{{\mathbf b}np}
2154 < (\hat{b}_p \cdot \hat{r})
2155 < \Bigr] v_{42}(r)  \\
2156 < % 4
2157 < &+ \frac{2}{4\pi \epsilon_0}
2158 < \sum_{lm} (\hat{r} \cdot \hat{a}_l)
2159 < Q_{{\mathbf a}lm}
2160 < (\hat{a}_m \cdot \hat{r})
2161 < \sum_{np}  (\hat{r} \times \hat{b}_n)
2162 < Q_{{\mathbf b}np}
2163 < (\hat{b}_p \cdot \hat{r}) v_{43}(r). \\
2164 < \end{split}
2165 < \end{equation}
2166 < %
2167 < \begin{table*}
2168 < \caption{\label{tab:tableFORCE2}Radial functions used in the force equations.}
2169 < \begin{ruledtabular}
2170 < \begin{tabular}{ccc}
2171 < Generic&Method 1&Method 2
2172 < \\ \hline
2173 < %
2174 < %
2175 < %
2176 < $w_a(r)$&
2177 < $g_0(r)$&
2178 < $g(r)-g(r_c)$ \\
2179 < %
2180 < %
2181 < $w_b(r)$ &
2182 < $\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ &
2183 < $h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\
2184 < %
2185 < $w_c(r)$ &
2186 < $\frac{g_1(r)}{r} $ &
2187 < $\frac{v_{11}(r)}{r}$ \\
2188 < %
2189 < %
2190 < $w_d(r)$&
2191 < $\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ &
2192 < $\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)
2193 < -\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $\\
2194 < %
2195 < $w_e(r)$ &
2196 < $\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ &
2197 < $\frac{v_{22}(r)}{r}$ \\
2198 < %
2199 < %
2200 < $w_f(r)$&
2201 < $\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ &
2202 < $\left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) - $ \\
2203 < &&$\left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)-\frac{2v_{22}(r)}{r}$\\
2204 < %
2205 < $w_g(r)$& $ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$&
2206 < $\frac{v_{31}(r)}{r}$\\
2207 < %
2208 < $w_h(r)$ &
2209 < $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
2210 < $\left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) - $\\
2211 < &&$\left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\
2212 < &&$-\frac{v_{31}(r)}{r}$\\
2213 < % 2
2214 < $w_i(r)$ &
2215 < $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
2216 < $\frac{v_{32}(r)}{r}$ \\
2217 < %
2218 < $w_j(r)$ &
2219 < $\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right)  $ &
2220 < $\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right) $  \\
2221 < &&$\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}$ \\
2222 < %
2223 < $w_k(r)$ &
2224 < $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
2225 < $\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}  \right)$  \\
2226 < &&$\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2}  \right)$ \\
2227 < %
2228 < $w_l(r)$ &
2229 < $\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)$ &
2230 < $\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\
2231 < &&$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)
2232 < -\frac{2v_{42}(r)}{r}$ \\
2233 < %
2234 < $w_m(r)$ &
2235 < $\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)$ &
2236 < $\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$ \\
2237 < &&$\left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}
2238 < +\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $ \\
2239 < &&$-\frac{4v_{43}(r)}{r}$ \\
2240 < %
2241 < $w_n(r)$ &
2242 < $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
2243 < $\frac{v_{42}(r)}{r}$ \\
2244 < %
2245 < $w_o(r)$ &
2246 < $\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)$ &
2247 < $\frac{v_{43}(r)}{r}$ \\
2248 < %
2249 < \end{tabular}
2250 < \end{ruledtabular}
2251 < \end{table*}
1643 > For undamped charges Table I lists these derivatives under the column
1644 > ``Bare Coulomb.''  Equations \ref{eq:b9} to \ref{eq:b13} are still
1645 > correct for GSF electrostatics if the subscript $n$ is eliminated.
1646  
1647   \newpage
1648  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines