Elsevier

Journal of Computational Physics

Volume 226, Issue 1, 10 September 2007, Pages 732-753
Journal of Computational Physics

Accelerated Cartesian expansions – A fast method for computing of potentials of the form Rν for all real ν

https://doi.org/10.1016/j.jcp.2007.04.033Get rights and content

Abstract

The need to compute potentials of the form Rν for ν  1 occur in a variety of areas ranging from electromagnetics to biophysics to molecular dynamics to astrophysics, etc. For instance, Coulomb, London, Lennard-Jones, H-bonds are of the form ν=1, 5, 6 (and 12), 10, respectively. For a systems with N source/observation points, the cost of computing these potentials scales proportional to O(N2). Methods to overcome this computational bottleneck have been a popular research topic for quite some time. For instance, the fast multipole method (FMM) and their cousins—tree codes—have revolutionized the computation of electrostatic potentials (ν = 1). These methods rely on a hierarchical decomposition of the computational domain, i.e, construction of a regular oct-tree decomposition, and exploits the principle of divide and conquer to compute potentials at each observation point. This paper presents two methods, in the vein of FMM, albeit based on Cartesian tensors. The salient features of the first are as follows: (i) it relies on totally symmetric tensors and (ii) the errors are independent of the height of the tree. The second method is presented specifically for ν = 1, and relies on traceless totally symmetric Cartesian tensors. Using the relationship between traceless Cartesian tensors and spherical harmonics, it will be shown that this technique has the same computational cost as the classical FMM. Generalization of the second method for ν  1 is trivial; however, one needs to use totally symmetric tensors instead. Finally, in the whole computation scheme, only the translation operator (that used to traverse across the tree) depends on ν. Convergence of the proposed method is proven for all νR. Numerical results that validate the cost and accuracy are presented for several potential functions; these include those typically encountered in the analysis of physical systems (Coulombic, Lennard-Jones, Lattice gas).

Introduction

Computation of pairwise interaction (for instance, Columbic interactions, London potentials, or Van der-Wall’s potential, etc.) is important in numerous research areas that are as diverse as biophysics, physics, computation chemistry, astrophysics, and electrical engineering to name a few. For example, ν = 1, 5, 6, 10 corresponds to Columbic interactions/gravitational potentials, London dispersion potentials, van der Waals potentials, and H-bonds, respectively. However, it is well known that computing these potentials requires prohibitive computation resources, both in terms of CPU cycles and memory. These costs are exacerbated when such computation is required either in a molecular dynamics or a Monte Carlo setting. It is well known that the CPU cost of computing mutual interactions between N particles distributed in R3 scales as O(N2). Consequently, this has engendered the need for computational methodologies that are efficient both in terms of memory and CPU time. Some of these include cutoff techniques [1], particle mesh algorithms [2], [3], [4], Ewald summation (based on an assumption of periodicity) [5], [6], tree-codes [7], [8], [9] and fast multipole methods (FMM) [10], [11], [12], [13]. Tree codes (and FMM) are based on subdividing the computational domain into hierarchical subdomains, and computing the influence between subdomains that are sufficiently separated using multipole/local expansions. The fundamental differences between FMM and tree codes notwithstanding, these methods have revolutionized analysis in application domains ranging from astrophysics to biophysics to engineering sciences.

At this point, we note that there is rampant confusion in terminology; in fact, the terms FMM and tree codes are used interchangeably in the literature that we have come across. This is not surprising as the two techniques are closely related. The differences between these two methods, albeit subtle, are significant. As was elucidated in [14], tree codes compute interactions between source pairs using one of three methods: (i) directly, (ii) evaluating fields at each observation point using multipole expansion due to a cluster of sources, or (iii) using local expansion at observation clusters to find fields. The decision on the operation used depends on which one is computationally efficient. On the other hand, the algorithmic structure of FMM enables the computation of potentials in an optimal manner [14]. Two addition operations that permit this are aggregation and disaggregation functions. These permit the computation of information at coarser (or finer) using information at finer (or coarser) levels. It so happens that for ν = 1, tree codes typically rely on Cartesian expansions, whereas FMM is based on spherical harmonics. In this paper, we develop theorems and operations necessary for constructing FMM methods (∀ν) using Cartesian tensors.

