ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/group/trunk/multipole/multipole1.tex

Comparing trunk/multipole/multipole1.tex (file contents):
Revision 3984 by gezelter, Sat Dec 28 18:41:19 2013 UTC vs.
Revision 3985 by gezelter, Tue Dec 31 01:33:31 2013 UTC

#	Line 34 | Line 34 | jcp]{revtex4-1}
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
#	Line 101 | Line 102 | and DSF is now used routinely in simulations of ionic
102		similar methods where the force vanishes at $R_\textrm{c}$ are
103		essentially quantitative.\cite{Izvekov:2008wo} The DSF and other
104		related methods have now been widely investigated,\cite{Hansen:2012uq}
105	<	and DSF is now used routinely in simulations of ionic
106	<	liquids,\cite{doi:10.1021/la400226g,McCann:2013fk} flow in carbon
107	<	nanotubes,\cite{kannam:094701} gas sorption in metal-organic framework
108	<	materials,\cite{Forrest:2012ly} thermal conductivity of methane
109	<	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.
105	>	and DSF is now used routinely in a diverse set of chemical
106	>	environments.\cite{doi:10.1021/la400226g,McCann:2013fk,kannam:094701,Forrest:2012ly,English:2008kx,Louden:2013ve,Tokumasu:2013zr}
107	>	DSF electrostatics provides a compromise between the computational
108	>	speed of real-space cutoffs and the accuracy of fully-periodic Ewald
109	>	treatments.
110
114	–	\subsection{Coarse Graining using Point Multipoles}
111		One common feature of many coarse-graining approaches, which treat
112		entire molecular subsystems as a single rigid body, is simplification
113		of the electrostatic interactions between these bodies so that fewer
114		site-site interactions are required to compute configurational
115	<	energies. Notably, the force matching approaches of Voth and coworkers
116	<	are an exciting development in their ability to represent realistic
117	<	(and {\it reactive}) chemical systems for very large length scales and
118	<	long times.  This approach utilized a coarse-graining in interaction
119	<	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}
115	>	energies.  The coarse-graining approaches of Ren \&
116	>	coworkers,\cite{Golubkov06} and Essex \&
117	>	coworkers,\cite{ISI:000276097500009,ISI:000298664400012} both utilize
118	>	point multipoles to model electrostatics for entire molecules or
119	>	functional groups.
120
127	–	The coarse-graining approaches of Ren \& coworkers,\cite{Golubkov06}
128	–	and Essex \&
129	–	coworkers,\cite{ISI:000276097500009,ISI:000298664400012}
130	–	both utilize Gay-Berne
131	–	ellipsoids~\cite{Berne72,Gay81,Luckhurst90,Cleaver96,Berardi98,Ravichandran:1999fk,Berardi99,Pasterny00}
132	–	to model dispersive interactions and point multipoles to model
133	–	electrostatics for entire molecules or functional groups.
134	–
121		Ichiye and coworkers have recently introduced a number of very fast
122		water models based on a ``sticky'' multipole model which are
123		qualitatively better at reproducing the behavior of real water than
124	<	the more common point-charge models (SPC/E, TIPnP).  The point charge
125	<	models are also substantially more computationally demanding than the
126	<	sticky multipole
127	<	approach.\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} The
128	<	SSDQO model requires the use of an approximate multipole expansion
129	<	(AME) as the exact multipole expansion is quite expensive
144	<	(particularly when handled via the Ewald sum).\cite{Ichiye:2006qy}
124	>	the more common point-charge models (SPC/E,
125	>	TIPnP).\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} The SSDQO
126	>	model requires the use of an approximate multipole expansion (AME) as
127	>	the exact multipole expansion is quite expensive (particularly when
128	>	handled via the Ewald sum).\cite{Ichiye:2006qy}
129	>
130		Another particularly important use of point multipoles (and multipole
131		polarizability) is in the very high-quality AMOEBA water model and
132		related force fields.\cite{Ponder:2010fk,schnieders:124114,Ren:2011uq}
#	Line 154 | Line 139 | kernel to point multipoles.  We have developed two dis
139		a different orientation can cause energy discontinuities.
140
141		This paper outlines an extension of the original DSF electrostatic
142	<	kernel to point multipoles.  We have developed two distinct real-space
142	>	kernel to point multipoles.  We describe two distinct real-space
143		interaction models for higher-order multipoles based on two truncated
144		Taylor expansions that are carried out at the cutoff radius.  We are
145		calling these models {\bf Taylor-shifted} and {\bf Gradient-shifted}
146		electrostatics.  Because of differences in the initial assumptions,
147		the two methods yield related, but different expressions for energies,
148	<	forces, and torques.
148	>	forces, and torques.
149
150		In this paper we outline the new methodology and give functional forms
151		for the energies, forces, and torques up to quadrupole-quadrupole
152		order.  We also compare the new methods to analytic energy constants
153		for periodic arrays of point multipoles.  In the following paper, we
154	<	provide extensive numerical comparisons to Ewald-based electrostatics
155	<	in common simulation enviornments.
154	>	provide numerical comparisons to Ewald-based electrostatics in common
155	>	simulation enviornments.
156
157		\section{Methodology}
158
#	Line 189 | Line 174 | + \frac{2\alpha}{\pi^{1/2}}
174		\right)\left(r_{ij}-R_\mathrm{c}\right)\ \Biggr{]}
175		\label{eq:DSFPot}
176		\end{equation*}
177	+	Note that in this potential and in all electrostatic quantities that
178	+	follow, the standard $4 \pi \epsilon_{0}$ has been omitted for
179	+	clarity.
180
181		To insure net charge neutrality within each cutoff sphere, an
182		additional ``self'' term is added to the potential. This term is
#	Line 199 | Line 187 | damping.\cite{deLeeuw80,Wolf99}
187		over the surface of the cutoff sphere.  A portion of the self term is
188		identical to the self term in the Ewald summation, and comes from the
189		utilization of the complimentary error function for electrostatic
190	<	damping.\cite{deLeeuw80,Wolf99}
190	>	damping.\cite{deLeeuw80,Wolf99}
191
192		There have been recent efforts to extend the Wolf self-neutralization
193		method to zero out the dipole and higher order multipoles contained
194		within the cutoff
195	<	sphere.\cite{Fukuda:2011jk,Fukuda:2012yu,Fukuda:2013qv}
195	>	sphere.\cite{Fukuda:2011jk,Fukuda:2012yu,Fukuda:2013qv}
196
197	<	In this work, we will be extending the idea of self-neutralization for
198	<	the point multipoles in a similar way.  In Figure
199	<	\ref{fig:shiftedMultipoles}, the central dipolar site $\mathbf{D}_i$
200	<	is interacting with point dipole $\mathbf{D}_j$ and point quadrupole,
201	<	$\mathbf{Q}_k$.  The self-neutralization scheme for point multipoles
202	<	involves projecting opposing multipoles for sites $j$ and $k$ on the
203	<	surface of the cutoff sphere.
197	>	In this work, we extend the idea of self-neutralization for the point
198	>	multipoles by insuring net charge-neutrality and net-zero moments
199	>	within each cutoff sphere.  In Figure \ref{fig:shiftedMultipoles}, the
200	>	central dipolar site $\mathbf{D}_i$ is interacting with point dipole
201	>	$\mathbf{D}_j$ and point quadrupole, $\mathbf{Q}_k$.  The
202	>	self-neutralization scheme for point multipoles involves projecting
203	>	opposing multipoles for sites $j$ and $k$ on the surface of the cutoff
204	>	sphere.  There are also significant modifications made to make the
205	>	forces and torques go smoothly to zero at the cutoff distance.
206
207		\begin{figure}
208		\includegraphics[width=3in]{SM}
#	Line 237 | Line 227 | V_a(\mathbf r) = \frac{1}{4\pi\epsilon_0}
227		a$.  Then the electrostatic potential of object $\bf a$ at
228		$\mathbf{r}$ is given by
229		\begin{equation}
230	<	V_a(\mathbf r) = \frac{1}{4\pi\epsilon_0}
230	>	V_a(\mathbf r) =
231		\sum_{k \, \text{in \bf a}} \frac{q_k}{\lvert \mathbf{r} - \mathbf{r}_k \rvert}.
232		\end{equation}
233		The Taylor expansion in $r$ can be written using an implied summation
#	Line 256 | Line 246 | V_{\bf a}(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\hat{M
246		can then be used to express the electrostatic potential on $\bf a$ in
247		terms of multipole operators,
248		\begin{equation}
249	<	V_{\bf a}(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\hat{M}_{\bf a} \frac{1}{r}
249	>	V_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r}
250		\end{equation}
251		where
252		\begin{equation}
#	Line 281 | Line 271 | U_{\bf{ab}}(r)
271		$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $), the interaction energy is given by
272		\begin{equation}
273		U_{\bf{ab}}(r)
274	<	= \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} .
274	>	= \hat{M}_a \sum_{j \, \text{in \bf b}} \frac {q_j}{\vert \bf{r}+\bf{r}_j \vert} .
275		\end{equation}
276		This can also be expanded as a Taylor series in $r$.  Using a notation
277		similar to before to define the multipoles on object {\bf b},
#	Line 292 | Line 282 | U_{\bf{ab}}(r)=\frac{\hat{M}_{\bf a} \hat{M}_{\bf b}}{
282		\end{equation}
283		we arrive at the multipole expression for the total interaction energy.
284		\begin{equation}
285	<	U_{\bf{ab}}(r)=\frac{\hat{M}_{\bf a} \hat{M}_{\bf b}}{4\pi \epsilon_0} \frac{1}{r}  \label{kernel}.
285	>	U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r}  \label{kernel}.
286		\end{equation}
287		This form has the benefit of separating out the energies of
288		interaction into contributions from the charge, dipole, and quadrupole
#	Line 316 | Line 306 | summarized in Appendix A.
306		these functions are known.  Smith's convenient functions $B_l(r)$ are
307		summarized in Appendix A.
308
319	–
309		The main goal of this work is to smoothly cut off the interaction
310		energy as well as forces and torques as $r\rightarrow r_c$.  To
311		describe how this goal may be met, we use two examples, charge-charge
#	Line 328 | Line 317 | U_{C_{\bf a}C_{\bf b}}(r)=\frac{1}{4\pi \epsilon_0} C_
317		charges $C_{\bf a}$ and $C_{\bf b}$ separated by a distance $r$ is
318		written:
319		\begin{equation}
320	<	U_{C_{\bf a}C_{\bf b}}(r)=\frac{1}{4\pi \epsilon_0} C_{\bf a} C_{\bf b}
320	>	U_{C_{\bf a}C_{\bf b}}(r)= C_{\bf a} C_{\bf b}
321		\left({ \frac{1}{r} - \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2}  }
322		\right) .
323		\end{equation}
#	Line 345 | Line 334 | There are a number of ways to generalize this derivati
334		\right) = \left(- \frac{1}{r^2} + \frac{1}{r_c^2}
335		\right) .
336		\end{equation}
337	<	There are a number of ways to generalize this derivative shift for
337	>	which clearly vanishes as the $r$ approaches the cutoff radius. There
338	>	are a number of ways to generalize this derivative shift for
339		higher-order multipoles.  Below, we present two methods, one based on
340		higher-order Taylor series for $r$ near $r_c$, and the other based on
341		linear shift of the kernel gradients at the cutoff itself.
#	Line 381 | Line 371 | U_{C_{\bf a}D_{\bf b}}(r)=
371		%
372		\begin{equation}
373		U_{C_{\bf a}D_{\bf b}}(r)=
374	<	\frac{C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}  \frac {\partial f_1(r) }{\partial r_\alpha}
375	<	=\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}
374	>	C_{\bf a} D_{{\bf b}\alpha}  \frac {\partial f_1(r) }{\partial r_\alpha}
375	>	= C_{\bf a} D_{{\bf b}\alpha}
376		\frac {r_\alpha}{r} \frac {\partial f_1(r)}{\partial r} .
377		\end{equation}
378		%
379		The force that dipole  $\bf b$ exerts on charge $\bf a$ is
380		%
381		\begin{equation}
382	<	F_{C_{\bf a}D_{\bf b}\beta} =\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}
382	>	F_{C_{\bf a}D_{\bf b}\beta} = C_{\bf a} D_{{\bf b}\alpha}
383		\left[ \frac{\delta_{\alpha\beta}}{r} \frac {\partial}{\partial r} +
384		\frac{r_\alpha r_\beta}{r^2}
385		\left( -\frac{1}{r} \frac {\partial} {\partial r}
#	Line 400 | Line 390 | F_{C_{\bf a}D_{\bf b}\beta} =
390		%
391		\begin{equation}
392		F_{C_{\bf a}D_{\bf b}\beta} =
393	<	\frac{C_{\bf a} D_{{\bf b}\beta} }{4\pi \epsilon_0r}
393	>	\frac{C_{\bf a} D_{{\bf b}\beta}}{r}
394		\left[  -\frac{1}{r^2}+\frac{1}{r_c^2}-\frac{2(r-r_c)}{r_c^3} \right]
395	<	+\frac{C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta }{4\pi \epsilon_0}
395	>	+C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta
396		\left[ \frac{3}{r^5}-\frac{3}{r^3r_c^2} \right] .
397		\end{equation}
398		%
#	Line 411 | Line 401 | U=\frac{1}{4\pi \epsilon_0} (\text{prefactor}) (\text{
401		In general, we can write
402		%
403		\begin{equation}
404	<	U=\frac{1}{4\pi \epsilon_0} (\text{prefactor}) (\text{derivatives}) f_n(r)
404	>	U= (\text{prefactor}) (\text{derivatives}) f_n(r)
405		\label{generic}
406		\end{equation}
407		%
408	<	where $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole
409	<	and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for quadrupole-quadrupole.
410	<	An example is the case of quadrupole-quadrupole for which the $\text{prefactor}$ is
411	<	$Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$ and the derivatives are
412	<	$\partial^4/\partial r_\alpha \partial r_\beta \partial r_\gamma \partial r_\delta$, with
413	<	implied summation combining the space indices.
408	>	with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for
409	>	charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and
410	>	$n=4$ for quadrupole-quadrupole.  For example, in
411	>	quadrupole-quadrupole interactions for which the $\text{prefactor}$ is
412	>	$Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$, the derivatives are
413	>	$\partial^4/\partial r_\alpha \partial r_\beta \partial
414	>	r_\gamma \partial r_\delta$, with implied summation combining the
415	>	space indices.
416
417		In the formulas presented in the tables below, the placeholder
418		function $f(r)$ is used to represent the electrostatic kernel (either
419		damped or undamped).  The main functions that go into the force and
420	<	torque terms, $f_n(r), g_n(r), h_n(r), s_n(r), \mathrm{~and~} t_n(r)$
421	<	are successive derivatives of the shifted electrostatic kernel of the
422	<	same index $n$.  The algebra required to evaluate energies, forces and
423	<	torques is somewhat tedious and are summarized in Appendices A and B.
420	>	torque terms, $g_n(r), h_n(r), s_n(r), \mathrm{~and~} t_n(r)$ are
421	>	successive derivatives of the shifted electrostatic kernel, $f_n(r)$
422	>	of the same index $n$.  The algebra required to evaluate energies,
423	>	forces and torques is somewhat tedious, so only the final forms are
424	>	presented in tables XX and YY.
425
426		\subsection{Gradient-shifted force (GSF) electrostatics}
427	<	Note the method used in the previous subsection to smoothly shift the
428	<	force to zero is a truncated Taylor Series in the radius $r$.  The
429	<	second method maintains only the linear $(r-r_c)$ term and has a
430	<	similar interaction energy for all multipole orders:
427	>	The second, and conceptually simpler approach to force-shifting
428	>	maintains only the linear $(r-r_c)$ term in the truncated Taylor
429	>	expansion, and has a similar interaction energy for all multipole
430	>	orders:
431		\begin{equation}
432	<	U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big \lvert  _{r_c} .
432	>	U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big
433	>	\lvert  _{r_c} .
434	>	\label{generic2}
435		\end{equation}
436	<	No higher order terms $(r-r_c)^n$ appear.  The primary difference
437	<	between the TSF and GSF methods is the stage at which the Taylor
438	<	Series is applied; in the Taylor-shifted approach, it is applied to
439	<	the kernel itself.  In the Gradient-shifted approach, it is applied to
440	<	individual radial interactions terms in the multipole expansion.
441	<	Terms from this method thus have the general form:
436	>	Here the gradient for force shifting is evaluated for an image
437	>	multipole on the surface of the cutoff sphere (see fig
438	>	\ref{fig:shiftedMultipoles}). No higher order terms $(r-r_c)^n$
439	>	appear.  The primary difference between the TSF and GSF methods is the
440	>	stage at which the Taylor Series is applied; in the Taylor-shifted
441	>	approach, it is applied to the kernel itself.  In the Gradient-shifted
442	>	approach, it is applied to individual radial interactions terms in the
443	>	multipole expansion.  Energies from this method thus have the general
444	>	form:
445		\begin{equation}
446	<	U=\frac{1}{4\pi \epsilon_0}\sum  (\text{angular factor}) (\text{radial factor}).
447	<	\label{generic2}
446	>	U= \sum  (\text{angular factor}) (\text{radial factor}).
447	>	\label{generic3}
448		\end{equation}
449
450	<	Results for both methods can be summarized using the form of
451	<	Eq.~(\ref{generic2}) and are listed in Table I below.
450	>	Functional forms for both methods (TSF and GSF) can be summarized
451	>	using the form of Eq.~(\ref{generic3}).  The basic forms for the
452	>	energy, force, and torque expressions are tabulated for both shifting
453	>	approaches below - for each separate orientational contribution, only
454	>	the radial factors differ between the two methods.
455
456		\subsection{\label{sec:level2}Body and space axes}
457
458	+	[XXX Do we need this section in the main paper? or should it go in the
459	+	extra materials?]
460	+
461		So far, all energies and forces have been written in terms of fixed
462	<	space coordinates $x$, $y$, $z$.  Interaction energies are computed
463	<	from the generic formulas Eq.~(\ref{generic}) and ~(\ref{generic2})
464	<	which combine prefactors with radial functions.  Because objects $\bf
462	>	space coordinates.  Interaction energies are computed from the generic
463	>	formulas Eq.~(\ref{generic}) and ~(\ref{generic2}) which combine
464	>	orientational prefactors with radial functions.  Because objects $\bf
465		a$ and $\bf b$ both translate and rotate during a molecular dynamics
466		(MD) simulation, it is desirable to contract all $r$-dependent terms
467		with dipole and quadrupole moments expressed in terms of their body
468	<	axes.  To do so, we follow the methodology of Allen and
469	<	Germano,\cite{Allen:2006fk} which was itself based on an earlier paper
470	<	by Price {\em et al.}\cite{Price:1984fk}
468	>	axes.  To do so, we have followed the methodology of Allen and
469	>	Germano,\cite{Allen:2006fk} which was itself based on earlier work by
470	>	Price {\em et al.}\cite{Price:1984fk}
471
472		We denote body axes for objects $\bf a$ and $\bf b$ by unit vectors
473		$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$
#	Line 475 | Line 479 | Allen and Germano define matrices $\hat{\mathbf {a}}$
479		\hat{a}_m= a_{mx}\hat{x} + a_{my}\hat{y} + a_{mz}\hat{z}  \\
480		\hat{b}_m= b_{mx}\hat{x} + b_{my}\hat{y} + b_{mz}\hat{z}  .
481		\end{eqnarray}
482	<	Allen and Germano define matrices $\hat{\mathbf {a}}$
483	<	and  $\hat{\mathbf {b}}$ using these unit vectors:
482	>	Rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be
483	>	expressed using these unit vectors:
484		\begin{eqnarray}
485		\hat{\mathbf {a}} =
486		\begin{pmatrix}
#	Line 498 | Line 502 | b_{1x}\quad  b_{1y} \quad b_{1z} \\
502		\end{pmatrix}
503		=
504		\begin{pmatrix}
505	<	b_{1x}\quad  b_{1y} \quad b_{1z} \\
505	>	b_{1x} \quad  b_{1y} \quad b_{1z} \\
506		b_{2x} \quad b_{2y} \quad b_{2z} \\
507		b_{3x} \quad b_{3y} \quad b_{3z}
508		\end{pmatrix}  .
509		\end{eqnarray}
510		%
511	<	These matrices convert from space-fixed $(xyz)$ to object-fixed $(123)$ coordinates.
512	<	All contractions of prefactors with derivatives of functions can be written in terms of these matrices.
513	<	It proves to be equally convenient to just write any contraction in terms of unit vectors
514	<	$\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$.
515	<	We have found it useful to write angular-dependent terms in three different fashions,
516	<	illustrated by the following
517	<	three examples from the interaction-energy expressions:
511	>	These matrices convert from space-fixed $(xyz)$ to body-fixed $(123)$
512	>	coordinates.  All contractions of prefactors with derivatives of
513	>	functions can be written in terms of these matrices. It proves to be
514	>	equally convenient to just write any contraction in terms of unit
515	>	vectors $\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$. In the torque
516	>	expressions, it is useful to have the angular-dependent terms
517	>	available in three different fashions, e.g. for the dipole-dipole
518	>	contraction:
519		%
520	<	\begin{eqnarray}
520	>	\begin{equation}
521		\mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}}
522	<	=D_{\bf {a}\alpha} D_{\bf {b}\alpha}=
523	<	\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}} \\
524	<	r^2 \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right)=
520	<	r_\alpha Q_{\bf b \alpha \beta} r_\beta = r^2
521	<	\sum_{mn}(\hat{r} \cdot \hat{b}_m)  Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \\
522	<	r ( \mathbf{D}_{\mathbf{a}} \cdot
523	<	\mathbf{Q}_{\mathbf{b}}  \cdot \hat{r})=
524	<	D_{\bf {a}\alpha}  Q_{\bf b \alpha \beta} r_\beta
525	<	=r \sum_{lmn} D_{\mathbf{a}l} (\hat{a}_l \cdot \hat{b}_m ) Q_{\mathbf{b}mn}
526	<	(\hat{b}_n \cdot \hat{r}) .
527	<	\end{eqnarray}
522	>	= D_{\bf {a}\alpha} D_{\bf {b}\alpha} =
523	>	\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}}
524	>	\end{equation}
525		%
526	<	[Dan, perhaps there are better examples to show here.]
526	>	The first two forms are written using space coordinates.  The first
527	>	form is standard in the chemistry literature, while the second is
528	>	expressed using implied summation notation.  The third form shows
529	>	explicit sums over body indices and the dot products now indicate
530	>	contractions using space indices.
531
531	–	In each line, the first two terms are written using space coordinates.  The first form is standard
532	–	in the chemistry literature, and the second is ``physicist style'' using implied summation notation.  The third
533	–	form shows explicitly sums over body indices and the dot products now indicate contractions using space indices.
534	–	We find the first form to be useful in writing equations prior to converting to computer code.  The second form is helpful
535	–	in derivations of the interaction energy expressions.  The third one is specifically helpful when deriving forces and torques, as will
536	–	be discussed below.
532
538	–
533		\subsection{The Self-Interaction \label{sec:selfTerm}}
534
535	<	The Wolf summation~\cite{Wolf99} and the later damped shifted force
536	<	(DSF) extension~\cite{Fennell06} included self-interactions that are
537	<	handled separately from the pairwise interactions between sites. The
538	<	self-term is normally calculated via a single loop over all sites in
539	<	the system, and is relatively cheap to evaluate. The self-interaction
540	<	has contributions from two sources:
541	<	\begin{itemize}
542	<	\item The neutralization procedure within the cutoff radius requires a
543	<	contribution from a charge opposite in sign, but equal in magnitude,
544	<	to the central charge, which has been spread out over the surface of
545	<	the cutoff sphere.  For a system of undamped charges, the total
546	<	self-term is
535	>	In addition to cutoff-sphere neutralization, the Wolf
536	>	summation~\cite{Wolf99} and the damped shifted force (DSF)
537	>	extension~\cite{Fennell:2006zl} also included self-interactions that
538	>	are handled separately from the pairwise interactions between
539	>	sites. The self-term is normally calculated via a single loop over all
540	>	sites in the system, and is relatively cheap to evaluate. The
541	>	self-interaction has contributions from two sources.
542	>
543	>	First, the neutralization procedure within the cutoff radius requires
544	>	a contribution from a charge opposite in sign, but equal in magnitude,
545	>	to the central charge, which has been spread out over the surface of
546	>	the cutoff sphere.  For a system of undamped charges, the total
547	>	self-term is
548		\begin{equation}
549		V_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2}
550		\label{eq:selfTerm}
551		\end{equation}
552	<	Note that in this potential and in all electrostatic quantities that
553	<	follow, the standard $4 \pi \epsilon_{0}$ has been omitted for
554	<	clarity.
555	<	\item Charge damping with the complementary error function is a
556	<	partial analogy to the Ewald procedure which splits the interaction
557	<	into real and reciprocal space sums.  The real space sum is retained
558	<	in the Wolf and DSF methods.  The reciprocal space sum is first
559	<	minimized by folding the largest contribution (the self-interaction)
560	<	into the self-interaction from charge neutralization of the damped
561	<	potential.  The remainder of the reciprocal space portion is then
562	<	discarded (as this contributes the largest computational cost and
568	<	complexity to the Ewald sum).  For a system containing only damped
569	<	charges, the complete self-interaction can be written as
552	>
553	>	Second, charge damping with the complementary error function is a
554	>	partial analogy to the Ewald procedure which splits the interaction
555	>	into real and reciprocal space sums.  The real space sum is retained
556	>	in the Wolf and DSF methods.  The reciprocal space sum is first
557	>	minimized by folding the largest contribution (the self-interaction)
558	>	into the self-interaction from charge neutralization of the damped
559	>	potential.  The remainder of the reciprocal space portion is then
560	>	discarded (as this contributes the largest computational cost and
561	>	complexity to the Ewald sum).  For a system containing only damped
562	>	charges, the complete self-interaction can be written as
563		\begin{equation}
564		V_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} +
565		\frac{\alpha}{\sqrt{\pi}} \right) \sum_{{\bf a}=1}^N
566		C_{\bf a}^{2}.
567		\label{eq:dampSelfTerm}
568		\end{equation}
576	–	\end{itemize}
569
570		The extension of DSF electrostatics to point multipoles requires
571		treatment of {\it both} the self-neutralization and reciprocal
#	Line 616 | Line 608 | terms, $f(r), g(r), h(r), s(r), \mathrm{~and~} t(r)$ a
608		multipole orders.  Symmetry prevents the charge-dipole and
609		dipole-quadrupole interactions from contributing to the
610		self-interaction.  The functions that go into the self-neutralization
611	<	terms, $f(r), g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive
612	<	derivatives of the electrostatic kernel (either the undamped $1/r$ or
613	<	the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that are
614	<	evaluated at the cutoff distance.  For undamped interactions, $f(r_c)
615	<	= 1/r_c$, $g(r_c) = -1/r_c^{2}$, and so on.  For damped interactions,
616	<	$f(r_c) = B_0(r_c)$, $g(r_c) = B_0'(r_c)$, and so on.  Appendix XX
617	<	contains recursion relations that allow rapid evaluation of these
618	<	derivatives.
611	>	terms, $g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive
612	>	derivatives of the electrostatic kernel, $f(r)$ (either the undamped
613	>	$1/r$ or the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that
614	>	have been evaluated at the cutoff distance.  For undamped
615	>	interactions, $f(r_c) = 1/r_c$, $g(r_c) = -1/r_c^{2}$, and so on.  For
616	>	damped interactions, $f(r_c) = B_0(r_c)$, $g(r_c) = B_0'(r_c)$, and so
617	>	on.  Appendix \ref{SmithFunc} contains recursion relations that allow
618	>	rapid evaluation of these derivatives.
619
620	<	\section{Energies, forces, and torques}
621	<	\subsection{Interaction energies}
620	>	\section{Interaction energies, forces, and torques}
621	>	The main result of this paper is a set of expressions for the
622	>	energies, forces and torques (up to quadrupole-quadrupole order) that
623	>	work for both the Taylor-shifted and Gradient-shifted approximations.
624	>	These expressions were derived using a set of generic radial
625	>	functions.  Without using the shifting approximations mentioned above,
626	>	some of these radial functions would be identical, and the expressions
627	>	coalesce into the familiar forms for unmodified multipole-multipole
628	>	interactions.  Table \ref{tab:tableenergy} maps between the generic
629	>	functions and the radial functions derived for both the Taylor-shifted
630	>	and Gradient-shifted methods.  The energy equations are written in
631	>	terms of lab-frame representations of the dipoles, quadrupoles, and
632	>	the unit vector connecting the two objects,
633
631	–	We now list multipole interaction energies using a set of generic
632	–	radial functions.  Table \ref{tab:tableenergy} maps between the
633	–	generic functions and the radial functions derived for both the
634	–	Taylor-shifted and Gradient-shifted methods.  This set of equations is
635	–	written in terms of space coordinates:
636	–
634		% Energy in space coordinate form ----------------------------------------------------------------------------------------------
635		%
636		%
#	Line 641 | Line 638 | U_{C_{\bf a}C_{\bf b}}(r)=&
638		%
639		\begin{align}
640		U_{C_{\bf a}C_{\bf b}}(r)=&
641	<	\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r)  \label{uchch}
641	>	C_{\bf a} C_{\bf b}  v_{01}(r)  \label{uchch}
642		\\
643		%
644		% u ca db
645		%
646		U_{C_{\bf a}D_{\bf b}}(r)=&
647	<	\frac{C_{\bf a}}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)  v_{11}(r)
647	>	C_{\bf a} \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)  v_{11}(r)
648		\label{uchdip}
649		\\
650		%
651		% u ca qb
652		%
653	<	U_{C_{\bf a}Q_{\bf b}}(r)=&
654	<	\frac{C_{\bf a }}{4\pi \epsilon_0} \Bigl[ \text{Tr}Q_{\bf b} v_{21}(r)
655	<	\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{22}(r) \Bigr]
653	>	U_{C_{\bf a}Q_{\bf b}}(r)=& C_{\bf a } \Bigl[ \text{Tr}Q_{\bf b}
654	>	v_{21}(r) + \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot
655	>	\hat{r} \right) v_{22}(r) \Bigr]
656		\label{uchquad}
657		\\
658		%
#	Line 669 | Line 666 | -\frac{1}{4\pi \epsilon_0} \Bigr[ \left( \mathbf{D}_{\
666		% u da db
667		%
668		U_{D_{\bf a}D_{\bf b}}(r)=&
669	<	-\frac{1}{4\pi \epsilon_0} \Bigr[ \left( \mathbf{D}_{\mathbf {a}} \cdot
669	>	-\Bigr[ \left( \mathbf{D}_{\mathbf {a}} \cdot
670		\mathbf{D}_{\mathbf{b}} \right)  v_{21}(r)
671		+\left( \mathbf{D}_{\mathbf {a}} \cdot \hat{r} \right)
672		\left( \mathbf{D}_{\mathbf {b}} \cdot \hat{r} \right)
#	Line 682 | Line 679 | -\frac{1}{4\pi \epsilon_0} \Bigl[
679		\begin{split}
680		% 1
681		U_{D_{\bf a}Q_{\bf b}}(r) =&
682	<	-\frac{1}{4\pi \epsilon_0} \Bigl[
682	>	-\Bigl[
683		\text{Tr}\mathbf{Q}_{\mathbf{b}}
684		\left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
685		+2 ( \mathbf{D}_{\mathbf{a}} \cdot
686		\mathbf{Q}_{\mathbf{b}} \cdot \hat{r} ) \Bigr] v_{31}(r) \\
687		% 2
688	<	&-\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
688	>	&- \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
689		\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{32}(r)
690		\label{udipquad}
691		\end{split}
#	Line 725 | Line 722 | U_{Q_{\bf a}Q_{\bf b}}(r)=&
722		\begin{split}
723		%1
724		U_{Q_{\bf a}Q_{\bf b}}(r)=&
725	<	\frac{1}{4\pi \epsilon_0} \Bigl[
725	>	\Bigl[
726		\text{Tr} \mathbf{Q}_{\mathbf{a}} \text{Tr} \mathbf{Q}_{\mathbf{b}}
727		+2 \text{Tr} \left(
728		\mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} \right) \Bigr] v_{41}(r)
729		\\
730		% 2
731	<	&+ \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
731	>	&+\Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
732		\left( \hat{r} \cdot
733		\mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right)
734		+\text{Tr}\mathbf{Q}_{\mathbf{b}}
#	Line 741 | Line 738 | +\text{Tr}\mathbf{Q}_{\mathbf{b}}
738		\Bigr] v_{42}(r)
739		\\
740		% 4
741	<	&+ \frac{1}{4\pi \epsilon_0}
741	>	&+
742		\left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right)
743		\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}}  \cdot \hat{r} \right) v_{43}(r).
744		\label{uquadquad}
745		\end{split}
746		\end{align}
747	<
747	>	%
748		Note that the energies of multipoles on site $\mathbf{b}$ interacting
749		with those on site $\mathbf{a}$ can be obtained by swapping indices
750		along with the sign of the intersite vector, $\hat{r}$.
#	Line 757 | Line 754 | along with the sign of the intersite vector, $\hat{r}$
754		% TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
755		%
756
757	<	\begin{table*}
758	<	\caption{\label{tab:tableenergy}Radial functions used in the energy and torque equations.  Functions
759	<	used in this table are defined in Appendices B and C.}
760	<	\begin{ruledtabular}
761	<	\begin{tabular}{|l|c|l|l}
762	<	Generic&Coulomb&Taylor-Shifted&Gradient-Shifted
757	>	\begin{sidewaystable}
758	>	\caption{\label{tab:tableenergy}Radial functions used in the energy
759	>	and torque equations.  The $f, g, h, s, t, \mathrm{and} u$
760	>	functions used in this table are defined in Appendices B and C.}
761	>	\begin{tabular}{|c|c|l|l|} \hline
762	>	Generic&Bare Coulomb&Taylor-Shifted&Gradient-Shifted
763		\\ \hline
764		%
765		%
#	Line 794 | Line 791 | -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right) $ \\
791		$\frac{3}{r^3}  $ &
792		$\left(-\frac{g_2(r)}{r} + h_2(r) \right)$ &
793		$\left(-\frac{g(r)}{r}+h(r) \right)
794	<	-\left(-\frac{g(r_c)}{r_c}+h(r_c) \right) $ \\
795	<	&&&$  -(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$
794	>	-\left(-\frac{g(r_c)}{r_c}+h(r_c) \right)$  \\
795	>	&&& $ ~~~-(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$
796		\\
797		%
798		%
#	Line 806 | Line 803 | -\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right
803		$\left(-\frac{g_3(r)}{r^2} + \frac{h_3(r)}{r} \right)$ &
804		$\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r} \right)
805		-\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right) $\\
806	<	&&&$  -(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)$
806	>	&&&$ ~~~ -(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)$
807		\\
808		%
809		$v_{32}(r)$ &
#	Line 814 | Line 811 | - \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} +
811		$\left( \frac{3g_3(r)}{r^2} - \frac{3h_3(r)}{r} + s_3(r) \right)$ &
812		$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)
813		- \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c) \right)$ \\
814	<	&&&$  -(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)$
814	>	&&&$ ~~~ -(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)$
815		\\
816		%
817		%
#	Line 825 | Line 822 | - \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2}
822		$\left(-\frac{g_4(r)}{r^3} +\frac{h_4(r)}{r^2} \right) $  &
823		$\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right)
824		- \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2} \right)$ \\
825	<	&&&$  -(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)$
825	>	&&&$ ~~~ -(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)$
826		\\
827		% 2
828		$v_{42}(r)$ &
#	Line 833 | Line 830 | -\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+
830		$\left( \frac{3g_4(r)}{r^3} - \frac{3h_4(r)}{r^2}+\frac{s_4(r)}{r} \right)$ &
831		$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right)
832		-\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ \\
833	<	&&&$  -(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}
833	>	&&&$ ~~~ -(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}
834		-\frac{4s(r_c)}{r_c^2} + \frac{t(r_c)}{r_c}\right)$
835		\\
836		% 3
#	Line 841 | Line 838 | $\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s
838		$ \frac{105}{r^5}  $ &
839		$\left(-\frac{15g_4(r)}{r^3}+\frac{15h_4(r)}{r^2}-\frac{6s_4(r)}{r} + t_4(r)\right) $ &
840		$\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s(r)}{r} + t(r)\right)$ \\
841	<	&&&$ -\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)$ \\
842	<	&&&$ -(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}
843	<	-\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\
841	>	&&&$~~~ -\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)$ \\
842	>	&&&$~~~ -(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}
843	>	-\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\ \hline
844		\end{tabular}
845	<	\end{ruledtabular}
849	<	\end{table*}
845	>	\end{sidewaystable}
846		%
847		%
848		% FORCE  TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
849		%
850
851	<	\begin{table}
851	>	\begin{sidewaystable}
852		\caption{\label{tab:tableFORCE}Radial functions used in the force equations.}
853	<	\begin{ruledtabular}
854	<	\begin{tabular}{cc}
859	<	Generic&Method 1 or Method 2
853	>	\begin{tabular}{|c|c|l|l|} \hline
854	>	Function&Definition&Taylor-Shifted&Gradient-Shifted
855		\\ \hline
856		%
857		%
858		%
859		$w_a(r)$&
860	<	$\frac{d v_{01}}{dr}$ \\
860	>	$\frac{d v_{01}}{dr}$&
861	>	$g_0(r)$&
862	>	$g(r)-g(r_c)$ \\
863		%
864		%
865		$w_b(r)$ &
866	<	$\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $ \\
866	>	$\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $&
867	>	$\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ &
868	>	$h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\
869		%
870		$w_c(r)$ &
871	<	$\frac{v_{11}(r)}{r}$ \\
871	>	$\frac{v_{11}(r)}{r}$ &
872	>	$\frac{g_1(r)}{r} $ &
873	>	$\frac{v_{11}(r)}{r}$\\
874		%
875		%
876		$w_d(r)$&
877	<	$\frac{d v_{21}}{dr}$ \\
877	>	$\frac{d v_{21}}{dr}$&
878	>	$\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ &
879	>	$\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)
880	>	-\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $ \\
881		%
882		$w_e(r)$ &
883	+	$\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ &
884	+	$\frac{v_{22}(r)}{r}$ &
885		$\frac{v_{22}(r)}{r}$ \\
886		%
887		%
888		$w_f(r)$&
889	<	$\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$\\
889	>	$\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$&
890	>	$\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ &
891	>	$ \left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) $  \\
892	>	&&& $ ~~~- \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c)
893	>	\right)-\frac{2v_{22}(r)}{r}$\\
894		%
895		$w_g(r)$&
896	+	$\frac{v_{31}(r)}{r}$&
897	+	$ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$&
898		$\frac{v_{31}(r)}{r}$\\
899		%
900		$w_h(r)$ &
901	<	$\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$\\
901	>	$\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$&
902	>	$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
903	>	$ \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) $ \\
904	>	&&& $ ~~~ -\frac{v_{31}(r)}{r}$ \\
905		% 2
906		$w_i(r)$ &
907	<	$\frac{v_{32}(r)}{r}$ \\
907	>	$\frac{v_{32}(r)}{r}$ &
908	>	$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
909	>	$\frac{v_{32}(r)}{r}$\\
910		%
911		$w_j(r)$ &
912	<	$\frac{d v_{32}}{dr}  - \frac{3v_{32}}{r}$ \\
912	>	$\frac{d v_{32}}{dr}  - \frac{3v_{32}}{r}$&
913	>	$\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right)  $ &
914	>	$\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right)$ \\
915	>	&&& $~~~-\left(\frac{-6g(_cr)}{r_c^3} +\frac{6h(r_c)}{r_c^2}
916	>	-\frac{3s(r_c)}{r_c} +t(r_c) \right) -\frac{3v_{32}}{r}$ \\
917		%
918		$w_k(r)$ &
919	<	$\frac{d v_{41}}{dr} $  \\
919	>	$\frac{d v_{41}}{dr} $ &
920	>	$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
921	>	$\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}  \right)
922	>	-\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2}  \right)$ \\
923		%
924		$w_l(r)$ &
925	<	$\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ \\
925	>	$\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ &
926	>	$\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)$ &
927	>	$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\
928	>	&&& $~~~ -\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)
929	>	-\frac{2v_{42}(r)}{r}$\\
930		%
931		$w_m(r)$ &
932	<	$\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$ \\
932	>	$\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$&
933	>	$\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)$ &
934	>	$\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$\\
935	>	&&& $~~~- \left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}
936	>	+\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $\\
937	>	&&& $~~~-\frac{4v_{43}(r)}{r}$  \\
938		%
939		$w_n(r)$ &
940	<	$\frac{v_{42}(r)}{r}$ \\
940	>	$\frac{v_{42}(r)}{r}$ &
941	>	$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
942	>	$\frac{v_{42}(r)}{r}$\\
943		%
944		$w_o(r)$ &
945	<	$\frac{v_{43}(r)}{r}$ \\
945	>	$\frac{v_{43}(r)}{r}$&
946	>	$\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)$ &
947	>	$\frac{v_{43}(r)}{r}$  \\ \hline
948		%
949
950		\end{tabular}
951	<	\end{ruledtabular}
915	<	\end{table}
951	>	\end{sidewaystable}
952		%
953		%
954		%
955
956		\subsection{Forces}
957	<
958	<	The force $\mathbf{F}_{\bf a}$ on $\bf{a}$ due to $\bf{b}$ is the negative of
959	<	the force  $\mathbf{F}_{\bf b}$ on $\bf{b}$ due to $\bf{a}$.  For a simple charge-charge
960	<	interaction, these forces will point along the $\pm \hat{r}$ directions, where
961	<	$\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $.  Thus
957	>	The force on object $\bf{a}$, $\mathbf{F}_{\bf a}$, due to object
958	>	$\bf{b}$ is the negative of the force on $\bf{b}$ due to $\bf{a}$. For
959	>	a simple charge-charge interaction, these forces will point along the
960	>	$\pm \hat{r}$ directions, where $\mathbf{r}=\mathbf{r}_b -
961	>	\mathbf{r}_a $.  Thus
962		%
963		\begin{equation}
964		F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}}{\partial r}
#	Line 930 | Line 966 | The concept of obtaining a force from an energy by tak
966		= - \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}} {\partial r}  .
967		\end{equation}
968		%
969	<	The concept of obtaining a force from an energy by taking a gradient is the same for
970	<	higher-order multipole interactions, the trick is to make sure that all
971	<	$r$-dependent derivatives are considered.
972	<	As is pointed out by Allen and Germano,\cite{Allen:2006fk} this is straightforward if the
973	<	interaction energies are written recognizing explicit
974	<	$\hat{r}$ and body axes ($\hat{a}_m$, $\hat{b}_n$) dependences:
969	>	Obtaining the force from the interaction energy expressions is the
970	>	same for higher-order multipole interactions -- the trick is to make
971	>	sure that all $r$-dependent derivatives are considered.  This is
972	>	straighforward if the interaction energies are written explicitly in
973	>	terms of $\hat{r}$ and the body axes ($\hat{a}_m$,
974	>	$\hat{b}_n$) :
975		%
976		\begin{equation}
977		U(r,\{\hat{a}_m \cdot \hat{r} \},
978	<	\{\hat{b}_n\cdot \hat{r} \}
978	>	\{\hat{b}_n\cdot \hat{r} \},
979		\{\hat{a}_m \cdot \hat{b}_n \}) .
980		\label{ugeneral}
981		\end{equation}
982		%
983	<	Then,
983	>	Allen and Germano,\cite{Allen:2006fk} showed that if the energy is
984	>	written in this form, the forces come out relatively cleanly,
985		%
986		\begin{equation}
987		\mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} =  \frac{\partial U}{\partial \mathbf{r}}
#	Line 957 | Line 994 | Note our definition of $\mathbf{r}=\mathbf{r}_b - \mat
994		\right] \label{forceequation}.
995		\end{equation}
996		%
997	<	Note our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ is opposite
998	<	that of Allen and Germano.\cite{Allen:2006fk}  In simplifying the algebra, we also use:
997	>	Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $
998	>	is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk}
999	>	In simplifying the algebra, we have also used:
1000		%
1001	<	\begin{eqnarray}
1001	>	\begin{align}
1002		\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}}
1003	<	= \frac{1}{r} \left(  \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r}
1003	>	=& \frac{1}{r} \left(  \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r}
1004		\right) \\
1005		\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}}
1006	<	= \frac{1}{r} \left(  \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r}
1006	>	=& \frac{1}{r} \left(  \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r}
1007		\right) .
1008	<	\end{eqnarray}
1008	>	\end{align}
1009		%
1010	<	We list below the force equations written in terms  of space coordinates.  The
1011	<	radial functions used in the two methods are listed in Table II.
1010	>	We list below the force equations written in terms of lab-frame
1011	>	coordinates.  The radial functions used in the two methods are listed
1012	>	in Table \ref{tab:tableFORCE}
1013		%
1014	<	%SPACE COORDINATES FORCE EQUTIONS
1014	>	%SPACE COORDINATES FORCE EQUATIONS
1015		%
1016		% **************************************************************************
1017		% f ca cb
1018		%
1019	<	\begin{equation}
1020	<	\mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =
1021	<	\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r}
983	<	\end{equation}
1019	>	\begin{align}
1020	>	\mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =&
1021	>	C_{\bf a} C_{\bf b}  w_a(r) \hat{r} \\
1022		%
1023		%
1024		%
1025	<	\begin{equation}
1026	<	\mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =
989	<	\frac{C_{\bf a}}{4\pi \epsilon_0} \Bigl[
1025	>	\mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =&
1026	>	C_{\bf a} \Bigl[
1027		\left( \hat{r} \cdot \mathbf{D}_{\mathbf{b}} \right)
1028		w_b(r) \hat{r}
1029	<	+ \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr]
993	<	\end{equation}
1029	>	+ \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr] \\
1030		%
1031		%
1032		%
1033	<	\begin{equation}
1034	<	\mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =
999	<	\frac{C_{\bf a }}{4\pi \epsilon_0} \Bigr[
1033	>	\mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =&
1034	>	C_{\bf a } \Bigr[
1035		\text{Tr}\mathbf{Q}_{\bf b} w_d(r) \hat{r}
1036		+ 2  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} w_e(r)
1037	<	+ \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
1038	<	\end{equation}
1037	>	+ \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}
1038	>	\right) w_f(r) \hat{r} \Bigr] \\
1039		%
1040		%
1041		%
1042	<	\begin{equation}
1043	<	\mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1044	<	-\frac{C_{\bf{b}}}{4\pi \epsilon_0} \Bigl[
1045	<	\left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) w_b(r) \hat{r}
1046	<	+ \mathbf{D}_{\mathbf{a}} w_c(r) \Bigr]
1047	<	\end{equation}
1042	>	% \begin{equation}
1043	>	% \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1044	>	% -C_{\bf{b}} \Bigl[
1045	>	% \left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) w_b(r) \hat{r}
1046	>	% + \mathbf{D}_{\mathbf{a}} w_c(r) \Bigr]
1047	>	% \end{equation}
1048		%
1049		%
1050		%
1051	<	\begin{equation}
1052	<	\mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} =
1018	<	\frac{1}{4\pi \epsilon_0} \Bigl[
1051	>	\begin{split}
1052	>	\mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} =&
1053		- \mathbf{D}_{\mathbf {a}} \cdot  \mathbf{D}_{\mathbf{b}} w_d(r) \hat{r}
1054		+ \left( \mathbf{D}_{\mathbf {a}}
1055		\left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)
1056	<	+ \mathbf{D}_{\mathbf {b}} \left( \mathbf{D}_{\mathbf{a}}  \cdot \hat{r} \right) \right) w_e(r)
1056	>	+ \mathbf{D}_{\mathbf {b}} \left( \mathbf{D}_{\mathbf{a}}  \cdot \hat{r} \right) \right) w_e(r)\\
1057		% 2
1058	<	- \left( \hat{r} \cdot \mathbf{D}_{\mathbf {a}} \right)
1059	<	\left( \hat{r} \cdot \mathbf{D}_{\mathbf {b}} \right) w_f(r) \hat{r} \Bigr]
1060	<	\end{equation}
1058	>	& - \left( \hat{r} \cdot \mathbf{D}_{\mathbf {a}} \right)
1059	>	\left( \hat{r} \cdot \mathbf{D}_{\mathbf {b}} \right) w_f(r) \hat{r}
1060	>	\end{split}\\
1061		%
1062		%
1063		%
1030	–	\begin{equation}
1064		\begin{split}
1065	<	\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1033	<	& - \frac{1}{4\pi \epsilon_0} \Bigl[
1065	>	\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =& - \Bigl[
1066		\text{Tr}\mathbf{Q}_{\mathbf{b}} \mathbf{ D}_{\mathbf{a}}
1067		+2 \mathbf{D}_{\mathbf{a}} \cdot
1068		\mathbf{Q}_{\mathbf{b}} \Bigr] w_g(r)
1069	<	- \frac{1}{4\pi \epsilon_0} \Bigl[
1069	>	- \Bigl[
1070		\text{Tr}\mathbf{Q}_{\mathbf{b}}
1071		\left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right)
1072		+2 ( \mathbf{D}_{\mathbf{a}} \cdot
1073		\mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
1074		% 3
1075	<	& - \frac{1}{4\pi \epsilon_0} \Bigl[\mathbf{ D}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
1075	>	& - \Bigl[\mathbf{ D}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
1076		+2  (\hat{r} \cdot \mathbf{D}_{\mathbf{a}} ) (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} )  \Bigr]
1077		w_i(r)
1078		% 4
1079	<	-\frac{1}{4\pi \epsilon_0}
1079	>	-
1080		(\hat{r} \cdot \mathbf{D}_{\mathbf{a}} )
1081	<	(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r}
1050	<	\end{split}
1051	<	\end{equation}
1081	>	(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r} \end{split} \\
1082		%
1083		%
1084	<	\begin{equation}
1085	<	\mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} =
1086	<	\frac{C_{\bf b }}{4\pi \epsilon_0} \Bigr[
1087	<	\text{Tr}\mathbf{Q}_{\bf a} w_d(r) \hat{r}
1088	<	+ 2  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} w_e(r)
1089	<	+ \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
1090	<	\end{equation}
1084	>	% \begin{equation}
1085	>	% \mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} =
1086	>	% \frac{C_{\bf b }}{4\pi \epsilon_0} \Bigr[
1087	>	% \text{Tr}\mathbf{Q}_{\bf a} w_d(r) \hat{r}
1088	>	% + 2  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} w_e(r)
1089	>	%  + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
1090	>	% \end{equation}
1091	>	% %
1092	>	% \begin{equation}
1093	>	% \begin{split}
1094	>	% \mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1095	>	% &\frac{1}{4\pi \epsilon_0} \Bigl[
1096	>	% \text{Tr}\mathbf{Q}_{\mathbf{a}} \mathbf{D}_{\mathbf{b}}
1097	>	% +2 \mathbf{D}_{\mathbf{b}} \cdot \mathbf{Q}_{\mathbf{a}}  \Bigr] w_g(r)
1098	>	% % 2
1099	>	% + \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
1100	>	% (\hat{r} \cdot  \mathbf{D}_{\mathbf{b}})
1101	>	% +2 (\mathbf{D}_{\mathbf{b}} \cdot
1102	>	% \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
1103	>	% % 3
1104	>	% &+ \frac{1}{4\pi \epsilon_0} \Bigl[ \mathbf{D}_{\mathbf{b}}
1105	>	% (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1106	>	% +2  (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1107	>	% (\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} ) \Bigr]   w_i(r)
1108	>	% % 4
1109	>	% +\frac{1}{4\pi \epsilon_0}
1110	>	% (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1111	>	% (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}}  \cdot \hat{r}) w_j(r) \hat{r}
1112	>	% \end{split}
1113	>	% \end{equation}
1114		%
1062	–	\begin{equation}
1063	–	\begin{split}
1064	–	\mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1065	–	&\frac{1}{4\pi \epsilon_0} \Bigl[
1066	–	\text{Tr}\mathbf{Q}_{\mathbf{a}} \mathbf{D}_{\mathbf{b}}
1067	–	+2 \mathbf{D}_{\mathbf{b}} \cdot \mathbf{Q}_{\mathbf{a}}  \Bigr] w_g(r)
1068	–	% 2
1069	–	+ \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
1070	–	(\hat{r} \cdot  \mathbf{D}_{\mathbf{b}})
1071	–	+2 (\mathbf{D}_{\mathbf{b}} \cdot
1072	–	\mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
1073	–	% 3
1074	–	&+ \frac{1}{4\pi \epsilon_0} \Bigl[ \mathbf{D}_{\mathbf{b}}
1075	–	(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1076	–	+2  (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1077	–	(\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} ) \Bigr]   w_i(r)
1078	–	% 4
1079	–	+\frac{1}{4\pi \epsilon_0}
1080	–	(\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1081	–	(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}}  \cdot \hat{r}) w_j(r) \hat{r}
1082	–	\end{split}
1083	–	\end{equation}
1115		%
1116		%
1086	–	%
1087	–	\begin{equation}
1117		\begin{split}
1118	<	\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1119	<	+\frac{1}{4\pi \epsilon_0} \Bigl[
1118	>	\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =&
1119	>	\Bigl[
1120		\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} \hat{r}
1121		+ 2 \text{Tr} ( \mathbf{Q}_{\mathbf{a}} \cdot  \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_k(r) \hat{r} \\
1122		% 2
1123	<	+\frac{1}{4\pi \epsilon_0} \Bigl[
1123	>	&+ \Bigl[
1124		2\text{Tr}\mathbf{Q}_{\mathbf{b}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )
1125		+ 2\text{Tr}\mathbf{Q}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} )
1126		% 3
1127		+4 (\mathbf{Q}_{\mathbf{a}}  \cdot  \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1128		+  4(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}}) \Bigr] w_n(r) \\
1129		% 4
1130	<	+ \frac{1}{4\pi \epsilon_0} \Bigl[
1130	>	&+  \Bigl[
1131		\text{Tr}\mathbf{Q}_{\mathbf{a}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1132		+ \text{Tr}\mathbf{Q}_{\mathbf{b}}
1133		(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
#	Line 1106 | Line 1135 | + \frac{1}{4\pi \epsilon_0} \Bigl[
1135		+4 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot
1136		\mathbf{Q}_{\mathbf{b}}   \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
1137		%
1138	<	+ \frac{1}{4\pi \epsilon_0} \Bigl[
1138	>	&+ \Bigl[
1139		+ 2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )
1140		(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1141		%6
1142		+2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \hat{r})
1143		(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_o(r) \\
1144		%  7
1145	<	+ \frac{1}{4\pi \epsilon_0}
1145	>	&+
1146		(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
1147	<	(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) w_m(r) \hat{r}
1148	<	\end{split}
1149	<	\end{equation}
1147	>	(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) w_m(r) \hat{r} \end{split}
1148	>	\end{align}
1149	>	Note that the forces for higher multipoles on site $\mathbf{a}$
1150	>	interacting with those of lower order on site $\mathbf{b}$ can be
1151	>	obtained by swapping indices in the expressions above.
1152	>
1153		%
1154	+	% Torques SECTION -----------------------------------------------------------------------------------------
1155		%
1123	–	% TORQUES SECTION -----------------------------------------------------------------------------------------
1124	–	%
1156		\subsection{Torques}
1157	<
1158	<	Following again Allen and Germano,\cite{Allen:2006fk} when energies are written in the form
1159	<	of Eq.~({\ref{ugeneral}), then torques can be expressed as:
1157	>	When energies are written in the form of Eq.~({\ref{ugeneral}), then
1158	>	torques can be found in a relatively straightforward
1159	>	manner,\cite{Allen:2006fk}
1160		%
1161		\begin{eqnarray}
1162		\mathbf{\tau}_{\bf a} =
#	Line 1146 | Line 1177 | Here we list the torque equations written in terms of
1177		\end{eqnarray}
1178		%
1179		%
1180	<	Here we list the torque equations written in terms of space coordinates.
1180	>	The torques for both the Taylor-Shifted as well as Gradient-Shifted
1181	>	methods are given in space-frame coordinates:
1182		%
1183		%
1184	+	\begin{align}
1185	+	\mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =&
1186	+	C_{\bf a}  (\hat{r} \times  \mathbf{D}_{\mathbf{b}}) v_{11}(r) \\
1187		%
1153	–	\begin{equation}
1154	–	\mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =
1155	–	\frac{C_{\bf a}}{4\pi \epsilon_0}  (\hat{r} \times  \mathbf{D}_{\mathbf{b}}) v_{11}(r)
1156	–	\end{equation}
1188		%
1189		%
1190	+	\mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =&
1191	+	2C_{\bf a}
1192	+	\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{22}(r) \\
1193		%
1160	–	\begin{equation}
1161	–	\mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =
1162	–	\frac{2C_{\bf a}}{4\pi \epsilon_0}
1163	–	\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{22}(r)
1164	–	\end{equation}
1194		%
1195		%
1196	+	% \begin{equation}
1197	+	% \mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =
1198	+	% -\frac{C_{\bf b}}{4\pi \epsilon_0}
1199	+	% (\hat{r} \times \mathbf{D}_{\mathbf{a}})  v_{11}(r)
1200	+	% \end{equation}
1201		%
1168	–	\begin{equation}
1169	–	\mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =
1170	–	-\frac{C_{\bf b}}{4\pi \epsilon_0}
1171	–	(\hat{r} \times \mathbf{D}_{\mathbf{a}})  v_{11}(r)
1172	–	\end{equation}
1202		%
1203		%
1204	<	%
1205	<	\begin{equation}
1177	<	\mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} =
1178	<	\frac{1}{4\pi \epsilon_0}  \mathbf{D}_{\mathbf {a}}  \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1204	>	\mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} =&
1205	>	\mathbf{D}_{\mathbf {a}}  \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1206		% 2
1207	<	-\frac{1}{4\pi \epsilon_0}
1207	>	-
1208		(\hat{r} \times \mathbf{D}_{\mathbf {a}} )
1209	<	(\hat{r} \cdot \mathbf{D}_{\mathbf {b}} )  v_{22}(r)
1183	<	\end{equation}
1209	>	(\hat{r} \cdot \mathbf{D}_{\mathbf {b}} )  v_{22}(r)\\
1210		%
1211		%
1212		%
1213	<	\begin{equation}
1214	<	\mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} =
1215	<	-\frac{1}{4\pi \epsilon_0} \mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1216	<	% 2
1217	<	+\frac{1}{4\pi \epsilon_0}
1218	<	(\hat{r} \cdot \mathbf{D}_{\mathbf {a}} )
1219	<	(\hat{r} \times \mathbf{D}_{\mathbf {b}} ) v_{22}(r)
1220	<	\end{equation}
1213	>	% \begin{equation}
1214	>	% \mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} =
1215	>	% -\frac{1}{4\pi \epsilon_0} \mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1216	>	% % 2
1217	>	% +\frac{1}{4\pi \epsilon_0}
1218	>	% (\hat{r} \cdot \mathbf{D}_{\mathbf {a}} )
1219	>	% (\hat{r} \times \mathbf{D}_{\mathbf {b}} ) v_{22}(r)
1220	>	% \end{equation}
1221		%
1222		%
1223		%
1224	<	\begin{equation}
1225	<	\mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1200	<	\frac{1}{4\pi \epsilon_0} \Bigl[
1224	>	\mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} =&
1225	>	\Bigl[
1226		-\text{Tr}\mathbf{Q}_{\mathbf{b}}
1227		(\hat{r} \times \mathbf{D}_{\mathbf{a}} )
1228		+2 \mathbf{D}_{\mathbf{a}}  \times
1229		(\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1230		\Bigr] v_{31}(r)
1231		% 3
1232	<	-\frac{1}{4\pi \epsilon_0}
1233	<	\ (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
1209	<	(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)
1210	<	\end{equation}
1232	>	- (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
1233	>	(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)\\
1234		%
1235		%
1236		%
1237	<	\begin{equation}
1238	<	\mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =
1216	<	\frac{1}{4\pi \epsilon_0} \Bigl[
1237	>	\mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =&
1238	>	\Bigl[
1239		+2 ( \mathbf{D}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} ) \times
1240		\hat{r}
1241		-2 \mathbf{D}_{\mathbf{a}}  \times
1242		(\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1243		\Bigr] v_{31}(r)
1244		% 2
1245	<	+\frac{2}{4\pi \epsilon_0}
1245	>	+
1246		(\hat{r} \cdot \mathbf{D}_{\mathbf{a}})
1247	<	(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)
1226	<	\end{equation}
1247	>	(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)\\
1248		%
1249		%
1250		%
1251	<	\begin{equation}
1252	<	\mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1253	<	\frac{1}{4\pi \epsilon_0} \Bigl[
1254	<	-2 (\mathbf{D}_{\mathbf{b}}  \cdot \mathbf{Q}_{\mathbf{a}} ) \times \hat{r}
1255	<	+2 \mathbf{D}_{\mathbf{b}}  \times
1256	<	(\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
1257	<	\Bigr] v_{31}(r)
1258	<	% 3
1259	<	- \frac{2}{4\pi \epsilon_0}
1260	<	(\hat{r} \cdot \mathbf{D}_{\mathbf{b}} )
1261	<	(\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}}) \times \hat{r} v_{32}(r)
1262	<	\end{equation}
1251	>	% \begin{equation}
1252	>	% \mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1253	>	% \frac{1}{4\pi \epsilon_0} \Bigl[
1254	>	% -2 (\mathbf{D}_{\mathbf{b}}  \cdot \mathbf{Q}_{\mathbf{a}} ) \times \hat{r}
1255	>	% +2 \mathbf{D}_{\mathbf{b}}  \times
1256	>	% (\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
1257	>	% \Bigr] v_{31}(r)
1258	>	% % 3
1259	>	% - \frac{2}{4\pi \epsilon_0}
1260	>	% (\hat{r} \cdot \mathbf{D}_{\mathbf{b}} )
1261	>	% (\hat{r} \cdot  \mathbf
1262	>	% {Q}_{{\mathbf a}}) \times \hat{r} v_{32}(r)
1263	>	% \end{equation}
1264		%
1265		%
1266		%
1267	<	\begin{equation}
1268	<	\mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} =
1269	<	\frac{1}{4\pi \epsilon_0} \Bigl[
1270	<	\text{Tr}\mathbf{Q}_{\mathbf{a}}
1271	<	(\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1272	<	+2 \mathbf{D}_{\mathbf{b}}  \times
1273	<	( \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r)
1274	<	% 2
1275	<	+\frac{1}{4\pi \epsilon_0}
1276	<	(\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1277	<	(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) v_{32}(r)
1278	<	\end{equation}
1267	>	% \begin{equation}
1268	>	% \mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} =
1269	>	% \frac{1}{4\pi \epsilon_0} \Bigl[
1270	>	% \text{Tr}\mathbf{Q}_{\mathbf{a}}
1271	>	% (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1272	>	% +2 \mathbf{D}_{\mathbf{b}}  \times
1273	>	% ( \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r)
1274	>	% % 2
1275	>	% +\frac{1}{4\pi \epsilon_0}
1276	>	% (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1277	>	% (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) v_{32}(r)
1278	>	% \end{equation}
1279		%
1280		%
1281		%
1260	–	\begin{equation}
1282		\begin{split}
1283	<	\mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1284	<	&-\frac{4}{4\pi \epsilon_0}
1283	>	\mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =&
1284	>	-4
1285		\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}}
1286		v_{41}(r) \\
1287		% 2
1288	<	&+ \frac{1}{4\pi \epsilon_0}
1288	>	&+
1289		\Bigl[-2\text{Tr}\mathbf{Q}_{\mathbf{b}}
1290		(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times \hat{r}
1291		+4 \hat{r} \times
#	Line 1273 | Line 1294 | -4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} )\times
1294		-4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} )\times
1295		( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} ) \Bigr] v_{42}(r) \\
1296		% 4
1297	<	&+ \frac{2}{4\pi \epsilon_0}
1297	>	&+ 2
1298		\hat{r} \times ( \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1299	<	(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
1279	<	\end{split}
1280	<	\end{equation}
1299	>	(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r) \end{split}\\
1300		%
1301		%
1302		%
1284	–	\begin{equation}
1303		\begin{split}
1304		\mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} =
1305	<	&\frac{4}{4\pi \epsilon_0}
1305	>	&4
1306		\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}} v_{41}(r) \\
1307		% 2
1308	<	&+ \frac{1}{4\pi \epsilon_0} \Bigl[- 2\text{Tr}\mathbf{Q}_{\mathbf{a}}
1308	>	&+  \Bigl[- 2\text{Tr}\mathbf{Q}_{\mathbf{a}}
1309		(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} ) \times \hat{r}
1310		-4  (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot
1311		\mathbf{Q}_{{\mathbf b}} ) \times
#	Line 1296 | Line 1314 | +4 ( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times
1314		( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
1315		\Bigr] v_{42}(r) \\
1316		% 4
1317	<	&+ \frac{2}{4\pi \epsilon_0}
1317	>	&+2
1318		(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1319	<	\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
1320	<	\end{split}
1303	<	\end{equation}
1319	>	\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)\end{split}
1320	>	\end{align}
1321		%
1322	+	Here, we have defined the matrix cross product in an identical form
1323	+	as in Ref. \onlinecite{Smith98}:
1324	+	\begin{equation}
1325	+	\left[\mathbf{A} \times \mathbf{B}\right]_\alpha = \sum_\beta
1326	+	\left[\mathbf{A}_{\alpha+1,\beta} \mathbf{B}_{\alpha+2,\beta}
1327	+	-\mathbf{A}_{\alpha+2,\beta} \mathbf{B}_{\alpha+2,\beta}
1328	+	\right]
1329	+	\end{equation}
1330	+	where $\alpha+1$ and $\alpha+2$ are regarded as cyclic
1331	+	permuations of the matrix indices.
1332
1333	+	All of the radial functions required for torques are identical with
1334	+	the radial functions previously computed for the interaction energies.
1335	+	These are tabulated for both shifted force methods in table
1336	+	\ref{tab:tableenergy}.  The torques for higher multipoles on site
1337	+	$\mathbf{a}$ interacting with those of lower order on site
1338	+	$\mathbf{b}$ can be obtained by swapping indices in the expressions
1339	+	above.
1340	+
1341		\section{Comparison to known multipolar energies}
1342
1343		To understand how these new real-space multipole methods behave in
#	Line 1394 | Line 1429 | where $Q$ is the quadrupole moment.
1429		\end{equation}
1430		where $Q$ is the quadrupole moment.
1431
1432	+	\section{Conclusion}
1433	+	We have presented two efficient real-space methods for computing the
1434	+	interactions between point multipoles.  These methods have the benefit
1435	+	of smoothly truncating the energies, forces, and torques at the cutoff
1436	+	radius, making them attractive for both molecular dynamics (MD) and
1437	+	Monte Carlo (MC) simulations.  We find that the Gradient-Shifted Force
1438	+	(GSF) and the Shifted-Potential (SP) methods converge rapidly to the
1439	+	correct lattice energies for ordered dipolar and quadrupolar arrays,
1440	+	while the Taylor-Shifted Force (TSF) is too severe an approximation to
1441	+	provide accurate convergence to lattice energies.
1442
1443	<
1444	<
1445	<
1446	<
1443	>	In most cases, GSF can obtain nearly quantitative agreement with the
1444	>	lattice energy constants with reasonably small cutoff radii.  The only
1445	>	exception we have observed is for crystals which exhibit a bulk
1446	>	macroscopic dipole moment (e.g. Luttinger \& Tisza's $Z_1$ lattice).
1447	>	In this particular case, the multipole neutralization scheme can
1448	>	interfere with the correct computation of the energies.  We note that
1449	>	the energies for these arrangements are typically much larger than for
1450	>	crystals with net-zero moments, so this is not expected to be an issue
1451	>	in most simulations.
1452
1453	+	In large systems, these new methods can be made to scale approximately
1454	+	linearly with system size, and detailed comparisons with the Ewald sum
1455	+	for a wide range of chemical environments follows in the second paper.
1456
1457		\begin{acknowledgments}
1458	<	Support for this project was provided by the National Science
1459	<	Foundation under grant CHE-0848243. Computational time was provided by
1460	<	the Center for Research Computing (CRC) at the University of Notre
1461	<	Dame.
1458	>	JDG acknowledges helpful discussions with Christopher
1459	>	Fennell. Support for this project was provided by the National
1460	>	Science Foundation under grant CHE-0848243. Computational time was
1461	>	provided by the Center for Research Computing (CRC) at the
1462	>	University of Notre Dame.
1463		\end{acknowledgments}
1464
1465		\newpage
1466		\appendix
1467
1468		\section{Smith's $B_l(r)$ functions for damped-charge distributions}
1469	<
1469	>	\label{SmithFunc}
1470		The following summarizes Smith's $B_l(r)$ functions and includes
1471		formulas given in his appendix.\cite{Smith98} The first function
1472		$B_0(r)$ is defined by
#	Line 1576 | Line 1630 | $n$ is eliminated.
1630		\ref{eq:b13} are still correct for GSF electrostatics if the subscript
1631		$n$ is eliminated.
1632
1633	<	\section{Extra Material}
1634	<	%
1635	<	%
1636	<	%Energy in body coordinate form ---------------------------------------------------------------
1637	<	%
1638	<	Here are the interaction energies written in terms of the body coordinates:
1585	<
1586	<	%
1587	<	% u ca cb
1588	<	%
1589	<	\begin{equation}
1590	<	U_{C_{\bf a}C_{\bf b}}(r)=
1591	<	\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r)
1592	<	\end{equation}
1593	<	%
1594	<	% u ca db
1595	<	%
1596	<	\begin{equation}
1597	<	U_{C_{\bf a}D_{\bf b}}(r)=
1598	<	\frac{C_{\bf a}}{4\pi \epsilon_0}
1599	<	\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} \,  v_{11}(r)
1600	<	\end{equation}
1601	<	%
1602	<	% u ca qb
1603	<	%
1604	<	\begin{equation}
1605	<	U_{C_{\bf a}Q_{\bf b}}(r)=
1606	<	\frac{C_{\bf a }\text{Tr}Q_{\bf b}}{4\pi \epsilon_0}
1607	<	v_{21}(r) \nonumber \\
1608	<	+\frac{C_{\bf a}}{4\pi \epsilon_0}
1609	<	\sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r})
1610	<	v_{22}(r)
1611	<	\end{equation}
1612	<	%
1613	<	% u da cb
1614	<	%
1615	<	\begin{equation}
1616	<	U_{D_{\bf a}C_{\bf b}}(r)=
1617	<	-\frac{C_{\bf b}}{4\pi \epsilon_0}
1618	<	\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} \,  v_{11}(r)
1619	<	\end{equation}
1620	<	%
1621	<	% u da db
1622	<	%
1623	<	\begin{equation}
1624	<	\begin{split}
1625	<	% 1
1626	<	U_{D_{\bf a}D_{\bf b}}(r)&=
1627	<	-\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
1628	<	(\hat{a}_m \cdot   \hat{b}_n)
1629	<	D_{\mathbf{b}n} v_{21}(r) \\
1630	<	% 2
1631	<	&-\frac{1}{4\pi \epsilon_0}
1632	<	\sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1633	<	\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n}
1634	<	v_{22}(r)
1635	<	\end{split}
1636	<	\end{equation}
1637	<	%
1638	<	% u da qb
1639	<	%
1640	<	\begin{equation}
1641	<	\begin{split}
1642	<	% 1
1643	<	U_{D_{\bf a}Q_{\bf b}}(r)&=
1644	<	-\frac{1}{4\pi \epsilon_0} \left(
1645	<	\text{Tr}Q_{\mathbf{b}}
1646	<	\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
1647	<	+2\sum_{lmn}D_{\mathbf{a}l}
1648	<	(\hat{a}_l \cdot \hat{b}_m)
1649	<	Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
1650	<	\right)  v_{31}(r) \\
1651	<	% 2
1652	<	&-\frac{1}{4\pi \epsilon_0}
1653	<	\sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
1654	<	\sum_{mn} (\hat{r} \cdot \hat{b}_m)
1655	<	Q_{{\mathbf b}mn}
1656	<	(\hat{b}_n \cdot \hat{r}) v_{32}(r)
1657	<	\end{split}
1658	<	\end{equation}
1659	<	%
1660	<	% u qa cb
1661	<	%
1662	<	\begin{equation}
1663	<	U_{Q_{\bf a}C_{\bf b}}(r)=
1664	<	\frac{C_{\bf b }\text{Tr}Q_{\bf a}}{4\pi \epsilon_0}  v_{21}(r)
1665	<	+\frac{C_{\bf b}}{4\pi \epsilon_0}
1666	<	\sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) v_{22}(r)
1667	<	\end{equation}
1668	<	%
1669	<	% u qa db
1670	<	%
1671	<	\begin{equation}
1672	<	\begin{split}
1673	<	%1
1674	<	U_{Q_{\bf a}D_{\bf b}}(r)&=
1675	<	\frac{1}{4\pi \epsilon_0} \left(
1676	<	\text{Tr}Q_{\mathbf{a}}
1677	<	\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n}
1678	<	+2\sum_{lmn}D_{\mathbf{b}l}
1679	<	(\hat{b}_l \cdot \hat{a}_m)
1680	<	Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
1681	<	\right) v_{31}(r)  \\
1682	<	% 2
1683	<	&+\frac{1}{4\pi \epsilon_0}
1684	<	\sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
1685	<	\sum_{mn} (\hat{r} \cdot \hat{a}_m)
1686	<	Q_{{\mathbf a}mn}
1687	<	(\hat{a}_n \cdot \hat{r}) v_{32}(r)
1688	<	\end{split}
1689	<	\end{equation}
1690	<	%
1691	<	% u qa qb
1692	<	%
1693	<	\begin{equation}
1694	<	\begin{split}
1695	<	%1
1696	<	U_{Q_{\bf a}Q_{\bf b}}(r)&=
1697	<	\frac{1}{4\pi \epsilon_0} \Bigl[
1698	<	\text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}}
1699	<	+2\sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
1700	<	Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
1701	<	(\hat{a}_m \cdot \hat{b}_n) \Bigr]
1702	<	v_{41}(r) \\
1703	<	% 2
1704	<	&+ \frac{1}{4\pi \epsilon_0}
1705	<	\Bigl[ \text{Tr}Q_{\mathbf{a}}
1706	<	\sum_{lm} (\hat{r} \cdot \hat{b}_l )
1707	<	Q_{{\mathbf b}lm}
1708	<	(\hat{b}_m \cdot \hat{r})
1709	<	+\text{Tr}Q_{\mathbf{b}}
1710	<	\sum_{lm} (\hat{r} \cdot \hat{a}_l )
1711	<	Q_{{\mathbf a}lm}
1712	<	(\hat{a}_m \cdot \hat{r}) \\
1713	<	% 3
1714	<	&+4 \sum_{lmnp}
1715	<	(\hat{r} \cdot \hat{a}_l )
1716	<	Q_{{\mathbf a}lm}
1717	<	(\hat{a}_m \cdot \hat{b}_n)
1718	<	Q_{{\mathbf b}np}
1719	<	(\hat{b}_p \cdot \hat{r})
1720	<	\Bigr] v_{42}(r)  \\
1721	<	% 4
1722	<	&+ \frac{1}{4\pi \epsilon_0}
1723	<	\sum_{lm} (\hat{r} \cdot \hat{a}_l)
1724	<	Q_{{\mathbf a}lm}
1725	<	(\hat{a}_m \cdot \hat{r})
1726	<	\sum_{np}  (\hat{r} \cdot \hat{b}_n)
1727	<	Q_{{\mathbf b}np}
1728	<	(\hat{b}_p \cdot \hat{r})  v_{43}(r).
1729	<	\end{split}
1730	<	\end{equation}
1731	<	%
1633	>	% \section{Extra Material}
1634	>	% %
1635	>	% %
1636	>	% %Energy in body coordinate form ---------------------------------------------------------------
1637	>	% %
1638	>	% Here are the interaction energies written in terms of the body coordinates:
1639
1640	<
1641	<	% BODY coordinates force equations --------------------------------------------
1642	<	%
1643	<	%
1644	<	Here are the force equations written in terms of body coordinates.
1645	<	%
1646	<	% f ca cb
1647	<	%
1648	<	\begin{equation}
1649	<	\mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =
1650	<	\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r}
1651	<	\end{equation}
1652	<	%
1653	<	% f ca db
1654	<	%
1655	<	\begin{equation}
1656	<	\mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =
1657	<	\frac{C_{\bf a}}{4\pi \epsilon_0}
1658	<	\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} w_b(r) \hat{r}
1659	<	+\frac{C_{\bf a}}{4\pi \epsilon_0}
1660	<	\sum_n  D_{\mathbf{b}n} \hat{b}_n w_c(r)
1661	<	\end{equation}
1662	<	%
1663	<	% f ca qb
1664	<	%
1665	<	\begin{equation}
1666	<	\begin{split}
1667	<	% 1
1668	<	\mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =
1669	<	\frac{1}{4\pi \epsilon_0}
1670	<	C_{\bf a }\text{Tr}Q_{\bf b} w_d(r) \hat{r}
1671	<	+ 2C_{\bf a } \sum_l  \hat{b}_l Q_{{\mathbf b}ln} (\hat{b}_n \cdot \hat{r}) w_e(r) \\
1672	<	% 2
1673	<	+\frac{C_{\bf a}}{4\pi \epsilon_0}
1674	<	\sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r}) w_f(r) \hat{r}
1675	<	\end{split}
1676	<	\end{equation}
1677	<	%
1678	<	% f da cb
1679	<	%
1680	<	\begin{equation}
1681	<	\mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1682	<	-\frac{C_{\bf{b}}}{4\pi \epsilon_0}
1683	<	\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} w_b(r)  \hat{r}
1684	<	-\frac{C_{\bf{b}}}{4\pi \epsilon_0}
1685	<	\sum_n  D_{\mathbf{a}n} \hat{a}_n w_c(r)
1686	<	\end{equation}
1687	<	%
1688	<	% f da db
1689	<	%
1690	<	\begin{equation}
1691	<	\begin{split}
1692	<	% 1
1693	<	\mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} &=
1694	<	-\frac{1}{4\pi \epsilon_0}
1695	<	\sum_{mn} D_{\mathbf {a}m}
1696	<	(\hat{a}_m \cdot   \hat{b}_n)
1697	<	D_{\mathbf{b}n}  w_d(r) \hat{r}
1698	<	-\frac{1}{4\pi \epsilon_0}
1699	<	\sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1700	<	\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} w_f(r) \hat{r} \\
1701	<	% 2
1702	<	& \quad + \frac{1}{4\pi \epsilon_0}
1703	<	\Bigl[ \sum_m D_{\mathbf {a}m}
1704	<	\hat{a}_m \sum_n D_{\mathbf{b}n}
1705	<	(\hat{b}_n \cdot \hat{r})
1706	<	+ \sum_m D_{\mathbf {b}m}
1707	<	\hat{b}_m \sum_n D_{\mathbf{a}n}
1708	<	(\hat{a}_n \cdot \hat{r}) \Bigr] w_e(r)  \\
1709	<	\end{split}
1710	<	\end{equation}
1711	<	%
1712	<	% f da qb
1713	<	%
1714	<	\begin{equation}
1715	<	\begin{split}
1716	<	% 1
1717	<	&\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1718	<	- \frac{1}{4\pi \epsilon_0} \Bigl[
1719	<	\text{Tr}Q_{\mathbf{b}}
1720	<	\sum_l  D_{\mathbf{a}l} \hat{a}_l
1721	<	+2\sum_{lmn} D_{\mathbf{a}l}
1722	<	(\hat{a}_l \cdot \hat{b}_m)
1723	<	Q_{\mathbf{b}mn} \hat{b}_n  \Bigr] w_g(r) \\
1724	<	% 3
1725	<	&  - \frac{1}{4\pi \epsilon_0} \Bigl[
1726	<	\text{Tr}Q_{\mathbf{b}}
1727	<	\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
1728	<	+2\sum_{lmn}D_{\mathbf{a}l}
1729	<	(\hat{a}_l \cdot \hat{b}_m)
1730	<	Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
1731	<	% 4
1732	<	&+ \frac{1}{4\pi \epsilon_0}
1733	<	\Bigl[\sum_l  D_{\mathbf{a}l} \hat{a}_l
1734	<	\sum_{mn} (\hat{r} \cdot \hat{b}_m)
1735	<	Q_{{\mathbf b}mn}
1736	<	(\hat{b}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{a}_l)
1737	<	D_{\mathbf{a}l}
1738	<	\sum_{mn} (\hat{r} \cdot \hat{b}_m)
1739	<	Q_{{\mathbf b}mn} \hat{b}_n \Bigr]   w_i(r)\\
1740	<	% 6
1741	<	&  -\frac{1}{4\pi \epsilon_0}
1742	<	\sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
1743	<	\sum_{mn} (\hat{r} \cdot \hat{b}_m)
1744	<	Q_{{\mathbf b}mn}
1745	<	(\hat{b}_n \cdot \hat{r})  w_j(r)  \hat{r}
1746	<	\end{split}
1747	<	\end{equation}
1748	<	%
1749	<	% force qa cb
1750	<	%
1751	<	\begin{equation}
1752	<	\begin{split}
1753	<	% 1
1754	<	\mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} &=
1755	<	\frac{1}{4\pi \epsilon_0}
1756	<	C_{\bf b }\text{Tr}Q_{\bf a} \hat{r} w_d(r)
1757	<	+ \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) \\
1758	<	% 2
1759	<	&  +\frac{C_{\bf b}}{4\pi \epsilon_0}
1760	<	\sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) w_f(r) \hat{r}
1761	<	\end{split}
1762	<	\end{equation}
1763	<	%
1764	<	% f qa db
1765	<	%
1766	<	\begin{equation}
1767	<	\begin{split}
1768	<	% 1
1769	<	&\mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1770	<	\frac{1}{4\pi \epsilon_0} \Bigl[
1771	<	\text{Tr}Q_{\mathbf{a}}
1772	<	\sum_l  D_{\mathbf{b}l} \hat{b}_l
1773	<	+2\sum_{lmn} D_{\mathbf{b}l}
1774	<	(\hat{b}_l \cdot \hat{a}_m)
1775	<	Q_{\mathbf{a}mn} \hat{a}_n  \Bigr]
1776	<	w_g(r)\\
1777	<	% 3
1778	<	&  + \frac{1}{4\pi \epsilon_0} \Bigl[
1779	<	\text{Tr}Q_{\mathbf{a}}
1780	<	\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n}
1781	<	+2\sum_{lmn}D_{\mathbf{b}l}
1782	<	(\hat{b}_l \cdot \hat{a}_m)
1783	<	Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
1784	<	% 4
1785	<	&  + \frac{1}{4\pi \epsilon_0} \Bigl[ \sum_l  D_{\mathbf{b}l} \hat{b}_l
1786	<	\sum_{mn} (\hat{r} \cdot \hat{a}_m)
1787	<	Q_{{\mathbf a}mn}
1788	<	(\hat{a}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{b}_l)
1789	<	D_{\mathbf{b}l}
1790	<	\sum_{mn} (\hat{r} \cdot \hat{a}_m)
1791	<	Q_{{\mathbf a}mn} \hat{a}_n \Bigr]   w_i(r) \\
1792	<	% 6
1793	<	&  +\frac{1}{4\pi \epsilon_0}
1794	<	\sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
1795	<	\sum_{mn} (\hat{r} \cdot \hat{a}_m)
1796	<	Q_{{\mathbf a}mn}
1797	<	(\hat{a}_n \cdot \hat{r})  w_j(r)  \hat{r}
1798	<	\end{split}
1799	<	\end{equation}
1800	<	%
1801	<	% f qa qb
1802	<	%
1803	<	\begin{equation}
1804	<	\begin{split}
1805	<	&\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1806	<	\frac{1}{4\pi \epsilon_0} \Bigl[
1807	<	\text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}}
1808	<	+ 2 \sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
1809	<	Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
1810	<	(\hat{a}_m \cdot \hat{b}_n) \Bigr] w_k(r) \hat{r}\\
1811	<	&+\frac{1}{4\pi \epsilon_0} \Bigl[
1812	<	2\text{Tr}Q_{\mathbf{b}} \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  \hat{a}_m
1813	<	+ 2\text{Tr}Q_{\mathbf{a}} \sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm}  \hat{b}_m \\
1814	<	&+ 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})
1815	<	+ 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
1816	<	\Bigr] w_n(r) \\
1817	<	&+ \frac{1}{4\pi \epsilon_0}
1818	<	\Bigl[ \text{Tr}Q_{\mathbf{a}}
1819	<	\sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm} (\hat{b}_m \cdot \hat{r})
1820	<	+ \text{Tr}Q_{\mathbf{b}}
1821	<	\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  (\hat{a}_m \cdot \hat{r}) \\
1822	<	&+4\sum_{lmnp} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{b}_n)
1823	<	Q_{\mathbf{b}np}  (\hat{b}_p \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
1824	<	%
1825	<	&+\frac{1}{4\pi \epsilon_0} \Bigl[
1826	<	2\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} \hat{a}_m
1827	<	\sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_n \cdot \hat{r}) \\
1828	<	&+2 \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
1829	<	\sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} \hat{b}_n \Bigr] w_o(r) \hat{r} \\
1830	<	&  + \frac{1}{4\pi \epsilon_0}
1831	<	\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
1832	<	\sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_p \cdot \hat{r}) w_m(r) \hat{r}
1833	<	\end{split}
1834	<	\end{equation}
1835	<	%
1836	<	Here we list the form of the non-zero damped shifted multipole torques showing
1837	<	explicitly dependences on body axes:
1838	<	%
1839	<	%  t ca db
1840	<	%
1841	<	\begin{equation}
1842	<	\mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =
1843	<	\frac{C_{\bf a}}{4\pi \epsilon_0}
1844	<	\sum_n  (\hat{r} \times \hat{b}_n)  D_{\mathbf{b}n} \,  v_{11}(r)
1845	<	\end{equation}
1846	<	%
1847	<	% t ca qb
1848	<	%
1849	<	\begin{equation}
1850	<	\mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =
1851	<	\frac{2C_{\bf a}}{4\pi \epsilon_0}
1852	<	\sum_{lm} (\hat{r} \times \hat{b}_l) Q_{{\mathbf b}lm} (\hat{b}_m \cdot \hat{r}) v_{22}(r)
1853	<	\end{equation}
1854	<	%
1855	<	%  t da cb
1856	<	%
1857	<	\begin{equation}
1858	<	\mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =
1859	<	-\frac{C_{\bf b}}{4\pi \epsilon_0}
1860	<	\sum_n  (\hat{r} \times \hat{a}_n)  D_{\mathbf{a}n} \,  v_{11}(r)
1861	<	\end{equation}%
1862	<	%
1863	<	%
1864	<	%  ta da db
1865	<	%
1866	<	\begin{equation}
1867	<	\begin{split}
1868	<	% 1
1869	<	\mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} &=
1870	<	\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
1871	<	(\hat{a}_m \times  \hat{b}_n)
1872	<	D_{\mathbf{b}n} v_{21}(r) \\
1873	<	% 2
1874	<	&-\frac{1}{4\pi \epsilon_0}
1875	<	\sum_m (\hat{r} \times \hat{a}_m) D_{\mathbf {a}m}
1876	<	\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
1877	<	\end{split}
1878	<	\end{equation}
1879	<	%
1880	<	%  tb da db
1881	<	%
1882	<	\begin{equation}
1883	<	\begin{split}
1884	<	% 1
1885	<	\mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} &=
1886	<	-\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
1887	<	(\hat{a}_m \times  \hat{b}_n)
1888	<	D_{\mathbf{b}n} v_{21}(r) \\
1889	<	% 2
1890	<	&+\frac{1}{4\pi \epsilon_0}
1891	<	\sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1892	<	\sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
1893	<	\end{split}
1894	<	\end{equation}
1895	<	%
1896	<	% ta da qb
1897	<	%
1898	<	\begin{equation}
1899	<	\begin{split}
1900	<	% 1
1901	<	\mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} &=
1902	<	\frac{1}{4\pi \epsilon_0} \left(
1903	<	-\text{Tr}Q_{\mathbf{b}}
1904	<	\sum_n (\hat{r} \times \hat{a}_n) D_{\mathbf{a}n}
1905	<	+2\sum_{lmn}D_{\mathbf{a}l}
1906	<	(\hat{a}_l \times \hat{b}_m)
1907	<	Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
1908	<	\right) v_{31}(r)\\
1909	<	% 2
1910	<	&-\frac{1}{4\pi \epsilon_0}
1911	<	\sum_l (\hat{r} \times \hat{a}_l) D_{\mathbf{a}l}
1912	<	\sum_{mn} (\hat{r} \cdot \hat{b}_m)
1913	<	Q_{{\mathbf b}mn}
1914	<	(\hat{b}_n \cdot \hat{r}) v_{32}(r)
1915	<	\end{split}
1916	<	\end{equation}
1917	<	%
1918	<	% tb da qb
1919	<	%
1920	<	\begin{equation}
1921	<	\begin{split}
1922	<	% 1
1923	<	\mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} &=
1924	<	\frac{1}{4\pi \epsilon_0} \left(
1925	<	-2\sum_{lmn}D_{\mathbf{a}l}
1926	<	(\hat{a}_l \cdot \hat{b}_m)
1927	<	Q_{\mathbf{b}mn} (\hat{r} \times \hat{b}_n)
1928	<	-2\sum_{lmn}D_{\mathbf{a}l}
1929	<	(\hat{a}_l \times \hat{b}_m)
1930	<	Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
1931	<	\right) v_{31}(r) \\
1932	<	% 2
1933	<	&-\frac{2}{4\pi \epsilon_0}
1934	<	\sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
1935	<	\sum_{mn} (\hat{r} \cdot \hat{b}_m)
1936	<	Q_{{\mathbf b}mn}
1937	<	(\hat{r}\times \hat{b}_n) v_{32}(r)
1938	<	\end{split}
1939	<	\end{equation}
1940	<	%
1941	<	% ta qa cb
1942	<	%
1943	<	\begin{equation}
1944	<	\mathbf{\tau}_{{\bf a}Q_{\bf a}C_{\bf b}} =
1945	<	\frac{2C_{\bf a}}{4\pi \epsilon_0}
1946	<	\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{{\mathbf a}lm} (\hat{r} \times \hat{a}_m) v_{22}(r)
1947	<	\end{equation}
1948	<	%
1949	<	% ta qa db
1950	<	%
1951	<	\begin{equation}
1952	<	\begin{split}
1953	<	% 1
1954	<	\mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} &=
1955	<	\frac{1}{4\pi \epsilon_0} \left(
1956	<	2\sum_{lmn}D_{\mathbf{b}l}
1957	<	(\hat{b}_l \cdot \hat{a}_m)
1958	<	Q_{\mathbf{a}mn} (\hat{r} \times \hat{a}_n)
1959	<	+2\sum_{lmn}D_{\mathbf{b}l}
1960	<	(\hat{a}_l \times \hat{b}_m)
1961	<	Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
1962	<	\right) v_{31}(r) \\
1963	<	% 2
1964	<	&+\frac{2}{4\pi \epsilon_0}
1965	<	\sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
1966	<	\sum_{mn} (\hat{r} \cdot \hat{a}_m)
1967	<	Q_{{\mathbf a}mn}
1968	<	(\hat{r}\times \hat{a}_n) v_{32}(r)
1969	<	\end{split}
1970	<	\end{equation}
1971	<	%
1972	<	% tb qa db
1973	<	%
1974	<	\begin{equation}
1975	<	\begin{split}
1976	<	% 1
1977	<	\mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} &=
1978	<	\frac{1}{4\pi \epsilon_0} \left(
1979	<	\text{Tr}Q_{\mathbf{a}}
1980	<	\sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf{b}n}
1981	<	+2\sum_{lmn}D_{\mathbf{b}l}
1982	<	(\hat{a}_l \times \hat{b}_m)
1983	<	Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
1984	<	\right) v_{31}(r)\\
1985	<	% 2
1986	<	&\frac{1}{4\pi \epsilon_0}
1987	<	\sum_l (\hat{r} \times \hat{b}_l) D_{\mathbf{b}l}
1988	<	\sum_{mn} (\hat{r} \cdot \hat{a}_m)
1989	<	Q_{{\mathbf a}mn}
1990	<	(\hat{a}_n \cdot \hat{r}) v_{32}(r)
1991	<	\end{split}
1992	<	\end{equation}
1993	<	%
1994	<	% ta qa qb
1995	<	%
1996	<	\begin{equation}
1997	<	\begin{split}
1998	<	% 1
1999	<	\mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} &=
2000	<	-\frac{4}{4\pi \epsilon_0}
2001	<	\sum_{lmnp} (\hat{a}_l \times \hat{b}_p)
2002	<	Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
2003	<	(\hat{a}_m \cdot \hat{b}_n) v_{41}(r) \\
2004	<	% 2
2005	<	&+ \frac{1}{4\pi \epsilon_0}
2006	<	\Bigl[
2007	<	2\text{Tr}Q_{\mathbf{b}}
2008	<	\sum_{lm} (\hat{r} \cdot \hat{a}_l )
2009	<	Q_{{\mathbf a}lm}
2010	<	(\hat{r} \times \hat{a}_m)
2011	<	+4 \sum_{lmnp}
2012	<	(\hat{r} \times \hat{a}_l )
2013	<	Q_{{\mathbf a}lm}
2014	<	(\hat{a}_m \cdot \hat{b}_n)
2015	<	Q_{{\mathbf b}np}
2016	<	(\hat{b}_p \cdot \hat{r}) \\
2017	<	% 3
2018	<	&-4 \sum_{lmnp}
2019	<	(\hat{r} \cdot \hat{a}_l )
2020	<	Q_{{\mathbf a}lm}
2021	<	(\hat{a}_m \times \hat{b}_n)
2022	<	Q_{{\mathbf b}np}
2023	<	(\hat{b}_p \cdot \hat{r})
2024	<	\Bigr] v_{42}(r) \\
2025	<	% 4
2026	<	&+ \frac{2}{4\pi \epsilon_0}
2027	<	\sum_{lm} (\hat{r} \times \hat{a}_l)
2028	<	Q_{{\mathbf a}lm}
2029	<	(\hat{a}_m \cdot \hat{r})
2030	<	\sum_{np}  (\hat{r} \cdot \hat{b}_n)
2031	<	Q_{{\mathbf b}np}
2032	<	(\hat{b}_p \cdot \hat{r})  v_{43}(r)\\
2033	<	\end{split}
2034	<	\end{equation}
2035	<	%
2036	<	% tb qa qb
2037	<	%
2038	<	\begin{equation}
2039	<	\begin{split}
2040	<	% 1
2041	<	\mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} &=
2042	<	\frac{4}{4\pi \epsilon_0}
2043	<	\sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
2044	<	Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
2045	<	(\hat{a}_m \times \hat{b}_n) v_{41}(r) \\
2046	<	% 2
2047	<	&+ \frac{1}{4\pi \epsilon_0}
2048	<	\Bigl[
2049	<	2\text{Tr}Q_{\mathbf{a}}
2050	<	\sum_{lm} (\hat{r} \cdot \hat{b}_l )
2051	<	Q_{{\mathbf b}lm}
2052	<	(\hat{r} \times \hat{b}_m)
2053	<	+4 \sum_{lmnp}
2054	<	(\hat{r} \cdot \hat{a}_l )
2055	<	Q_{{\mathbf a}lm}
2056	<	(\hat{a}_m \cdot \hat{b}_n)
2057	<	Q_{{\mathbf b}np}
2058	<	(\hat{r} \times \hat{b}_p) \\
2059	<	% 3
2060	<	&+4 \sum_{lmnp}
2061	<	(\hat{r} \cdot \hat{a}_l )
2062	<	Q_{{\mathbf a}lm}
2063	<	(\hat{a}_m \times \hat{b}_n)
2064	<	Q_{{\mathbf b}np}
2065	<	(\hat{b}_p \cdot \hat{r})
2066	<	\Bigr] v_{42}(r)  \\
2067	<	% 4
2068	<	&+ \frac{2}{4\pi \epsilon_0}
2069	<	\sum_{lm} (\hat{r} \cdot \hat{a}_l)
2070	<	Q_{{\mathbf a}lm}
2071	<	(\hat{a}_m \cdot \hat{r})
2072	<	\sum_{np}  (\hat{r} \times \hat{b}_n)
2073	<	Q_{{\mathbf b}np}
2074	<	(\hat{b}_p \cdot \hat{r}) v_{43}(r). \\
2075	<	\end{split}
2076	<	\end{equation}
2077	<	%
2078	<	\begin{table*}
2079	<	\caption{\label{tab:tableFORCE2}Radial functions used in the force equations.}
2080	<	\begin{ruledtabular}
2081	<	\begin{tabular}{|l|l|l|}
2082	<	Generic&Taylor-shifted Force&Gradient-shifted Force
2083	<	\\ \hline
2084	<	%
2085	<	%
2086	<	%
2087	<	$w_a(r)$&
2088	<	$g_0(r)$&
2089	<	$g(r)-g(r_c)$ \\
2090	<	%
2091	<	%
2092	<	$w_b(r)$ &
2093	<	$\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ &
2094	<	$h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\
2095	<	%
2096	<	$w_c(r)$ &
2097	<	$\frac{g_1(r)}{r} $ &
2098	<	$\frac{v_{11}(r)}{r}$ \\
2099	<	%
2100	<	%
2101	<	$w_d(r)$&
2102	<	$\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ &
2103	<	$\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)
2104	<	-\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $\\
2105	<	%
2106	<	$w_e(r)$ &
2107	<	$\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ &
2108	<	$\frac{v_{22}(r)}{r}$ \\
2109	<	%
2110	<	%
2111	<	$w_f(r)$&
2112	<	$\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ &
2113	<	$\left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) - $ \\
2114	<	&&$\left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)-\frac{2v_{22}(r)}{r}$\\
2115	<	%
2116	<	$w_g(r)$& $ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$&
2117	<	$\frac{v_{31}(r)}{r}$\\
2118	<	%
2119	<	$w_h(r)$ &
2120	<	$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
2121	<	$\left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) - $\\
2122	<	&&$\left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\
2123	<	&&$-\frac{v_{31}(r)}{r}$\\
2124	<	% 2
2125	<	$w_i(r)$ &
2126	<	$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
2127	<	$\frac{v_{32}(r)}{r}$ \\
2128	<	%
2129	<	$w_j(r)$ &
2130	<	$\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right)  $ &
2131	<	$\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right) $  \\
2132	<	&&$\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}$ \\
2133	<	%
2134	<	$w_k(r)$ &
2135	<	$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
2136	<	$\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}  \right)$  \\
2137	<	&&$\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2}  \right)$ \\
2138	<	%
2139	<	$w_l(r)$ &
2140	<	$\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)$ &
2141	<	$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\
2142	<	&&$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)
2143	<	-\frac{2v_{42}(r)}{r}$ \\
2144	<	%
2145	<	$w_m(r)$ &
2146	<	$\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)$ &
2147	<	$\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$ \\
2148	<	&&$\left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}
2149	<	+\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $ \\
2150	<	&&$-\frac{4v_{43}(r)}{r}$ \\
2151	<	%
2152	<	$w_n(r)$ &
2153	<	$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
2154	<	$\frac{v_{42}(r)}{r}$ \\
2155	<	%
2156	<	$w_o(r)$ &
2157	<	$\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)$ &
2158	<	$\frac{v_{43}(r)}{r}$ \\
1640	>	% %
1641	>	% % u ca cb
1642	>	% %
1643	>	% \begin{equation}
1644	>	% U_{C_{\bf a}C_{\bf b}}(r)=
1645	>	% \frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r)
1646	>	% \end{equation}
1647	>	% %
1648	>	% % u ca db
1649	>	% %
1650	>	% \begin{equation}
1651	>	% U_{C_{\bf a}D_{\bf b}}(r)=
1652	>	% \frac{C_{\bf a}}{4\pi \epsilon_0}
1653	>	% \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} \,  v_{11}(r)
1654	>	% \end{equation}
1655	>	% %
1656	>	% % u ca qb
1657	>	% %
1658	>	% \begin{equation}
1659	>	% U_{C_{\bf a}Q_{\bf b}}(r)=
1660	>	% \frac{C_{\bf a }\text{Tr}Q_{\bf b}}{4\pi \epsilon_0}
1661	>	% v_{21}(r) \nonumber \\
1662	>	% +\frac{C_{\bf a}}{4\pi \epsilon_0}
1663	>	% \sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r})
1664	>	% v_{22}(r)
1665	>	% \end{equation}
1666	>	% %
1667	>	% % u da cb
1668	>	% %
1669	>	% \begin{equation}
1670	>	% U_{D_{\bf a}C_{\bf b}}(r)=
1671	>	% -\frac{C_{\bf b}}{4\pi \epsilon_0}
1672	>	% \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} \,  v_{11}(r)
1673	>	% \end{equation}
1674	>	% %
1675	>	% % u da db
1676	>	% %
1677	>	% \begin{equation}
1678	>	% \begin{split}
1679	>	% % 1
1680	>	% U_{D_{\bf a}D_{\bf b}}(r)&=
1681	>	% -\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
1682	>	% (\hat{a}_m \cdot   \hat{b}_n)
1683	>	% D_{\mathbf{b}n} v_{21}(r) \\
1684	>	% % 2
1685	>	% &-\frac{1}{4\pi \epsilon_0}
1686	>	% \sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1687	>	% \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n}
1688	>	% v_{22}(r)
1689	>	% \end{split}
1690	>	% \end{equation}
1691	>	% %
1692	>	% % u da qb
1693	>	% %
1694	>	% \begin{equation}
1695	>	% \begin{split}
1696	>	% % 1
1697	>	% U_{D_{\bf a}Q_{\bf b}}(r)&=
1698	>	% -\frac{1}{4\pi \epsilon_0} \left(
1699	>	% \text{Tr}Q_{\mathbf{b}}
1700	>	% \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
1701	>	% +2\sum_{lmn}D_{\mathbf{a}l}
1702	>	% (\hat{a}_l \cdot \hat{b}_m)
1703	>	% Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
1704	>	% \right)  v_{31}(r) \\
1705	>	% % 2
1706	>	% &-\frac{1}{4\pi \epsilon_0}
1707	>	% \sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
1708	>	% \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1709	>	% Q_{{\mathbf b}mn}
1710	>	% (\hat{b}_n \cdot \hat{r}) v_{32}(r)
1711	>	% \end{split}
1712	>	% \end{equation}
1713	>	% %
1714	>	% % u qa cb
1715	>	% %
1716	>	% \begin{equation}
1717	>	% U_{Q_{\bf a}C_{\bf b}}(r)=
1718	>	% \frac{C_{\bf b }\text{Tr}Q_{\bf a}}{4\pi \epsilon_0}  v_{21}(r)
1719	>	% +\frac{C_{\bf b}}{4\pi \epsilon_0}
1720	>	% \sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) v_{22}(r)
1721	>	% \end{equation}
1722	>	% %
1723	>	% % u qa db
1724	>	% %
1725	>	% \begin{equation}
1726	>	% \begin{split}
1727	>	% %1
1728	>	% U_{Q_{\bf a}D_{\bf b}}(r)&=
1729	>	% \frac{1}{4\pi \epsilon_0} \left(
1730	>	% \text{Tr}Q_{\mathbf{a}}
1731	>	% \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n}
1732	>	% +2\sum_{lmn}D_{\mathbf{b}l}
1733	>	% (\hat{b}_l \cdot \hat{a}_m)
1734	>	% Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
1735	>	% \right) v_{31}(r)  \\
1736	>	% % 2
1737	>	% &+\frac{1}{4\pi \epsilon_0}
1738	>	% \sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
1739	>	% \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1740	>	% Q_{{\mathbf a}mn}
1741	>	% (\hat{a}_n \cdot \hat{r}) v_{32}(r)
1742	>	% \end{split}
1743	>	% \end{equation}
1744	>	% %
1745	>	% % u qa qb
1746	>	% %
1747	>	% \begin{equation}
1748	>	% \begin{split}
1749	>	% %1
1750	>	% U_{Q_{\bf a}Q_{\bf b}}(r)&=
1751	>	% \frac{1}{4\pi \epsilon_0} \Bigl[
1752	>	% \text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}}
1753	>	% +2\sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
1754	>	% Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
1755	>	% (\hat{a}_m \cdot \hat{b}_n) \Bigr]
1756	>	% v_{41}(r) \\
1757	>	% % 2
1758	>	% &+ \frac{1}{4\pi \epsilon_0}
1759	>	% \Bigl[ \text{Tr}Q_{\mathbf{a}}
1760	>	% \sum_{lm} (\hat{r} \cdot \hat{b}_l )
1761	>	% Q_{{\mathbf b}lm}
1762	>	% (\hat{b}_m \cdot \hat{r})
1763	>	% +\text{Tr}Q_{\mathbf{b}}
1764	>	% \sum_{lm} (\hat{r} \cdot \hat{a}_l )
1765	>	% Q_{{\mathbf a}lm}
1766	>	% (\hat{a}_m \cdot \hat{r}) \\
1767	>	% % 3
1768	>	% &+4 \sum_{lmnp}
1769	>	% (\hat{r} \cdot \hat{a}_l )
1770	>	% Q_{{\mathbf a}lm}
1771	>	% (\hat{a}_m \cdot \hat{b}_n)
1772	>	% Q_{{\mathbf b}np}
1773	>	% (\hat{b}_p \cdot \hat{r})
1774	>	% \Bigr] v_{42}(r)  \\
1775	>	% % 4
1776	>	% &+ \frac{1}{4\pi \epsilon_0}
1777	>	% \sum_{lm} (\hat{r} \cdot \hat{a}_l)
1778	>	% Q_{{\mathbf a}lm}
1779	>	% (\hat{a}_m \cdot \hat{r})
1780	>	% \sum_{np}  (\hat{r} \cdot \hat{b}_n)
1781	>	% Q_{{\mathbf b}np}
1782	>	% (\hat{b}_p \cdot \hat{r})  v_{43}(r).
1783	>	% \end{split}
1784	>	% \end{equation}
1785	>	% %
1786	>
1787	>
1788	>	% % BODY coordinates force equations --------------------------------------------
1789	>	% %
1790	>	% %
1791	>	% Here are the force equations written in terms of body coordinates.
1792	>	% %
1793	>	% % f ca cb
1794	>	% %
1795	>	% \begin{equation}
1796	>	% \mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =
1797	>	% \frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r}
1798	>	% \end{equation}
1799	>	% %
1800	>	% % f ca db
1801	>	% %
1802	>	% \begin{equation}
1803	>	% \mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =
1804	>	% \frac{C_{\bf a}}{4\pi \epsilon_0}
1805	>	% \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} w_b(r) \hat{r}
1806	>	% +\frac{C_{\bf a}}{4\pi \epsilon_0}
1807	>	% \sum_n  D_{\mathbf{b}n} \hat{b}_n w_c(r)
1808	>	% \end{equation}
1809	>	% %
1810	>	% % f ca qb
1811	>	% %
1812	>	% \begin{equation}
1813	>	% \begin{split}
1814	>	% % 1
1815	>	% \mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =
1816	>	% \frac{1}{4\pi \epsilon_0}
1817	>	% C_{\bf a }\text{Tr}Q_{\bf b} w_d(r) \hat{r}
1818	>	% + 2C_{\bf a } \sum_l  \hat{b}_l Q_{{\mathbf b}ln} (\hat{b}_n \cdot \hat{r}) w_e(r) \\
1819	>	% % 2
1820	>	% +\frac{C_{\bf a}}{4\pi \epsilon_0}
1821	>	% \sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r}) w_f(r) \hat{r}
1822	>	% \end{split}
1823	>	% \end{equation}
1824	>	% %
1825	>	% % f da cb
1826	>	% %
1827	>	% \begin{equation}
1828	>	% \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1829	>	% -\frac{C_{\bf{b}}}{4\pi \epsilon_0}
1830	>	% \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} w_b(r)  \hat{r}
1831	>	% -\frac{C_{\bf{b}}}{4\pi \epsilon_0}
1832	>	% \sum_n  D_{\mathbf{a}n} \hat{a}_n w_c(r)
1833	>	% \end{equation}
1834	>	% %
1835	>	% % f da db
1836	>	% %
1837	>	% \begin{equation}
1838	>	% \begin{split}
1839	>	% % 1
1840	>	% \mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} &=
1841	>	% -\frac{1}{4\pi \epsilon_0}
1842	>	%  \sum_{mn} D_{\mathbf {a}m}
1843	>	% (\hat{a}_m \cdot   \hat{b}_n)
1844	>	% D_{\mathbf{b}n}  w_d(r) \hat{r}
1845	>	% -\frac{1}{4\pi \epsilon_0}
1846	>	% \sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1847	>	% \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} w_f(r) \hat{r} \\
1848	>	% % 2
1849	>	% & \quad + \frac{1}{4\pi \epsilon_0}
1850	>	% \Bigl[ \sum_m D_{\mathbf {a}m}
1851	>	% \hat{a}_m \sum_n D_{\mathbf{b}n}
1852	>	% (\hat{b}_n \cdot \hat{r})
1853	>	% + \sum_m D_{\mathbf {b}m}
1854	>	% \hat{b}_m \sum_n D_{\mathbf{a}n}
1855	>	% (\hat{a}_n \cdot \hat{r}) \Bigr] w_e(r)  \\
1856	>	% \end{split}
1857	>	% \end{equation}
1858	>	% %
1859	>	% % f da qb
1860	>	% %
1861	>	% \begin{equation}
1862	>	% \begin{split}
1863	>	% % 1
1864	>	% &\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1865	>	%  - \frac{1}{4\pi \epsilon_0} \Bigl[
1866	>	% \text{Tr}Q_{\mathbf{b}}
1867	>	% \sum_l  D_{\mathbf{a}l} \hat{a}_l
1868	>	% +2\sum_{lmn} D_{\mathbf{a}l}
1869	>	% (\hat{a}_l \cdot \hat{b}_m)
1870	>	% Q_{\mathbf{b}mn} \hat{b}_n  \Bigr] w_g(r) \\
1871	>	% % 3
1872	>	% &  - \frac{1}{4\pi \epsilon_0} \Bigl[
1873	>	% \text{Tr}Q_{\mathbf{b}}
1874	>	% \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
1875	>	% +2\sum_{lmn}D_{\mathbf{a}l}
1876	>	% (\hat{a}_l \cdot \hat{b}_m)
1877	>	% Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
1878	>	% % 4
1879	>	% &+ \frac{1}{4\pi \epsilon_0}
1880	>	% \Bigl[\sum_l  D_{\mathbf{a}l} \hat{a}_l
1881	>	% \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1882	>	% Q_{{\mathbf b}mn}
1883	>	% (\hat{b}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{a}_l)
1884	>	%  D_{\mathbf{a}l}
1885	>	% \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1886	>	% Q_{{\mathbf b}mn} \hat{b}_n \Bigr]   w_i(r)\\
1887	>	% % 6
1888	>	% &  -\frac{1}{4\pi \epsilon_0}
1889	>	% \sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
1890	>	% \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1891	>	% Q_{{\mathbf b}mn}
1892	>	% (\hat{b}_n \cdot \hat{r})  w_j(r)  \hat{r}
1893	>	% \end{split}
1894	>	% \end{equation}
1895	>	% %
1896	>	% % force qa cb
1897	>	% %
1898	>	% \begin{equation}
1899	>	% \begin{split}
1900	>	% % 1
1901	>	% \mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} &=
1902	>	% \frac{1}{4\pi \epsilon_0}
1903	>	% C_{\bf b }\text{Tr}Q_{\bf a} \hat{r} w_d(r)
1904	>	% + \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) \\
1905	>	% % 2
1906	>	% &  +\frac{C_{\bf b}}{4\pi \epsilon_0}
1907	>	% \sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) w_f(r) \hat{r}
1908	>	% \end{split}
1909	>	% \end{equation}
1910	>	% %
1911	>	% % f qa db
1912	>	% %
1913	>	% \begin{equation}
1914	>	% \begin{split}
1915	>	% % 1
1916	>	% &\mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1917	>	% \frac{1}{4\pi \epsilon_0} \Bigl[
1918	>	% \text{Tr}Q_{\mathbf{a}}
1919	>	% \sum_l  D_{\mathbf{b}l} \hat{b}_l
1920	>	% +2\sum_{lmn} D_{\mathbf{b}l}
1921	>	% (\hat{b}_l \cdot \hat{a}_m)
1922	>	% Q_{\mathbf{a}mn} \hat{a}_n  \Bigr]
1923	>	% w_g(r)\\
1924	>	% % 3
1925	>	% &  + \frac{1}{4\pi \epsilon_0} \Bigl[
1926	>	% \text{Tr}Q_{\mathbf{a}}
1927	>	% \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n}
1928	>	% +2\sum_{lmn}D_{\mathbf{b}l}
1929	>	% (\hat{b}_l \cdot \hat{a}_m)
1930	>	% Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
1931	>	% % 4
1932	>	% &  + \frac{1}{4\pi \epsilon_0} \Bigl[ \sum_l  D_{\mathbf{b}l} \hat{b}_l
1933	>	% \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1934	>	% Q_{{\mathbf a}mn}
1935	>	% (\hat{a}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{b}_l)
1936	>	%  D_{\mathbf{b}l}
1937	>	% \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1938	>	% Q_{{\mathbf a}mn} \hat{a}_n \Bigr]   w_i(r) \\
1939	>	% % 6
1940	>	% &  +\frac{1}{4\pi \epsilon_0}
1941	>	% \sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
1942	>	% \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1943	>	% Q_{{\mathbf a}mn}
1944	>	% (\hat{a}_n \cdot \hat{r})  w_j(r)  \hat{r}
1945	>	% \end{split}
1946	>	% \end{equation}
1947	>	% %
1948	>	% % f qa qb
1949	>	% %
1950	>	% \begin{equation}
1951	>	% \begin{split}
1952	>	% &\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1953	>	% \frac{1}{4\pi \epsilon_0} \Bigl[
1954	>	% \text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}}
1955	>	% + 2 \sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
1956	>	% Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
1957	>	% (\hat{a}_m \cdot \hat{b}_n) \Bigr] w_k(r) \hat{r}\\
1958	>	% &+\frac{1}{4\pi \epsilon_0} \Bigl[
1959	>	% 2\text{Tr}Q_{\mathbf{b}} \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  \hat{a}_m
1960	>	% + 2\text{Tr}Q_{\mathbf{a}} \sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm}  \hat{b}_m \\
1961	>	% &+ 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})
1962	>	% + 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
1963	>	% \Bigr] w_n(r) \\
1964	>	% &+ \frac{1}{4\pi \epsilon_0}
1965	>	% \Bigl[ \text{Tr}Q_{\mathbf{a}}
1966	>	% \sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm} (\hat{b}_m \cdot \hat{r})
1967	>	% + \text{Tr}Q_{\mathbf{b}}
1968	>	% \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  (\hat{a}_m \cdot \hat{r}) \\
1969	>	% &+4\sum_{lmnp} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{b}_n)
1970	>	% Q_{\mathbf{b}np}  (\hat{b}_p \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
1971	>	% %
1972	>	% &+\frac{1}{4\pi \epsilon_0} \Bigl[
1973	>	% 2\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} \hat{a}_m
1974	>	% \sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_n \cdot \hat{r}) \\
1975	>	% &+2 \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
1976	>	% \sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} \hat{b}_n \Bigr] w_o(r) \hat{r} \\
1977	>	% &  + \frac{1}{4\pi \epsilon_0}
1978	>	% \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
1979	>	% \sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_p \cdot \hat{r}) w_m(r) \hat{r}
1980	>	% \end{split}
1981	>	% \end{equation}
1982	>	% %
1983	>	% Here we list the form of the non-zero damped shifted multipole torques showing
1984	>	% explicitly dependences on body axes:
1985	>	% %
1986	>	% %  t ca db
1987	>	% %
1988	>	% \begin{equation}
1989	>	% \mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =
1990	>	% \frac{C_{\bf a}}{4\pi \epsilon_0}
1991	>	% \sum_n  (\hat{r} \times \hat{b}_n)  D_{\mathbf{b}n} \,  v_{11}(r)
1992	>	% \end{equation}
1993	>	% %
1994	>	% % t ca qb
1995	>	% %
1996	>	% \begin{equation}
1997	>	% \mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =
1998	>	% \frac{2C_{\bf a}}{4\pi \epsilon_0}
1999	>	% \sum_{lm} (\hat{r} \times \hat{b}_l) Q_{{\mathbf b}lm} (\hat{b}_m \cdot \hat{r}) v_{22}(r)
2000	>	% \end{equation}
2001	>	% %
2002	>	% %  t da cb
2003	>	% %
2004	>	% \begin{equation}
2005	>	% \mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =
2006	>	% -\frac{C_{\bf b}}{4\pi \epsilon_0}
2007	>	% \sum_n  (\hat{r} \times \hat{a}_n)  D_{\mathbf{a}n} \,  v_{11}(r)
2008	>	% \end{equation}%
2009	>	% %
2010	>	% %
2011	>	% %  ta da db
2012	>	% %
2013	>	% \begin{equation}
2014	>	% \begin{split}
2015	>	% % 1
2016	>	% \mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} &=
2017	>	% \frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
2018	>	% (\hat{a}_m \times  \hat{b}_n)
2019	>	% D_{\mathbf{b}n} v_{21}(r) \\
2020	>	% % 2
2021	>	% &-\frac{1}{4\pi \epsilon_0}
2022	>	% \sum_m (\hat{r} \times \hat{a}_m) D_{\mathbf {a}m}
2023	>	% \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
2024	>	% \end{split}
2025	>	% \end{equation}
2026	>	% %
2027	>	% %  tb da db
2028	>	% %
2029	>	% \begin{equation}
2030	>	% \begin{split}
2031	>	% % 1
2032	>	% \mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} &=
2033	>	% -\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
2034	>	% (\hat{a}_m \times  \hat{b}_n)
2035	>	% D_{\mathbf{b}n} v_{21}(r) \\
2036	>	% % 2
2037	>	% &+\frac{1}{4\pi \epsilon_0}
2038	>	% \sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
2039	>	% \sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
2040	>	% \end{split}
2041	>	% \end{equation}
2042	>	% %
2043	>	% % ta da qb
2044	>	% %
2045	>	% \begin{equation}
2046	>	% \begin{split}
2047	>	% % 1
2048	>	% \mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} &=
2049	>	% \frac{1}{4\pi \epsilon_0} \left(
2050	>	% -\text{Tr}Q_{\mathbf{b}}
2051	>	% \sum_n (\hat{r} \times \hat{a}_n) D_{\mathbf{a}n}
2052	>	% +2\sum_{lmn}D_{\mathbf{a}l}
2053	>	% (\hat{a}_l \times \hat{b}_m)
2054	>	% Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
2055	>	% \right) v_{31}(r)\\
2056	>	% % 2
2057	>	% &-\frac{1}{4\pi \epsilon_0}
2058	>	% \sum_l (\hat{r} \times \hat{a}_l) D_{\mathbf{a}l}
2059	>	% \sum_{mn} (\hat{r} \cdot \hat{b}_m)
2060	>	% Q_{{\mathbf b}mn}
2061	>	% (\hat{b}_n \cdot \hat{r}) v_{32}(r)
2062	>	% \end{split}
2063	>	% \end{equation}
2064	>	% %
2065	>	% % tb da qb
2066	>	% %
2067	>	% \begin{equation}
2068	>	% \begin{split}
2069	>	% % 1
2070	>	% \mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} &=
2071	>	% \frac{1}{4\pi \epsilon_0} \left(
2072	>	% -2\sum_{lmn}D_{\mathbf{a}l}
2073	>	% (\hat{a}_l \cdot \hat{b}_m)
2074	>	% Q_{\mathbf{b}mn} (\hat{r} \times \hat{b}_n)
2075	>	% -2\sum_{lmn}D_{\mathbf{a}l}
2076	>	% (\hat{a}_l \times \hat{b}_m)
2077	>	% Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
2078	>	% \right) v_{31}(r) \\
2079	>	% % 2
2080	>	% &-\frac{2}{4\pi \epsilon_0}
2081	>	% \sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
2082	>	% \sum_{mn} (\hat{r} \cdot \hat{b}_m)
2083	>	% Q_{{\mathbf b}mn}
2084	>	% (\hat{r}\times \hat{b}_n) v_{32}(r)
2085	>	% \end{split}
2086	>	% \end{equation}
2087	>	% %
2088	>	% % ta qa cb
2089	>	% %
2090	>	% \begin{equation}
2091	>	% \mathbf{\tau}_{{\bf a}Q_{\bf a}C_{\bf b}} =
2092	>	% \frac{2C_{\bf a}}{4\pi \epsilon_0}
2093	>	% \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{{\mathbf a}lm} (\hat{r} \times \hat{a}_m) v_{22}(r)
2094	>	% \end{equation}
2095	>	% %
2096	>	% % ta qa db
2097	>	% %
2098	>	% \begin{equation}
2099	>	% \begin{split}
2100	>	% % 1
2101	>	% \mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} &=
2102	>	% \frac{1}{4\pi \epsilon_0} \left(
2103	>	% 2\sum_{lmn}D_{\mathbf{b}l}
2104	>	% (\hat{b}_l \cdot \hat{a}_m)
2105	>	% Q_{\mathbf{a}mn} (\hat{r} \times \hat{a}_n)
2106	>	% +2\sum_{lmn}D_{\mathbf{b}l}
2107	>	% (\hat{a}_l \times \hat{b}_m)
2108	>	% Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
2109	>	% \right) v_{31}(r) \\
2110	>	% % 2
2111	>	% &+\frac{2}{4\pi \epsilon_0}
2112	>	% \sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
2113	>	% \sum_{mn} (\hat{r} \cdot \hat{a}_m)
2114	>	% Q_{{\mathbf a}mn}
2115	>	% (\hat{r}\times \hat{a}_n) v_{32}(r)
2116	>	% \end{split}
2117	>	% \end{equation}
2118	>	% %
2119	>	% % tb qa db
2120	>	% %
2121	>	% \begin{equation}
2122	>	% \begin{split}
2123	>	% % 1
2124	>	% \mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} &=
2125	>	% \frac{1}{4\pi \epsilon_0} \left(
2126	>	% \text{Tr}Q_{\mathbf{a}}
2127	>	% \sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf{b}n}
2128	>	% +2\sum_{lmn}D_{\mathbf{b}l}
2129	>	% (\hat{a}_l \times \hat{b}_m)
2130	>	% Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
2131	>	% \right) v_{31}(r)\\
2132	>	% % 2
2133	>	% &\frac{1}{4\pi \epsilon_0}
2134	>	% \sum_l (\hat{r} \times \hat{b}_l) D_{\mathbf{b}l}
2135	>	% \sum_{mn} (\hat{r} \cdot \hat{a}_m)
2136	>	% Q_{{\mathbf a}mn}
2137	>	% (\hat{a}_n \cdot \hat{r}) v_{32}(r)
2138	>	% \end{split}
2139	>	% \end{equation}
2140	>	% %
2141	>	% % ta qa qb
2142	>	% %
2143	>	% \begin{equation}
2144	>	% \begin{split}
2145	>	% % 1
2146	>	% \mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} &=
2147	>	% -\frac{4}{4\pi \epsilon_0}
2148	>	% \sum_{lmnp} (\hat{a}_l \times \hat{b}_p)
2149	>	% Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
2150	>	% (\hat{a}_m \cdot \hat{b}_n) v_{41}(r) \\
2151	>	% % 2
2152	>	% &+ \frac{1}{4\pi \epsilon_0}
2153	>	% \Bigl[
2154	>	% 2\text{Tr}Q_{\mathbf{b}}
2155	>	% \sum_{lm} (\hat{r} \cdot \hat{a}_l )
2156	>	% Q_{{\mathbf a}lm}
2157	>	% (\hat{r} \times \hat{a}_m)
2158	>	% +4 \sum_{lmnp}
2159	>	% (\hat{r} \times \hat{a}_l )
2160	>	% Q_{{\mathbf a}lm}
2161	>	% (\hat{a}_m \cdot \hat{b}_n)
2162	>	% Q_{{\mathbf b}np}
2163	>	% (\hat{b}_p \cdot \hat{r}) \\
2164	>	% % 3
2165	>	% &-4 \sum_{lmnp}
2166	>	% (\hat{r} \cdot \hat{a}_l )
2167	>	% Q_{{\mathbf a}lm}
2168	>	% (\hat{a}_m \times \hat{b}_n)
2169	>	% Q_{{\mathbf b}np}
2170	>	% (\hat{b}_p \cdot \hat{r})
2171	>	% \Bigr] v_{42}(r) \\
2172	>	% % 4
2173	>	% &+ \frac{2}{4\pi \epsilon_0}
2174	>	% \sum_{lm} (\hat{r} \times \hat{a}_l)
2175	>	% Q_{{\mathbf a}lm}
2176	>	% (\hat{a}_m \cdot \hat{r})
2177	>	% \sum_{np}  (\hat{r} \cdot \hat{b}_n)
2178	>	% Q_{{\mathbf b}np}
2179	>	% (\hat{b}_p \cdot \hat{r})  v_{43}(r)\\
2180	>	% \end{split}
2181	>	% \end{equation}
2182	>	% %
2183	>	% % tb qa qb
2184	>	% %
2185	>	% \begin{equation}
2186	>	% \begin{split}
2187	>	% % 1
2188	>	% \mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} &=
2189	>	% \frac{4}{4\pi \epsilon_0}
2190	>	% \sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
2191	>	% Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
2192	>	% (\hat{a}_m \times \hat{b}_n) v_{41}(r) \\
2193	>	% % 2
2194	>	% &+ \frac{1}{4\pi \epsilon_0}
2195	>	% \Bigl[
2196	>	% 2\text{Tr}Q_{\mathbf{a}}
2197	>	% \sum_{lm} (\hat{r} \cdot \hat{b}_l )
2198	>	% Q_{{\mathbf b}lm}
2199	>	% (\hat{r} \times \hat{b}_m)
2200	>	% +4 \sum_{lmnp}
2201	>	% (\hat{r} \cdot \hat{a}_l )
2202	>	% Q_{{\mathbf a}lm}
2203	>	% (\hat{a}_m \cdot \hat{b}_n)
2204	>	% Q_{{\mathbf b}np}
2205	>	% (\hat{r} \times \hat{b}_p) \\
2206	>	% % 3
2207	>	% &+4 \sum_{lmnp}
2208	>	% (\hat{r} \cdot \hat{a}_l )
2209	>	% Q_{{\mathbf a}lm}
2210	>	% (\hat{a}_m \times \hat{b}_n)
2211	>	% Q_{{\mathbf b}np}
2212	>	% (\hat{b}_p \cdot \hat{r})
2213	>	% \Bigr] v_{42}(r)  \\
2214	>	% % 4
2215	>	% &+ \frac{2}{4\pi \epsilon_0}
2216	>	% \sum_{lm} (\hat{r} \cdot \hat{a}_l)
2217	>	% Q_{{\mathbf a}lm}
2218	>	% (\hat{a}_m \cdot \hat{r})
2219	>	% \sum_{np}  (\hat{r} \times \hat{b}_n)
2220	>	% Q_{{\mathbf b}np}
2221	>	% (\hat{b}_p \cdot \hat{r}) v_{43}(r). \\
2222	>	% \end{split}
2223	>	% \end{equation}
2224		%
2225	<	\end{tabular}
2226	<	\end{ruledtabular}
2227	<	\end{table*}
2225	>	% \begin{table*}
2226	>	% \caption{\label{tab:tableFORCE2}Radial functions used in the force equations.}
2227	>	% \begin{ruledtabular}
2228	>	% \begin{tabular}{|l|l|l|}
2229	>	% Generic&Taylor-shifted Force&Gradient-shifted Force
2230	>	% \\ \hline
2231	>	% %
2232	>	% %
2233	>	% %
2234	>	% $w_a(r)$&
2235	>	% $g_0(r)$&
2236	>	% $g(r)-g(r_c)$ \\
2237	>	% %
2238	>	% %
2239	>	% $w_b(r)$ &
2240	>	% $\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ &
2241	>	% $h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\
2242	>	% %
2243	>	% $w_c(r)$ &
2244	>	% $\frac{g_1(r)}{r} $ &
2245	>	% $\frac{v_{11}(r)}{r}$ \\
2246	>	% %
2247	>	% %
2248	>	% $w_d(r)$&
2249	>	% $\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ &
2250	>	% $\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)
2251	>	% -\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $\\
2252	>	% %
2253	>	% $w_e(r)$ &
2254	>	% $\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ &
2255	>	% $\frac{v_{22}(r)}{r}$ \\
2256	>	% %
2257	>	% %
2258	>	% $w_f(r)$&
2259	>	% $\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ &
2260	>	% $\left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) - $ \\
2261	>	% &&$\left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)-\frac{2v_{22}(r)}{r}$\\
2262	>	% %
2263	>	% $w_g(r)$& $ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$&
2264	>	% $\frac{v_{31}(r)}{r}$\\
2265	>	% %
2266	>	% $w_h(r)$ &
2267	>	% $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
2268	>	% $\left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) - $\\
2269	>	% &&$\left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\
2270	>	% &&$-\frac{v_{31}(r)}{r}$\\
2271	>	% % 2
2272	>	% $w_i(r)$ &
2273	>	% $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
2274	>	% $\frac{v_{32}(r)}{r}$ \\
2275	>	% %
2276	>	% $w_j(r)$ &
2277	>	% $\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right)  $ &
2278	>	% $\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right) $  \\
2279	>	% &&$\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}$ \\
2280	>	% %
2281	>	% $w_k(r)$ &
2282	>	% $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
2283	>	% $\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}  \right)$  \\
2284	>	% &&$\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2}  \right)$ \\
2285	>	% %
2286	>	% $w_l(r)$ &
2287	>	% $\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)$ &
2288	>	% $\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\
2289	>	% &&$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)
2290	>	% -\frac{2v_{42}(r)}{r}$ \\
2291	>	% %
2292	>	% $w_m(r)$ &
2293	>	% $\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)$ &
2294	>	% $\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$ \\
2295	>	% &&$\left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}
2296	>	% +\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $ \\
2297	>	% &&$-\frac{4v_{43}(r)}{r}$ \\
2298	>	% %
2299	>	% $w_n(r)$ &
2300	>	% $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
2301	>	% $\frac{v_{42}(r)}{r}$ \\
2302	>	% %
2303	>	% $w_o(r)$ &
2304	>	% $\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)$ &
2305	>	% $\frac{v_{43}(r)}{r}$ \\
2306	>	% %
2307	>	% \end{tabular}
2308	>	% \end{ruledtabular}
2309	>	% \end{table*}
2310
2311		\newpage
2312

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines