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
.
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:
Container
- type of particle container; pulled from
the Interface
template argument
Calculator( size_t order, Container& particles)
-
constructor; order
is the order of multipole expansion.
void set_max_number_per_box( size_t n)
- sets
the maximum number of particles in a bottom-level box. The
default is 20.
size_t max_number_per_box()
void set_distance_ratio( double)
- Two boxes
are considered near neighbors if the ratio of the
maximum dimension of Box 1 or Box 2 (whichever is larger),
and the center-to-center distance, is larger than a certain
value. This value can be changed by this function; the
default is 0.5.
double distance_ratio()
void set_padding( double)
- if there is no
padding, then the outermost box has its bounds set at the
extremes of particle positions. The padding value increases
the box size by a certain amount in each direction. The
default value is 0.0. Adding padding might be a good idea
if you think your particles might drift outside the box.
void setup_boxes()
- sets up the containing boxes and
precomputes many of the multipole coefficients. Must be done
at the beginning of a simulation, and again only if the distribution
of the particles changes significantly. You probably don't want
to call this at each time step of a simulation, since it can
take significantly more time than the other functions described below.
The three
setting functions above, set_padding
, set_max_number_per_box
,
and set_distance_ratio
, must be called before this
function if the user doesn't want to use the default values.
void compute_expansion()
- performs the computation
steps: the upward combining of multipoles, the translation to
far boxes, and the downward combining of Taylor expansions.
Must be done each time the particles change positions or
interaction strength (e.g., charge). Basically, this computes
the far-field stuff.
void compute_interactions_with_potential()
-
finishes the pairwise interactions using the value of
the potential field. Loops over the particles in each box
and near-neighbor boxes, and includes the far field computed
in compute_expansion
. Makes use of the engine function
apply_potential
(see below). The
compute_expansion
function must be called beforehand
to compute the far-field stuff.
void compute_interactions_with_field()
-
like previous, but uses the gradient of the
field. Makes use of the engine function
apply_field
(see below).
reload_contents()
- if the particles have
changed position enough that many of them are no longer in
their original boxes, then this function redistributes them.
It does not change the boxes themselves. In an MD simulation,
it is probably a good idea to call this function frequently,
since it doesn't appear to consume much time.
double potential_at_point( const double* r)
-
given a point position (x,y,z) in r, returns the value of
the field at the point, whether inside or outside the
outer box. Only valid once compute_expansion
has
been called.
void get_field_at_point( const double* r, double* F)
-
just like the previous function, except it put the gradient
of the field at that point into F
.
typedef Container
- type of particle container. This
actually contains the particles. At no place does the multipole
code ever address either individual particles or their references,
so this gives lots of flexibility.
typedef Refs
- Each bottom-level box maintains references
to the particles, and Refs
is the type of the container.
For example, it could be an array of pointers, or it could be
an array of integers that address the Container object
.
size_t size( const Container& container, const Refs& refs)
-
number of contained references. Note that in this function and
several others below, the reference container refs
is
always accompanied by the original particle container; some implementations
might require it (e.g., using integers as references would require
a starting memory address to address a particle).
void copy_references( Container& container, Refs& refs)
-
copies references from container
to refs
void copy_references( Container& container, Refs& refs0, Refs& refs1)
-
copies references from refs0
to refs1
.
void empty_references( Container& container, Refs& refs)
- empties
refs
void split_container( Container& cont, Refs& refs,
Refs* const refs_cube[][2][2]
)
- splits refs
and distributes the particles into the eight subcubes of the current cube.
Pointers to eight Refs
data structures are provided by the
multipole code and presented to the user in a 2x2x2 array. If 0 and 1
are indices describing coordinates of the eight cubes, then the
first index of refs_cube
is the x-coordinate.
void get_bounds( const Container& cont,
double low[3], double high[3])
-
places the lowest and
highest x,y, and z coordinate of refs
into
low
and high
.
void process_bodies( Container& cont, Refs& refs)
-
directly computes the interactions among all particles
referenced in refs
. This and the following function are
used to compute the near-field interactions.
void process_bodies( Container& cont, Refs& refs0, Refs& refs1)
-
directly computes the interactions between the particles in
refs0
and those in refs1
. Does not compute interactions
among particles in the same container.
template< class Function> void apply_potential( Container& cont,
Refs& refs, Function& f)
- the function object f
, which
is provided by and
passed to this function by the multipole code, has
overloaded operator double operator()( const double* r)
.
This function object f
returns the field value at point r
.
The apply_potential
function must loop through each particle in refs
,
get its position, feed the position to f
, get the
field value, and do something with it (e.g., multiply it by the
particle charge and add to the total energy). This is used to
compute the far-field interactions on the individual particles.
template< class Function> void apply_field( Container& cont,
Refs& refs, Function& f)
- like the previous function, except
that f
has void operator( const double* r, double* F)
.
This takes particle position r
and places the field
gradient at that point into F
, which is provided by the
user.
For example, this function would be
used to compute forces in a Coulombic system by getting F
for
each particle, multiplying it by the charge and -1 to get the far-field
contribution to the force on that particle.
template< class Function> void apply_multipole_function(
const Container& cont, const Refs& refs, Function& f
)
- this function is used to get multipoles from individual
particles in refs
. The function object f
has
the () operator overloaded as
void operator()( size_t order, const double* r, const double* mpole)
. The apply_multipole_function
must loop through each particle and call the function object f
with
the particle's multipole order in order
,
the position in r
, and its multipole moments (more below) in mpole
.
For example, if we were
doing a Coulombic calculation with point charges, order
would be 0, and the particle charge would be placed in mpole[0]
.
If the particles were dipoles, order
would be 1,
the charge would be placed
into mpole[0]
, and the x,y and z components of the
dipole would be placed
into mpole[1]
, mpole[2]
and mpole[3]
.
See the note below on how the arguments to mpole
are defined.
void get_r_derivatives( size_t order, double r, double* derivs)
-
This places the ith derivative of the spherically symmetric
potential (or Green's function) at r
into derivs[i]
.
Derivatives up to and including order
must be included.
double near_potential( const Container& cont, const Refs& refs,
const double* r)
- given a position r
, computes the
field value at that point from the particles in refs
.
Required only if the potential_at_point
function is
being used.
double get_near_field( const Container& cont, const Refs& refs,
const double* r, double* F)
- like the previous, except
computes the field gradient at
and puts it into
F
.
Required only if the get_field_at_point
function is
being used.
size_t n_sources( const Container& cont, const Refs& refs)
- number of source particles in refs
.
If the user is not
distinguishing between sources and receivers, this should return the same number
as size
.
size_t n_receivers( const Container& cont, const Refs& refs)
- number of receiver particles in refs
.
If the user is not
distinguishing between sources and receivers, this should return the same number
as size
.
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.
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.