Tree codes for computing the Lennard-Jones potential was developed as early as 1996 [15]. More refined methods for computing the same were proposed by [16], [17], [18]. All three methods essentially use Taylor’s expansion in the Cartesian coordinates, and some of these use Gegenbauer polynomials based recursion techniques to accelerate translation of multipole expansion to fields at the observer. In 2005, [19] proposed a variation of these schemes using Taylor’s series expansion in a different coordinate systems. Other techniques that have been proposed rely on precorrected fast Fourier transform and using a singular value decomposition. Developing FMM-like techniques using special functions has proven difficult [19] as Gegenbauer polynomials in the spherical coordinate system are not separable. However, it is well known that Gegenbauer polynomials can be written in terms of Legendre polynomials. Sarin et. al. used this fact to develop a tree code; since Legendre polynomials are used, an FMM scheme may be readily derived from these expressions as well. Recently, Chowdry and Jandhyala [20] developed operators necessary to extend their scheme to a multilevel setting. However, using Gegenbauer polynomials for either recursion or developing tree codes has a singular disadvantage; these polynomials are defined for ν   1/2. Consequently, methods that rely explicitly on these cannot be generalized to non-oscillatory potentials (like the lattice gas potential) nor can one prove convergence νR.

The motivation for writing this paper are fourfold: (i) To formulate a fast method for 1/Rν in terms of totally symmetric tensors, and exploit these to reduce the costs. (ii) To introduce new theorems that enable exact traversal up and down the tree. This implies that one does not accrue error as the height of the tree increases. (iii) To prove convergence νR. (iv) To demonstrate the intimate connection between the classical FMM introduced by Greengard [10] and its Cartesian counterparts using traceless Maxwell Cartesian tensors. This connection will show that properly constructed Cartesian FMM schemes have the same computational complexity as classical FMM. Another salient feature of the proposed method is that it can readily form the framework for rapidly computing potentials that are non-oscillatory (e.g. Lattice gas, Yukawa, Gauss Transform, Lattice Sums, etc.) without the use of special functions (application of this methodology to some of these potentials will be presented elsewhere). Furthermore, this algorithm is capable of aggregating different potentials within the same simulation tree.

Development in fast methods for computing Columbic interactions (ν = 1) have progressed along two fronts: exploiting the representation of the potential using either spherical harmonics or Taylor’s expansions. The latter was introduced at approximately the same time as the Greengard’s first paper on three-dimensional FMM [21], and has found extensive application in tree-codes. Recently, FMM codes based on Cartesian expansions have used recurrence relations to avoid derivatives [22]. Typically, asymptotic computational cost of schemes based on Cartesian expansions is higher. This is because spherical harmonics are optimal for representing harmonic functions in three dimensions. For a system with N mutually interacting particles, it can be readily proven that for a one-level implementation, the computational complexity scales as O(P4N) and O(P6N) for FMMs based on spherical harmonics and Cartesian expansions, respectively. Here, P is the number of harmonics used in the computation. This cost can be reduced by choosing the number of particles at the lowest level in the tree in an optimal manner. Our interest in revisiting these schemes is motivated by the following observations: (i) Taylor’s series expansions provides a natural framework for developing addition theorems [23]; (ii) Taylor’s expansion involves representing the fields in terms of Cartesian tensors; (iii) there is an intimate connection between Cartesian tensors and the spherical harmonics. These connections are well known, and have been explored extensively (as early as Maxwell!); see [24], [25], [26] and references therein. The following statements hold true: (i) components of a traceless tensor of rank n serve as constant coefficients in a spherical harmonic of degree n, and (ii) there is a class of traceless tensors of rank n whose components are n-degree spherical harmonics functions of x,y,z. These connections imply that there should be an intimate connection between the two seemingly disparate methodologies, and one should be able to obtain a similar cost structure for both methods. As we will show, the recurrence relationship that were conjectured for translating multipole expansions [22] can be rigorously derived using traceless tensors. This implies that one should be able to derive a computational scheme using Cartesian tensors that are optimal in the sense of FMM. The method presented herein can be readily generalized to analyze potentials for all ν or, for that matter, to any potential function whose power series converges rapidly [27] (for instance, the Yukawa potential), and more importantly, without the need to use special functions!

