Generic Cartesian-Based Adaptive Multipole Calculator

Multipole methods allow the rapid computation of pairwise interactions among a large number of particles. If there are N particles, then the time required is proportional to N, whereas the simple method of doing a double loop and computing the interaction of each pair requires time proportional to N2.

This particular code has the following features. First of all, it uses a Cartesian formulation of the multipole method, rather than spherical harmonics. Rather than using explicit formulas, it uses a sparse matrix approach to automatically generate expansion terms, using the method of Visscher and Apalkov (J. of Magnetism and Magnetic Materials 322(2) 275-281 (2010)) This allows us to specify any order of expansion at the time of computation. This approach also allows the code to use any spherically symmetric potential of interaction. So, in addition to the usual Coulombic interaction, one can just as easily use a screened Coulombic, or any other. The partitioning of space in this code is adaptive; the particles are bounded by a cubes which are recursively subdivided into octrees. Particles can also be divided into "sources" and "receivers"; the sources provide the field and the field is evaluated at the receivers. (The sources and receivers can be the same particles.) Providing this information can help the software optimize the division of space into cubic cells. Finally, the code itself is generic; it can be used with any container of particles; one just has to provide the necessary types and functions to an Interface class, as discussed in previous progress reports.

In tests I've run with Coulombic particles evenly distributed in a box, using order 5 gives forces with relative errors on the average less than 0.01 percent. At this order, the break-even point, above which the multipole method is faster than the simple N2 double loop, is a bit less than 1000 particles.

My motivation for writing this is the fact that multipole methods are needed if current methods for doing protein-protein Brownian dynamics are to scale to larger things like ribosomes and other organelles. For example, boundary element methods would need to use Coulombic and screened Coulombic fields. Also, in computing desolvation forces in my Brownian dynamics code, I need to compute interactions that scale as 1/r4. So, it's nice that one piece of code can handle all three types of interactions.

This code is found in two files in browndye/src. The header file cartesian_multipole.hh is included in the user code, and it is linked to the object code of cartesian_multipole.cc.

Class and Member Functions

The templated C++ class used for this is the Calculator class, found in namespace Cartesian_Multipole. The Calculator takes an Interface class as its template argument. It has the following member functions and types:

Interface Class

The Interface class must have the following types and functions defined:

Further Comments

As can be seen, all individual particles are accessed only by the user-supplied code. This allows the use of containers with more that one type of objects in each container. For example, boundary element calculations use both boundary elements and permanent charges, but the multipole code does not need to know that it is dealing with two different kinds of objects.

Regarding the mpole argument to the function object given to apply_multipole_function, the multipoles are laid out in a one-dimensional array, even though they represent higher-order tensors. The multipole (i,j,k) is defined simply as the integral over the charge distribution of the Green's function f times (xi yj zk). This is different than the conventional definition of a quadrupole, for example. The three indices are arranged in order by the following generation algorithm: 

i = previous_i; // order among indices 
for( int kz = 0; kz <= order, kz++)
  for( int ky = kz; ky <= order; ky++){
    indices[i][0] = order - ky; // x power
    indices[i][1] = ky - kz;    // y power
    indices[i][2] = kz;         // z power
    ++i;
}
It can be seen that the 3 indices always sum to order. This double loop is run for all orders, from 0 up to the actual order of the computation, and indices of the next higher order are placed right after the previous one. For example, the six quadrupole indices would be (2,0,0), (1,1,0), (0,2,0), (1,0,1), (0,1,1), (0,0,2), and they would be placed in slots 4 through 9, because the 0 slot has the 0th-order pole, and slots 1 through 3 have the dipoles.

An actual simulation loop would look something like this: 

size_t order = 5;
Particles particles;
size_t nsteps = ...;

// initialize the particles
...

Cartesian_Multipole::Calculator< Interface> ms( order, particles);
ms.set_max_number_per_cell( 10);
ms.setup_cells();

for( size_t istep = 0; istep < nsteps; istep++){
  ms.compute_expansion();
  ms.compute_interactions_with_field(); // get forces

  // do other stuff for the simulation, including updating particle
  // positions
  ...
  ms.reload_contents();
}
It is the loop in the code above that is the basis for my determination of the breakeven point.

Future Work

I would like to include code to compute interactions between two such multipole boxes that have been rotated and translated relative to each other. To do that, I need to figure out how to efficiently apply rotation operators to these multipoles. Fortunately, this isn't immediately necessary for what I'm working on, but it would be nice further down the road.

It would also be nice to further genericize this code to allow people to add their own methods for doing multipole expansions. There are some very fast and sophisticated algorithms based on spherical harmonics and plane waves, and allowing others to plug these into a generic framework could be very useful.