trunk/multipole/multipole1.tex

% ****** Start of file aipsamp.tex ******
%
%   This file is part of the AIP files in the AIP distribution for REVTeX 4.
%   Version 4.1 of REVTeX, October 2009
%
%   Copyright (c) 2009 American Institute of Physics.
%
%   See the AIP README file for restrictions and more information.
%
% TeX'ing this file requires that you have AMS-LaTeX 2.0 installed
% as well as the rest of the prerequisites for REVTeX 4.1
%
% It also requires running BibTeX. The commands are as follows:
%
%  1)  latex  aipsamp
%  2)  bibtex aipsamp
%  3)  latex  aipsamp
%  4)  latex  aipsamp
%
% Use this file as a source of example code for your aip document.
% Use the file aiptemplate.tex as a template for your document.
\documentclass[%
 aip,
 jmp,
 amsmath,amssymb,
preprint,%
% reprint,%
%author-year,%
%author-numerical,%
]{revtex4-1}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math
%\usepackage[mathlines]{lineno}% Enable numbering of text and display math
%\linenumbers\relax % Commence numbering lines

\begin{document}

\preprint{AIP/123-QED}

\title[Generalization of the Shifted-Force Potential to Higher-Order Potentials]
{Generalization of the Shifted-Force Potential to Higher-Order Potentials}

\author{Madan Lamichhane}
\affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}

\author{J. Daniel Gezelter}
 \email{gezelter@nd.edu.}
\affiliation{Department of Chemistry and Biochemistry, University
of Notre Dame, Notre Dame, IN 46556}

\author{Kathie E. Newman}
\affiliation{Department of Physics, University
of Notre Dame, Notre Dame, IN 46556}


\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
Over the past several years, there has been increasing interest
in real space methods for calculating electrostatic interactions
in computer simulations of condensed molecular systems. We
have extended our original damped-shifted force (DSF)
electrostatic kernel and have been able to derive a set of
interaction models for higher-order multipoles based on
truncated Taylor expansions around the cutoff. For multipolemultipole
interactions, we find that each of the distinct
orientational contributions has a separate radial function to
ensure that the overall forces and torques vanish at the cutoff
radius.
\end{abstract}

\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
                             % Classification Scheme.
\keywords{Suggested keywords}%Use showkeys class option if keyword
                              %display desired
\maketitle

\section{Introduction}

The Coulomb electrostatic interaction is of importance in a number of physical chemistry problems 
[background needed, do we mention gases, liquids, solids?].

[...]

The methods that we develop in this paper are meant specifically for problems involving interacting rigid molecules which will be described
in terms of classical mechanics and electrodynamics.  From mechanics, the molecule's mass distribution
determines its total mass and moment of inertia tensor.
From electrostatics, its charge distribution is conveniently described using multipoles.
Our goal is to advance methods for handling inter-molecular interactions in molecular dynamics simulations.
To do this, we must develop consistent approximate equations for 
interaction energies, forces, and torques. 

[...]

This paper extends the shifted-force potential method 
of Fennel and Gezelter to higher-order multipole interactions.  [describe?]
Extending an idea from Wolf, multipole images are placed on the surface of a 
``cutoff'' sphere of radius $r_c$. These images are used to make all  interaction energies, forces, and torques
be zero for $r < r_c$.
Two such methods have been developed, both based on Taylor-series expansions.  
The first is applied to the Coulomb kernel of the multipole expansion.  The second is 
applied to individual terms for interaction energies in the multipole expansion.
Because of differences in the initial assumptions, the two methods yield different results.

Also explored here is the effect of replacing the bare Coulomb kernel with that of a smeared
charge distribution.  Thus four methods are compared in this paper:  
(1) Shifted force, Coulomb, method 1; 
(2) Sihfted force, Coulomb, method 2;
(3) Shifted force, smeared charge, method 1; and
(4) Shifted force, smeared charge, method 2.
The last of these methods is our preferred method and is called the Extended Shifted Force Method.  
Subsequent papers will apply this method to various problems of physical and chemical interest.

[...]


\section{Development of the Methods}

\subsection{Multipole Expansion}

Consider two discrete rigid collections of atoms and ions, denoted as objects $\bf a$ and $\bf b$.  
In the following, we assume
that the two objects only interact via electrostatics and describe those interactions in terms of 
a multipole expansion.  Putting the origin of the coordinate system at the center of mass of $\bf a$, we use
vectors $\mathbf{r}_k$ to denote the positions of all charges $q_k$ in $\bf a$. 
Then the electrostatic potential of object $\bf a$ at $\mathbf{r}$ is given by
\begin{equation}
V_a(\mathbf r) = \frac{1}{4\pi\epsilon_0}
 \sum_{k \, \text{in \bf a}} \frac{q_k}{\lvert \mathbf{r} - \mathbf{r}_k \rvert}.
\end{equation}
We write the Taylor series expansion in $r$ using an implied summation notation,
Greek indices are used to indicate space coordinates $x$, $y$, $z$ and the subscripts
$k$ and $j$ are reserved for labelling specific charges in $\bf a$ and $\bf b$ respectively, and find:
\begin{equation}
 \frac{1}{\lvert \mathbf{r} - \mathbf{r}_k \rvert} = 
\left( 1
-  r_{k\alpha} \frac{\partial}{\partial r_{\alpha}} 
+ \frac{1}{2}  r_{k\alpha} r_{k\beta} \frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} +\dots 
\right)
\frac{1}{r}  .
\end{equation}
We then follow Smith in defining an operator for the expansion:
\begin{equation}
V_{\bf a}(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\hat{M}_{\bf a} \frac{1}{r} 
\end{equation}
where
\begin{equation}
\hat{M}_{\bf a} = C_{\bf a} - D_{{\bf a}\alpha} \frac{\partial}{\partial r_{\alpha}} 
+  Q_{{\bf a}\alpha\beta}
 \frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots
\end{equation}
and the charge $C_{\bf a}$, dipole $D_{{\bf a}\alpha}$, 
and quadrupole $Q_{{\bf a}\alpha\beta}$ are defined by
\begin{equation}
C_{\bf a}=\sum_{k \, \text{in \bf a}} q_k ,
\end{equation}
\begin{equation}
D_{{\bf a}\alpha} = \sum_{k \, \text{in \bf a}} q_k r_{k\alpha} ,
\end{equation}
\begin{equation}
Q_{{\bf a}\alpha\beta} = \frac{1}{2} \sum_{k \, \text{in \bf a}} q_k r_{k\alpha}  r_{k\beta} .
\end{equation}

It is convenient to locate charges $q_j$ relative to the center of mass of  $\bf b$.  Then with $\bf{r}$ pointing from
$\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $), the interaction energy is given by
\begin{eqnarray}
U_{\bf{ab}}(r) =&& \frac{1}{4\pi \epsilon_0}
\sum_{k \, \text{in \bf a}} \sum_{j \, \text{in \bf b}}
\frac{q_k q_j}{\vert \bf{r}_k - (\bf{r}+\bf{r}_j) \vert} \nonumber\\
=&& \frac{1}{4\pi \epsilon_0}
\sum_{k \, \text{in \bf a}} \sum_{j \, \text{in \bf b}}
\frac{q_k q_j}{\vert \bf{r}+ (\bf{r}_j-\bf{r}_k) \vert} \nonumber\\
=&&\frac{1}{4\pi \epsilon_0}  \sum_{j \, \text{in \bf b}}  q_j V_a(\bf{r}+\bf{r}_j) \nonumber\\
=&&\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} .
\end{eqnarray}
The last expression can also be expanded as a Taylor series in $r$.  Using a notation similar to before, we define
\begin{equation}
\hat{M}_{\bf b} = C_{\bf b} + D_{{\bf b}\alpha} \frac{\partial}{\partial r_{\alpha}} 
+   Q_{{\bf b}\alpha\beta}
 \frac{\partial^2}{\partial r_{\alpha} \partial r_{\beta}} + \dots
\end{equation}
and
\begin{equation}
U_{\bf{ab}}(r)=\frac{\hat{M}_{\bf a} \hat{M}_{\bf b}}{4\pi \epsilon_0} \frac{1}{r}  \label{kernel}.
\end{equation}
Note the ease of separting out the respective energies of interaction of the charge, dipole, and quadrupole of $\bf a$  from those of $\bf b$.

\subsection{Bare Coulomb versus smeared charge}

With the four types of methods being considered here, it is desirable to list the approximations in as transparent a form
as possible.  First, one may use the bare Coulomb potential, with radial dependence $1/r$, 
as shown in Eq.~(\ref{kernel}).  Alternatively, one may use
a smeared charge distribution, then the``kernel'' $1/r$ of the expansion is replaced with a function:
\begin{equation}
B_0(r)=\frac{\text{erfc}(\alpha r)}{r} = \frac{2}{\sqrt{\pi}r}
\int_{\alpha r}^{\infty} \text{e}^{-s^2} ds .
\end{equation}
We develop equations using a function $f(r)$ to represent either $1/r$ or $B_0(r)$, dependent on the the type of approach being considered.
Smith's convenient functions $B_l(r)$ are summarized in Appendix A.

\subsection{Shifting the force, method 1}

As discussed in the introduction, it is desirable to cutoff the electrostatic energy at a radius
$r_c$.  For consistency in approximation, we want the interaction energy as well as the force and
torque to go to zero at $r=r_c$.  
To describe how this goal may be met using a radial approximation, we use two examples, charge-charge
and charge-dipole, using the bare Coulomb kernel $f(r)=1/r$ to explain the idea.

In the shifted-force approximation, the interaction energy $U_{\bf{ab}}(r_c)=0$ 
for two charges $C_{\bf a}$ and $C_{\bf b}$ separated by a distance $r$ is written:
\begin{equation}
U_{C_{\bf a}C_{\bf b}}(r)=\frac{1}{4\pi \epsilon_0} C_{\bf a} C_{\bf b}
\left({ \frac{1}{r} - \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2}  }
\right) .
\end{equation}
Two shifting terms appear in this equations because we want the force to
also be shifted due to an image charge located at a distance $r_c$.  
Since one derivative of the interaction energy is needed for the force, we want a term
linear in $(r-r_c)$ in the interaction energy, that is:
\begin{equation}
\frac{d\,}{dr} 
\left( {\frac{1}{r} - \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2}  }
\right) = \left(- \frac{1}{r^2} + \frac{1}{r_c^2} 
\right) .
\end{equation}
This demonstrates the need of the third term in the brackets of the energy expression, but  leads to the question, how does this idea generalize for higher-order multipole energies?

In method 1, the procedure that we follow is based on the number of derivatives need for each energy, force, or torque.  That is,
a quadrupole-quadrupole interaction energy will have four derivatives,
$\partial^4/\partial r_\alpha \partial r_\beta \partial r_\gamma \partial r_\delta$,
and the force or torque will bring in yet another derivative.  
We thus want shifted energy expressions to include terms so that all energies, forces, and torques
are zero at $r=r_c$.  In each case, we will subtract off a function $f_n^{\text{shift}}(r)$ from the
kernel $f(r)=1/r$.  The index $n$ indicates the number of derivatives to be taken when
deriving a given multipole energy.
We choose a function with guaranteed smooth derivatives --- a truncated Taylor series of the function
$f(r)$, e.g., 
%
\begin{equation}
f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)} \Big \lvert  _{r_c}  .
\end{equation}
%
The combination of $f(r)$ with the shifted function is denoted $f_n(r)=f(r)-f_n^{\text{shift}}(r)$.
Thus, for $f(r)=1/r$, we find
%
\begin{equation}
f_1(r)=\frac{1}{r}- \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2} - \frac{(r-r_c)^2}{r_c^3} .
\end{equation}
%
Continuing with the example of a charge  $\bf a$  interacting with a dipole  $\bf b$, we write
%
\begin{equation}
U_{C_{\bf a}D_{\bf b}}(r)=
\frac{C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}  \frac {\partial f_1(r) }{\partial r_\alpha} 
=\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}   
\frac {r_\alpha}{r} \frac {\partial f_1(r)}{\partial r} .
\end{equation}
%
The force that dipole  $\bf b$ puts on charge $\bf a$ is
%
\begin{equation}
F_{C_{\bf a}D_{\bf b}\beta} =\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}
\left[ \frac{\delta_{\alpha\beta}}{r} \frac {\partial}{\partial r} + 
\frac{r_\alpha r_\beta}{r^2}
\left( -\frac{1}{r} \frac {\partial} {\partial r} 
+ \frac {\partial ^2} {\partial r^2} \right) \right] f_1(r) .
\end{equation}
%
For $f(r)=1/r$, we find
%
\begin{equation}
F_{C_{\bf a}D_{\bf b}\beta} =
\frac{C_{\bf a} D_{{\bf b}\beta} }{4\pi \epsilon_0r}
\left[  -\frac{1}{r^2}+\frac{1}{r_c^2}-\frac{2(r-r_c)}{r_c^3} \right] 
+\frac{C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta }{4\pi \epsilon_0}
\left[ \frac{3}{r^5}-\frac{3}{r^3r_c^2} \right] .
\end{equation}
%
This expansion shows the expected $1/r^3$ dependence of the force.  

In general, we write
%
\begin{equation}
U=\frac{1}{4\pi \epsilon_0} (\text{prefactor}) (\text{derivatives}) f_n(r)
\label{generic}
\end{equation}
%
where $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole
and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for quadrupole-quadrupole.
An example is the case of quadrupole-quadrupole for which the $\text{prefactor}$ is 
$Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$ and the derivatives are
$\partial^4/\partial r_\alpha \partial r_\beta \partial r_\gamma \partial r_\delta$, with
implied summation combining the space indices.

To apply this method  to the smeared-charge approach,  
we write $f(r)=\text{erfc}(\alpha r)/r$.  By using one function $f(r)$ for both
approaches, we simplify the tabulation of equations used.  Because
of the many derivatives that are taken, the algebra is tedious and are summarized
in Appendices A and B.  

\subsection{Shifting the force, method 2}

Note the method used in the previous subsection to shift a force is basically that of using
a truncated Taylor Series in the radius $r$.  An alternate method exists, best explained by
writing one shifted formula for all interaction energies $U(r)$:
\begin{equation}
U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big \lvert  _{r_c} .
\end{equation}
Note that this method uses only the linear term, $(r-r_c)$ in the Taylor series, no higher order terms
$(r-r_c)^n$ appear.   The primary difference between methods 1 and 2 originates
with the stage in the derivation where the Taylor Series is applied; in method 1, it is applied to the 
kernel.  In method 2, it is applied to individual interaction energies of the multipole expansion.
Terms from this method thus have the general form:
\begin{equation}
U=\frac{1}{4\pi \epsilon_0}\sum  (\text{angular factor}) (\text{radial factor}).
\label{generic2}
\end{equation}

Results for both methods can be summarized using the form of Eq.~(\ref{generic2})
and are listed in Table I below.

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

Up to this point, all energies and forces have been written in terms of fixed space
coordinates $x$, $y$, $z$.  Interaction energies are computed from the generic formulas Eq.~(\ref{generic}) and ~(\ref{generic2}) which 
combine prefactors with radial functions.  But because objects
 $\bf a$ and  $\bf b$ both translate and rotate as part of a MD simulation,
it is desirable to contract all $r$-dependent terms with dipole and quadrupole
moments expressed in terms of their body axes.  
Since the interaction energy expressions will be used to derive both forces and torques,
we follow here the development of Allen and Germano, which was itself based on an 
earlier paper by Price \em et al.\em

Denote body axes for objects  $\bf a$ and  $\bf b$ by unit vectors
$\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$ referring to a convenient
set of inertial body axes.  (Note, these body axes are generally not the same as those for which the
quadrupole moment is diagonal.)  Then,
%
\begin{eqnarray}
\hat{a}_m= a_{mx}\hat{x} + a_{my}\hat{y} + a_{mz}\hat{z}  \\
\hat{b}_m= b_{mx}\hat{x} + b_{my}\hat{y} + b_{mz}\hat{z}  .
\end{eqnarray}
Allen and Germano define matrices $\hat{\mathbf {a}}$ 
and  $\hat{\mathbf {b}}$ using these unit vectors:
\begin{eqnarray}
\hat{\mathbf {a}} = 
\begin{pmatrix}
\hat{a}_1 \\
\hat{a}_2 \\
\hat{a}_3
\end{pmatrix}
=
\begin{pmatrix}
a_{1x} \quad a_{1y} \quad a_{1z} \\
a_{2x} \quad a_{2y} \quad a_{2z} \\
a_{3x} \quad a_{3y} \quad a_{3z}
\end{pmatrix}\\
\hat{\mathbf {b}} = 
\begin{pmatrix}
\hat{b}_1 \\
\hat{b}_2 \\
\hat{b}_3
\end{pmatrix}
=
\begin{pmatrix}
b_{1x}\quad  b_{1y} \quad b_{1z} \\
b_{2x} \quad b_{2y} \quad b_{2z} \\
b_{3x} \quad b_{3y} \quad b_{3z}
\end{pmatrix}  .
\end{eqnarray}
%
These matrices convert from space-fixed $(xyz)$ to object-fixed $(123)$ coordinates.
All contractions of prefactors with derivatives of functions can be written in terms of these matrices.
It proves to be equally convenient to just write any contraction in terms of unit vectors
$\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$.  
We have found it useful to write angular-dependent terms in three different fashions,
illustrated by the following 
three examples from the interaction-energy expressions:
%
\begin{eqnarray}
 \mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}}
=D_{\bf {a}\alpha} D_{\bf {b}\alpha}=
\sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}} \\
r^2 \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right)=
r_\alpha Q_{\bf b \alpha \beta} r_\beta = r^2
\sum_{mn}(\hat{r} \cdot \hat{b}_m)  Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \\
r ( \mathbf{D}_{\mathbf{a}} \cdot 
\mathbf{Q}_{\mathbf{b}}  \cdot \hat{r})=
D_{\bf {a}\alpha}  Q_{\bf b \alpha \beta} r_\beta
=r \sum_{lmn} D_{\mathbf{a}l} (\hat{a}_l \cdot \hat{b}_m ) Q_{\mathbf{b}mn}
(\hat{b}_n \cdot \hat{r}) .
\end{eqnarray}
%
[Dan, perhaps there are better examples to show here.]

In each line, the first two terms are written using space coordinates.  The first form is standard
in the chemistry literature, and the second is ``physicist style'' using implied summation notation.  The third 
form shows explicitly sums over body indices and the dot products now indicate contractions using space indices.
We find the first form to be useful in writing equations prior to converting to computer code.  The second form is helpful
in derivations of the interaction energy expressions.  The third one is specifically helpful when deriving forces and torques, as will
be discussed below.  

\section{Energies, forces, and torques}
\subsection{Interaction energies}

We now list multipole interaction energies for the four types of approximation.  
A ``generic'' set of radial functions is introduced so to be able to present the results in Table I.  This set of
equations is written in terms of space coordinates:

% Energy in space coordinate form ----------------------------------------------------------------------------------------------
%
%
% u ca cb
%
\begin{equation}
U_{C_{\bf a}C_{\bf b}}(r)=
\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r)  \label{uchch}
\end{equation}
%
% u ca db
%
\begin{equation}
U_{C_{\bf a}D_{\bf b}}(r)=
\frac{C_{\bf a}}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)  v_{11}(r)  
 \label{uchdip}
\end{equation}
%
% u ca qb
%
\begin{equation}
U_{C_{\bf a}Q_{\bf b}}(r)=
\frac{C_{\bf a }}{4\pi \epsilon_0} \Bigl[ \text{Tr}Q_{\bf b} v_{21}(r) 
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{22}(r) \Bigr]
\label{uchquad}
\end{equation}
%
% u da cb
%
\begin{equation}
U_{D_{\bf a}C_{\bf b}}(r)=
-\frac{C_{\bf b}}{4\pi \epsilon_0}  
\left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)   v_{11}(r) \label{udipch}
\end{equation}
%
% u da db
%
\begin{equation}
U_{D_{\bf a}D_{\bf b}}(r)=
-\frac{1}{4\pi \epsilon_0} \Bigr[ \left( \mathbf{D}_{\mathbf {a}} \cdot
\mathbf{D}_{\mathbf{b}} \right)  v_{21}(r)
+\left( \mathbf{D}_{\mathbf {a}} \cdot \hat{r} \right)
\left( \mathbf{D}_{\mathbf {b}} \cdot \hat{r} \right)  
v_{22}(r) \Bigr]
\label{udipdip}
\end{equation}
%
% u da qb
%
\begin{equation}
\begin{split}
% 1
U_{D_{\bf a}Q_{\bf b}}(r)&=
-\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{b}}
\left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
+2 ( \mathbf{D}_{\mathbf{a}} \cdot
\mathbf{Q}_{\mathbf{b}} \cdot \hat{r} ) \Bigr] v_{31}(r) \\
% 2
&-\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
 \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{32}(r)
\label{udipquad}
\end{split}
\end{equation}
%
% u qa cb
%
\begin{equation}
U_{Q_{\bf a}C_{\bf b}}(r)=
\frac{C_{\bf b }}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\bf a}  v_{21}(r)
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right)  v_{22}(r)  \Bigr]
\label{uquadch}
\end{equation}
%
% u qa db
%
\begin{equation}
\begin{split}
%1
U_{Q_{\bf a}D_{\bf b}}(r)&=
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{a}}
\left(  \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)
+2 ( \mathbf{D}_{\mathbf{b}} \cdot 
\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r}) \Bigr] v_{31}(r)
% 2
+\frac{1}{4\pi \epsilon_0}
\left(  \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) v_{32}(r)
\label{uquaddip}
\end{split}
\end{equation}
%
% u qa qb
%
\begin{equation}
\begin{split}
%1
U_{Q_{\bf a}Q_{\bf b}}(r)&=
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr} \mathbf{Q}_{\mathbf{a}} \text{Tr} \mathbf{Q}_{\mathbf{b}}
+2 \text{Tr} \left(
\mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} \right) \Bigr] v_{41}(r)
\\
% 2
&+ \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
 \left( \hat{r} \cdot 
\mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right)
+\text{Tr}\mathbf{Q}_{\mathbf{b}}
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}}
 \cdot \hat{r} \right)  +4 (\hat{r}  \cdot
\mathbf{Q}_{{\mathbf a}}\cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
\Bigr] v_{42}(r)
 \\
% 4
&+ \frac{1}{4\pi \epsilon_0}
\left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right)
\left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}}  \cdot \hat{r} \right) v_{43}(r).
\label{uquadquad}
\end{split}
\end{equation}


%
%
% TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
%