Thus, this paper will focus on the use of Cartesian tensors to derive fast computational schemes for all ν. Similar to classical FMMs, the methodologies developed herein will rely on a divide and conquer computational strategy that is facilitated by a hierarchical partitioning of the computational domain through the construction of an oct-tree data structure. The underlying mathematics for two different computational methods will be derived; in the first, operators will be derived for traversing up, down and across the tree. This technique will rely on using totally symmetric tensors. The salient feature of this method is that the traversal up and down the tree (or shifting the origin of the multipole/local expansion) is exact. The second method produces optimal technique, in the sense of FMM, for ν = 1. This optimality is achieved using traceless totally symmetric tensors. For ν  1, it yields the same computational complexity as the first, albeit a different error bound. One of the most interesting features of both algorithms, which separates the proposed technique from [28], is the fact that the operators derived for traversing up and down the tree are independent of ν, and only traversing across the tree depends on ν. This is a powerful feature, as one can use a single traversal up and down the tree to compute the combined effects of different potentials! Finally, the algorithm presented does not involve any explicit (or numerical computation of) derivatives, and quantities are expressed in terms of (products of) traceless (or totally symmetric) tensors.

This paper is organized as follows: in Section 2 we will present relevant details and theorems regarding tensors, detracer, Maxwell Cartesian tensors, and homogeneous polynomials. Here, we also demonstrate the intimate connection between traceless tensors and Legendre polynomials. In Section 3, we present the two FMM-like algorithms, implementation details, and computational cost. In Section 4, we present numerous results to validate the accuracy and to demonstrate the efficiency of the proposed method. Finally, in Section 5, we summarize the contribution of this paper.

Section snippets

Preliminaries

This section introduces basic notation and theorems that will be used in the rest of this paper. The material presented builds upon some of the earlier work by Applequist [26], [29], [30]. For completeness, we have described some of the theorems given in his papers (without proofs) as well as added some of our own (with the necessary proofs).

Divide and conquer strategy

Typically, potentials are evaluated between source and observation pairs that are randomly distributed in a domain ΩR3. The computational scheme developed here will follow those typically used in FMM. To this end, the entire domain is embedded in a fictitious cube that is then divided into eight sub-cubes, and so on. This process continues recursively until the desired level of refinement is reached; an Nl-level scheme implies Nl  1 recursive divisions of the domain. At any level, the

Results

In this section, we will demonstrate the validity of the numerical method presented via numerous examples. The overarching goals of this section are as follows: (i) numerically show that the traversal up and down the tree can be performed exactly, (ii) demonstrate that the proposed method produces accurate results for different values of ν, (iii) demonstrate that this scheme can be used seamlessly for computing potentials that are superposition of potential of the form Rν for multiple values

Summary

In this paper, we have developed two methods for rapidly computing potentials of the form Rν. Both these methods are founded on addition theorems based on Taylor expansions. Taylor’s series has a couple of inherent advantages: (i) it forms a natural framework for developing addition theorem based computational schemes for a range of potentials; (ii) only Cartesian tensors (or products of Cartesian quantities) are used as opposed to special functions. This makes creating a fast scheme possible

Acknowledgements

We are grateful to NSF for support under CCR-0306436 and to the MSU-IRGP program. We are also grateful to the three anonymous reviewers who went through this document with painstaking thoroughness; the paper has improved significantly thanks in large part to their comments.

References (35)

  • A. Appel

    An efficient program for many-body simulations

    SIAM J. Sci. Comput.

    (1985)
  • J. Barnes et al.

    A hierarchical o(n log n) force calculation algorithm

    Nature

    (1986)
  • J.A. Board, W.J. Blamke, D.C. Gray, Z.S. Hukara, W. Elliott, J. Leathrum, Scalable implementations of...
  • L. Greengard et al.

    A fast algorithm for particle simulations

    J. Comput. Phys.

    (1987)
  • L. Greengard

    The Rapid Evaluation of Potential Fields in Particle Systems

    (1988)
  • L. Greengard et al.

    A new version of the fast multipole method for the laplace equation in three dimensions

    Acta Numer.

    (1997)
  • G.L. Xue et al.

    Am. Math. Soc.

    (1996)
  • Cited by (0)

    View full text