Multibody Multipole Methods
Abstract
A threebody potential function can account for interactions among triples of particles which are uncaptured by pairwise interaction functions such as Coulombic or LennardJones potentials. Likewise, a multibody potential of order can account for interactions among tuples of particles uncaptured by interaction functions of lower orders. To date, the computation of multibody potential functions for a large number of particles has not been possible due to its scaling cost. In this paper we describe a fast treecode for efficiently approximating multibody potentials that can be factorized as products of functions of pairwise distances. For the first time, we show how to derive a BarnesHut type algorithm for handling interactions among more than two particles. Our algorithm uses two approximation schemes: 1) a deterministic series expansionbased method; 2) a Monte Carlobased approximation based on the central limit theorem. Our approach guarantees a userspecified bound on the absolute or relative error in the computed potential with an asymptotic probability guarantee. We provide speedup results on a threebody dispersion potential, the AxilrodTeller potential.
keywords:
Fast multipole methods; Data structures; kdtrees; AxilrodTeller potential; Multitree algorithms1 Introduction
In this paper, we generalize previous algorithmic frameworks for rapidly computing pairwise summations to include higherorder summations. Suppose we are given a set of particles in dimensional space.
For and a tuple function , we are interested in computing the following form^{1}^{1}1In computing , we fix one of the arguments of as and choose a subset from which does not contain .:
(1) 
Sums of the form Equation (1) occur in molecular dynamics, protein structure prediction, and other similar contexts. Biomolecular simulations usually break down the interactions in complex chemical systems into ballsandsprings mechanical models augmented by torsional terms, pairwise point charge electrostatic terms, and simple pairwise dispersion (van der Waals) interactions, etc. However, such pairwise () interactions often fail to capture important, complex nonadditive interactions found in real systems. Though many researchers have argued that multibody potentials enable more accurate and realistic molecular modeling, the evaluation of body forces for in systems beyond tiny sizes (less than 10,000 particles) has not been possible due to the unavailability of an efficient way to realize the computation.
In this paper we focus on computing multibody potentials of the third order (), but frame our presentation so that the methods can easily be generalized to handle higherorder potentials. For concreteness, we consider the AxilrodTeller potential (dispersion potential):
(2) 
where , , are the angles at the vertices of the triangle and is the Euclidean distance metric. This potential axilrod1943interaction () describes induced dipole interactions between triples of atoms, and is known to be important for the accurate computation of the physical properties of certain noble gases.
This Paper. For the first time, we introduce a fast algorithm for efficiently computing multibody potentials for a large number of particles. We restrict the class of multibody potentials to those that can be factorized as products of functions of pairwise Euclidean distances. That is,
(3) 
Our algorithm achieves speedup by utilizing two approximation methods: a deterministic and a probabilistic one. The deterministic approximation is based on the analytic seriesexpansionbased approach in cheng1999fam (); greengard2002nvf (); greengard1991fgt () to handle potential functions that describe body interactions with . The probabilistic approach uses a Monte Carlobased approximation based on the central limit theorem. Our algorithm can compute multibody potentials within userspecified bounds for relative or absolute error with an asymptotic probability guarantee.
However, we would like to point out the following limitations in our algorithm. First of all, we do not present a fullfledged derivation of all three translation operators (namely the fartofar, the fartolocal, and the localtolocal translation operators) for the general multibody case. While we define the farfield expansion for a restricted class of multibody potentials, defining the local expansion for this same class is harder (see Section 3.4). We would also like to point out that the hybrid deterministic/probabilistic approximation heuristic works under some partial distributions but not all. Indeed, there are configurations for which the speedup factor over the naive bruteforce method is minimal. The MonteCarlo based approximation relies on two theorems: 1) the central limit theorem from which we determine the number of required samples; 2) the BerryEsseen theorem which characterizes the the rate at which the sample average converges to the true average. Both theorems provide only asymptotic guarantees.
Our work utilizes and extends a framework for efficient algorithms for socalled generalized Body Problems gray2001nbp (), which introduced multitree methods. The framework was originally developed to accelerate common bottleneck statistical computations based on distances; it utilizes multiple trees and other spatial data structures to reduce computation times both asymptotically and practically by multiple orders of magnitude. This work extends the framework with higherorder hierarchical series approximation techniques, demonstrating a fast multipoletype method for higherorder interactions for the first time, effectively creating a Multibody Multipole Method.
Section 3 introduces the generalized body framework and describes a partial extension of fast multipoletype methods to handle higherorder interactions; we will discuss the technical difficulties for deriving all of the necessary tools for the general multibody case. As a result, we utilize only a simple but effective approximation using the centerofmass approximations. Section 4 focuses on threebody interactions and introduces methods to do potential computations under both deterministic and probabilistic error criteria; the section also provides a description of the fast algorithm for the threebody case. Section 5 proves that our proposed algorithms can approximate potentials within userspecified error bounds. Section 6 shows experimental scalability results for our proposed algorithms against the naive algorithm under different error parameter settings.
Notations. Throughout this paper, we use these common sets of notations:

(Normal Distribution). This is denoted by where and are the mean and the covariance respectively.

(Vector Component). For a given vector , we access its th component by where (i.e. based index).

(Multiindex Notation). Throughout this paper, we will be using the multiindex notation. A dimensional multiindex is a tuple of nonnegative integers and will be denoted using a bold lowercase Greek alphabet. For any dimensional multiindices , and any ,






for .
where is a th directional partial derivative. Define if , and for if for (and similarly for ).


(Size of a Point Set). Given a set , it size is denoted by .

(Probability Guarantee). We use the unbold Greek alphabet .

(A Tree Node). A tree node represents a subset of a point set represented by the root node. Hence, we use the same notation as the previous.

(Representative Point of a Tree Node). Usually a geometric center is used but any point inside the bounding primitive of a tree node is chosen as well. For the tree node , this is denoted as .

(Child Nodes of an Internal Tree Node). Given a node , denote its left and right child nodes by and respectively.
2 Related Work
2.1 Error Bounds
Due to its expensive computational cost, many algorithms approximate sums at the expense of reduced precision. The following error bounding criteria are used in the literature:
Definition 2.1.
absolute error bound: For each for , it computes such that .
Definition 2.2.
relative error bound: For each for , compute such that .
Bounding the relative error is much harder because the error bound criterion is in terms of the initially unknown exact quantity. As a result, many previous methods greengard1991fgt (); yang2003improved () have focused on bounding the absolute error. The relative error bound criterion is preferred to the absolute error bound criterion in statistical applications in which high accuracy is desired. Our framework can enforce the following error form:
Definition 2.3.
probabilistic relative/ absolute error: For each for , compute , such that with at least probability , .
2.2 Series Expansion
A series of papers first laid the foundations for efficiently computing sums of pairwise potentials such as Coulombic and Yukawa potentials cheng1999fam (); greengard2002nvf (); greengard1991fgt (). The common approach in these papers is to derive analytical series expansions of the given potential function in either Cartesian or spherical coordinate systems. The series expansion is then truncated after taking a fixed number of terms. The associated error bounds are derived from summing the truncated terms in an appropriate infinite geometric sum or bounding the remainder term using Taylor’s theorem. A recent line of work on efficient computation of pairwise function has focused on developing numerical representations of the potential matrix , rather than relying on analytical expansion of the potential function. martinsson2008accelerated () and kapur1998iee () use singular value decomposition and the QR decomposition to compute the compressed forms of the potential function and the three translation operators. anderson1992ifm (); ying2004kia () take the “pseudoparticle” approach by placing equivalent artificial charges on the bounding surface of the actual particles by solving appropriate integral equations. All of these works have been limited to pairwise potential functions, and the approach does not naturally suggest a generalization to body potentials with . To our knowledge, no research has been performed on the problem of evaluating multibody potentials using a method more sophisticated than the bruteforce algorithm with an adhoc cutoff distance. marcelli2001tro (); marcelli2009beyond ().
3 Generalized body Framework
We use a variant of trees bentley1975multidimensional () to form hierarchical groupings of points based on their locations using the recursive procedure shown in Algorithm 1. Initially, the algorithm starts with (the entire point set). We split a given set of points along the widest dimension of the bounding hyperrectangle into two equal halves at the splitting coordinate. We continue splitting until the number of points is below some userdefined threshold called the leaf threshold. If the number of points owned by a node exceeds the leaf threshold, then it is called an internal node. Otherwise it is called a leaf node. Assuming that each split on a level results in the equal number of points on the left subset and the right subset and respectively, the runtime cost is . We note that the cost of building a tree is negligible compared to the actual multibody potential computation (see Section 6). See Figure 2.
The general framework for computing Equation (1) is formalized in gray2001nbp (); gray2003nde (); gray2003vfm (); gray2003rem (). This approach consists of the following steps:

Build a spatial tree (such as kdtrees) for the set of particles and build farfield moments on each node of the tree (Bottomup phase).

Perform a multitree traversal over tuples of nodes (Approximation phase).

Preorder traverse the tree and propagate unincorporated bound changes downward (Topdown phase).
Step 2 utilizes the procedure shown in Algorithm 2 (called by setting each for ), a recursive function that allows us to consider the tuples formed by choosing each from ; we can gain efficiency over the naive enumeration of the tuples by using the bounding box and the moment information stored in each . One such information is the distance bound computed using the bounding box (see Figure 3).
CanSummarize function first eliminates redundant recursive calls for the list of node tuples that satisfy the following condition: if there exists a pair of nodes and () among the node list , such that the maximum depthfirst rank of is less than the minimum depthfirst rank of . In this case, the function returns true. See Figure 2 and moore2001fast (). In addition, if any one of the nodes in the list includes one of the other nodes (i.e. there exists nodes and such that the minimum depthfirst rank of the minimum depthfirst rank of the maximum depthfirst rank of the maximum depthfirst rank of ), CanSummarize returns false. We do this because it is a bit tricky to count the number of tuples for each point in this case (see Figure 4).
Otherwise, CanSummarize function tests whether each potential sum for can be approximated within the error tolerance determined by the algorithm. For example, if , we test for each , , , , the following exact quantities can be approximated:
If the approximation is not possible, then the algorithm continues to consider the data at a finer granularity; it chooses an internal node (typically the one with the largest diameter) to split among . Before recursing to two subcalls in Line 9 and Line 10 of Algorithm 2, the algorithm can optionally push quantities from a node that is being split to its child nodes (Line 8). After returning from the recursive calls, the node that was just split can refine summary statistics based on the results accumulated on its child nodes. The details of these operations are available in earlier papers gray2001nbp (); gray2003nde (); gray2003vfm (); gray2003rem (); 1102.2878 ().
The basic idea is to terminate the recursion as soon as possible, i.e. by considering a tuple of large subsets and avoiding the number of exhaustive leafleafleaf computations. We note that the CanSummarize and Summarize functions effectively replace unwieldy interaction lists used in FMM algorithms. Interaction lists in tuple interaction, if naively enumerated, can be large depending on the potential function and the dimensionality of the problem, whereas the generalized body approach can handle a wide spectrum of problems without this drawback.
3.1 Algorithm for Pairwise Potentials ()
The general algorithmic strategy for pairwise potentials is described in gray2001nbp (); gray2003nde (); gray2003vfm (); gray2003rem (), and consists of the following three main phases (see Figure 5). Suppose we are given a set of “source” points (denoted as reference points) and a set of “target” points (denoted as query points)^{2}^{2}2The terms “reference/query” have been used in the general framework we are applying.

Bottomup phase: Compute farfield moments of order in every leaf node of the reference tree. The resulting farfield expansion of each reference node is given by:
reads as “the potential sum on due to the contribution of ” and as “the th farfield coefficient of centered at .” Because it is impossible to store an infinite number of farfield moments , we truncate the Taylor expansion up to the order (determined either arbitrarily or by an appropriate error criterion):
(4) such that is sufficiently small. reads as “the approximated potential sum on due to the points owned by using up to the th order farfield expansion of centered at .”
For internal reference nodes, perform the fartofar (F2F) translation to convert the farfield moments owned by the child nodes to form the farfield moments for their common parent node . For example, the farfield moments of centered at is shifted to by:
(5) where
(6) Note that there is no error incurred in each F2F translation, i.e. for any query point from the intersection of the domains of for and ; the domain for which the farfield expansion remains valid depends on the error bound criterion for each potential. The farfield moments of the parent node is the sum of the translated moments of its child nodes:

Approximation phase: For a given pair of the query and the reference nodes, determine the order of approximation and either (1) translate the farfield moments of the reference node to the local moments of the query node (2) or recurse to their subsets, if the F2L translation is more costly than the direct exhaustive method.
Let us rewrite the exact contribution of to a point :
(7) where reads as “the exact local moments ^{3}^{3}3We use to denote the local moments because a “nearfield” expansion is another widely used term for a local expansion. It avoids the potential notational confusion in the later parts of the paper. contributed by the points in centered at .” Truncating Equation (7) at for some yields a direct local accumulation of order .
From the bottomup phase, we know that . Similarly, we can store only a finite number of local moments up to the order and thus . We get the local expansion for formed due to translated farfield moments of :
(8) where reads as “approximation to the exact local moments ” and as “the approximated potential sum on due to the points in using up to the th order inexact local moments centered at ”. The F2L translation is applied only if is sufficiently small.

Topdown phase: Propagate the local moments of each query node (i.e. pruned quantities) to its child nodes using the localtolocal (L2L) operator. Suppose we have the following local expansion for :
where is the maximum approximation order among (1) the F2L translations performed for and all of the ancestor nodes of (denoted by ); and (2) the direct local accumulations of and those passed down from all of the ancestors of (denoted by ). Shifting the expansion to another center is given by:
(9) (10) This shifted moments are added to the local moments of each child of , in effect transmitting the pruned contributions downward. At each query leaf, we evaluate the resulting local expansion at each query point.
3.2 Farfield Expansion for Threebody Potentials ()
In this section, we define farfield expansions for a threebody potential that is a product of functions of pairwise distances (see Equation (1)):
(11) 
We define the farfield moments of a node the same way defined for the pairwise potential case. Suppose we are given three nodes from the tree. The following nested sum expresses the contribution for due to the other nodes and :
(12) 
The basic goal here is to decompose Equation (12) into sums of products of the farfield moments of each node. A farfield expansion for induced by the farfield moments of and is given by (see Figure 6):
By setting and pushing the summations over and inside, we get:
(13) 
Truncating at th order yields:
(14) 
where reads as “the th order farfield expansion at due to the moments of centered at and the moments of centered at .”
Computational Cost of Evaluating the Farfield Expansion. The first three summations over , , collectively contribute terms, and the inner summation contributing at most terms. Thus, evaluating the th order farfield expansion for a threebody potential on a single point takes time.
3.3 Farfield Expansion for General Multibody Potentials
For a general multibody potential that can be expressed as products of pairwise functions (see Equation (3)), the farfield expansion induced by the points in for is:
Focus on grouping and multiplying monomial powers of for each :
Let and . Then,
(15) 
Equation (15) is a convolution of farfield moments of . We can truncate the expansion above for terms for for some . Note that Equation (15) includes the and cases.