\begin{table*}
\caption{\label{tab:tableenergy}Radial functions used in the energy and torque equations.  Functions
used in this table are defined in Appendices B and C.}
\begin{ruledtabular}
\begin{tabular}{cccc}
Generic&Coulomb&Method 1&Method 2
\\ \hline
%
%
%
%Ch-Ch&
$v_{01}(r)$ &
$\frac{1}{r}$ &
$f_0(r)$ &
$f(r)-f(r_c)-(r-r_c)g(r_c)$
\\
%
%
%
%Ch-Di&
$v_{11}(r)$ &
$-\frac{1}{r^2}$ &
$g_1(r)$ &
$g(r)-g(r_c)-(r-r_c)h(r_c)$ \\
%
%
%
%Ch-Qu/Di-Di&
$v_{21}(r)$ &
$-\frac{1}{r^3}  $ & 
$\frac{g_2(r)}{r} $ &
$\frac{g(r)}{r}-\frac{g(r_c)}{r_c} -(r-r_c)
\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right)$ \\
$v_{22}(r)$ &
$\frac{3}{r^3}  $ &
$\left(-\frac{g_2(r)}{r} + h_2(r) \right)$ &
$\left(-\frac{g(r)}{r}+h(r) \right)
-\left(-\frac{g(r_c)}{r_c}+h(r_c) \right) $ \\
&&&$  -(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$ 
\\
%
%
%
%Di-Qu &
$v_{31}(r)$ &
$\frac{3}{r^4}  $ &
$\left(-\frac{g_3(r)}{r^2} + \frac{h_3(r)}{r} \right)$ &
$\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r} \right)
-\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right) $\\
&&&$  -(r-r_c) \left(\frac{2g(r_c)}{r_c^3}-\frac{2h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ 
\\
% 
$v_{32}(r)$ &
$-\frac{15}{r^4}  $ &
$\left( \frac{3g_3(r)}{r^2} - \frac{3h_3(r)}{r} + s_3(r) \right)$ &
$\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)
- \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c) \right)$ \\
&&&$  -(r-r_c) \left( \frac{-6g(r_c)}{r_c^3}+\frac{6h(r_c)}{r_c^2}-\frac{3s(r_c)}{r_c}+t(r_c) \right)$ 
\\
%
%
%
%Qu-Qu&
$v_{41}(r)$ &
$\frac{3}{r^5} $ &
$\left(-\frac{g_4(r)}{r^3} +\frac{h_4(r)}{r^2} \right) $  &
$\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right)
- \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2} \right)$ \\
&&&$  -(r-r_c) \left( \frac{3g(r_c)}{r_c^4}-\frac{3h(r_c)}{r_c^3}+\frac{s(r_c)}{r_c^2} \right)$ 
\\
% 2
$v_{42}(r)$ &
$- \frac{15}{r^5}   $ &
$\left( \frac{3g_4(r)}{r^3} - \frac{3h_4(r)}{r^2}+\frac{s_4(r)}{r} \right)$ &
$\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right)
-\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ \\
&&&$  -(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}
-\frac{4s(r_c)}{r_c^2} + \frac{t(r_c)}{r_c}\right)$ 
\\
% 3
$v_{43}(r)$ &
$ \frac{105}{r^5}  $ &
$\left(-\frac{15g_4(r)}{r^3}+\frac{15h_4(r)}{r^2}-\frac{6s_4(r)}{r} + t_4(r)\right) $ &
$\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s(r)}{r} + t(r)\right)$ \\
&&&$ -\left(-\frac{15g(r_c)}{r_c^3}+\frac{15h(r_c)}{r_c^2}-\frac{6s(r_c)}{r_c} + t(r_c)\right)$ \\
&&&$ -(r-r_c)\left(\frac{45g(r_c)}{r_c^4}-\frac{45h(r_c)}{r_c^3}+\frac{21s(r_c)}{r_c^2}
-\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\
\end{tabular}
\end{ruledtabular}
\end{table*}
%
%
% FORCE  TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
%

\begin{table}
\caption{\label{tab:tableFORCE}Radial functions used in the force equations.}
\begin{ruledtabular}
\begin{tabular}{cc}
Generic&Method 1 or Method 2
\\ \hline
%
%
%
$w_a(r)$&
$\frac{d v_{01}}{dr}$ \\
%
%
$w_b(r)$ &
$\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $ \\
%
$w_c(r)$ &
$\frac{v_{11}(r)}{r}$ \\
%
%
$w_d(r)$&
$\frac{d v_{21}}{dr}$ \\
%
$w_e(r)$ &
$\frac{v_{22}(r)}{r}$ \\
% 
%
$w_f(r)$&
$\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$\\
%
$w_g(r)$& 
$\frac{v_{31}(r)}{r}$\\
%
$w_h(r)$ &
$\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$\\
% 2
$w_i(r)$ &
$\frac{v_{32}(r)}{r}$ \\
% 
$w_j(r)$ &
$\frac{d v_{32}}{dr}  - \frac{3v_{32}}{r}$ \\
%
$w_k(r)$ &
$\frac{d v_{41}}{dr} $  \\
%
$w_l(r)$ &
$\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ \\
%
$w_m(r)$ &
$\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$ \\
%
$w_n(r)$ &
$\frac{v_{42}(r)}{r}$ \\
%
$w_o(r)$ &
$\frac{v_{43}(r)}{r}$ \\
%

\end{tabular}
\end{ruledtabular}
\end{table}
%
%
%

\subsection{Forces}

The force $\mathbf{F}_{\bf a}$ on $\bf{a}$ due to $\bf{b}$ is the negative of
the force  $\mathbf{F}_{\bf b}$ on $\bf{b}$ due to $\bf{a}$.  For a simple charge-charge
interaction, these forces will point along the $\pm \hat{r}$ directions, where 
$\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $.  Thus
%
\begin{equation}
F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}}{\partial r} 
\quad \text{and} \quad  F_{\bf b \alpha} 
= - \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}} {\partial r}  .
\end{equation}
%
The concept of obtaining a force from an energy by taking a gradient is the same for
higher-order multipole interactions, the trick is to make sure that all
$r$-dependent derivatives are considered.
As is pointed out by Allen and Germano, this is straightforward if the
interaction energies are written recognizing explicit 
$\hat{r}$ and body axes ($\hat{a}_m$, $\hat{b}_n$) dependences:
%
\begin{equation}
U(r,\{\hat{a}_m \cdot \hat{r} \}, 
\{\hat{b}_n\cdot \hat{r} \} 
\{\hat{a}_m \cdot \hat{b}_n \}) .
\label{ugeneral}
\end{equation}
%
Then,
%
\begin{equation}
\mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} =  \frac{\partial U}{\partial \mathbf{r}} 
= \frac{\partial U}{\partial r} \hat{r} 
 + \sum_m \left[ 
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} 
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} 
+ \frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} 
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} 
\right] \label{forceequation}.
\end{equation}
%
Note our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ is opposite
that of Allen and Germano.  In simplifying the algebra, we also use:
%
\begin{eqnarray}
\frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}} 
= \frac{1}{r} \left(  \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r}
\right) \\
\frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}} 
= \frac{1}{r} \left(  \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r} 
\right) .
\end{eqnarray}
%
We list below the force equations written in terms  of space coordinates.  The
radial functions used in the two methods are listed in Table II.
%
%SPACE COORDINATES FORCE EQUTIONS
%
% **************************************************************************
% f ca cb
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} = 
\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r} 
\end{equation}
%
%
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} = 
\frac{C_{\bf a}}{4\pi \epsilon_0} \Bigl[
 \left( \hat{r} \cdot \mathbf{D}_{\mathbf{b}} \right) 
w_b(r) \hat{r}  
+ \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr]
\end{equation}
%
%
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} = 
\frac{C_{\bf a }}{4\pi \epsilon_0} \Bigr[
\text{Tr}\mathbf{Q}_{\bf b} w_d(r) \hat{r}
+ 2  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} w_e(r)
 + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
\end{equation}
%
%
%
\begin{equation}
\mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} = 
-\frac{C_{\bf{b}}}{4\pi \epsilon_0} \Bigl[
\left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) w_b(r) \hat{r}
+ \mathbf{D}_{\mathbf{a}} w_c(r) \Bigr]
\end{equation}
%
%
%
\begin{equation}
\mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
- \mathbf{D}_{\mathbf {a}} \cdot  \mathbf{D}_{\mathbf{b}} w_d(r) \hat{r}
+ \left( \mathbf{D}_{\mathbf {a}} 
\left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)
+ \mathbf{D}_{\mathbf {b}} \left( \mathbf{D}_{\mathbf{a}}  \cdot \hat{r} \right) \right) w_e(r)
% 2
 - \left( \hat{r} \cdot \mathbf{D}_{\mathbf {a}} \right) 
\left( \hat{r} \cdot \mathbf{D}_{\mathbf {b}} \right) w_f(r) \hat{r} \Bigr]
\end{equation}
%
%
%
\begin{equation}
\begin{split}
\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} = 
& - \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{b}} \mathbf{ D}_{\mathbf{a}}
+2 \mathbf{D}_{\mathbf{a}} \cdot 
\mathbf{Q}_{\mathbf{b}} \Bigr] w_g(r)
 - \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{b}}
\left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) 
+2 ( \mathbf{D}_{\mathbf{a}} \cdot 
\mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
% 3
& - \frac{1}{4\pi \epsilon_0} \Bigl[\mathbf{ D}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) 
+2  (\hat{r} \cdot \mathbf{D}_{\mathbf{a}} ) (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} )  \Bigr]
w_i(r)
% 4
 -\frac{1}{4\pi \epsilon_0} 
(\hat{r} \cdot \mathbf{D}_{\mathbf{a}} )
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r}  
\end{split}
\end{equation}
%
%
\begin{equation}
\mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} = 
\frac{C_{\bf b }}{4\pi \epsilon_0} \Bigr[
\text{Tr}\mathbf{Q}_{\bf a} w_d(r) \hat{r}
+ 2  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} w_e(r)
 + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
\end{equation}
%
\begin{equation}
\begin{split}
\mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} = 
&\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{a}} \mathbf{D}_{\mathbf{b}} 
+2 \mathbf{D}_{\mathbf{b}} \cdot \mathbf{Q}_{\mathbf{a}}  \Bigr] w_g(r)
% 2
+ \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
(\hat{r} \cdot  \mathbf{D}_{\mathbf{b}})
+2 (\mathbf{D}_{\mathbf{b}} \cdot 
\mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
% 3
&+ \frac{1}{4\pi \epsilon_0} \Bigl[ \mathbf{D}_{\mathbf{b}} 
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) 
+2  (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
(\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} ) \Bigr]   w_i(r)
% 4
+\frac{1}{4\pi \epsilon_0} 
(\hat{r} \cdot \mathbf{D}_{\mathbf{b}}) 
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}}  \cdot \hat{r}) w_j(r) \hat{r} 
\end{split}
\end{equation}
%
%
%
\begin{equation}
\begin{split}
\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} = 
+\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} \hat{r}
+ 2 \text{Tr} ( \mathbf{Q}_{\mathbf{a}} \cdot  \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_k(r) \hat{r} \\
% 2
+\frac{1}{4\pi \epsilon_0} \Bigl[
2\text{Tr}\mathbf{Q}_{\mathbf{b}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )  
+ 2\text{Tr}\mathbf{Q}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} ) 
% 3
+4 (\mathbf{Q}_{\mathbf{a}}  \cdot  \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})  
+  4(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}}) \Bigr] w_n(r) \\
% 4
+ \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{a}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) 
+ \text{Tr}\mathbf{Q}_{\mathbf{b}}
(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})  
% 5
+4 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot  
\mathbf{Q}_{\mathbf{b}}   \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
%
+ \frac{1}{4\pi \epsilon_0} \Bigl[
+ 2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
%6
+2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \hat{r})
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_o(r) \\
%  7
+ \frac{1}{4\pi \epsilon_0} 
(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r}) 
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) w_m(r) \hat{r} 
\end{split}
\end{equation}
%
%
% TORQUES SECTION -----------------------------------------------------------------------------------------
%
\subsection{Torques}

