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

Comparing trunk/multipole/multipole1.tex (file contents):
Revision 3984 by gezelter, Sat Dec 28 18:41:19 2013 UTC vs.
Revision 3988 by gezelter, Thu Jan 2 15:46:58 2014 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
41  
42   \begin{document}
43  
44 < \preprint{AIP/123-QED}
44 > %\preprint{AIP/123-QED}
45  
46 < \title[Taylor-shifted and Gradient-shifted electrostatics for multipoles]
46 < {Real space alternatives to the Ewald
46 > \title{Real space alternatives to the Ewald
47   Sum. I. Taylor-shifted and Gradient-shifted electrostatics for multipoles}
48  
49   \author{Madan Lamichhane}
# Line 76 | Line 76 | of Notre Dame, Notre Dame, IN 46556}
76    results for ordered arrays of multipoles.
77   \end{abstract}
78  
79 < \pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
79 > %\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
80                               % Classification Scheme.
81 < \keywords{Suggested keywords}%Use showkeys class option if keyword
81 > %\keywords{Suggested keywords}%Use showkeys class option if keyword
82                                %display desired
83   \maketitle
84  
# Line 98 | Line 98 | similar methods where the force vanishes at $R_\textrm
98   The computational efficiency and the accuracy of the DSF method are
99   surprisingly good, particularly for systems with uniform charge
100   density. Additionally, dielectric constants obtained using DSF and
101 < similar methods where the force vanishes at $R_\textrm{c}$ are
101 > similar methods where the force vanishes at $r_{c}$ are
102   essentially quantitative.\cite{Izvekov:2008wo} The DSF and other
103   related methods have now been widely investigated,\cite{Hansen:2012uq}
104 < and DSF is now used routinely in simulations of ionic
105 < liquids,\cite{doi:10.1021/la400226g,McCann:2013fk} flow in carbon
106 < nanotubes,\cite{kannam:094701} gas sorption in metal-organic framework
107 < materials,\cite{Forrest:2012ly} thermal conductivity of methane
108 < hydrates,\cite{English:2008kx} condensation coefficients of
109 < water,\cite{Louden:2013ve} and in tribology at solid-liquid-solid
110 < interfaces.\cite{Tokumasu:2013zr} DSF electrostatics provides a
111 < compromise between the computational speed of real-space cutoffs and
112 < the accuracy of fully-periodic Ewald treatments.
104 > and DSF is now used routinely in a diverse set of chemical
105 > environments.\cite{doi:10.1021/la400226g,McCann:2013fk,kannam:094701,Forrest:2012ly,English:2008kx,Louden:2013ve,Tokumasu:2013zr}
106 > DSF electrostatics provides a compromise between the computational
107 > speed of real-space cutoffs and the accuracy of fully-periodic Ewald
108 > treatments.
109  
114 \subsection{Coarse Graining using Point Multipoles}
110   One common feature of many coarse-graining approaches, which treat
111   entire molecular subsystems as a single rigid body, is simplification
112   of the electrostatic interactions between these bodies so that fewer
113   site-site interactions are required to compute configurational
114 < energies. Notably, the force matching approaches of Voth and coworkers
115 < are an exciting development in their ability to represent realistic
116 < (and {\it reactive}) chemical systems for very large length scales and
122 < long times.  This approach utilized a coarse-graining in interaction
123 < space (CGIS) which fits an effective force for the coarse grained
124 < system using a variational force-matching method to a fine-grained
125 < simulation.\cite{Izvekov:2008wo}
114 > energies. To do this, the interactions between coarse-grained sites
115 > are typically taken to be point
116 > multipoles.\cite{Golubkov06,ISI:000276097500009,ISI:000298664400012}
117  
118 < The coarse-graining approaches of Ren \& coworkers,\cite{Golubkov06}
119 < and Essex \&
120 < coworkers,\cite{ISI:000276097500009,ISI:000298664400012}
121 < both utilize Gay-Berne
122 < ellipsoids~\cite{Berne72,Gay81,Luckhurst90,Cleaver96,Berardi98,Ravichandran:1999fk,Berardi99,Pasterny00}
123 < to model dispersive interactions and point multipoles to model
124 < electrostatics for entire molecules or functional groups.
125 <
135 < Ichiye and coworkers have recently introduced a number of very fast
136 < water models based on a ``sticky'' multipole model which are
137 < qualitatively better at reproducing the behavior of real water than
138 < the more common point-charge models (SPC/E, TIPnP).  The point charge
139 < models are also substantially more computationally demanding than the
140 < sticky multipole
141 < approach.\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} The
142 < SSDQO model requires the use of an approximate multipole expansion
143 < (AME) as the exact multipole expansion is quite expensive
144 < (particularly when handled via the Ewald sum).\cite{Ichiye:2006qy}
145 < Another particularly important use of point multipoles (and multipole
146 < polarizability) is in the very high-quality AMOEBA water model and
118 > Water, in particular, has been modeled recently with point multipoles
119 > up to octupolar
120 > order.\cite{Chowdhuri:2006lr,Te:2010rt,Te:2010ys,Te:2010vn} For
121 > maximum efficiency, these models require the use of an approximate
122 > multipole expansion as the exact multipole expansion can become quite
123 > expensive (particularly when handled via the Ewald
124 > sum).\cite{Ichiye:2006qy} Point multipoles and multipole
125 > polarizability have also been utilized in the AMOEBA water model and
126   related force fields.\cite{Ponder:2010fk,schnieders:124114,Ren:2011uq}
127  
128   Higher-order multipoles present a peculiar issue for molecular
# Line 154 | Line 133 | kernel to point multipoles.  We have developed two dis
133   a different orientation can cause energy discontinuities.
134  
135   This paper outlines an extension of the original DSF electrostatic
136 < kernel to point multipoles.  We have developed two distinct real-space
136 > kernel to point multipoles.  We describe two distinct real-space
137   interaction models for higher-order multipoles based on two truncated
138   Taylor expansions that are carried out at the cutoff radius.  We are
139   calling these models {\bf Taylor-shifted} and {\bf Gradient-shifted}
140   electrostatics.  Because of differences in the initial assumptions,
141 < the two methods yield related, but different expressions for energies,
142 < forces, and torques.  
141 > the two methods yield related, but somewhat different expressions for
142 > energies, forces, and torques.
143  
144   In this paper we outline the new methodology and give functional forms
145   for the energies, forces, and torques up to quadrupole-quadrupole
146   order.  We also compare the new methods to analytic energy constants
147 < for periodic arrays of point multipoles.  In the following paper, we
148 < provide extensive numerical comparisons to Ewald-based electrostatics
149 < in common simulation enviornments.  
147 > for periodic arrays of point multipoles. In the following paper, we
148 > provide numerical comparisons to Ewald-based electrostatics in common
149 > simulation enviornments.
150  
151   \section{Methodology}
152 + An efficient real-space electrostatic method involves the use of a
153 + pair-wise functional form,
154 + \begin{equation}
155 + V = \sum_i \sum_{j>i} V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) +
156 + \sum_i V_i^\mathrm{correction}
157 + \end{equation}
158 + that is short-ranged and easily truncated at a cutoff radius,
159 + \begin{equation}
160 +  V_\mathrm{pair}(r_{ij}, \Omega_i, \Omega_j) = \left\{
161 + \begin{array}{ll}
162 + V_\mathrm{approx} (r_{ij}, \Omega_i, \Omega_j) & \quad r \le r_c \\
163 + 0 & \quad r > r_c ,
164 + \end{array}
165 + \right.
166 + \end{equation}
167 + along with an easily computed correction term ($\sum_i
168 + V_i^\mathrm{correction}$) which has linear-scaling with the number of
169 + particles.  Here $\Omega_i$ and $\Omega_j$ represent orientational
170 + coordinates of the two sites.  The computational efficiency, energy
171 + conservation, and even some physical properties of a simulation can
172 + depend dramatically on how the $V_\mathrm{approx}$ function behaves at
173 + the cutoff radius. The goal of any approximation method should be to
174 + mimic the real behavior of the electrostatic interactions as closely
175 + as possible without sacrificing the near-linear scaling of a cutoff
176 + method.
177  
178   \subsection{Self-neutralization, damping, and force-shifting}
179   The DSF and Wolf methods operate by neutralizing the total charge
180   contained within the cutoff sphere surrounding each particle.  This is
181   accomplished by shifting the potential functions to generate image
182   charges on the surface of the cutoff sphere for each pair interaction
183 < computed within $R_\textrm{c}$. Damping using a complementary error
183 > computed within $r_c$. Damping using a complementary error
184   function is applied to the potential to accelerate convergence. The
185   potential for the DSF method, where $\alpha$ is the adjustable damping
186   parameter, is given by
187   \begin{equation*}
188 < V_\mathrm{DSF}(r) = C_a C_b \Biggr{[}
188 > V_\mathrm{DSF}(r) = C_i C_j \Biggr{[}
189   \frac{\mathrm{erfc}\left(\alpha r_{ij}\right)}{r_{ij}}
190 < - \frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}} + \left(\frac{\mathrm{erfc}\left(\alpha R_\mathrm{c}\right)}{R_\mathrm{c}^2}
190 > - \frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c} + \left(\frac{\mathrm{erfc}\left(\alpha r_c\right)}{r_c^2}
191   + \frac{2\alpha}{\pi^{1/2}}
192 < \frac{\exp\left(-\alpha^2R_\mathrm{c}^2\right)}{R_\mathrm{c}}
193 < \right)\left(r_{ij}-R_\mathrm{c}\right)\ \Biggr{]}
192 > \frac{\exp\left(-\alpha^2r_c^2\right)}{r_c}
193 > \right)\left(r_{ij}-r_c\right)\ \Biggr{]}
194   \label{eq:DSFPot}
195   \end{equation*}
196 + Note that in this potential and in all electrostatic quantities that
197 + follow, the standard $1/4 \pi \epsilon_{0}$ has been omitted for
198 + clarity.
199  
200   To insure net charge neutrality within each cutoff sphere, an
201   additional ``self'' term is added to the potential. This term is
# Line 199 | Line 206 | damping.\cite{deLeeuw80,Wolf99}
206   over the surface of the cutoff sphere.  A portion of the self term is
207   identical to the self term in the Ewald summation, and comes from the
208   utilization of the complimentary error function for electrostatic
209 < damping.\cite{deLeeuw80,Wolf99}
209 > damping.\cite{deLeeuw80,Wolf99} There have also been recent efforts to
210 > extend the Wolf self-neutralization method to zero out the dipole and
211 > higher order multipoles contained within the cutoff
212 > sphere.\cite{Fukuda:2011jk,Fukuda:2012yu,Fukuda:2013qv}
213  
214 < There have been recent efforts to extend the Wolf self-neutralization
215 < method to zero out the dipole and higher order multipoles contained
216 < within the cutoff
217 < sphere.\cite{Fukuda:2011jk,Fukuda:2012yu,Fukuda:2013qv}  
214 > In this work, we extend the idea of self-neutralization for the point
215 > multipoles by insuring net charge-neutrality and net-zero moments
216 > within each cutoff sphere.  In Figure \ref{fig:shiftedMultipoles}, the
217 > central dipolar site $\mathbf{D}_i$ is interacting with point dipole
218 > $\mathbf{D}_j$ and point quadrupole, $\mathbf{Q}_k$.  The
219 > self-neutralization scheme for point multipoles involves projecting
220 > opposing multipoles for sites $j$ and $k$ on the surface of the cutoff
221 > sphere.  There are also significant modifications made to make the
222 > forces and torques go smoothly to zero at the cutoff distance.
223  
209 In this work, we will be extending the idea of self-neutralization for
210 the point multipoles in a similar way.  In Figure
211 \ref{fig:shiftedMultipoles}, the central dipolar site $\mathbf{D}_i$
212 is interacting with point dipole $\mathbf{D}_j$ and point quadrupole,
213 $\mathbf{Q}_k$.  The self-neutralization scheme for point multipoles
214 involves projecting opposing multipoles for sites $j$ and $k$ on the
215 surface of the cutoff sphere.
216
224   \begin{figure}
225   \includegraphics[width=3in]{SM}
226   \caption{Reversed multipoles are projected onto the surface of the
# Line 222 | Line 229 | As in the point-charge approach, there is a contributi
229   \label{fig:shiftedMultipoles}
230   \end{figure}
231  
232 < As in the point-charge approach, there is a contribution from
233 < self-neutralization of site $i$.  The self term for multipoles is
232 > As in the point-charge approach, there is an additional contribution
233 > from self-neutralization of site $i$.  The self term for multipoles is
234   described in section \ref{sec:selfTerm}.
235  
236   \subsection{The multipole expansion}
# Line 237 | Line 244 | V_a(\mathbf r) = \frac{1}{4\pi\epsilon_0}
244   a$.  Then the electrostatic potential of object $\bf a$ at
245   $\mathbf{r}$ is given by
246   \begin{equation}
247 < V_a(\mathbf r) = \frac{1}{4\pi\epsilon_0}
247 > V_a(\mathbf r) =
248   \sum_{k \, \text{in \bf a}} \frac{q_k}{\lvert \mathbf{r} - \mathbf{r}_k \rvert}.
249   \end{equation}
250   The Taylor expansion in $r$ can be written using an implied summation
# Line 256 | Line 263 | V_{\bf a}(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\hat{M
263   can then be used to express the electrostatic potential on $\bf a$ in
264   terms of multipole operators,
265   \begin{equation}
266 < V_{\bf a}(\mathbf{r}) = \frac{1}{4\pi\epsilon_0}\hat{M}_{\bf a} \frac{1}{r}
266 > V_{\bf a}(\mathbf{r}) =\hat{M}_{\bf a} \frac{1}{r}
267   \end{equation}
268   where
269   \begin{equation}
# Line 281 | Line 288 | U_{\bf{ab}}(r)
288   $\bf a$ to $\bf b$ ($\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $), the interaction energy is given by
289   \begin{equation}
290   U_{\bf{ab}}(r)
291 < = \frac{1}{4\pi \epsilon_0} \hat{M}_a \sum_{j \, \text{in \bf b}} \frac {q_j}{\vert \bf{r}+\bf{r}_j \vert} .
291 > = \hat{M}_a \sum_{j \, \text{in \bf b}} \frac {q_j}{\vert \bf{r}+\bf{r}_j \vert} .
292   \end{equation}
293   This can also be expanded as a Taylor series in $r$.  Using a notation
294   similar to before to define the multipoles on object {\bf b},
# Line 292 | Line 299 | U_{\bf{ab}}(r)=\frac{\hat{M}_{\bf a} \hat{M}_{\bf b}}{
299   \end{equation}
300   we arrive at the multipole expression for the total interaction energy.
301   \begin{equation}
302 < U_{\bf{ab}}(r)=\frac{\hat{M}_{\bf a} \hat{M}_{\bf b}}{4\pi \epsilon_0} \frac{1}{r}  \label{kernel}.
302 > U_{\bf{ab}}(r)=\hat{M}_{\bf a} \hat{M}_{\bf b} \frac{1}{r}  \label{kernel}.
303   \end{equation}
304   This form has the benefit of separating out the energies of
305   interaction into contributions from the charge, dipole, and quadrupole
306 < of $\bf a$ interacting with the same multipoles $\bf b$.
306 > of $\bf a$ interacting with the same multipoles on $\bf b$.
307  
308   \subsection{Damped Coulomb interactions}
309   In the standard multipole expansion, one typically uses the bare
# Line 311 | Line 318 | either $1/r$ or $B_0(r)$, and all of the techniques ca
318   \int_{\alpha r}^{\infty} \text{e}^{-s^2} ds .
319   \end{equation}
320   We develop equations below using the function $f(r)$ to represent
321 < either $1/r$ or $B_0(r)$, and all of the techniques can be applied
322 < either to bare or damped Coulomb kernels as long as derivatives of
323 < these functions are known.  Smith's convenient functions $B_l(r)$ are
324 < summarized in Appendix A.
321 > either $1/r$ or $B_0(r)$, and all of the techniques can be applied to
322 > bare or damped Coulomb kernels (or any other function) as long as
323 > derivatives of these functions are known.  Smith's convenient
324 > functions $B_l(r)$ are summarized in Appendix A.
325  
319
326   The main goal of this work is to smoothly cut off the interaction
327   energy as well as forces and torques as $r\rightarrow r_c$.  To
328   describe how this goal may be met, we use two examples, charge-charge
329 < and charge-dipole, using the bare Coulomb kernel $f(r)=1/r$ to explain
330 < the idea.
329 > and charge-dipole, using the bare Coulomb kernel, $f(r)=1/r$, to
330 > explain the idea.
331  
332   \subsection{Shifted-force methods}
333   In the shifted-force approximation, the interaction energy for two
334   charges $C_{\bf a}$ and $C_{\bf b}$ separated by a distance $r$ is
335   written:
336   \begin{equation}
337 < U_{C_{\bf a}C_{\bf b}}(r)=\frac{1}{4\pi \epsilon_0} C_{\bf a} C_{\bf b}
337 > U_{C_{\bf a}C_{\bf b}}(r)= C_{\bf a} C_{\bf b}
338   \left({ \frac{1}{r} - \frac{1}{r_c} + (r - r_c) \frac{1}{r_c^2}  }
339   \right) .
340   \end{equation}
# Line 345 | Line 351 | There are a number of ways to generalize this derivati
351   \right) = \left(- \frac{1}{r^2} + \frac{1}{r_c^2}
352   \right) .
353   \end{equation}
354 < There are a number of ways to generalize this derivative shift for
354 > which clearly vanishes as the $r$ approaches the cutoff radius. There
355 > are a number of ways to generalize this derivative shift for
356   higher-order multipoles.  Below, we present two methods, one based on
357   higher-order Taylor series for $r$ near $r_c$, and the other based on
358   linear shift of the kernel gradients at the cutoff itself.
# Line 356 | Line 363 | derivative. We therefore require shifted energy expres
363   derivatives required for each force term to vanish at the cutoff.  For
364   example, the quadrupole-quadrupole interaction energy requires four
365   derivatives of the kernel, and the force requires one additional
366 < derivative. We therefore require shifted energy expressions that
367 < include enough terms so that all energies, forces, and torques are
368 < zero as $r \rightarrow r_c$.  In each case, we will subtract off a
369 < function $f_n^{\text{shift}}(r)$ from the kernel $f(r)=1/r$.  The
370 < subscript $n$ indicates the number of derivatives to be taken when
371 < deriving a given multipole energy.  We choose a function with
372 < guaranteed smooth derivatives --- a truncated Taylor series of the
373 < function $f(r)$, e.g.,
366 > derivative. For quadrupole-quadrupole interactions, we therefore
367 > require shifted energy expressions that include up to $(r-r_c)^5$ so
368 > that all energies, forces, and torques are zero as $r \rightarrow
369 > r_c$. In each case, we subtract off a function $f_n^{\text{shift}}(r)$
370 > from the kernel $f(r)=1/r$.  The subscript $n$ indicates the number of
371 > derivatives to be taken when deriving a given multipole energy.  We
372 > choose a function with guaranteed smooth derivatives -- a truncated
373 > Taylor series of the function $f(r)$, e.g.,
374   %
375   \begin{equation}
376   f_n^{\text{shift}}(r)=\sum_{m=0}^{n+1} \frac {(r-r_c)^m}{m!} f^{(m)}(r_c)  .
# Line 381 | Line 388 | U_{C_{\bf a}D_{\bf b}}(r)=
388   %
389   \begin{equation}
390   U_{C_{\bf a}D_{\bf b}}(r)=
391 < \frac{C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}  \frac {\partial f_1(r) }{\partial r_\alpha}
392 < =\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}  
391 > C_{\bf a} D_{{\bf b}\alpha}  \frac {\partial f_1(r) }{\partial r_\alpha}
392 > = C_{\bf a} D_{{\bf b}\alpha}
393   \frac {r_\alpha}{r} \frac {\partial f_1(r)}{\partial r} .
394   \end{equation}
395   %
396   The force that dipole  $\bf b$ exerts on charge $\bf a$ is
397   %
398   \begin{equation}
399 < F_{C_{\bf a}D_{\bf b}\beta} =\frac{ C_{\bf a} D_{{\bf b}\alpha}}{4\pi \epsilon_0}
399 > F_{C_{\bf a}D_{\bf b}\beta} = C_{\bf a} D_{{\bf b}\alpha}
400   \left[ \frac{\delta_{\alpha\beta}}{r} \frac {\partial}{\partial r} +
401   \frac{r_\alpha r_\beta}{r^2}
402   \left( -\frac{1}{r} \frac {\partial} {\partial r}
# Line 400 | Line 407 | F_{C_{\bf a}D_{\bf b}\beta} =
407   %
408   \begin{equation}
409   F_{C_{\bf a}D_{\bf b}\beta} =
410 < \frac{C_{\bf a} D_{{\bf b}\beta} }{4\pi \epsilon_0r}
410 > \frac{C_{\bf a} D_{{\bf b}\beta}}{r}
411   \left[  -\frac{1}{r^2}+\frac{1}{r_c^2}-\frac{2(r-r_c)}{r_c^3} \right]
412 < +\frac{C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta }{4\pi \epsilon_0}
412 > +C_{\bf a} D_{{\bf b}\alpha}r_\alpha r_\beta
413   \left[ \frac{3}{r^5}-\frac{3}{r^3r_c^2} \right] .
414   \end{equation}
415   %
# Line 411 | Line 418 | U=\frac{1}{4\pi \epsilon_0} (\text{prefactor}) (\text{
418   In general, we can write
419   %
420   \begin{equation}
421 < U=\frac{1}{4\pi \epsilon_0} (\text{prefactor}) (\text{derivatives}) f_n(r)
421 > U= (\text{prefactor}) (\text{derivatives}) f_n(r)
422   \label{generic}
423   \end{equation}
424   %
425 < where $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for charge-quadrupole
426 < and dipole-dipole, $n=3$ for dipole-quadrupole, and $n=4$ for quadrupole-quadrupole.
427 < An example is the case of quadrupole-quadrupole for which the $\text{prefactor}$ is
428 < $Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$ and the derivatives are
429 < $\partial^4/\partial r_\alpha \partial r_\beta \partial r_\gamma \partial r_\delta$, with
430 < implied summation combining the space indices.
425 > with $n=0$ for charge-charge, $n=1$ for charge-dipole, $n=2$ for
426 > charge-quadrupole and dipole-dipole, $n=3$ for dipole-quadrupole, and
427 > $n=4$ for quadrupole-quadrupole.  For example, in
428 > quadrupole-quadrupole interactions for which the $\text{prefactor}$ is
429 > $Q_{{\bf a}\alpha\beta}Q_{{\bf b}\gamma\delta}$, the derivatives are
430 > $\partial^4/\partial r_\alpha \partial r_\beta \partial
431 > r_\gamma \partial r_\delta$, with implied summation combining the
432 > space indices.
433  
434   In the formulas presented in the tables below, the placeholder
435   function $f(r)$ is used to represent the electrostatic kernel (either
436   damped or undamped).  The main functions that go into the force and
437 < torque terms, $f_n(r), g_n(r), h_n(r), s_n(r), \mathrm{~and~} t_n(r)$
438 < are successive derivatives of the shifted electrostatic kernel of the
439 < same index $n$.  The algebra required to evaluate energies, forces and
440 < torques is somewhat tedious and are summarized in Appendices A and B.
437 > torque terms, $g_n(r), h_n(r), s_n(r), \mathrm{~and~} t_n(r)$ are
438 > successive derivatives of the shifted electrostatic kernel, $f_n(r)$
439 > of the same index $n$.  The algebra required to evaluate energies,
440 > forces and torques is somewhat tedious, so only the final forms are
441 > presented in tables \ref{tab:tableenergy} and \ref{tab:tableFORCE}.
442  
443   \subsection{Gradient-shifted force (GSF) electrostatics}
444 < Note the method used in the previous subsection to smoothly shift the
445 < force to zero is a truncated Taylor Series in the radius $r$.  The
446 < second method maintains only the linear $(r-r_c)$ term and has a
447 < similar interaction energy for all multipole orders:
438 < \begin{equation}
439 < U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big \lvert  _{r_c} .
440 < \end{equation}
441 < No higher order terms $(r-r_c)^n$ appear.  The primary difference
442 < between the TSF and GSF methods is the stage at which the Taylor
443 < Series is applied; in the Taylor-shifted approach, it is applied to
444 < the kernel itself.  In the Gradient-shifted approach, it is applied to
445 < individual radial interactions terms in the multipole expansion.
446 < Terms from this method thus have the general form:
444 > The second, and conceptually simpler approach to force-shifting
445 > maintains only the linear $(r-r_c)$ term in the truncated Taylor
446 > expansion, and has a similar interaction energy for all multipole
447 > orders:
448   \begin{equation}
449 < U=\frac{1}{4\pi \epsilon_0}\sum  (\text{angular factor}) (\text{radial factor}).
449 > U^{\text{shift}}(r)=U(r)-U(r_c)-(r-r_c)\hat{r}\cdot \nabla U(r) \Big
450 > \lvert  _{r_c} .
451   \label{generic2}
452   \end{equation}
453 + Here the gradient for force shifting is evaluated for an image
454 + multipole projected onto the surface of the cutoff sphere (see fig
455 + \ref{fig:shiftedMultipoles}). No higher order terms $(r-r_c)^n$
456 + appear.  The primary difference between the TSF and GSF methods is the
457 + stage at which the Taylor Series is applied; in the Taylor-shifted
458 + approach, it is applied to the kernel itself.  In the Gradient-shifted
459 + approach, it is applied to individual radial interactions terms in the
460 + multipole expansion.  Energies from this method thus have the general
461 + form:
462 + \begin{equation}
463 + U= \sum  (\text{angular factor}) (\text{radial factor}).
464 + \label{generic3}
465 + \end{equation}
466  
467 < Results for both methods can be summarized using the form of
468 < Eq.~(\ref{generic2}) and are listed in Table I below.
467 > Functional forms for both methods (TSF and GSF) can both be summarized
468 > using the form of Eq.~(\ref{generic3}).  The basic forms for the
469 > energy, force, and torque expressions are tabulated for both shifting
470 > approaches below -- for each separate orientational contribution, only
471 > the radial factors differ between the two methods.
472  
473   \subsection{\label{sec:level2}Body and space axes}
474  
475 + [XXX Do we need this section in the main paper? or should it go in the
476 + extra materials?]
477 +
478   So far, all energies and forces have been written in terms of fixed
479 < space coordinates $x$, $y$, $z$.  Interaction energies are computed
480 < from the generic formulas Eq.~(\ref{generic}) and ~(\ref{generic2})
481 < which combine prefactors with radial functions.  Because objects $\bf
479 > space coordinates.  Interaction energies are computed from the generic
480 > formulas Eq.~(\ref{generic}) and ~(\ref{generic2}) which combine
481 > orientational prefactors with radial functions.  Because objects $\bf
482   a$ and $\bf b$ both translate and rotate during a molecular dynamics
483   (MD) simulation, it is desirable to contract all $r$-dependent terms
484   with dipole and quadrupole moments expressed in terms of their body
485 < axes.  To do so, we follow the methodology of Allen and
486 < Germano,\cite{Allen:2006fk} which was itself based on an earlier paper
487 < by Price {\em et al.}\cite{Price:1984fk}
485 > axes.  To do so, we have followed the methodology of Allen and
486 > Germano,\cite{Allen:2006fk} which was itself based on earlier work by
487 > Price {\em et al.}\cite{Price:1984fk}
488  
489   We denote body axes for objects $\bf a$ and $\bf b$ by unit vectors
490   $\hat{a}_m$ and $\hat{b}_m$, respectively, with the index $m=(123)$
# Line 475 | Line 496 | Allen and Germano define matrices $\hat{\mathbf {a}}$
496   \hat{a}_m= a_{mx}\hat{x} + a_{my}\hat{y} + a_{mz}\hat{z}  \\
497   \hat{b}_m= b_{mx}\hat{x} + b_{my}\hat{y} + b_{mz}\hat{z}  .
498   \end{eqnarray}
499 < Allen and Germano define matrices $\hat{\mathbf {a}}$
500 < and  $\hat{\mathbf {b}}$ using these unit vectors:
499 > Rotation matrices $\hat{\mathbf {a}}$ and $\hat{\mathbf {b}}$ can be
500 > expressed using these unit vectors:
501   \begin{eqnarray}
502   \hat{\mathbf {a}} =
503   \begin{pmatrix}
# Line 498 | Line 519 | b_{1x}\quad  b_{1y} \quad b_{1z} \\
519   \end{pmatrix}
520   =
521   \begin{pmatrix}
522 < b_{1x}\quad  b_{1y} \quad b_{1z} \\
522 > b_{1x} \quad  b_{1y} \quad b_{1z} \\
523   b_{2x} \quad b_{2y} \quad b_{2z} \\
524   b_{3x} \quad b_{3y} \quad b_{3z}
525   \end{pmatrix}  .
526   \end{eqnarray}
527   %
528 < These matrices convert from space-fixed $(xyz)$ to object-fixed $(123)$ coordinates.
529 < All contractions of prefactors with derivatives of functions can be written in terms of these matrices.
530 < It proves to be equally convenient to just write any contraction in terms of unit vectors
531 < $\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$.  
532 < We have found it useful to write angular-dependent terms in three different fashions,
533 < illustrated by the following
534 < three examples from the interaction-energy expressions:
528 > These matrices convert from space-fixed $(xyz)$ to body-fixed $(123)$
529 > coordinates.  All contractions of prefactors with derivatives of
530 > functions can be written in terms of these matrices. It proves to be
531 > equally convenient to just write any contraction in terms of unit
532 > vectors $\hat{r}$, $\hat{a}_m$, and $\hat{b}_n$. In the torque
533 > expressions, it is useful to have the angular-dependent terms
534 > available in three different fashions, e.g. for the dipole-dipole
535 > contraction:
536   %
537 < \begin{eqnarray}
537 > \begin{equation}
538   \mathbf{D}_{\mathbf {a}} \cdot \mathbf{D}_{\mathbf{b}}
539 < =D_{\bf {a}\alpha} D_{\bf {b}\alpha}=
540 < \sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}} \\
541 < 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}
539 > = D_{\bf {a}\alpha} D_{\bf {b}\alpha} =
540 > \sum_{mn} {D_{\mathbf{a}m} \hat{a}_m \cdot \hat{b}_n D_{\mathbf{b}n}}
541 > \end{equation}
542   %
543 < [Dan, perhaps there are better examples to show here.]
543 > The first two forms are written using space coordinates.  The first
544 > form is standard in the chemistry literature, while the second is
545 > expressed using implied summation notation.  The third form shows
546 > explicit sums over body indices and the dot products now indicate
547 > contractions using space indices.
548  
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.  
549  
538
550   \subsection{The Self-Interaction \label{sec:selfTerm}}
551  
552 < The Wolf summation~\cite{Wolf99} and the later damped shifted force
553 < (DSF) extension~\cite{Fennell06} included self-interactions that are
554 < handled separately from the pairwise interactions between sites. The
555 < self-term is normally calculated via a single loop over all sites in
556 < the system, and is relatively cheap to evaluate. The self-interaction
557 < has contributions from two sources:
558 < \begin{itemize}
559 < \item The neutralization procedure within the cutoff radius requires a
560 <  contribution from a charge opposite in sign, but equal in magnitude,
561 <  to the central charge, which has been spread out over the surface of
562 <  the cutoff sphere.  For a system of undamped charges, the total
563 <  self-term is
552 > In addition to cutoff-sphere neutralization, the Wolf
553 > summation~\cite{Wolf99} and the damped shifted force (DSF)
554 > extension~\cite{Fennell:2006zl} also included self-interactions that
555 > are handled separately from the pairwise interactions between
556 > sites. The self-term is normally calculated via a single loop over all
557 > sites in the system, and is relatively cheap to evaluate. The
558 > self-interaction has contributions from two sources.
559 >
560 > First, the neutralization procedure within the cutoff radius requires
561 > a contribution from a charge opposite in sign, but equal in magnitude,
562 > to the central charge, which has been spread out over the surface of
563 > the cutoff sphere.  For a system of undamped charges, the total
564 > self-term is
565   \begin{equation}
566   V_\textrm{self} = - \frac{1}{r_c} \sum_{{\bf a}=1}^N C_{\bf a}^{2}
567   \label{eq:selfTerm}
568   \end{equation}
569 < Note that in this potential and in all electrostatic quantities that
570 < follow, the standard $4 \pi \epsilon_{0}$ has been omitted for
571 < clarity.
572 < \item Charge damping with the complementary error function is a
573 <  partial analogy to the Ewald procedure which splits the interaction
574 <  into real and reciprocal space sums.  The real space sum is retained
575 <  in the Wolf and DSF methods.  The reciprocal space sum is first
576 <  minimized by folding the largest contribution (the self-interaction)
577 <  into the self-interaction from charge neutralization of the damped
578 <  potential.  The remainder of the reciprocal space portion is then
579 <  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
569 >
570 > Second, charge damping with the complementary error function is a
571 > partial analogy to the Ewald procedure which splits the interaction
572 > into real and reciprocal space sums.  The real space sum is retained
573 > in the Wolf and DSF methods.  The reciprocal space sum is first
574 > minimized by folding the largest contribution (the self-interaction)
575 > into the self-interaction from charge neutralization of the damped
576 > potential.  The remainder of the reciprocal space portion is then
577 > discarded (as this contributes the largest computational cost and
578 > complexity to the Ewald sum).  For a system containing only damped
579 > charges, the complete self-interaction can be written as
580   \begin{equation}
581   V_\textrm{self} = - \left(\frac{\textrm{erfc}(\alpha r_c)}{r_c} +
582    \frac{\alpha}{\sqrt{\pi}} \right) \sum_{{\bf a}=1}^N
583        C_{\bf a}^{2}.
584   \label{eq:dampSelfTerm}
585   \end{equation}
576 \end{itemize}
586  
587   The extension of DSF electrostatics to point multipoles requires
588   treatment of {\it both} the self-neutralization and reciprocal
# Line 616 | Line 625 | terms, $f(r), g(r), h(r), s(r), \mathrm{~and~} t(r)$ a
625   multipole orders.  Symmetry prevents the charge-dipole and
626   dipole-quadrupole interactions from contributing to the
627   self-interaction.  The functions that go into the self-neutralization
628 < terms, $f(r), g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive
629 < derivatives of the electrostatic kernel (either the undamped $1/r$ or
630 < the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that are
631 < evaluated at the cutoff distance.  For undamped interactions, $f(r_c)
632 < = 1/r_c$, $g(r_c) = -1/r_c^{2}$, and so on.  For damped interactions,
633 < $f(r_c) = B_0(r_c)$, $g(r_c) = B_0'(r_c)$, and so on.  Appendix XX
634 < contains recursion relations that allow rapid evaluation of these
635 < derivatives.
628 > terms, $g(r), h(r), s(r), \mathrm{~and~} t(r)$ are successive
629 > derivatives of the electrostatic kernel, $f(r)$ (either the undamped
630 > $1/r$ or the damped $B_0(r)=\mathrm{erfc}(\alpha r)/r$ function) that
631 > have been evaluated at the cutoff distance.  For undamped
632 > interactions, $f(r_c) = 1/r_c$, $g(r_c) = -1/r_c^{2}$, and so on.  For
633 > damped interactions, $f(r_c) = B_0(r_c)$, $g(r_c) = B_0'(r_c)$, and so
634 > on.  Appendix \ref{SmithFunc} contains recursion relations that allow
635 > rapid evaluation of these derivatives.
636  
637 < \section{Energies, forces, and torques}
638 < \subsection{Interaction energies}
637 > \section{Interaction energies, forces, and torques}
638 > The main result of this paper is a set of expressions for the
639 > energies, forces and torques (up to quadrupole-quadrupole order) that
640 > work for both the Taylor-shifted and Gradient-shifted approximations.
641 > These expressions were derived using a set of generic radial
642 > functions.  Without using the shifting approximations mentioned above,
643 > some of these radial functions would be identical, and the expressions
644 > coalesce into the familiar forms for unmodified multipole-multipole
645 > interactions.  Table \ref{tab:tableenergy} maps between the generic
646 > functions and the radial functions derived for both the Taylor-shifted
647 > and Gradient-shifted methods.  The energy equations are written in
648 > terms of lab-frame representations of the dipoles, quadrupoles, and
649 > the unit vector connecting the two objects,
650  
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
651   % Energy in space coordinate form ----------------------------------------------------------------------------------------------
652   %
653   %
# Line 641 | Line 655 | U_{C_{\bf a}C_{\bf b}}(r)=&
655   %
656   \begin{align}
657   U_{C_{\bf a}C_{\bf b}}(r)=&
658 < \frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  v_{01}(r)  \label{uchch}
658 > C_{\bf a} C_{\bf b}  v_{01}(r)  \label{uchch}
659   \\
660   %
661   % u ca db
662   %
663   U_{C_{\bf a}D_{\bf b}}(r)=&
664 < \frac{C_{\bf a}}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)  v_{11}(r)  
664 > C_{\bf a} \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)  v_{11}(r)  
665   \label{uchdip}
666   \\
667   %
668   % u ca qb
669   %
670 < U_{C_{\bf a}Q_{\bf b}}(r)=&
671 < \frac{C_{\bf a }}{4\pi \epsilon_0} \Bigl[ \text{Tr}Q_{\bf b} v_{21}(r)
672 < \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{22}(r) \Bigr]
670 > U_{C_{\bf a}Q_{\bf b}}(r)=& C_{\bf a } \Bigl[ \text{Tr}Q_{\bf b}
671 > v_{21}(r) + \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot
672 >  \hat{r} \right) v_{22}(r) \Bigr]
673   \label{uchquad}
674   \\
675   %
# Line 669 | Line 683 | -\frac{1}{4\pi \epsilon_0} \Bigr[ \left( \mathbf{D}_{\
683   % u da db
684   %
685   U_{D_{\bf a}D_{\bf b}}(r)=&
686 < -\frac{1}{4\pi \epsilon_0} \Bigr[ \left( \mathbf{D}_{\mathbf {a}} \cdot
686 > -\Bigr[ \left( \mathbf{D}_{\mathbf {a}} \cdot
687   \mathbf{D}_{\mathbf{b}} \right)  v_{21}(r)
688   +\left( \mathbf{D}_{\mathbf {a}} \cdot \hat{r} \right)
689   \left( \mathbf{D}_{\mathbf {b}} \cdot \hat{r} \right)  
# Line 682 | Line 696 | -\frac{1}{4\pi \epsilon_0} \Bigl[
696   \begin{split}
697   % 1
698   U_{D_{\bf a}Q_{\bf b}}(r) =&
699 < -\frac{1}{4\pi \epsilon_0} \Bigl[
699 > -\Bigl[
700   \text{Tr}\mathbf{Q}_{\mathbf{b}}
701   \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
702   +2 ( \mathbf{D}_{\mathbf{a}} \cdot
703   \mathbf{Q}_{\mathbf{b}} \cdot \hat{r} ) \Bigr] v_{31}(r) \\
704   % 2
705 < &-\frac{1}{4\pi \epsilon_0} \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
705 > &- \left( \mathbf{D}_{\mathbf{a}} \cdot \hat{r} \right)
706   \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) v_{32}(r)
707   \label{udipquad}
708   \end{split}
# Line 725 | Line 739 | U_{Q_{\bf a}Q_{\bf b}}(r)=&
739   \begin{split}
740   %1
741   U_{Q_{\bf a}Q_{\bf b}}(r)=&
742 < \frac{1}{4\pi \epsilon_0} \Bigl[
742 > \Bigl[
743   \text{Tr} \mathbf{Q}_{\mathbf{a}} \text{Tr} \mathbf{Q}_{\mathbf{b}}
744   +2 \text{Tr} \left(
745   \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} \right) \Bigr] v_{41}(r)
746   \\
747   % 2
748 < &+ \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
748 > &+\Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
749   \left( \hat{r} \cdot
750   \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right)
751   +\text{Tr}\mathbf{Q}_{\mathbf{b}}
# Line 741 | Line 755 | +\text{Tr}\mathbf{Q}_{\mathbf{b}}
755   \Bigr] v_{42}(r)
756   \\
757   % 4
758 < &+ \frac{1}{4\pi \epsilon_0}
758 > &+
759   \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right)
760   \left( \hat{r} \cdot \mathbf{Q}_{{\mathbf b}}  \cdot \hat{r} \right) v_{43}(r).
761   \label{uquadquad}
762   \end{split}
763   \end{align}
764 <
764 > %
765   Note that the energies of multipoles on site $\mathbf{b}$ interacting
766   with those on site $\mathbf{a}$ can be obtained by swapping indices
767   along with the sign of the intersite vector, $\hat{r}$.
# Line 757 | Line 771 | along with the sign of the intersite vector, $\hat{r}$
771   % TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
772   %
773  
774 < \begin{table*}
775 < \caption{\label{tab:tableenergy}Radial functions used in the energy and torque equations.  Functions
776 < used in this table are defined in Appendices B and C.}
777 < \begin{ruledtabular}
778 < \begin{tabular}{|l|c|l|l}
779 < Generic&Coulomb&Taylor-Shifted&Gradient-Shifted
774 > \begin{sidewaystable}
775 >  \caption{\label{tab:tableenergy}Radial functions used in the energy
776 >    and torque equations.  The $f, g, h, s, t, \mathrm{and} u$
777 >    functions used in this table are defined in Appendices B and C.}
778 > \begin{tabular}{|c|c|l|l|} \hline
779 > Generic&Bare Coulomb&Taylor-Shifted&Gradient-Shifted
780   \\ \hline
781   %
782   %
# Line 794 | Line 808 | -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right) $ \\
808   $\frac{3}{r^3}  $ &
809   $\left(-\frac{g_2(r)}{r} + h_2(r) \right)$ &
810   $\left(-\frac{g(r)}{r}+h(r) \right)
811 < -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right) $ \\
812 < &&&$  -(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$
811 > -\left(-\frac{g(r_c)}{r_c}+h(r_c) \right)$  \\
812 > &&& $ ~~~-(r-r_c) \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)$
813   \\
814   %
815   %
# Line 806 | Line 820 | -\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right
820   $\left(-\frac{g_3(r)}{r^2} + \frac{h_3(r)}{r} \right)$ &
821   $\left( -\frac{g(r)}{r^2}+\frac{h(r)}{r} \right)
822   -\left(-\frac{g(r_c)}{r_c^2}+\frac{h(r_c)}{r_c} \right) $\\
823 < &&&$  -(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)$
823 > &&&$ ~~~ -(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)$
824   \\
825   %
826   $v_{32}(r)$ &
# Line 814 | Line 828 | - \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} +
828   $\left( \frac{3g_3(r)}{r^2} - \frac{3h_3(r)}{r} + s_3(r) \right)$ &
829   $\left( \frac{3g(r)}{r^2} - \frac{3h(r)}{r} + s(r) \right)
830   - \left( \frac{3g(r_c)}{r_c^2} - \frac{3h(r_c)}{r_c} + s(r_c) \right)$ \\
831 < &&&$  -(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)$
831 > &&&$ ~~~ -(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)$
832   \\
833   %
834   %
# Line 825 | Line 839 | - \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2}
839   $\left(-\frac{g_4(r)}{r^3} +\frac{h_4(r)}{r^2} \right) $  &
840   $\left( -\frac{g(r)}{r^3} + \frac{h(r)}{r^2} \right)
841   - \left( -\frac{g(r_c)}{r_c^3} + \frac{h(r_c)}{r_c^2} \right)$ \\
842 < &&&$  -(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)$
842 > &&&$ ~~~ -(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)$
843   \\
844   % 2
845   $v_{42}(r)$ &
# Line 833 | Line 847 | -\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+
847   $\left( \frac{3g_4(r)}{r^3} - \frac{3h_4(r)}{r^2}+\frac{s_4(r)}{r} \right)$ &
848   $\left( \frac{3g(r)}{r^3} - \frac{3h(r)}{r^2}+\frac{s(r)}{r} \right)
849   -\left( \frac{3g(r_c)}{r_c^3} - \frac{3h(r_c)}{r_c^2}+\frac{s(r_c)}{r_c} \right)$ \\
850 < &&&$  -(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}
850 > &&&$ ~~~ -(r-r_c) \left(- \frac{9g(r_c)}{r_c^4}+\frac{9h(r_c)}{r_c^3}
851   -\frac{4s(r_c)}{r_c^2} + \frac{t(r_c)}{r_c}\right)$
852   \\
853   % 3
# Line 841 | Line 855 | $\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s
855   $ \frac{105}{r^5}  $ &
856   $\left(-\frac{15g_4(r)}{r^3}+\frac{15h_4(r)}{r^2}-\frac{6s_4(r)}{r} + t_4(r)\right) $ &
857   $\left(-\frac{15g(r)}{r^3}+\frac{15h(r)}{r^2}-\frac{6s(r)}{r} + t(r)\right)$ \\
858 < &&&$ -\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)$ \\
859 < &&&$ -(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}
860 < -\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\
858 > &&&$~~~ -\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)$ \\
859 > &&&$~~~ -(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}
860 > -\frac{6t(r_c)}{r_c}+u(r_c) \right)$ \\ \hline
861   \end{tabular}
862 < \end{ruledtabular}
849 < \end{table*}
862 > \end{sidewaystable}
863   %
864   %
865   % FORCE  TABLE of radial functions  ----------------------------------------------------------------------------------------------------------------
866   %
867  
868 < \begin{table}
868 > \begin{sidewaystable}
869   \caption{\label{tab:tableFORCE}Radial functions used in the force equations.}
870 < \begin{ruledtabular}
871 < \begin{tabular}{cc}
859 < Generic&Method 1 or Method 2
870 > \begin{tabular}{|c|c|l|l|} \hline
871 > Function&Definition&Taylor-Shifted&Gradient-Shifted
872   \\ \hline
873   %
874   %
875   %
876   $w_a(r)$&
877 < $\frac{d v_{01}}{dr}$ \\
877 > $\frac{d v_{01}}{dr}$&
878 > $g_0(r)$&
879 > $g(r)-g(r_c)$ \\
880   %
881   %
882   $w_b(r)$ &
883 < $\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $ \\
883 > $\frac{d v_{11}}{dr} - \frac{v_{11}(r)}{r} $&
884 > $\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ &
885 > $h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\
886   %
887   $w_c(r)$ &
888 < $\frac{v_{11}(r)}{r}$ \\
888 > $\frac{v_{11}(r)}{r}$ &
889 > $\frac{g_1(r)}{r} $ &
890 > $\frac{v_{11}(r)}{r}$\\
891   %
892   %
893   $w_d(r)$&
894 < $\frac{d v_{21}}{dr}$ \\
894 > $\frac{d v_{21}}{dr}$&
895 > $\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ &
896 > $\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)
897 > -\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $ \\
898   %
899   $w_e(r)$ &
900 + $\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ &
901 + $\frac{v_{22}(r)}{r}$ &
902   $\frac{v_{22}(r)}{r}$ \\
903   %
904   %
905   $w_f(r)$&
906 < $\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$\\
906 > $\frac{d v_{22}}{dr} - \frac{2v_{22}(r)}{r}$&
907 > $\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ &
908 >  $ \left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) $  \\
909 >  &&& $ ~~~- \left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c)
910 >     \right)-\frac{2v_{22}(r)}{r}$\\
911   %
912   $w_g(r)$&
913 + $\frac{v_{31}(r)}{r}$&
914 + $ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$&
915   $\frac{v_{31}(r)}{r}$\\
916   %
917   $w_h(r)$ &
918 < $\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$\\
918 > $\frac{d v_{31}}{dr} -\frac{v_{31}(r)}{r}$&
919 > $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
920 > $ \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) $ \\
921 > &&& $ ~~~ -\frac{v_{31}(r)}{r}$ \\
922   % 2
923   $w_i(r)$ &
924 < $\frac{v_{32}(r)}{r}$ \\
924 > $\frac{v_{32}(r)}{r}$ &
925 > $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
926 > $\frac{v_{32}(r)}{r}$\\
927   %
928   $w_j(r)$ &
929 < $\frac{d v_{32}}{dr}  - \frac{3v_{32}}{r}$ \\
929 > $\frac{d v_{32}}{dr}  - \frac{3v_{32}}{r}$&
930 > $\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right)  $ &
931 > $\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right)$ \\
932 > &&& $~~~-\left(\frac{-6g(_cr)}{r_c^3} +\frac{6h(r_c)}{r_c^2}
933 >  -\frac{3s(r_c)}{r_c} +t(r_c) \right) -\frac{3v_{32}}{r}$ \\
934   %
935   $w_k(r)$ &
936 < $\frac{d v_{41}}{dr} $  \\
936 > $\frac{d v_{41}}{dr} $ &
937 > $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
938 > $\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}  \right)
939 > -\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2}  \right)$ \\
940   %
941   $w_l(r)$ &
942 < $\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ \\
942 > $\frac{d v_{42}}{dr} -\frac{2v_{42}(r)}{r}$ &
943 > $\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)$ &
944 > $\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\
945 > &&& $~~~ -\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)
946 > -\frac{2v_{42}(r)}{r}$\\
947   %
948   $w_m(r)$ &
949 < $\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$ \\
949 > $\frac{d v_{43}}{dr} -\frac{4v_{43}(r)}{r}$&
950 > $\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)$ &
951 > $\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$\\
952 > &&& $~~~- \left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}
953 > +\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $\\
954 > &&& $~~~-\frac{4v_{43}(r)}{r}$  \\
955   %
956   $w_n(r)$ &
957 < $\frac{v_{42}(r)}{r}$ \\
957 > $\frac{v_{42}(r)}{r}$ &
958 > $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
959 > $\frac{v_{42}(r)}{r}$\\
960   %
961   $w_o(r)$ &
962 < $\frac{v_{43}(r)}{r}$ \\
962 > $\frac{v_{43}(r)}{r}$&
963 > $\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)$ &
964 > $\frac{v_{43}(r)}{r}$  \\ \hline
965   %
966  
967   \end{tabular}
968 < \end{ruledtabular}
915 < \end{table}
968 > \end{sidewaystable}
969   %
970   %
971   %
972  
973   \subsection{Forces}
974 <
975 < The force $\mathbf{F}_{\bf a}$ on $\bf{a}$ due to $\bf{b}$ is the negative of
976 < the force  $\mathbf{F}_{\bf b}$ on $\bf{b}$ due to $\bf{a}$.  For a simple charge-charge
977 < interaction, these forces will point along the $\pm \hat{r}$ directions, where
978 < $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_a $.  Thus
974 > The force on object $\bf{a}$, $\mathbf{F}_{\bf a}$, due to object
975 > $\bf{b}$ is the negative of the force on $\bf{b}$ due to $\bf{a}$. For
976 > a simple charge-charge interaction, these forces will point along the
977 > $\pm \hat{r}$ directions, where $\mathbf{r}=\mathbf{r}_b -
978 > \mathbf{r}_a $.  Thus
979   %
980   \begin{equation}
981   F_{\bf a \alpha} = \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}}{\partial r}
# Line 930 | Line 983 | The concept of obtaining a force from an energy by tak
983   = - \hat{r}_\alpha \frac{\partial U_{C_{\bf a}C_{\bf b}}} {\partial r}  .
984   \end{equation}
985   %
986 < The concept of obtaining a force from an energy by taking a gradient is the same for
987 < higher-order multipole interactions, the trick is to make sure that all
988 < $r$-dependent derivatives are considered.
989 < As is pointed out by Allen and Germano,\cite{Allen:2006fk} this is straightforward if the
990 < interaction energies are written recognizing explicit
991 < $\hat{r}$ and body axes ($\hat{a}_m$, $\hat{b}_n$) dependences:
986 > Obtaining the force from the interaction energy expressions is the
987 > same for higher-order multipole interactions -- the trick is to make
988 > sure that all $r$-dependent derivatives are considered.  This is
989 > straighforward if the interaction energies are written explicitly in
990 > terms of $\hat{r}$ and the body axes ($\hat{a}_m$,
991 > $\hat{b}_n$) :
992   %
993   \begin{equation}
994   U(r,\{\hat{a}_m \cdot \hat{r} \},
995 < \{\hat{b}_n\cdot \hat{r} \}
995 > \{\hat{b}_n\cdot \hat{r} \},
996   \{\hat{a}_m \cdot \hat{b}_n \}) .
997   \label{ugeneral}
998   \end{equation}
999   %
1000 < Then,
1000 > Allen and Germano,\cite{Allen:2006fk} showed that if the energy is
1001 > written in this form, the forces come out relatively cleanly,
1002   %
1003   \begin{equation}
1004   \mathbf{F}_{\bf a}=-\mathbf{F}_{\bf b} =  \frac{\partial U}{\partial \mathbf{r}}
# Line 957 | Line 1011 | Note our definition of $\mathbf{r}=\mathbf{r}_b - \mat
1011   \right] \label{forceequation}.
1012   \end{equation}
1013   %
1014 < Note our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $ is opposite
1015 < that of Allen and Germano.\cite{Allen:2006fk}  In simplifying the algebra, we also use:
1014 > Note that our definition of $\mathbf{r}=\mathbf{r}_b - \mathbf{r}_b $
1015 > is opposite in sign to that of Allen and Germano.\cite{Allen:2006fk}
1016 > In simplifying the algebra, we have also used:
1017   %
1018 < \begin{eqnarray}
1018 > \begin{align}
1019   \frac { \partial (\hat{a}_m \cdot \hat{r})}{\partial \mathbf{r}}
1020 < = \frac{1}{r} \left(  \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r}
1020 > =& \frac{1}{r} \left(  \hat{a}_m - (\hat{a}_m \cdot \hat{r})\hat{r}
1021   \right) \\
1022   \frac { \partial (\hat{b}_m \cdot \hat{r})}{\partial \mathbf{r}}
1023 < = \frac{1}{r} \left(  \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r}
1023 > =& \frac{1}{r} \left(  \hat{b}_m - (\hat{b}_m \cdot \hat{r})\hat{r}
1024   \right) .
1025 < \end{eqnarray}
1025 > \end{align}
1026   %
1027 < We list below the force equations written in terms  of space coordinates.  The
1028 < radial functions used in the two methods are listed in Table II.
1027 > We list below the force equations written in terms of lab-frame
1028 > coordinates.  The radial functions used in the two methods are listed
1029 > in Table \ref{tab:tableFORCE}
1030   %
1031 < %SPACE COORDINATES FORCE EQUTIONS
1031 > %SPACE COORDINATES FORCE EQUATIONS
1032   %
1033   % **************************************************************************
1034   % f ca cb
1035   %
1036 < \begin{equation}
1037 < \mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =
1038 < \frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r}
983 < \end{equation}
1036 > \begin{align}
1037 > \mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =&
1038 > C_{\bf a} C_{\bf b}  w_a(r) \hat{r} \\
1039   %
1040   %
1041   %
1042 < \begin{equation}
1043 < \mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =
989 < \frac{C_{\bf a}}{4\pi \epsilon_0} \Bigl[
1042 > \mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =&
1043 > C_{\bf a} \Bigl[
1044   \left( \hat{r} \cdot \mathbf{D}_{\mathbf{b}} \right)
1045   w_b(r) \hat{r}  
1046 < + \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr]
993 < \end{equation}
1046 > + \mathbf{D}_{\mathbf{b}} w_c(r) \Bigr] \\
1047   %
1048   %
1049   %
1050 < \begin{equation}
1051 < \mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =
999 < \frac{C_{\bf a }}{4\pi \epsilon_0} \Bigr[
1050 > \mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =&
1051 > C_{\bf a } \Bigr[
1052   \text{Tr}\mathbf{Q}_{\bf b} w_d(r) \hat{r}
1053   + 2  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} w_e(r)
1054 < + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
1055 < \end{equation}
1054 > + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}
1055 > \right) w_f(r) \hat{r} \Bigr] \\
1056   %
1057   %
1058   %
1059 < \begin{equation}
1060 < \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1061 < -\frac{C_{\bf{b}}}{4\pi \epsilon_0} \Bigl[
1062 < \left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) w_b(r) \hat{r}
1063 < + \mathbf{D}_{\mathbf{a}} w_c(r) \Bigr]
1064 < \end{equation}
1059 > % \begin{equation}
1060 > % \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1061 > % -C_{\bf{b}} \Bigl[
1062 > % \left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right) w_b(r) \hat{r}
1063 > % + \mathbf{D}_{\mathbf{a}} w_c(r) \Bigr]
1064 > % \end{equation}
1065   %
1066   %
1067   %
1068 < \begin{equation}
1069 < \mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} =
1018 < \frac{1}{4\pi \epsilon_0} \Bigl[
1068 > \begin{split}
1069 > \mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} =&
1070   - \mathbf{D}_{\mathbf {a}} \cdot  \mathbf{D}_{\mathbf{b}} w_d(r) \hat{r}
1071   + \left( \mathbf{D}_{\mathbf {a}}
1072   \left( \mathbf{D}_{\mathbf{b}} \cdot \hat{r} \right)
1073 < + \mathbf{D}_{\mathbf {b}} \left( \mathbf{D}_{\mathbf{a}}  \cdot \hat{r} \right) \right) w_e(r)
1073 > + \mathbf{D}_{\mathbf {b}} \left( \mathbf{D}_{\mathbf{a}}  \cdot \hat{r} \right) \right) w_e(r)\\
1074   % 2
1075 < - \left( \hat{r} \cdot \mathbf{D}_{\mathbf {a}} \right)
1076 < \left( \hat{r} \cdot \mathbf{D}_{\mathbf {b}} \right) w_f(r) \hat{r} \Bigr]
1077 < \end{equation}
1075 > & - \left( \hat{r} \cdot \mathbf{D}_{\mathbf {a}} \right)
1076 > \left( \hat{r} \cdot \mathbf{D}_{\mathbf {b}} \right) w_f(r) \hat{r}
1077 > \end{split}\\
1078   %
1079   %
1080   %
1030 \begin{equation}
1081   \begin{split}
1082 < \mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1033 < & - \frac{1}{4\pi \epsilon_0} \Bigl[
1082 > \mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =& - \Bigl[
1083   \text{Tr}\mathbf{Q}_{\mathbf{b}} \mathbf{ D}_{\mathbf{a}}
1084   +2 \mathbf{D}_{\mathbf{a}} \cdot
1085   \mathbf{Q}_{\mathbf{b}} \Bigr] w_g(r)
1086 < - \frac{1}{4\pi \epsilon_0} \Bigl[
1086 > - \Bigl[
1087   \text{Tr}\mathbf{Q}_{\mathbf{b}}
1088   \left( \hat{r} \cdot  \mathbf{D}_{\mathbf{a}} \right)
1089   +2 ( \mathbf{D}_{\mathbf{a}} \cdot
1090   \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
1091   % 3
1092 < & - \frac{1}{4\pi \epsilon_0} \Bigl[\mathbf{ D}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
1092 > & - \Bigl[\mathbf{ D}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
1093   +2  (\hat{r} \cdot \mathbf{D}_{\mathbf{a}} ) (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} )  \Bigr]
1094   w_i(r)
1095   % 4
1096 < -\frac{1}{4\pi \epsilon_0}
1096 > -
1097   (\hat{r} \cdot \mathbf{D}_{\mathbf{a}} )
1098 < (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r}  
1050 < \end{split}
1051 < \end{equation}
1098 > (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) w_j(r) \hat{r} \end{split} \\
1099   %
1100   %
1101 < \begin{equation}
1102 < \mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} =
1103 < \frac{C_{\bf b }}{4\pi \epsilon_0} \Bigr[
1104 < \text{Tr}\mathbf{Q}_{\bf a} w_d(r) \hat{r}
1105 < + 2  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} w_e(r)
1106 < + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
1107 < \end{equation}
1101 > % \begin{equation}
1102 > % \mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} =
1103 > % \frac{C_{\bf b }}{4\pi \epsilon_0} \Bigr[
1104 > % \text{Tr}\mathbf{Q}_{\bf a} w_d(r) \hat{r}
1105 > % + 2  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} w_e(r)
1106 > %  + \left( \hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} \cdot \hat{r} \right) w_f(r) \hat{r} \Bigr]
1107 > % \end{equation}
1108 > % %
1109 > % \begin{equation}
1110 > % \begin{split}
1111 > % \mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1112 > % &\frac{1}{4\pi \epsilon_0} \Bigl[
1113 > % \text{Tr}\mathbf{Q}_{\mathbf{a}} \mathbf{D}_{\mathbf{b}}
1114 > % +2 \mathbf{D}_{\mathbf{b}} \cdot \mathbf{Q}_{\mathbf{a}}  \Bigr] w_g(r)
1115 > % % 2
1116 > % + \frac{1}{4\pi \epsilon_0} \Bigl[ \text{Tr}\mathbf{Q}_{\mathbf{a}}
1117 > % (\hat{r} \cdot  \mathbf{D}_{\mathbf{b}})
1118 > % +2 (\mathbf{D}_{\mathbf{b}} \cdot
1119 > % \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] w_h(r) \hat{r}  \\
1120 > % % 3
1121 > % &+ \frac{1}{4\pi \epsilon_0} \Bigl[ \mathbf{D}_{\mathbf{b}}
1122 > % (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1123 > % +2  (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1124 > % (\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}} ) \Bigr]   w_i(r)
1125 > % % 4
1126 > % +\frac{1}{4\pi \epsilon_0}
1127 > % (\hat{r} \cdot \mathbf{D}_{\mathbf{b}})
1128 > % (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}}  \cdot \hat{r}) w_j(r) \hat{r}
1129 > % \end{split}
1130 > % \end{equation}
1131   %
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}
1132   %
1133   %
1086 %
1087 \begin{equation}
1134   \begin{split}
1135 < \mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1136 < +\frac{1}{4\pi \epsilon_0} \Bigl[
1135 > \mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =&
1136 > \Bigl[
1137   \text{Tr}\mathbf{Q}_{\mathbf{a}} \text{Tr}\mathbf{Q}_{\mathbf{b}} \hat{r}
1138   + 2 \text{Tr} ( \mathbf{Q}_{\mathbf{a}} \cdot  \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_k(r) \hat{r} \\
1139   % 2
1140 < +\frac{1}{4\pi \epsilon_0} \Bigl[
1140 > &+ \Bigl[
1141   2\text{Tr}\mathbf{Q}_{\mathbf{b}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )  
1142   + 2\text{Tr}\mathbf{Q}_{\mathbf{a}}  (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} )
1143   % 3
1144   +4 (\mathbf{Q}_{\mathbf{a}}  \cdot  \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})  
1145   +  4(\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}}) \Bigr] w_n(r) \\
1146   % 4
1147 < + \frac{1}{4\pi \epsilon_0} \Bigl[
1147 > &+  \Bigl[
1148   \text{Tr}\mathbf{Q}_{\mathbf{a}} (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1149   + \text{Tr}\mathbf{Q}_{\mathbf{b}}
1150   (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})  
# Line 1106 | Line 1152 | + \frac{1}{4\pi \epsilon_0} \Bigl[
1152   +4 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot  
1153   \mathbf{Q}_{\mathbf{b}}   \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
1154   %
1155 < + \frac{1}{4\pi \epsilon_0} \Bigl[
1155 > &+ \Bigl[
1156   + 2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} )
1157   (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1158   %6
1159   +2 (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}} \cdot \hat{r})
1160   (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} ) \Bigr] w_o(r) \\
1161   %  7
1162 < + \frac{1}{4\pi \epsilon_0}
1162 > &+
1163   (\hat{r} \cdot \mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
1164 < (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) w_m(r) \hat{r}
1165 < \end{split}
1166 < \end{equation}
1164 > (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}} \cdot \hat{r}) w_m(r) \hat{r} \end{split}
1165 > \end{align}
1166 > Note that the forces for higher multipoles on site $\mathbf{a}$
1167 > interacting with those of lower order on site $\mathbf{b}$ can be
1168 > obtained by swapping indices in the expressions above.
1169 >
1170   %
1171 + % Torques SECTION -----------------------------------------------------------------------------------------
1172   %
1123 % TORQUES SECTION -----------------------------------------------------------------------------------------
1124 %
1173   \subsection{Torques}
1174 <
1175 < Following again Allen and Germano,\cite{Allen:2006fk} when energies are written in the form
1176 < of Eq.~({\ref{ugeneral}), then torques can be expressed as:
1174 > When energies are written in the form of Eq.~({\ref{ugeneral}), then
1175 >  torques can be found in a relatively straightforward
1176 >  manner,\cite{Allen:2006fk}
1177   %
1178   \begin{eqnarray}
1179   \mathbf{\tau}_{\bf a} =
# Line 1146 | Line 1194 | Here we list the torque equations written in terms of
1194   \end{eqnarray}
1195   %
1196   %
1197 < Here we list the torque equations written in terms of space coordinates.
1197 > The torques for both the Taylor-Shifted as well as Gradient-Shifted
1198 > methods are given in space-frame coordinates:
1199   %
1200   %
1201 + \begin{align}
1202 + \mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =&
1203 + C_{\bf a}  (\hat{r} \times  \mathbf{D}_{\mathbf{b}}) v_{11}(r) \\
1204   %
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}
1205   %
1206   %
1207 + \mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =&
1208 + 2C_{\bf a}
1209 + \hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{22}(r) \\
1210   %
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}
1211   %
1212   %
1213 + % \begin{equation}
1214 + % \mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =  
1215 + % -\frac{C_{\bf b}}{4\pi \epsilon_0}  
1216 + % (\hat{r} \times \mathbf{D}_{\mathbf{a}})  v_{11}(r)
1217 + % \end{equation}
1218   %
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}
1219   %
1220   %
1221 < %
1222 < \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)
1221 > \mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} =&
1222 > \mathbf{D}_{\mathbf {a}}  \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1223   % 2
1224 < -\frac{1}{4\pi \epsilon_0}
1224 > -
1225   (\hat{r} \times \mathbf{D}_{\mathbf {a}} )
1226 < (\hat{r} \cdot \mathbf{D}_{\mathbf {b}} )  v_{22}(r)
1183 < \end{equation}
1226 > (\hat{r} \cdot \mathbf{D}_{\mathbf {b}} )  v_{22}(r)\\
1227   %
1228   %
1229   %
1230 < \begin{equation}
1231 < \mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} =
1232 < -\frac{1}{4\pi \epsilon_0} \mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1233 < % 2
1234 < +\frac{1}{4\pi \epsilon_0}
1235 < (\hat{r} \cdot \mathbf{D}_{\mathbf {a}} )
1236 < (\hat{r} \times \mathbf{D}_{\mathbf {b}} ) v_{22}(r)
1237 < \end{equation}
1230 > % \begin{equation}
1231 > % \mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} =
1232 > % -\frac{1}{4\pi \epsilon_0} \mathbf{D}_{\mathbf {a}} \times \mathbf{D}_{\mathbf{b}} v_{21}(r)
1233 > % % 2
1234 > % +\frac{1}{4\pi \epsilon_0}
1235 > % (\hat{r} \cdot \mathbf{D}_{\mathbf {a}} )
1236 > % (\hat{r} \times \mathbf{D}_{\mathbf {b}} ) v_{22}(r)
1237 > % \end{equation}
1238   %
1239   %
1240   %
1241 < \begin{equation}
1242 < \mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1200 < \frac{1}{4\pi \epsilon_0} \Bigl[
1241 > \mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} =&
1242 > \Bigl[
1243   -\text{Tr}\mathbf{Q}_{\mathbf{b}}
1244   (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
1245   +2 \mathbf{D}_{\mathbf{a}}  \times
1246   (\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1247   \Bigr] v_{31}(r)
1248   % 3
1249 < -\frac{1}{4\pi \epsilon_0}
1250 < \ (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
1209 < (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)
1210 < \end{equation}
1249 > - (\hat{r} \times \mathbf{D}_{\mathbf{a}} )
1250 > (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{32}(r)\\
1251   %
1252   %
1253   %
1254 < \begin{equation}
1255 < \mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =
1216 < \frac{1}{4\pi \epsilon_0} \Bigl[
1254 > \mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} =&
1255 > \Bigl[
1256   +2 ( \mathbf{D}_{\mathbf{a}} \cdot \mathbf{Q}_{\mathbf{b}} ) \times
1257   \hat{r}
1258   -2 \mathbf{D}_{\mathbf{a}}  \times
1259   (\mathbf{Q}_{\mathbf{b}} \cdot \hat{r})
1260   \Bigr] v_{31}(r)
1261   % 2
1262 < +\frac{2}{4\pi \epsilon_0}
1262 > +
1263   (\hat{r} \cdot \mathbf{D}_{\mathbf{a}})
1264 < (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)
1226 < \end{equation}
1264 > (\hat{r} \cdot \mathbf{Q}_{\mathbf{b}}) \times \hat{r} v_{32}(r)\\
1265   %
1266   %
1267   %
1268 < \begin{equation}
1269 < \mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1270 < \frac{1}{4\pi \epsilon_0} \Bigl[
1271 < -2 (\mathbf{D}_{\mathbf{b}}  \cdot \mathbf{Q}_{\mathbf{a}} ) \times \hat{r}
1272 < +2 \mathbf{D}_{\mathbf{b}}  \times
1273 < (\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
1274 < \Bigr] v_{31}(r)
1275 < % 3
1276 < - \frac{2}{4\pi \epsilon_0}
1277 < (\hat{r} \cdot \mathbf{D}_{\mathbf{b}} )
1278 < (\hat{r} \cdot  \mathbf{Q}_{{\mathbf a}}) \times \hat{r} v_{32}(r)
1279 < \end{equation}
1268 > % \begin{equation}
1269 > % \mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1270 > % \frac{1}{4\pi \epsilon_0} \Bigl[
1271 > % -2 (\mathbf{D}_{\mathbf{b}}  \cdot \mathbf{Q}_{\mathbf{a}} ) \times \hat{r}
1272 > % +2 \mathbf{D}_{\mathbf{b}}  \times
1273 > % (\mathbf{Q}_{\mathbf{a}}  \cdot \hat{r})
1274 > % \Bigr] v_{31}(r)
1275 > % % 3
1276 > % - \frac{2}{4\pi \epsilon_0}
1277 > % (\hat{r} \cdot \mathbf{D}_{\mathbf{b}} )
1278 > % (\hat{r} \cdot  \mathbf
1279 > % {Q}_{{\mathbf a}}) \times \hat{r} v_{32}(r)
1280 > % \end{equation}
1281   %
1282   %
1283   %
1284 < \begin{equation}
1285 < \mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} =
1286 < \frac{1}{4\pi \epsilon_0} \Bigl[
1287 < \text{Tr}\mathbf{Q}_{\mathbf{a}}
1288 < (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1289 < +2 \mathbf{D}_{\mathbf{b}}  \times
1290 < ( \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r)
1291 < % 2
1292 < +\frac{1}{4\pi \epsilon_0}
1293 < (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1294 < (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) v_{32}(r)
1295 < \end{equation}
1284 > % \begin{equation}
1285 > % \mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} =
1286 > % \frac{1}{4\pi \epsilon_0} \Bigl[
1287 > % \text{Tr}\mathbf{Q}_{\mathbf{a}}
1288 > % (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1289 > % +2 \mathbf{D}_{\mathbf{b}}  \times
1290 > % ( \mathbf{Q}_{\mathbf{a}} \cdot \hat{r}) \Bigr] v_{31}(r)
1291 > % % 2
1292 > % +\frac{1}{4\pi \epsilon_0}
1293 > % (\hat{r} \times \mathbf{D}_{\mathbf{b}} )
1294 > % (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r}) v_{32}(r)
1295 > % \end{equation}
1296   %
1297   %
1298   %
1260 \begin{equation}
1299   \begin{split}
1300 < \mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1301 < &-\frac{4}{4\pi \epsilon_0}
1300 > \mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} =&
1301 > -4
1302   \mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}}
1303   v_{41}(r) \\
1304   % 2
1305 < &+ \frac{1}{4\pi \epsilon_0}
1305 > &+
1306   \Bigl[-2\text{Tr}\mathbf{Q}_{\mathbf{b}}
1307   (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times \hat{r}
1308   +4 \hat{r} \times
# Line 1273 | Line 1311 | -4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} )\times
1311   -4 (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} )\times
1312   ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r} ) \Bigr] v_{42}(r) \\
1313   % 4
1314 < &+ \frac{2}{4\pi \epsilon_0}
1314 > &+ 2
1315   \hat{r} \times ( \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1316 < (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
1279 < \end{split}
1280 < \end{equation}
1316 > (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r) \end{split}\\
1317   %
1318   %
1319   %
1284 \begin{equation}
1320   \begin{split}
1321   \mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} =  
1322 < &\frac{4}{4\pi \epsilon_0}
1322 > &4
1323   \mathbf{Q}_{{\mathbf a}} \times \mathbf{Q}_{{\mathbf b}} v_{41}(r) \\
1324   % 2
1325 < &+ \frac{1}{4\pi \epsilon_0} \Bigl[- 2\text{Tr}\mathbf{Q}_{\mathbf{a}}
1325 > &+  \Bigl[- 2\text{Tr}\mathbf{Q}_{\mathbf{a}}
1326   (\hat{r} \cdot \mathbf{Q}_{{\mathbf b}} ) \times \hat{r}
1327   -4  (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot
1328   \mathbf{Q}_{{\mathbf b}} ) \times
# Line 1296 | Line 1331 | +4 ( \hat{r} \cdot \mathbf{Q}_{{\mathbf a}} ) \times
1331   ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r})
1332   \Bigr] v_{42}(r) \\
1333   % 4
1334 < &+ \frac{2}{4\pi \epsilon_0}
1334 > &+2
1335   (\hat{r} \cdot \mathbf{Q}_{{\mathbf a}} \cdot \hat{r})
1336 < \hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)
1337 < \end{split}
1303 < \end{equation}
1336 > \hat{r} \times ( \mathbf{Q}_{{\mathbf b}} \cdot \hat{r}) v_{43}(r)\end{split}
1337 > \end{align}
1338   %
1339 + Here, we have defined the matrix cross product in an identical form
1340 + as in Ref. \onlinecite{Smith98}:
1341 + \begin{equation}
1342 + \left[\mathbf{A} \times \mathbf{B}\right]_\alpha = \sum_\beta
1343 + \left[\mathbf{A}_{\alpha+1,\beta} \mathbf{B}_{\alpha+2,\beta}
1344 +  -\mathbf{A}_{\alpha+2,\beta} \mathbf{B}_{\alpha+2,\beta}
1345 + \right]
1346 + \end{equation}
1347 + where $\alpha+1$ and $\alpha+2$ are regarded as cyclic
1348 + permuations of the matrix indices.
1349  
1350 + All of the radial functions required for torques are identical with
1351 + the radial functions previously computed for the interaction energies.
1352 + These are tabulated for both shifted force methods in table
1353 + \ref{tab:tableenergy}.  The torques for higher multipoles on site
1354 + $\mathbf{a}$ interacting with those of lower order on site
1355 + $\mathbf{b}$ can be obtained by swapping indices in the expressions
1356 + above.
1357 +
1358   \section{Comparison to known multipolar energies}
1359  
1360   To understand how these new real-space multipole methods behave in
# Line 1323 | Line 1375 | Tisza\cite{LT,LT2} who tabulated energy constants for
1375   who computed the energies of ordered dipole arrays of zero
1376   magnetization and obtained a number of these constants.\cite{Sauer}
1377   This theory was developed more completely by Luttinger and
1378 < Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays and
1379 < other periodic structures.  We have repeated the Luttinger \& Tisza
1380 < series summations to much higher order and obtained the following
1381 < energy constants (converged to one part in $10^9$):
1382 < \begin{table*}
1378 > Tisza\cite{LT,LT2} who tabulated energy constants for the Sauer arrays
1379 > and other periodic structures.  We have repeated the Luttinger \&
1380 > Tisza series summations to much higher order and obtained the energy
1381 > constants (converged to one part in $10^9$) in table \ref{tab:LT}.
1382 >
1383 > \begin{table*}[h]
1384   \centering{
1385    \caption{Luttinger \& Tisza arrays and their associated
1386      energy constants. Type "A" arrays have nearest neighbor strings of
# Line 1345 | Line 1398 | Array Type &  Lattice &   Dipole Direction &    Energy
1398     A     &      BCC      &      111         &      -1.770078733 \\
1399     A     &      FCC      &      001         &       2.166932835 \\
1400     A     &      FCC      &      011         &      -1.083466417 \\
1401 <
1402 <   *     &      BCC      &    minimum       &      -1.985920929 \\
1403 <
1404 <   B     &      SC        &     001         &      -2.676788684 \\
1405 <   B     &      BCC       &     001         &      -1.338394342 \\
1406 <   B     &      BCC       &     111         &     -1.770078733 \\
1354 <   B     &      FCC       &     001         &      -1.083466417 \\
1355 <   B     &      FCC       &     011         &      -1.807573634
1401 >   B     &      SC       &      001         &      -2.676788684 \\
1402 >   B     &      BCC      &      001         &      -1.338394342 \\
1403 >   B     &      BCC      &      111         &      -1.770078733 \\
1404 >   B     &      FCC      &      001         &      -1.083466417 \\
1405 >   B     &      FCC      &      011         &      -1.807573634 \\
1406 >  --     &      BCC      &    minimum       &      -1.985920929 \\
1407   \end{tabular}
1408   \end{ruledtabular}
1409   \end{table*}
1410  
1411   In addition to the A and B arrays, there is an additional minimum
1412   energy structure for the BCC lattice that was found by Luttinger \&
1413 < Tisza.  The total electrostatic energy for an array is given by:
1413 > Tisza.  The total electrostatic energy for any of the arrays is given
1414 > by:
1415   \begin{equation}
1416    E = C N^2 \mu^2
1417   \end{equation}
1418 < where $C$ is the energy constant given above, $N$ is the number of
1419 < dipoles per unit volume, and $\mu$ is the strength of the dipole.
1418 > where $C$ is the energy constant given in table \ref{tab:LT}, $N$ is
1419 > the number of dipoles per unit volume, and $\mu$ is the strength of
1420 > the dipole.
1421  
1422 < {\it Quadrupolar} analogues to the Madelung constants were first worked out by Nagai and Nakamura who
1423 < computed the energies of selected quadrupole arrays based on
1424 < extensions to the Luttinger and Tisza
1425 < approach.\cite{Nagai01081960,Nagai01091963}  We have compared the
1422 > To test the new electrostatic methods, we have constructed very large,
1423 > $N$ = 8,000~(sc), 16,000~(bcc), or 32,000~(fcc) arrays of dipoles in
1424 > the orientations described in table \ref{tab:LT}.  For the purposes of
1425 > testing the energy expressions and the self-neutralization schemes,
1426 > the primary quantity of interest is the analytic energy constant for
1427 > the perfect arrays.  Convergence to these constants are shown as a
1428 > function of both the cutoff radius, $r_c$, and the damping parameter,
1429 > $\alpha$ in Figs.  \ref{fig:energyConstVsCutoff} and XXX. We have
1430 > simultaneously tested a hard cutoff (where the kernel is simply
1431 > truncated at the cutoff radius), as well as a shifted potential (SP)
1432 > form which includes a potential-shifting and self-interaction term,
1433 > but does not shift the forces and torques smoothly at the cutoff
1434 > radius.
1435 >
1436 > \begin{figure}
1437 > \includegraphics[width=4.5in]{energyConstVsCutoff}
1438 > \caption{Convergence to the analytic energy constants as a function of
1439 >  cutoff radius (normalized by the lattice constant) for the different
1440 >  real-space methods. The two crystals shown here are the ``B'' array
1441 >  for bcc crystals with the dipoles along the 001 direction (upper),
1442 >  as well as the minimum energy bcc lattice (lower).  The analytic
1443 >  energy constants are shown as a grey dashed line.  The left panel
1444 >  shows results for the undamped kernel ($1/r$), while the damped
1445 >  error function kernel, $B_0(r)$ was used in the right panel. }
1446 > \label{fig:energyConstVsCutoff}
1447 > \end{figure}
1448 >
1449 > The Hard cutoff exhibits oscillations around the analytic energy
1450 > constants, and converges to incorrect energies when the complementary
1451 > error function damping kernel is used.  The shifted potential (SP) and
1452 > gradient-shifted force (GSF) approximations converge to the correct
1453 > energy smoothly by $r_c / 6 a$ even for the undamped case.  This
1454 > indicates that the correction provided by the self term is required
1455 > for obtaining accurate energies.  The Taylor-shifted force (TSF)
1456 > approximation appears to perturb the potential too much inside the
1457 > cutoff region to provide accurate measures of the energy constants.
1458 >
1459 >
1460 > {\it Quadrupolar} analogues to the Madelung constants were first
1461 > worked out by Nagai and Nakamura who computed the energies of selected
1462 > quadrupole arrays based on extensions to the Luttinger and Tisza
1463 > approach.\cite{Nagai01081960,Nagai01091963} We have compared the
1464   energy constants for the lowest energy configurations for linear
1465   quadrupoles shown in table \ref{tab:NNQ}
1466  
# Line 1394 | Line 1485 | where $Q$ is the quadrupole moment.
1485   \end{equation}
1486   where $Q$ is the quadrupole moment.
1487  
1488 + \section{Conclusion}
1489 + We have presented two efficient real-space methods for computing the
1490 + interactions between point multipoles.  These methods have the benefit
1491 + of smoothly truncating the energies, forces, and torques at the cutoff
1492 + radius, making them attractive for both molecular dynamics (MD) and
1493 + Monte Carlo (MC) simulations.  We find that the Gradient-Shifted Force
1494 + (GSF) and the Shifted-Potential (SP) methods converge rapidly to the
1495 + correct lattice energies for ordered dipolar and quadrupolar arrays,
1496 + while the Taylor-Shifted Force (TSF) is too severe an approximation to
1497 + provide accurate convergence to lattice energies.  
1498  
1499 + In most cases, GSF can obtain nearly quantitative agreement with the
1500 + lattice energy constants with reasonably small cutoff radii.  The only
1501 + exception we have observed is for crystals which exhibit a bulk
1502 + macroscopic dipole moment (e.g. Luttinger \& Tisza's $Z_1$ lattice).
1503 + In this particular case, the multipole neutralization scheme can
1504 + interfere with the correct computation of the energies.  We note that
1505 + the energies for these arrangements are typically much larger than for
1506 + crystals with net-zero moments, so this is not expected to be an issue
1507 + in most simulations.
1508  
1509 + In large systems, these new methods can be made to scale approximately
1510 + linearly with system size, and detailed comparisons with the Ewald sum
1511 + for a wide range of chemical environments follows in the second paper.
1512  
1400
1401
1402
1403
1513   \begin{acknowledgments}
1514 < Support for this project was provided by the National Science
1515 < Foundation under grant CHE-0848243. Computational time was provided by
1516 < the Center for Research Computing (CRC) at the University of Notre
1517 < Dame.
1514 >  JDG acknowledges helpful discussions with Christopher
1515 >  Fennell. Support for this project was provided by the National
1516 >  Science Foundation under grant CHE-0848243. Computational time was
1517 >  provided by the Center for Research Computing (CRC) at the
1518 >  University of Notre Dame.
1519   \end{acknowledgments}
1520  
1521   \newpage
1522   \appendix
1523  
1524   \section{Smith's $B_l(r)$ functions for damped-charge distributions}
1525 <
1525 > \label{SmithFunc}
1526   The following summarizes Smith's $B_l(r)$ functions and includes
1527   formulas given in his appendix.\cite{Smith98} The first function
1528   $B_0(r)$ is defined by
# Line 1575 | Line 1685 | $n$ is eliminated.
1685   under the column ``Bare Coulomb.''  Equations \ref{eq:b9} to
1686   \ref{eq:b13} are still correct for GSF electrostatics if the subscript
1687   $n$ is eliminated.
1578
1579 \section{Extra Material}
1580 %
1581 %
1582 %Energy in body coordinate form ---------------------------------------------------------------
1583 %
1584 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 %
1732
1733
1734 % BODY coordinates force equations --------------------------------------------
1735 %
1736 %
1737 Here are the force equations written in terms of body coordinates.
1738 %
1739 % f ca cb
1740 %
1741 \begin{equation}
1742 \mathbf{F}_{{\bf a}C_{\bf a}C_{\bf b}} =
1743 \frac{C_{\bf a} C_{\bf b}}{4\pi \epsilon_0}  w_a(r) \hat{r}
1744 \end{equation}
1745 %
1746 % f ca db
1747 %
1748 \begin{equation}
1749 \mathbf{F}_{{\bf a}C_{\bf a}D_{\bf b}} =
1750 \frac{C_{\bf a}}{4\pi \epsilon_0}  
1751 \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n} w_b(r) \hat{r}
1752 +\frac{C_{\bf a}}{4\pi \epsilon_0}  
1753 \sum_n  D_{\mathbf{b}n} \hat{b}_n w_c(r)
1754 \end{equation}
1755 %
1756 % f ca qb
1757 %
1758 \begin{equation}
1759 \begin{split}
1760 % 1
1761 \mathbf{F}_{{\bf a}C_{\bf a}Q_{\bf b}} =
1762 \frac{1}{4\pi \epsilon_0}  
1763 C_{\bf a }\text{Tr}Q_{\bf b} w_d(r) \hat{r}
1764 + 2C_{\bf a } \sum_l  \hat{b}_l Q_{{\mathbf b}ln} (\hat{b}_n \cdot \hat{r}) w_e(r) \\
1765 % 2
1766 +\frac{C_{\bf a}}{4\pi \epsilon_0}
1767 \sum_{mn} (\hat{r} \cdot \hat{b}_m) Q_{{\mathbf b}mn} (\hat{b}_n \cdot \hat{r}) w_f(r) \hat{r}
1768 \end{split}
1769 \end{equation}
1770 %
1771 % f da cb
1772 %
1773 \begin{equation}
1774 \mathbf{F}_{{\bf a}D_{\bf a}C_{\bf b}} =
1775 -\frac{C_{\bf{b}}}{4\pi \epsilon_0}
1776 \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n} w_b(r)  \hat{r}
1777 -\frac{C_{\bf{b}}}{4\pi \epsilon_0}
1778 \sum_n  D_{\mathbf{a}n} \hat{a}_n w_c(r)
1779 \end{equation}
1780 %
1781 % f da db
1782 %
1783 \begin{equation}
1784 \begin{split}
1785 % 1
1786 \mathbf{F}_{{\bf a}D_{\bf a}D_{\bf b}} &=
1787 -\frac{1}{4\pi \epsilon_0}
1788 \sum_{mn} D_{\mathbf {a}m}
1789 (\hat{a}_m \cdot   \hat{b}_n)
1790 D_{\mathbf{b}n}  w_d(r) \hat{r}
1791 -\frac{1}{4\pi \epsilon_0}
1792 \sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1793 \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} w_f(r) \hat{r} \\
1794 % 2
1795 & \quad + \frac{1}{4\pi \epsilon_0}
1796 \Bigl[ \sum_m D_{\mathbf {a}m}
1797 \hat{a}_m \sum_n D_{\mathbf{b}n}
1798 (\hat{b}_n \cdot \hat{r})
1799 + \sum_m D_{\mathbf {b}m}
1800 \hat{b}_m \sum_n D_{\mathbf{a}n}
1801 (\hat{a}_n \cdot \hat{r}) \Bigr] w_e(r)  \\
1802 \end{split}
1803 \end{equation}
1804 %
1805 % f da qb
1806 %
1807 \begin{equation}
1808 \begin{split}
1809 % 1
1810 &\mathbf{F}_{{\bf a}D_{\bf a}Q_{\bf b}} =
1811 - \frac{1}{4\pi \epsilon_0} \Bigl[
1812 \text{Tr}Q_{\mathbf{b}}
1813 \sum_l  D_{\mathbf{a}l} \hat{a}_l
1814 +2\sum_{lmn} D_{\mathbf{a}l}
1815 (\hat{a}_l \cdot \hat{b}_m)
1816 Q_{\mathbf{b}mn} \hat{b}_n  \Bigr] w_g(r) \\
1817 % 3
1818 &  - \frac{1}{4\pi \epsilon_0} \Bigl[
1819 \text{Tr}Q_{\mathbf{b}}
1820 \sum_n (\hat{r} \cdot \hat{a}_n) D_{\mathbf{a}n}
1821 +2\sum_{lmn}D_{\mathbf{a}l}
1822 (\hat{a}_l \cdot \hat{b}_m)
1823 Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
1824 % 4
1825 &+ \frac{1}{4\pi \epsilon_0}
1826 \Bigl[\sum_l  D_{\mathbf{a}l} \hat{a}_l
1827 \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1828 Q_{{\mathbf b}mn}
1829 (\hat{b}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{a}_l)
1830 D_{\mathbf{a}l}
1831 \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1832 Q_{{\mathbf b}mn} \hat{b}_n \Bigr]   w_i(r)\\
1833 % 6
1834 &  -\frac{1}{4\pi \epsilon_0}
1835 \sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
1836 \sum_{mn} (\hat{r} \cdot \hat{b}_m)
1837 Q_{{\mathbf b}mn}
1838 (\hat{b}_n \cdot \hat{r})  w_j(r)  \hat{r}
1839 \end{split}
1840 \end{equation}
1841 %
1842 % force qa cb
1843 %
1844 \begin{equation}
1845 \begin{split}
1846 % 1
1847 \mathbf{F}_{{\bf a}Q_{\bf a}C_{\bf b}} &=
1848 \frac{1}{4\pi \epsilon_0}  
1849 C_{\bf b }\text{Tr}Q_{\bf a} \hat{r} w_d(r)
1850 + \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) \\
1851 % 2
1852 &  +\frac{C_{\bf b}}{4\pi \epsilon_0}
1853 \sum_{mn} (\hat{r} \cdot \hat{a}_m) Q_{{\mathbf a}mn} (\hat{a}_n \cdot \hat{r}) w_f(r) \hat{r}
1854 \end{split}
1855 \end{equation}
1856 %
1857 % f qa db
1858 %
1859 \begin{equation}
1860 \begin{split}
1861 % 1
1862 &\mathbf{F}_{{\bf a}Q_{\bf a}D_{\bf b}} =
1863 \frac{1}{4\pi \epsilon_0} \Bigl[
1864 \text{Tr}Q_{\mathbf{a}}
1865 \sum_l  D_{\mathbf{b}l} \hat{b}_l
1866 +2\sum_{lmn} D_{\mathbf{b}l}
1867 (\hat{b}_l \cdot \hat{a}_m)
1868 Q_{\mathbf{a}mn} \hat{a}_n  \Bigr]
1869 w_g(r)\\
1870 % 3
1871 &  + \frac{1}{4\pi \epsilon_0} \Bigl[
1872 \text{Tr}Q_{\mathbf{a}}
1873 \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf{b}n}
1874 +2\sum_{lmn}D_{\mathbf{b}l}
1875 (\hat{b}_l \cdot \hat{a}_m)
1876 Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r}) \Bigr] w_h(r) \hat{r} \\
1877 % 4
1878 &  + \frac{1}{4\pi \epsilon_0} \Bigl[ \sum_l  D_{\mathbf{b}l} \hat{b}_l
1879 \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1880 Q_{{\mathbf a}mn}
1881 (\hat{a}_n \cdot \hat{r}) +2 \sum_l (\hat{r} \cdot \hat{b}_l)
1882 D_{\mathbf{b}l}
1883 \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1884 Q_{{\mathbf a}mn} \hat{a}_n \Bigr]   w_i(r) \\
1885 % 6
1886 &  +\frac{1}{4\pi \epsilon_0}
1887 \sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
1888 \sum_{mn} (\hat{r} \cdot \hat{a}_m)
1889 Q_{{\mathbf a}mn}
1890 (\hat{a}_n \cdot \hat{r})  w_j(r)  \hat{r}
1891 \end{split}
1892 \end{equation}
1893 %
1894 % f qa qb
1895 %
1896 \begin{equation}
1897 \begin{split}
1898 &\mathbf{F}_{{\bf a}Q_{\bf a}Q_{\bf b}} =
1899 \frac{1}{4\pi \epsilon_0} \Bigl[
1900 \text{Tr}Q_{\mathbf{a}} \text{Tr}Q_{\mathbf{b}}
1901 + 2 \sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
1902 Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
1903 (\hat{a}_m \cdot \hat{b}_n) \Bigr] w_k(r) \hat{r}\\
1904 &+\frac{1}{4\pi \epsilon_0} \Bigl[
1905 2\text{Tr}Q_{\mathbf{b}} \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  \hat{a}_m  
1906 + 2\text{Tr}Q_{\mathbf{a}} \sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm}  \hat{b}_m \\
1907 &+ 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})  
1908 + 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
1909 \Bigr] w_n(r) \\
1910 &+ \frac{1}{4\pi \epsilon_0}
1911 \Bigl[ \text{Tr}Q_{\mathbf{a}}
1912 \sum_{lm} (\hat{r} \cdot \hat{b}_l) Q_{\mathbf{b}lm} (\hat{b}_m \cdot \hat{r})
1913 + \text{Tr}Q_{\mathbf{b}}
1914 \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm}  (\hat{a}_m \cdot \hat{r}) \\
1915 &+4\sum_{lmnp} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{b}_n)
1916 Q_{\mathbf{b}np}  (\hat{b}_p \cdot \hat{r}) \Bigr] w_l(r) \hat{r} \\
1917 %
1918 &+\frac{1}{4\pi \epsilon_0} \Bigl[
1919 2\sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} \hat{a}_m
1920 \sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_n \cdot \hat{r}) \\
1921 &+2 \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
1922 \sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} \hat{b}_n \Bigr] w_o(r) \hat{r} \\
1923 &  + \frac{1}{4\pi \epsilon_0}
1924 \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{\mathbf{a}lm} (\hat{a}_m \cdot \hat{r})
1925 \sum_{np} (\hat{r} \cdot \hat{b}_n) Q_{\mathbf{b}np} (\hat{b}_p \cdot \hat{r}) w_m(r) \hat{r}
1926 \end{split}
1927 \end{equation}
1928 %
1929 Here we list the form of the non-zero damped shifted multipole torques showing
1930 explicitly dependences on body axes:
1931 %
1932 %  t ca db
1933 %
1934 \begin{equation}
1935 \mathbf{\tau}_{{\bf b}C_{\bf a}D_{\bf b}} =
1936 \frac{C_{\bf a}}{4\pi \epsilon_0}  
1937 \sum_n  (\hat{r} \times \hat{b}_n)  D_{\mathbf{b}n} \,  v_{11}(r)
1938 \end{equation}
1939 %
1940 % t ca qb
1941 %
1942 \begin{equation}
1943 \mathbf{\tau}_{{\bf b}C_{\bf a}Q_{\bf b}} =
1944 \frac{2C_{\bf a}}{4\pi \epsilon_0}
1945 \sum_{lm} (\hat{r} \times \hat{b}_l) Q_{{\mathbf b}lm} (\hat{b}_m \cdot \hat{r}) v_{22}(r)
1946 \end{equation}
1947 %
1948 %  t da cb
1949 %
1950 \begin{equation}
1951 \mathbf{\tau}_{{\bf a}D_{\bf a}C_{\bf b}} =
1952 -\frac{C_{\bf b}}{4\pi \epsilon_0}  
1953 \sum_n  (\hat{r} \times \hat{a}_n)  D_{\mathbf{a}n} \,  v_{11}(r)
1954 \end{equation}%
1955 %
1956 %
1957 %  ta da db
1958 %
1959 \begin{equation}
1960 \begin{split}
1961 % 1
1962 \mathbf{\tau}_{{\bf a}D_{\bf a}D_{\bf b}} &=
1963 \frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
1964 (\hat{a}_m \times  \hat{b}_n)
1965 D_{\mathbf{b}n} v_{21}(r) \\
1966 % 2
1967 &-\frac{1}{4\pi \epsilon_0}
1968 \sum_m (\hat{r} \times \hat{a}_m) D_{\mathbf {a}m}
1969 \sum_n (\hat{r} \cdot \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
1970 \end{split}
1971 \end{equation}
1972 %
1973 %  tb da db
1974 %
1975 \begin{equation}
1976 \begin{split}
1977 % 1
1978 \mathbf{\tau}_{{\bf b}D_{\bf a}D_{\bf b}} &=
1979 -\frac{1}{4\pi \epsilon_0}  \sum_{mn} D_{\mathbf {a}m}
1980 (\hat{a}_m \times  \hat{b}_n)
1981 D_{\mathbf{b}n} v_{21}(r) \\
1982 % 2
1983 &+\frac{1}{4\pi \epsilon_0}
1984 \sum_m (\hat{r} \cdot \hat{a}_m) D_{\mathbf {a}m}
1985 \sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf {b}n} v_{22}(r)
1986 \end{split}
1987 \end{equation}
1988 %
1989 % ta da qb
1990 %
1991 \begin{equation}
1992 \begin{split}
1993 % 1
1994 \mathbf{\tau}_{{\bf a}D_{\bf a}Q_{\bf b}} &=
1995 \frac{1}{4\pi \epsilon_0} \left(
1996 -\text{Tr}Q_{\mathbf{b}}
1997 \sum_n (\hat{r} \times \hat{a}_n) D_{\mathbf{a}n}
1998 +2\sum_{lmn}D_{\mathbf{a}l}
1999 (\hat{a}_l \times \hat{b}_m)
2000 Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
2001 \right) v_{31}(r)\\
2002 % 2
2003 &-\frac{1}{4\pi \epsilon_0}
2004 \sum_l (\hat{r} \times \hat{a}_l) D_{\mathbf{a}l}
2005 \sum_{mn} (\hat{r} \cdot \hat{b}_m)
2006 Q_{{\mathbf b}mn}
2007 (\hat{b}_n \cdot \hat{r}) v_{32}(r)
2008 \end{split}
2009 \end{equation}
2010 %
2011 % tb da qb
2012 %
2013 \begin{equation}
2014 \begin{split}
2015 % 1
2016 \mathbf{\tau}_{{\bf b}D_{\bf a}Q_{\bf b}} &=
2017 \frac{1}{4\pi \epsilon_0} \left(
2018 -2\sum_{lmn}D_{\mathbf{a}l}
2019 (\hat{a}_l \cdot \hat{b}_m)
2020 Q_{\mathbf{b}mn} (\hat{r} \times \hat{b}_n)
2021 -2\sum_{lmn}D_{\mathbf{a}l}
2022 (\hat{a}_l \times \hat{b}_m)
2023 Q_{\mathbf{b}mn} (\hat{b}_n \cdot \hat{r})
2024 \right) v_{31}(r) \\
2025 % 2
2026 &-\frac{2}{4\pi \epsilon_0}
2027 \sum_l (\hat{r} \cdot \hat{a}_l) D_{\mathbf{a}l}
2028 \sum_{mn} (\hat{r} \cdot \hat{b}_m)
2029 Q_{{\mathbf b}mn}
2030 (\hat{r}\times \hat{b}_n) v_{32}(r)
2031 \end{split}
2032 \end{equation}
2033 %
2034 % ta qa cb
2035 %
2036 \begin{equation}
2037 \mathbf{\tau}_{{\bf a}Q_{\bf a}C_{\bf b}} =
2038 \frac{2C_{\bf a}}{4\pi \epsilon_0}
2039 \sum_{lm} (\hat{r} \cdot \hat{a}_l) Q_{{\mathbf a}lm} (\hat{r} \times \hat{a}_m) v_{22}(r)
2040 \end{equation}
2041 %
2042 % ta qa db
2043 %
2044 \begin{equation}
2045 \begin{split}
2046 % 1
2047 \mathbf{\tau}_{{\bf a}Q_{\bf a}D_{\bf b}} &=
2048 \frac{1}{4\pi \epsilon_0} \left(
2049 2\sum_{lmn}D_{\mathbf{b}l}
2050 (\hat{b}_l \cdot \hat{a}_m)
2051 Q_{\mathbf{a}mn} (\hat{r} \times \hat{a}_n)
2052 +2\sum_{lmn}D_{\mathbf{b}l}
2053 (\hat{a}_l \times \hat{b}_m)
2054 Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
2055 \right) v_{31}(r) \\
2056 % 2
2057 &+\frac{2}{4\pi \epsilon_0}
2058 \sum_l (\hat{r} \cdot \hat{b}_l) D_{\mathbf{b}l}
2059 \sum_{mn} (\hat{r} \cdot \hat{a}_m)
2060 Q_{{\mathbf a}mn}
2061 (\hat{r}\times \hat{a}_n) v_{32}(r)
2062 \end{split}
2063 \end{equation}
2064 %
2065 % tb qa db
2066 %
2067 \begin{equation}
2068 \begin{split}
2069 % 1
2070 \mathbf{\tau}_{{\bf b}Q_{\bf a}D_{\bf b}} &=
2071 \frac{1}{4\pi \epsilon_0} \left(
2072 \text{Tr}Q_{\mathbf{a}}
2073 \sum_n (\hat{r} \times \hat{b}_n) D_{\mathbf{b}n}
2074 +2\sum_{lmn}D_{\mathbf{b}l}
2075 (\hat{a}_l \times \hat{b}_m)
2076 Q_{\mathbf{a}mn} (\hat{a}_n \cdot \hat{r})
2077 \right) v_{31}(r)\\
2078 % 2
2079 &\frac{1}{4\pi \epsilon_0}
2080 \sum_l (\hat{r} \times \hat{b}_l) D_{\mathbf{b}l}
2081 \sum_{mn} (\hat{r} \cdot \hat{a}_m)
2082 Q_{{\mathbf a}mn}
2083 (\hat{a}_n \cdot \hat{r}) v_{32}(r)
2084 \end{split}
2085 \end{equation}
2086 %
2087 % ta qa qb
2088 %
2089 \begin{equation}
2090 \begin{split}
2091 % 1
2092 \mathbf{\tau}_{{\bf a}Q_{\bf a}Q_{\bf b}} &=
2093 -\frac{4}{4\pi \epsilon_0}
2094 \sum_{lmnp} (\hat{a}_l \times \hat{b}_p)
2095 Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
2096 (\hat{a}_m \cdot \hat{b}_n) v_{41}(r) \\
2097 % 2
2098 &+ \frac{1}{4\pi \epsilon_0}
2099 \Bigl[
2100 2\text{Tr}Q_{\mathbf{b}}
2101 \sum_{lm} (\hat{r} \cdot \hat{a}_l )
2102 Q_{{\mathbf a}lm}
2103 (\hat{r} \times \hat{a}_m)
2104 +4 \sum_{lmnp}
2105 (\hat{r} \times \hat{a}_l )
2106 Q_{{\mathbf a}lm}
2107 (\hat{a}_m \cdot \hat{b}_n)
2108 Q_{{\mathbf b}np}
2109 (\hat{b}_p \cdot \hat{r}) \\
2110 % 3
2111 &-4 \sum_{lmnp}
2112 (\hat{r} \cdot \hat{a}_l )
2113 Q_{{\mathbf a}lm}
2114 (\hat{a}_m \times \hat{b}_n)
2115 Q_{{\mathbf b}np}
2116 (\hat{b}_p \cdot \hat{r})
2117 \Bigr] v_{42}(r) \\
2118 % 4
2119 &+ \frac{2}{4\pi \epsilon_0}
2120 \sum_{lm} (\hat{r} \times \hat{a}_l)
2121 Q_{{\mathbf a}lm}
2122 (\hat{a}_m \cdot \hat{r})
2123 \sum_{np}  (\hat{r} \cdot \hat{b}_n)
2124 Q_{{\mathbf b}np}
2125 (\hat{b}_p \cdot \hat{r})  v_{43}(r)\\
2126 \end{split}
2127 \end{equation}
2128 %
2129 % tb qa qb
2130 %
2131 \begin{equation}
2132 \begin{split}
2133 % 1
2134 \mathbf{\tau}_{{\bf b}Q_{\bf a}Q_{\bf b}} &=
2135 \frac{4}{4\pi \epsilon_0}
2136 \sum_{lmnp} (\hat{a}_l \cdot \hat{b}_p)
2137 Q_{\mathbf{a}lm}  Q_{\mathbf{b}np}
2138 (\hat{a}_m \times \hat{b}_n) v_{41}(r) \\
2139 % 2
2140 &+ \frac{1}{4\pi \epsilon_0}
2141 \Bigl[
2142 2\text{Tr}Q_{\mathbf{a}}
2143 \sum_{lm} (\hat{r} \cdot \hat{b}_l )
2144 Q_{{\mathbf b}lm}
2145 (\hat{r} \times \hat{b}_m)
2146 +4 \sum_{lmnp}
2147 (\hat{r} \cdot \hat{a}_l )
2148 Q_{{\mathbf a}lm}
2149 (\hat{a}_m \cdot \hat{b}_n)
2150 Q_{{\mathbf b}np}
2151 (\hat{r} \times \hat{b}_p) \\
2152 % 3
2153 &+4 \sum_{lmnp}
2154 (\hat{r} \cdot \hat{a}_l )
2155 Q_{{\mathbf a}lm}
2156 (\hat{a}_m \times \hat{b}_n)
2157 Q_{{\mathbf b}np}
2158 (\hat{b}_p \cdot \hat{r})
2159 \Bigr] v_{42}(r)  \\
2160 % 4
2161 &+ \frac{2}{4\pi \epsilon_0}
2162 \sum_{lm} (\hat{r} \cdot \hat{a}_l)
2163 Q_{{\mathbf a}lm}
2164 (\hat{a}_m \cdot \hat{r})
2165 \sum_{np}  (\hat{r} \times \hat{b}_n)
2166 Q_{{\mathbf b}np}
2167 (\hat{b}_p \cdot \hat{r}) v_{43}(r). \\
2168 \end{split}
2169 \end{equation}
2170 %
2171 \begin{table*}
2172 \caption{\label{tab:tableFORCE2}Radial functions used in the force equations.}
2173 \begin{ruledtabular}
2174 \begin{tabular}{|l|l|l|}
2175 Generic&Taylor-shifted Force&Gradient-shifted Force
2176 \\ \hline
2177 %
2178 %
2179 %
2180 $w_a(r)$&
2181 $g_0(r)$&
2182 $g(r)-g(r_c)$ \\
2183 %
2184 %
2185 $w_b(r)$ &
2186 $\left( -\frac{g_1(r)}{r}+h_1(r) \right)$ &
2187 $h(r)- h(r_c) - \frac{v_{11}(r)}{r} $ \\
2188 %
2189 $w_c(r)$ &
2190 $\frac{g_1(r)}{r} $ &
2191 $\frac{v_{11}(r)}{r}$ \\
2192 %
2193 %
2194 $w_d(r)$&
2195 $\left( -\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right) $ &
2196 $\left( -\frac{g(r)}{r^2} + \frac{h(r)}{r} \right)
2197 -\left( -\frac{g(r_c)}{r_c^2} + \frac{h(r_c)}{r_c} \right) $\\
2198 %
2199 $w_e(r)$ &
2200 $\left(-\frac{g_2(r)}{r^2} + \frac{h_2(r)}{r} \right)$ &
2201 $\frac{v_{22}(r)}{r}$ \\
2202 %
2203 %
2204 $w_f(r)$&
2205 $\left( \frac{3g_2(r)}{r^2}-\frac{3h_2(r)}{r}+s_2(r) \right)$ &
2206 $\left( \frac{g(r)}{r^2}-\frac{h(r)}{r}+s(r) \right) - $ \\
2207 &&$\left( \frac{g(r_c)}{r_c^2}-\frac{h(r_c)}{r_c}+s(r_c) \right)-\frac{2v_{22}(r)}{r}$\\
2208 %
2209 $w_g(r)$& $ \left( -\frac{g_3(r)}{r^3}+\frac{h_3(r)}{r^2} \right)$&
2210 $\frac{v_{31}(r)}{r}$\\
2211 %
2212 $w_h(r)$ &
2213 $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
2214 $\left(\frac{2g(r)}{r^3} -\frac{2h(r)}{r^2} +\frac{s(r)}{r} \right) - $\\
2215 &&$\left(\frac{2g(r_c)}{r_c^3} -\frac{2h(r_c)}{r_c^2} +\frac{s(r_c)}{r_c} \right) $ \\
2216 &&$-\frac{v_{31}(r)}{r}$\\
2217 % 2
2218 $w_i(r)$ &
2219 $\left(\frac{3g_3(r)}{r^3} -\frac{3h_3(r)}{r^2} +\frac{s_3(r)}{r} \right) $  &
2220 $\frac{v_{32}(r)}{r}$ \\
2221 %
2222 $w_j(r)$ &
2223 $\left(\frac{-15g_3(r)}{r^3} + \frac{15h_3(r)}{r^2} - \frac{6s_3(r)}{r} + t_3(r) \right)  $ &
2224 $\left(\frac{-6g(r)}{r^3} +\frac{6h(r)}{r^2} -\frac{3s(r)}{r} +t(r) \right) $  \\
2225 &&$\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}$ \\
2226 %
2227 $w_k(r)$ &
2228 $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
2229 $\left(\frac{3g(r)}{r^4} -\frac{3h(r)}{r^3} +\frac{s(r)}{r^2}  \right)$  \\
2230 &&$\left(\frac{3g(r_c)}{r_c^4} -\frac{3h(r_c)}{r_c^3} +\frac{s(r_c)}{r_c^2}  \right)$ \\
2231 %
2232 $w_l(r)$ &
2233 $\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)$ &
2234 $\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)$ \\
2235 &&$\left(-\frac{9g(r)}{r^4} +\frac{9h(r)}{r^3} -\frac{4s(r)}{r^2} +\frac{t(r)}{r} \right)
2236 -\frac{2v_{42}(r)}{r}$ \\
2237 %
2238 $w_m(r)$ &
2239 $\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)$ &
2240 $\left(\frac{45g(r)}{r^4} -\frac{45h(r)}{r^3} +\frac{21s(r)}{r^2} -\frac{6t(r)}{r} +u(r) \right)$ \\
2241 &&$\left(\frac{45g(r_c)}{r_c^4} -\frac{45h(r_c)}{r_c^3}
2242 +\frac{21s(r_c)}{r_c^2} -\frac{6t(r_c)}{r_c} +u(r_c) \right) $ \\
2243 &&$-\frac{4v_{43}(r)}{r}$ \\
2244 %
2245 $w_n(r)$ &
2246 $\left(\frac{3g_4(r)}{r^4} -\frac{3h_4(r)}{r^3} +\frac{s_4(r)}{r^2}  \right)$ &
2247 $\frac{v_{42}(r)}{r}$ \\
2248 %
2249 $w_o(r)$ &
2250 $\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)$ &
2251 $\frac{v_{43}(r)}{r}$ \\
2252 %
2253 \end{tabular}
2254 \end{ruledtabular}
2255 \end{table*}
1688  
1689   \newpage
1690  

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines