Multibody Multipole Methods

# Multibody Multipole Methods

Dongryeol Lee Computational Science and Engineering
Georgia Institute of Technology
266 Ferst Drive
Atlanta, GA 30332
Arkadas Ozakin Georgia Tech Research Institute
Georgia Institute of Technology
266 Ferst Drive
Atlanta, GA 30332
Alexander G. Gray
###### Abstract

A three-body potential function can account for interactions among triples of particles which are uncaptured by pairwise interaction functions such as Coulombic or Lennard-Jones 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 tree-code 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 Barnes-Hut type algorithm for handling interactions among more than two particles. Our algorithm uses two approximation schemes: 1) a deterministic series expansion-based method; 2) a Monte Carlo-based approximation based on the central limit theorem. Our approach guarantees a user-specified bound on the absolute or relative error in the computed potential with an asymptotic probability guarantee. We provide speedup results on a three-body dispersion potential, the Axilrod-Teller potential.

###### keywords:
Fast multipole methods; Data structures; kd-trees; Axilrod-Teller potential; Multi-tree algorithms
journal: Journal of Computational Physics

## 1 Introduction

In this paper, we generalize previous algorithmic frameworks for rapidly computing pair-wise summations to include higher-order summations. Suppose we are given a set of particles in -dimensional space.

For and a -tuple function , we are interested in computing the following form111In computing , we fix one of the arguments of as and choose a -subset from which does not contain .:

 Φ(x;X×⋯×X(n−1) copies)=∑xi2∈X∖{x} ∑xi3∈X∖{x}i2

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 balls-and-springs 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 non-additive 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 higher-order potentials. For concreteness, we consider the Axilrod-Teller potential (dispersion potential):

 ϕ(xi,xj,xk)=1+3cosθicosθjcosθk||xi−xj||3||xi−xk||3||xj−xk||3 (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,

 ϕ(xi1,⋯,xin)=∏1≤p

Our algorithm achieves speedup by utilizing two approximation methods: a deterministic and a probabilistic one. The deterministic approximation is based on the analytic series-expansion-based approach in cheng1999fam (); greengard2002nvf (); greengard1991fgt () to handle potential functions that describe -body interactions with . The probabilistic approach uses a Monte Carlo-based approximation based on the central limit theorem. Our algorithm can compute multibody potentials within user-specified 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 full-fledged derivation of all three translation operators (namely the far-to-far, the far-to-local, and the local-to-local translation operators) for the general multibody case. While we define the far-field 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 brute-force method is minimal. The Monte-Carlo based approximation relies on two theorems: 1) the central limit theorem from which we determine the number of required samples; 2) the Berry-Esseen 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 so-called generalized -Body Problems gray2001nbp (), which introduced multi-tree 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 higher-order hierarchical series approximation techniques, demonstrating a fast multipole-type method for higher-order 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 multipole-type methods to handle higher-order 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 center-of-mass approximations. Section 4 focuses on three-body 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 three-body case. Section 5 proves that our proposed algorithms can approximate potentials within user-specified 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).

• (Multi-index Notation). Throughout this paper, we will be using the multi-index notation. A -dimensional multi-index is a -tuple of non-negative integers and will be denoted using a bold lowercase Greek alphabet. For any -dimensional multi-indices , 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 “pseudo-particle” 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 brute-force algorithm with an ad-hoc cut-off distance. marcelli2001tro (); marcelli2009beyond ().

## 3 Generalized N-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 hyper-rectangle into two equal halves at the splitting coordinate. We continue splitting until the number of points is below some user-defined 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:

1. Build a spatial tree (such as kd-trees) for the set of particles and build far-field moments on each node of the tree (Bottom-up phase).

2. Perform a multi-tree traversal over -tuples of nodes (Approximation phase).