Following again Allen and Germano, when energies are written in the form
of Eq.~({\ref{ugeneral}), then torques can be expressed as:
%
\begin{eqnarray}
\mathbf{\tau}_{\bf a} =
 \sum_m 
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{r})} 
( \hat{r} \times \hat{a}_m )
-\sum_{mn}
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} 
(\hat{a}_m \times \hat{b}_n) \\
%
\mathbf{\tau}_{\bf b} =
 \sum_m 
\frac{\partial U}{\partial (\hat{b}_m \cdot \hat{r})} 
( \hat{r} \times \hat{b}_m)
+\sum_{mn}
\frac{\partial U}{\partial (\hat{a}_m \cdot \hat{b}_n)} 
(\hat{a}_m \times \hat{b}_n) .
\end{eqnarray}
%
%
Here we list the torque equations written in terms of space coordinates.
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} = 
\frac{C_{\bf a}}{4\pi \epsilon_0}  (\hat{r} \times  \mathbf{D}_{\mathbf{b}}) v_{11}(r) 
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} = 
\frac{2C_{\bf a}}{4\pi \epsilon_0}
\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{22}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =  
-\frac{C_{\bf b}}{4\pi \epsilon_0}  
(\hat{r} \times \mathbf{D}_{\mathbf{a}})  v_{11}(r) 
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} = 
\frac{1}{4\pi \epsilon_0}  \mathbf{D}_{\mathbf {a}}  \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
% 2
-\frac{1}{4\pi \epsilon_0} 
(\hat{r} \times \mathbf{D}_{\mathbf {a}} )
(\hat{r} \cdot \mathbf{D}_{\mathbf {b}} )  v_{22}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} = 
-\frac{1}{4\pi \epsilon_0} \mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
% 2
+\frac{1}{4\pi \epsilon_0} 
(\hat{r} \cdot \mathbf{D}_{\mathbf {a}} )
(\hat{r} \times \mathbf{D}_{\mathbf {b}} ) v_{22}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
-\text{Tr}\mathbf{Q}_{\mathbf{b}}
(\hat{r} \times \mathbf{D}_{\mathbf{a}} )
+2 \mathbf{D}_{\mathbf{a}}  \times 
(\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
\Bigr] v_{31}(r)
% 3
-\frac{1}{4\pi \epsilon_0}
\ (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =
 \frac{1}{4\pi \epsilon_0} \Bigl[
+2 ( \mathbf{D}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} ) \times
\hat{r} 
-2 \mathbf{D}_{\mathbf{a}}  \times 
(\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
\Bigr] v_{31}(r)
% 2
+\frac{2}{4\pi \epsilon_0}
(\hat{r} \cdot \mathbf{D}_{\mathbf{a}})
(\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} =
\frac{1}{4\pi \epsilon_0} \Bigl[
-2 (\mathbf{D}_{\mathbf{b}}  \cdot \mathbf{Q}_{\mathbf{a}} ) \times \hat{r} 
+2 \mathbf{D}_{\mathbf{b}}  \times 
(\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
\Bigr] v_{31}(r)
% 3
- \frac{2}{4\pi \epsilon_0}
(\hat{r} \cdot \mathbf{D}_{\mathbf{b}} )
(\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}}) \times \hat{r} v_{32}(r)
\end{equation}
%
%
%
\begin{equation}
\mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}\mathbf{Q}_{\mathbf{a}}
(\hat{r} \times \mathbf{D}_{\mathbf{b}} )
+2 \mathbf{D}_{\mathbf{b}}  \times
( \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r)
% 2
+\frac{1}{4\pi \epsilon_0}
(\hat{r} \times \mathbf{D}_{\mathbf{b}} )
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) v_{32}(r)
\end{equation}
%
%
%
\begin{equation}
\begin{split}
\mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
&-\frac{4}{4\pi \epsilon_0} 
\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}}
v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[-2\text{Tr}\mathbf{Q}_{\mathbf{b}}
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times \hat{r}
+4 \hat{r} \times 
( \mathbf{Q}_{{\mathbf a}} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) 
% 3
-4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} )\times 
( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} ) \Bigr] v_{42}(r) \\
% 4
&+ \frac{2}{4\pi \epsilon_0}
\hat{r} \times ( \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
\end{split}
\end{equation}
%
%
%
\begin{equation}
\begin{split}
\mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} =  
&\frac{4}{4\pi \epsilon_0} 
\mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}} v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} \Bigl[- 2\text{Tr}\mathbf{Q}_{\mathbf{a}}
(\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} ) \times \hat{r}
-4  (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot 
\mathbf{Q}_{{\mathbf b}} ) \times
\hat{r} 
+4 ( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times 
( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
\Bigr] v_{42}(r) \\
% 4
&+ \frac{2}{4\pi \epsilon_0}
(\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
\hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
\end{split}
\end{equation}
%
%
%
\begin{acknowledgments}
We wish to acknowledge the support of the author community in using
REV\TeX{}, offering suggestions and encouragement, testing new versions,
\dots.
\end{acknowledgments}

\appendix

\section{Smith's $B_l(r)$ functions for smeared-charge distributions}

The following summarizes Smith's $B_l(r)$ functions and
includes formulas given in his appendix.  

The first function $B_0(r)$ is defined by
%
\begin{equation}
B_0(r)=\frac{\text{erfc}(\alpha r)}{r} = \frac{2}{\sqrt{\pi}r}=
\int_{\alpha r}^{\infty} \text{e}^{-s^2} ds .
\end{equation}
%
The first derivative of this function is
%
\begin{equation}
\frac{dB_0(r)}{dr}=-\frac{1}{r^2}\text{erfc}(\alpha r)
-\frac{2\alpha}{r\sqrt{\pi}}\text{e}^{-{\alpha}^2r^2}
\end{equation}
%
and can be rewritten in terms of a function $B_1(r)$:
%
\begin{equation}
B_1(r)=-\frac{1}{r}\frac{dB_0(r)}{dr}
\end{equation}
%
In general,
\begin{equation}
B_l(r)=-\frac{1}{r}\frac{dB_{l-1}(r)}{dr} 
= \frac{1}{r^2} \left[ (2l-1)B_{l-1}(r) + \frac {(2\alpha^2)^l}{\alpha \sqrt{\pi}}
\text{e}^{-{\alpha}^2r^2}
\right] .
\end{equation}
%
Using these formulas, we find
%
\begin{eqnarray}
\frac{dB_0}{dr}=-rB_1(r) \\
\frac{d^2B_0}{dr^2}=-B_1(r) + r^2B_2(r) \\
\frac{d^3B_0}{dr^3}=3rB_2(r) - r^3B_3(r) \\
\frac{d^4B_0}{dr^4}=3B_2(r) - 6r^2B_3(r)+r^4B_4(r) \\
\frac{d^5B_0}{dr^5}=-15rB_3(r) + 10r^3B_4(r) -r^5B_5(r) .
\end{eqnarray}
%
As noted by Smith,
%
\begin{equation}
B_l(r)=\frac{(2l)!}{l!2^lr^{2l+1}} - \frac {(2\alpha^2)^{l+1}}{(2l+1)\alpha \sqrt{\pi}}
+\text{O}(r) .
\end{equation}

\section{Method 1, the $r$-dependent factors}

Using the shifted damped functions $f_n(r)$ defined by:
%
\begin{equation}
f_n(r)= B_0 \Big \lvert  _r -\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} B_0^{(m)} \Big \lvert  _{r_c}  ,
\end{equation}
%
we first provide formulas for successive derivatives of this function.  (If there is 
no damping, then $B_0(r)$ is replaced by $1/r$, as discussed in Section~\ref{damped???}.)  First, we find:
%
\begin{equation}
\frac{\partial f_n}{\partial r_\alpha}=\hat{r}_\alpha \frac{d f_n}{d r} .
\end{equation}
%
This formula clearly brings in derivatives of Smith's $B_0(r)$ function, motivating us to
define higher-order derivatives as  follows:
%
\begin{eqnarray}
g_n(r)= \frac{d f_n}{d r} =
B_0^{(1)} \Big \lvert  _r -\sum_{m=0}^{n} \frac {(r-r_c)^m}{m!} B_0^{(m+1)} \Big \lvert  _{r_c} \\
h_n(r)= \frac{d^2f_n}{d r^2} =
B_0^{(2)} \Big \lvert  _r -\sum_{m=0}^{n-1} \frac {(r-r_c)^m}{m!} B_0^{(m+2)} \Big \lvert  _{r_c} \\
s_n(r)= \frac{d^3f_n}{d r^3} =
B_0^{(3)} \Big \lvert  _r -\sum_{m=0}^{n-2} \frac {(r-r_c)^m}{m!} B_0^{(m+3)} \Big \lvert  _{r_c} \\
t_n(r)= \frac{d^4f_n}{d r^4} =
B_0^{(4)} \Big \lvert  _r -\sum_{m=0}^{n-3} \frac {(r-r_c)^m}{m!} B_0^{(m+4)} \Big \lvert  _{r_c} \\
u_n(r)= \frac{d^5f_n}{d r^5} =
B_0^{(5)} \Big \lvert  _r -\sum_{m=0}^{n-4} \frac {(r-r_c)^m}{m!} B_0^{(m+5)} \Big \lvert  _{r_c}  .
\end{eqnarray}
%
We note that the last function needed (for quadrupole-quadrupole) is
%
\begin{equation}
u_4(r)=B_0^{(5)} \Big \lvert  _r - B_0^{(5)} \Big \lvert  _{r_c} .
\end{equation}

The functions $f_n(r)$ to $u_n(r)$ are recursively computed and stored for values of $r$
from $0$ to $r_c$.  The functions needed are listed schematically below:
%
\begin{eqnarray}
f_0 \quad f_1 \qquad \qquad \quad & \nonumber \\
g_0 \quad g_1 \quad g_2 \quad g_3 \quad &g_4 \nonumber \\
h_1 \quad h_2 \quad h_3 \quad &h_4 \nonumber \\
s_2 \quad s_3 \quad &s_4 \nonumber \\
t_3 \quad &t_4 \nonumber \\
&u_4 \nonumber .
\end{eqnarray}

Using these functions, we find
%
\begin{equation}
\frac{\partial f_n}{\partial r_\alpha} =r_\alpha \frac {g_n}{r} 
\end{equation}
%
\begin{equation}
\frac{\partial^2 f_n}{\partial r_\alpha \partial r_\beta} =\delta_{\alpha \beta}\frac {g_n}{r} 
+r_\alpha r_\beta \left( -\frac{g_n}{r^3} +\frac{h_n}{r^2}\right) 
\end{equation}
%
\begin{equation}
\frac{\partial^3 f_n}{\partial r_\alpha \partial r_\beta r_\gamma} = 
\left( \delta_{\alpha \beta} r_\gamma + \delta_{\alpha \gamma} r_\beta + 
\delta_{ \beta \gamma} r_\alpha \right)  
\left(  -\frac{g_n}{r^3} +\frac{h_n}{r^2} \right)
+ r_\alpha r_\beta r_\gamma 
\left(  \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right)
\end{equation}
%
\begin{eqnarray}
\frac{\partial^4 f_n}{\partial r_\alpha \partial r_\beta r_\gamma r_\delta} = 
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} 
+ \delta_{\alpha \gamma} \delta_{\beta \delta}
 +\delta_{ \beta \gamma} \delta_{\alpha \delta} \right)
\left( - \frac{g_n}{r^3} + \frac{h_n}{r^2} \right)  \nonumber \\ 
+ \left( \delta_{\alpha \beta} r_\gamma r_\delta
+ \text{5 perm}
\right) \left( \frac{3 g_n}{r^5} - \frac{3h_n}{r^4} + \frac{s_n}{r^3} 
\right) \nonumber \\
+ r_\alpha r_\beta r_\gamma r_\delta
\left(  -\frac{15g_n}{r^7} + \frac{15h_n}{r^6} - \frac{6s_n}{r^5}
+ \frac{t_n}{r^4} \right)
\end{eqnarray}
%
\begin{eqnarray}
\frac{\partial^5 f_n}
{\partial r_\alpha \partial r_\beta r_\gamma r_\delta r_\epsilon} = 
\left( \delta_{\alpha \beta} \delta_{\gamma \delta} r_\epsilon
+ \text{14 perm} \right) 
\left(  \frac{3g_n}{r^5}-\frac{3h_n}{r^4} +\frac{s_n}{r^3} \right) \nonumber \\
+ \left( \delta_{\alpha \beta} r_\gamma r_\delta r_\epsilon
+ \text{9 perm}
\right) \left(- \frac{15g_n}{r^7}+\frac{15h_n}{r^7} -\frac{6s_n}{r^5} +\frac{t_n}{r^4} 
\right) \nonumber \\
+ r_\alpha r_\beta r_\gamma r_\delta r_\epsilon
\left(  \frac{105g_n}{r^9} - \frac{105h_n}{r^8} + \frac{45s_n}{r^7}
- \frac{10t_n}{r^6} +\frac{u_n}{r^5} \right)
\end{eqnarray}
%
%
%
\section{Method 2, the $r$-dependent factors}

In method 2, the kernel is not expanded, rather the individual terms in the multipole interaction energies,
see Eq. (20?).  For a smeared-charge distribution, this still brings into the algebra multiple derivatives
of the kernel $B_0(r)$.  To denote these terms, we generalize the notation of the previous appendix.
For $f(r)=1/r$ (bare Coulomb) or $f(r)=B_0(r)$ (smeared charge)
%
\begin{eqnarray}
g(r)= \frac{df}{d r}\\
h(r)= \frac{dg}{d r} = \frac{d^2f}{d r^2} \\
s(r)= \frac{dh}{d r} = \frac{d^3f}{d r^3} \\
t(r)= \frac{ds}{d r} = \frac{d^4f}{d r^4} \\
u(r)= \frac{dt}{d r} =\frac{d^5f}{d r^5} .
\end{eqnarray}
%
For $f(r)=1/r$, Table I lists these derivatives under the column ``Bare Coulomb.''  Checks of algebra can be made by using limiting forms
of equations, e.g., the leading term in the function $g_n(r)$ has $r$ dependence given by $g(r)$.  Equations (B9) to B(13)
are correct for method 2 if one just eliminates the subscript $n$.

\section{Extra Material}
%
%
%Energy in body coordinate form ---------------------------------------------------------------
%
Here are the interaction energies written in terms of the body coordinates:

%
% u ca cb
%
\begin{equation}
U_{C_{\bf a}C_{\bf b}}(r)=
\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r) 
\end{equation}
%
% u ca db
%
\begin{equation}
U_{C_{\bf a}D_{\bf b}}(r)=
\frac{C_{\bf a}}{4\pi \epsilon_0} 
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} \,  v_{11}(r) 
\end{equation}
%
% u ca qb
%
\begin{equation}
U_{C_{\bf a}Q_{\bf b}}(r)=
\frac{C_{\bf a }\text{Tr}Q_{\bf b}}{4\pi \epsilon_0}  
v_{21}(r) \nonumber \\
+\frac{C_{\bf a}}{4\pi \epsilon_0}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r})
v_{22}(r) 
\end{equation}
%
% u da cb
%
\begin{equation}
U_{D_{\bf a}C_{\bf b}}(r)=
-\frac{C_{\bf b}}{4\pi \epsilon_0}  
\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} \,  v_{11}(r) 
\end{equation}
%
% u da db
%
\begin{equation}
\begin{split}
% 1
U_{D_{\bf a}D_{\bf b}}(r)&=
-\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m} 
(\hat{a}_m \cdot   \hat{b}_n)
D_{\mathbf{b}n} v_{21}(r) \\
% 2
&-\frac{1}{4\pi \epsilon_0} 
\sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m} 
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} 
v_{22}(r)
\end{split}
\end{equation}
%
% u da qb
%
\begin{equation}
\begin{split}
% 1
U_{D_{\bf a}Q_{\bf b}}(r)&=
-\frac{1}{4\pi \epsilon_0} \left(
\text{Tr}Q_{\mathbf{b}}
\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
+2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \cdot \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
\right)  v_{31}(r) \\
% 2
&-\frac{1}{4\pi \epsilon_0}
\sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{b}_n \cdot \hat{r}) v_{32}(r)
\end{split}
\end{equation}
%
% u qa cb
%
\begin{equation}
U_{Q_{\bf a}C_{\bf b}}(r)=
\frac{C_{\bf b }\text{Tr}Q_{\bf a}}{4\pi \epsilon_0}  v_{21}(r)
+\frac{C_{\bf b}}{4\pi \epsilon_0}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) v_{22}(r)
\end{equation}
%
% u qa db
%
\begin{equation}
\begin{split}
%1
U_{Q_{\bf a}D_{\bf b}}(r)&=
\frac{1}{4\pi \epsilon_0} \left(
\text{Tr}Q_{\mathbf{a}}
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n}
+2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{b}_l \cdot \hat{a}_m)
Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
\right) v_{31}(r)  \\
% 2
&+\frac{1}{4\pi \epsilon_0}
\sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{a}_n \cdot \hat{r}) v_{32}(r)
\end{split}
\end{equation}
%
% u qa qb
%
\begin{equation}
\begin{split}
%1
U_{Q_{\bf a}Q_{\bf b}}(r)&=
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}}
+2\sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p) 
Q_{\mathbf{a}lm}  Q_{\mathbf{b}np} 
(\hat{a}_m \cdot \hat{b}_n) \Bigr] 
v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[ \text{Tr}Q_{\mathbf{a}}
\sum_{lm} (\hat{r} \cdot \hat{b}_l )
Q_{{\mathbf b}lm}
(\hat{b}_m \cdot \hat{r})
+\text{Tr}Q_{\mathbf{b}}
\sum_{lm} (\hat{r} \cdot \hat{a}_l )
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{r}) \\
% 3
&+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 \cdot \hat{r})
\Bigr] v_{42}(r)  \\
% 4
&+ \frac{1}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \cdot \hat{a}_l) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{r})
\sum_{np}  (\hat{r} \cdot \hat{b}_n) 
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r})  v_{43}(r).
\end{split}
\end{equation}
%