3. Pre-order traverse the tree and propagate unincorporated bound changes downward (Top-down 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 depth-first rank of is less than the minimum depth-first 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 depth-first rank of the minimum depth-first rank of the maximum depth-first rank of the maximum depth-first 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:

 Φ(x1;P2×P3×P4) =∑xi2∈P2∖{x1} ∑xi3∈P3∖{x1}i2
 Φ(x3;P1×P2×P4) =∑xi1∈P1∖{x3} ∑xi2∈P2∖{x3}i1

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 sub-calls 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 leaf-leaf-leaf 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 (n=2)

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)222The terms “reference/query” have been used in the general framework we are applying.

1. Bottom-up phase: Compute far-field moments of order in every leaf node of the reference tree. The resulting far-field expansion of each reference node is given by:

 Φ(x;P2)=∑α≥0⎡⎢⎣∑xi2∈P2(−1)αα!(xi2−cP2)α⎤⎥⎦Dαϕ(x−cP2)=∑α≥0Mα(P2,cP2)Dαϕ(x−cP2)

reads as “the potential sum on due to the contribution of ” and as “the -th far-field coefficient of centered at .” Because it is impossible to store an infinite number of far-field moments , we truncate the Taylor expansion up to the order (determined either arbitrarily or by an appropriate error criterion):

 ˜Φ(x;P2;F(cP2,p))=∑|α|≤pMα(P2,cP2)Dαϕ(x−cP2) (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 far-field expansion of centered at .”

For internal reference nodes, perform the far-to-far (F2F) translation to convert the far-field moments owned by the child nodes to form the far-field moments for their common parent node . For example, the far-field moments of centered at is shifted to by:

 ˜Φ(x;PL2;F(cP2,p))=∑γ≤pMγ(PL2,cP2)(−1)γDγϕ(x−cP2) (5)

where

 Mγ(PL2,cP2)=∑α≤γMα(PL2,cPL2)(cPL2−cP2)γ−α(γ−α)! (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 far-field expansion remains valid depends on the error bound criterion for each potential. The far-field moments of the parent node is the sum of the translated moments of its child nodes:

2. Approximation phase: For a given pair of the query and the reference nodes, determine the order of approximation and either (1) translate the far-field 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 re-write the exact contribution of to a point :

 Φ(x;P2)=∑β≥01β!∑α≥0Mα(P2,cP2)Dα+βϕ(cP1−cP2)(x−cP1)β = ∑β≥0⎡⎢⎣∑xi2∈P21β!Dβϕ(cP1−xi2)⎤⎥⎦(x−cP1)β=∑β≥0Nβ(P2,cP1)(x−cP1)β (7)

where reads as “the exact local moments 333We use to denote the local moments because a “near-field” 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 bottom-up 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 far-field moments of :

 ˜Φ(x;P2;˜N(cP1,p′)) =∑|β|≤p′⎡⎣1β!∑|α|≤p′Mα(P2,cP2)Dα+βϕ1,2(cP1−cP2)⎤⎦(x−cP1)β =∑|β|≤p′˜Nβ(P2,cP1)(x−cP1)β (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.

3. Top-down phase: Propagate the local moments of each query node (i.e. pruned quantities) to its child nodes using the local-to-local (L2L) operator. Suppose we have the following local expansion for :

 ˜Φ(x;F2L(P1)∪DL(P1);˜N(cP1,puP1))=∑|α|≤puP1˜Nα(F2L(P1)∪DL(P1),cP1)(xi1−cP1)α

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:

 ˜Φ(x;F2L(P1)∪DL(P1);˜N(c∗P1,puP1)) (9) = ∑|α|≤puP1⎡⎣∑β≥α(βα)˜Nβ(F2L(P1)∪DL(P1),cP1)(c∗P1−cP1)β−α⎤⎦(x−c∗P1)α = ∑|α|≤puP1˜Nα(F2L(P1)∪DL(P1),c∗P1)(x−c∗P1)α (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 Far-field Expansion for Three-body Potentials (n=3)

In this section, we define far-field expansions for a three-body potential that is a product of functions of pairwise distances (see Equation (1)):

 ϕ(xi1,xi2,xi3)=ϕ1,2(xi1,xi2)⋅ϕ1,3(xi1,xi3)⋅ϕ2,3(xi2,xi3) (11)

We define the far-field 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 :

 Φ(x;P2×P3)=∑xi2∈P2∑xi3∈P3ϕ(x,xi2,xi3) (12)

The basic goal here is to decompose Equation (12) into sums of products of the far-field moments of each node. A far-field expansion for induced by the far-field moments of and is given by (see Figure 6):

 Φ(x;P2×P3) = ∑xi2∈P2∑xi3∈P3∑α1,2≥0(xi2−cP2)α1,2α1,2!(−1)α1,2Dα1,2ϕ1,2(x−cP2) ∑α1,3≥0(xi3−cP3)α1,3α1,3!(−1)α1,3Dα1,3ϕ1,3(x−cP3) ∑α2,3≥0∑β2,3≤α2,3(xi2−cP2)β2,3β2,3!(xi3−cP3)α2,3−β2,3(α2,3−β2,3)!(−1)α2,3−β2,3Dα2,3ϕ2,3(cP2−cP3)

By setting and pushing the summations over and inside, we get:

 Φ(x;P2×P3)= ∑α≥0∑α1,2≤α∑α1,3≤α−α1,2∑β2,3≤α−α1,2−α1,3(α1,2+β2,3α1,2)(α−α1,2−β2,3α1,3) Mα1,2+β2,3(P2,cP2)Mα−α1,2−β2,3(P3,cP3)(−1)β2,3 Dα1,2ϕ1,2(xi1−cP2)Dα1,3ϕ1,3(xi1−cP3)Dα−α1,2−α1,3ϕ2,3(cP2−cP3) (13)

Truncating at -th order yields:

 ˜Φ(x;P2×P3;F(cP2×cP3,p)) = ∑|α|≤p∑α1,2≤α∑α1,3≤α−α1,2∑β2,3≤α−α1,2−α1,3(α1,2+β2,3α1,2)(α−α1,2−β2,3α1,3) Mα1,2+β2,3(P2,cP2)Mα−α1,2−β2,3(P3,cP3)(−1)β2,3 Dα1,2ϕ1,2(xi1−cP2)Dα1,3ϕ1,3(xi1−cP3)Dα−α1,2−α1,3ϕ2,3(cP2−cP3) (14)

where reads as “the -th order far-field expansion at due to the moments of centered at and the moments of centered at .”

Computational Cost of Evaluating the Far-field Expansion. The first three summations over , , collectively contribute terms, and the inner summation contributing at most terms. Thus, evaluating the -th order far-field expansion for a three-body potential on a single point takes time.

### 3.3 Far-field Expansion for General Multibody Potentials (n≥2)

For a general multibody potential that can be expressed as products of pairwise functions (see Equation (3)), the far-field expansion induced by the points in for is:

 Φ(x;P2×⋯×Pn) = ∏2≤k≤n∑xik∈Pk∑α1,k≥0(xik−cPk)α1,kα1,k!(−1)α1,kDα1,kϕ1,k(x−cPk) ∏2≤s

Focus on grouping and multiplying monomial powers of for each :

 (xik−cPk)α1,k+k−1∑u=2(αu,k−βu,k)+n∑v=k+1βk,vα1,k!k−1∏u=2(αu,k−βu,k)!n∏v=k+1βk,v!

Let and . Then,

 Φ(x;P2×⋯×Pn) = ∏2≤s

Equation (15) is a convolution of far-field moments of . We can truncate the expansion above for terms for for some . Note that Equation (15) includes the and cases.

 ˜Φ(x;P2×⋯×Pn;F(cP2×⋯×cPn,p)) = ∏2≤s