% BODY coordinates force equations --------------------------------------------
%
%
Here are the force equations written in terms of body coordinates.
%
% f ca cb
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} = 
\frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r}
\end{equation}
%
% f ca db
%
\begin{equation}
\mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} = 
\frac{C_{\bf a}}{4\pi \epsilon_0}  
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} w_b(r) \hat{r}
+\frac{C_{\bf a}}{4\pi \epsilon_0}  
\sum_n  D_{\mathbf{b}n} \hat{b}_n w_c(r)
\end{equation}
%
% f ca qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} = 
\frac{1}{4\pi \epsilon_0}  
C_{\bf a }\text{Tr}Q_{\bf b} w_d(r) \hat{r}
+ 2C_{\bf a } \sum_l  \hat{b}_l Q_{{\mathbf b}ln} (\hat{b}_n \cdot \hat{r}) w_e(r) \\
% 2
+\frac{C_{\bf a}}{4\pi \epsilon_0}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r}) w_f(r) \hat{r} 
\end{split}
\end{equation}
%
% f da cb
%
\begin{equation}
\mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} = 
-\frac{C_{\bf{b}}}{4\pi \epsilon_0} 
\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} w_b(r)  \hat{r}
-\frac{C_{\bf{b}}}{4\pi \epsilon_0} 
\sum_n  D_{\mathbf{a}n} \hat{a}_n w_c(r)
\end{equation}
% 
% f da db
%
\begin{equation}
\begin{split}
% 1
\mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} &= 
-\frac{1}{4\pi \epsilon_0}
 \sum_{mn} D_{\mathbf {a}m} 
(\hat{a}_m \cdot   \hat{b}_n)
D_{\mathbf{b}n}  w_d(r) \hat{r}
-\frac{1}{4\pi \epsilon_0} 
\sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m} 
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} w_f(r) \hat{r} \\
% 2
& \quad + \frac{1}{4\pi \epsilon_0} 
\Bigl[ \sum_m D_{\mathbf {a}m} 
\hat{a}_m \sum_n D_{\mathbf{b}n} 
(\hat{b}_n \cdot \hat{r}) 
+ \sum_m D_{\mathbf {b}m} 
\hat{b}_m \sum_n D_{\mathbf{a}n} 
(\hat{a}_n \cdot \hat{r}) \Bigr] w_e(r)  \\
\end{split}
\end{equation}
%
% f da qb
%
\begin{equation}
\begin{split}
% 1
&\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} = 
 - \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{b}}
\sum_l  D_{\mathbf{a}l} \hat{a}_l
+2\sum_{lmn} D_{\mathbf{a}l} 
(\hat{a}_l \cdot \hat{b}_m)
Q_{\mathbf{b}mn} \hat{b}_n  \Bigr] w_g(r) \\
% 3
&  - \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{b}}
\sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
+2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \cdot \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
% 4
&+ \frac{1}{4\pi \epsilon_0}
\Bigl[\sum_l  D_{\mathbf{a}l} \hat{a}_l
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{b}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{a}_l)
 D_{\mathbf{a}l}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn} \hat{b}_n \Bigr]   w_i(r)\\
% 6
&  -\frac{1}{4\pi \epsilon_0} 
\sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l} 
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{b}_n \cdot \hat{r})  w_j(r)  \hat{r}
\end{split}
\end{equation}
%
% force qa cb
%
\begin{equation}
\begin{split}
% 1
\mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0}  
C_{\bf b }\text{Tr}Q_{\bf a} \hat{r} w_d(r)
+ \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) \\
% 2
&  +\frac{C_{\bf b}}{4\pi \epsilon_0}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) w_f(r) \hat{r} 
\end{split}
\end{equation}
%
% f qa db 
%
\begin{equation}
\begin{split}
% 1
&\mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{a}}
\sum_l  D_{\mathbf{b}l} \hat{b}_l
+2\sum_{lmn} D_{\mathbf{b}l} 
(\hat{b}_l \cdot \hat{a}_m)
Q_{\mathbf{a}mn} \hat{a}_n  \Bigr] 
w_g(r)\\
% 3
&  + \frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{a}}
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} 
+2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{b}_l \cdot \hat{a}_m)
Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
% 4
&  + \frac{1}{4\pi \epsilon_0} \Bigl[ \sum_l  D_{\mathbf{b}l} \hat{b}_l
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{a}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{b}_l)
 D_{\mathbf{b}l}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn} \hat{a}_n \Bigr]   w_i(r) \\
% 6
&  +\frac{1}{4\pi \epsilon_0} 
\sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l} 
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{a}_n \cdot \hat{r})  w_j(r)  \hat{r}
\end{split}
\end{equation}
%
% f qa qb
%
\begin{equation}
\begin{split}
&\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} = 
\frac{1}{4\pi \epsilon_0} \Bigl[
\text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}} 
+ 2 \sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p) 
Q_{\mathbf{a}lm}  Q_{\mathbf{b}np} 
(\hat{a}_m \cdot \hat{b}_n) \Bigr] w_k(r) \hat{r}\\
&+\frac{1}{4\pi \epsilon_0} \Bigl[
2\text{Tr}Q_{\mathbf{b}} \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  \hat{a}_m  
+ 2\text{Tr}Q_{\mathbf{a}} \sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm}  \hat{b}_m \\
&+ 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})   
+ 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
\Bigr] w_n(r) \\
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[ \text{Tr}Q_{\mathbf{a}}
\sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm} (\hat{b}_m \cdot \hat{r}) 
+ \text{Tr}Q_{\mathbf{b}}
\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  (\hat{a}_m \cdot \hat{r}) \\
&+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 \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
%
&+\frac{1}{4\pi \epsilon_0} \Bigl[
2\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} \hat{a}_m
\sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_n \cdot \hat{r}) \\
&+2 \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
\sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} \hat{b}_n \Bigr] w_o(r) \hat{r} \\
&  + \frac{1}{4\pi \epsilon_0} 
\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
\sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_p \cdot \hat{r}) w_m(r) \hat{r} 
\end{split}
\end{equation}
%
Here we list the form of the non-zero damped shifted multipole torques showing
explicitly dependences on body axes:
%
%  t ca db
%
\begin{equation}
\mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} = 
\frac{C_{\bf a}}{4\pi \epsilon_0}  
\sum_n  (\hat{r} \times \hat{b}_n)  D_{\mathbf{b}n} \,  v_{11}(r) 
\end{equation}
%
% t ca qb
%
\begin{equation}
\mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} = 
\frac{2C_{\bf a}}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \times \hat{b}_l) Q_{{\mathbf b}lm} (\hat{b}_m \cdot \hat{r}) v_{22}(r)
\end{equation}
%
%  t da cb
%
\begin{equation}
\mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} = 
-\frac{C_{\bf b}}{4\pi \epsilon_0}  
\sum_n  (\hat{r} \times \hat{a}_n)  D_{\mathbf{a}n} \,  v_{11}(r) 
\end{equation}%
%
%
%  ta da db
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m} 
(\hat{a}_m \times  \hat{b}_n)
D_{\mathbf{b}n} v_{21}(r) \\
% 2
&-\frac{1}{4\pi \epsilon_0} 
\sum_m (\hat{r} \times \hat{a}_m) D_{\mathbf {a}m} 
\sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
\end{split}
\end{equation}
%
%  tb da db
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} &= 
-\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m} 
(\hat{a}_m \times  \hat{b}_n)
D_{\mathbf{b}n} v_{21}(r) \\
% 2
&+\frac{1}{4\pi \epsilon_0} 
\sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m} 
\sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf {b}n} v_{22}(r) 
\end{split}
\end{equation}
%
% ta da qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0} \left(
-\text{Tr}Q_{\mathbf{b}}
\sum_n (\hat{r} \times \hat{a}_n) D_{\mathbf{a}n}
+2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \times \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
\right) v_{31}(r)\\
% 2
&-\frac{1}{4\pi \epsilon_0}
\sum_l (\hat{r} \times \hat{a}_l) D_{\mathbf{a}l}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{b}_n \cdot \hat{r}) v_{32}(r)
\end{split}
\end{equation}
%
% tb da qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0} \left(
-2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \cdot \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{r} \times \hat{b}_n)
-2\sum_{lmn}D_{\mathbf{a}l} 
(\hat{a}_l \times \hat{b}_m)
Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
\right) v_{31}(r) \\
% 2
&-\frac{2}{4\pi \epsilon_0}
\sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
\sum_{mn} (\hat{r} \cdot \hat{b}_m) 
Q_{{\mathbf b}mn}
(\hat{r}\times \hat{b}_n) v_{32}(r)
\end{split}
\end{equation}
%
% ta qa cb
%
\begin{equation}
\mathbf{\tau}_{{\bf a}Q_{\bf a}C_{\bf b}} = 
\frac{2C_{\bf a}}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{{\mathbf a}lm} (\hat{r} \times \hat{a}_m) v_{22}(r)
\end{equation}
%
% ta qa db
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0} \left(
2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{b}_l \cdot \hat{a}_m)
Q_{\mathbf{a}mn} (\hat{r} \times \hat{a}_n)
+2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{a}_l \times \hat{b}_m)
Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
\right) v_{31}(r) \\
% 2
&+\frac{2}{4\pi \epsilon_0}
\sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{r}\times \hat{a}_n) v_{32}(r)
\end{split}
\end{equation}
%
% tb qa db
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} &= 
\frac{1}{4\pi \epsilon_0} \left(
\text{Tr}Q_{\mathbf{a}}
\sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf{b}n}
+2\sum_{lmn}D_{\mathbf{b}l} 
(\hat{a}_l \times \hat{b}_m)
Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
\right) v_{31}(r)\\
% 2
&\frac{1}{4\pi \epsilon_0}
\sum_l (\hat{r} \times \hat{b}_l) D_{\mathbf{b}l}
\sum_{mn} (\hat{r} \cdot \hat{a}_m) 
Q_{{\mathbf a}mn}
(\hat{a}_n \cdot \hat{r}) v_{32}(r)
\end{split}
\end{equation}
%
% ta qa qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} &= 
-\frac{4}{4\pi \epsilon_0} 
\sum_{lmnp} (\hat{a}_l \times \hat{b}_p)
Q_{\mathbf{a}lm}  Q_{\mathbf{b}np} 
(\hat{a}_m \cdot \hat{b}_n) v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[
2\text{Tr}Q_{\mathbf{b}}
\sum_{lm} (\hat{r} \cdot \hat{a}_l )
Q_{{\mathbf a}lm}
(\hat{r} \times \hat{a}_m)
+4 \sum_{lmnp}
(\hat{r} \times \hat{a}_l ) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{b}_n)
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r}) \\
% 3
&-4 \sum_{lmnp}
(\hat{r} \cdot \hat{a}_l ) 
Q_{{\mathbf a}lm}
(\hat{a}_m \times \hat{b}_n)
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r})
\Bigr] v_{42}(r) \\
% 4
&+ \frac{2}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \times \hat{a}_l) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{r})
\sum_{np}  (\hat{r} \cdot \hat{b}_n) 
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r})  v_{43}(r)\\
\end{split}
\end{equation}
%
% tb qa qb
%
\begin{equation}
\begin{split}
% 1
\mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} &= 
\frac{4}{4\pi \epsilon_0} 
\sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
Q_{\mathbf{a}lm}  Q_{\mathbf{b}np} 
(\hat{a}_m \times \hat{b}_n) v_{41}(r) \\
% 2
&+ \frac{1}{4\pi \epsilon_0} 
\Bigl[
2\text{Tr}Q_{\mathbf{a}}
\sum_{lm} (\hat{r} \cdot \hat{b}_l )
Q_{{\mathbf b}lm}
(\hat{r} \times \hat{b}_m)
+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{r} \times \hat{b}_p) \\
% 3
&+4 \sum_{lmnp}
(\hat{r} \cdot \hat{a}_l ) 
Q_{{\mathbf a}lm}
(\hat{a}_m \times \hat{b}_n)
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r})
\Bigr] v_{42}(r)  \\
% 4
&+ \frac{2}{4\pi \epsilon_0}
\sum_{lm} (\hat{r} \cdot \hat{a}_l) 
Q_{{\mathbf a}lm}
(\hat{a}_m \cdot \hat{r})
\sum_{np}  (\hat{r} \times \hat{b}_n) 
Q_{{\mathbf b}np}
(\hat{b}_p \cdot \hat{r}) v_{43}(r). \\
\end{split}
\end{equation}
%
\begin{table*}
\caption{\label{tab:tableFORCE2}Radial functions used in the force equations.}
\begin{ruledtabular}
\begin{tabular}{ccc}
Generic&Method 1&Method 2
\\ \hline
%
%
%
$w_a(r)$&
$g_0(r)$&
$g(r)-g(r_c)$ \\
%
%
$w_b(r)$ &
$\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ &
$h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\
%
$w_c(r)$ &
$\frac{g_1(r)}{r} $ &
$\frac{v_{11}(r)}{r}$ \\
%
%
$w_d(r)$&
$\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ &
$\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)
-\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $\\
%
$w_e(r)$ &
$\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ &
$\frac{v_{22}(r)}{r}$ \\
% 
%
$w_f(r)$&
$\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ &
$\left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) - $ \\
&&$\left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)-\frac{2v_{22}(r)}{r}$\\
%
$w_g(r)$& $ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$& 
$\frac{v_{31}(r)}{r}$\\
%
$w_h(r)$ &
$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
$\left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) - $\\
&&$\left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\
&&$-\frac{v_{31}(r)}{r}$\\
% 2
$w_i(r)$ &
$\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
$\frac{v_{32}(r)}{r}$ \\
% 
$w_j(r)$ &
$\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right)  $ &
$\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right) $  \\
&&$\left(\frac{-6g(_cr)}{r_c^3} +\frac{6h(r_c)}{r_c^2} -\frac{3s(r_c)}{r_c} +t(r_c) \right) -\frac{3v_{32}}{r}$ \\
%
$w_k(r)$ &
$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
$\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}  \right)$  \\
&&$\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2}  \right)$ \\
%
$w_l(r)$ &
$\left(-\frac{15g_4(r)}{r^4} +\frac{15h_4(r)}{r^3} -\frac{6s_4(r)}{r^2} +\frac{t_4(r)}{r} \right)$ &
$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\
&&$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)
-\frac{2v_{42}(r)}{r}$ \\
%
$w_m(r)$ &
$\left(\frac{105g_4(r)}{r^4} - \frac{105h_4(r)}{r^3} + \frac{45s_4(r)}{r^2} - \frac{10t_4(r)}{r} +u_4(r) \right)$ &
$\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$ \\
&&$\left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3} 
+\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $ \\ 
&&$-\frac{4v_{43}(r)}{r}$ \\
%
$w_n(r)$ &
$\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
$\frac{v_{42}(r)}{r}$ \\
%
$w_o(r)$ &
$\left(-\frac{15g_4(r)}{r^4} +\frac{15h_4(r)}{r^3} -\frac{6s_4(r)}{r^2} +\frac{t_4(r)}{r} \right)$ &
$\frac{v_{43}(r)}{r}$ \\
%
\end{tabular}
\end{ruledtabular}
\end{table*}
\end{document}
%
% ****** End of file multipole.tex ******
Revision:	3906
Committed:	Tue Jul 9 21:45:05 2013 UTC (12 years, 4 months ago) by gezelter
Content type:	application/x-tex
Original Path:	*trunk/multipole/multipole3.tex*
File size:	63023 byte(s)
Log Message:	first import