Chapter 3. Cell Multipole Method
The most expensive computation in standard molecular mechanics calculations is
the evaluation of the nonbonded energy. Exact computation requires
O(N^2) operations, which is infeasible for large-scale
systems. We must thus find a way to represent groups of atoms in an
approximate way, reducing the operation count while maintaining as much
accuracy as possible. Traditionally, cutoff methods have represented all atoms
sufficiently far away from an atom of interest as not being present at all,
assuming that their effects would be negligible. Unfortunately, for many
interesting systems this assumption may not be valid and accuracy suffers, even
though the operation count is reduced to O(N).
The Cell Multipole Method (CMM)  was developed to overcome these
limitations in handling long-range power-law forces in molecular systems. In
particular, it can be used to handle the R^-1 (or
R^-2 if screened) Coulomb interaction and the
R^-6 attractive portion of the usual Lennard-Jones 12-6 or the
exponential-6 van der Waals potentials. It is a true
O(N)-operation algorithm , with substantially better accuracy
than cutoff methods of the same speed. It is thus the most suitable method for
handling large-scale molecular mechanics problems.
The key feature of the CMM and other multipole methods [3,4,5,6] is that
they replace the effects of atoms with multipole expansions representing the
fields due to those atoms. This replacement reduces the number of computations
that must be performed, while not ignoring these effects altogether as in
cutoff methods. In order to maintain accuracy, stronger, nearby interactions
are represented more accurately than weaker, more distant interactions.
Computing and using these multipole expansions is the major computational task
in all multipole methods. The CMM is a particularly regular, easily
parallelizable multipole method that enables this task to be performed
efficiently. It uses Cartesian coordinates only, unlike the other fast
multipole formulations that use spherical harmonics [5,6,7], further
simplifying the method. We can divide the CMM into four parts:
In the CMM, the system of interest is surrounded by a bounding box, typically a
cube, but in general a parallelepiped specified by three side lengths and three
corner angles. This box, the level 0 cell, is then subdivided into eight
octants by bisecting each side. The eight "child" cells are then themselves
further subdivided into octants. This process continues recursively, forming
an octree decomposition of the original box. The maximum level of the tree is
a parameter of the method and is chosen to best maximize computational speed
and accuracy. Increasing the maximum level will tend to increase computational
speed but decrease accuracy; decreasing the maximum level will have the
opposite effect. The cells at the maximum tree level ("lowest" level) will be
referred to as "leaf" cells. Note that the decomposition is carried out to the
same level across the entire system, unlike other adaptive multipole methods.
This regularity makes it easier to determine the neighbors of any given cell
and may increase the number of timesteps that can be executed before rebuilding
the cell tree.
- Octree decomposition.
- The space occupied by the system of interest is
divided into cells that form a tree. Each cell's effects will be represented
by multipole expansions.
- Multipole expansion computation.
- The multipole expansions are computed for
each cell at each level within the tree, starting from the leaves, or smallest
cells, and working upwards in the tree to the root, which represents the entire
- Taylor series expansion computation.
- The multipoles represent the effects of
atoms within a given cell. What we require, however, is the effect on an atom
within a cell of all the other atoms in the system. We represent the
long-range component of this effect with a Taylor series expansion that applies
to all atoms within a given cell. The coefficients of this series are computed
from the multipoles from the previous step, starting at the root of the tree
and working downwards to the leaves.
- Farfield and nearfield computation.
- Once the Taylor series expansion
coefficients have been computed for each leaf cell, the force on each atom can
be computed as a combination of the effects represented by the Taylor series,
evaluated at the atom's location, and the effects of nearby atoms, which are
Figure 3-1 shows a two-dimensional projection of some of the cells in the
octree. The bounding box has been divided into level 1 cells, one of which has
been further subdivided into level 2 cells. A level 2 cell has in turn been
split into level 3 cells. All of the level 1 and level 2 cells would be
subdivided in this fashion.
Within the leaf cells, the effects of the atoms are replaced by multipole
expansions. The level of multipole expansion (quadrupole, octupole, etc.) is
another free parameter. Increasing the level of expansion will increase
accuracy at the cost of extra computation time. The equations for the cell
charges, dipoles, and quadrupoles are as follows:
Figure 3-1. Two-dimensional representation of octree decomposition.
are labelled with their level.
are the charge, dipole, and quadrupole moments, respectively;
is the position of the atom;
is the position of the center of the expansion; the power-law for the energy is
The centers of the multipole expansions can be chosen in several different
ways. The simplest choice is to use the geometric centers of the cells. If
the atoms within a cell are not distributed evenly, however, expanding about
the centroid (or average) of the atom locations produces improved accuracy for
a given level of multipole expansion.
Another possible alternative is to use a weighted average of the atom
locations, with the weights being, for example, the absolute values of the
charges on the atoms. This tends to place the center close to atoms with large
charges, which will again tend to reduce the higher-order multipole
The centroids or weighted averages can easily be computed in a hierarchical
fashion using the existing cell tree, as the centroid of a higher level cell is
equal to the weighted average of the centroids of its children, where the
weights are the number of atoms in each child cell.
We have found that the increased accuracy from the centroid formulation allows
the use of more highly truncated expansions than would otherwise be required.
Further details are given in chapter 7.
Once the leaf cell multipoles have been computed, the expansions may be
translated to the centers of the next higher-level ("parent") cells and
combined. In this way, the multipole expansion representing the parent cell is
This process continues up the tree until the root (level 0 cell) is reached.
At this point, every cell in the system has associated with it a multipole
expansion representing the field due to all of the atoms contained within the
In order to compute energies and forces on a given atom within the system, we
need to obtain the field due to all of the other atoms in the system. This can
be broken into two components: the "nearfield" interaction due to all of the
atoms in the same cell or neighboring cells, and the "farfield" interaction due
to all of the rest of the atoms. In the CMM, the nearfield interactions are
evaluated explicitly to ensure high accuracy while the farfield interactions
are evaluated using effective fields.
A Taylor series expansion of the farfield for a given cell is used to enable
the computation of the farfield at any point within the cell. The coefficients
of this expansion are determined from the multipole expansions computed in the
We define the Taylor series expansion of the field at the position of an atom
are the zeroth, first, and second order expansion coefficients, respectively.
The centers used for the Taylor series expansion are the same as those used for
the multipole expansion.
If we assume that a cell's parent's Taylor expansion has been computed and
properly represents the field due to all atoms in the system outside the
parent's neighbors, the remaining contribution we need to add to obtain the
cell's farfield expansion is just the fields from the children of the parent's
neighbors and from the other children of the parent, minus the fields from the
cell itself and its immediate neighbors at its level. This is 27 cells (parent
and its neighbors) times 8 (children per cell) minus 27 cells (cell and its
neighbors), or 189 cells. Note that all of the children of the same parent are
immediate neighbors of each other and can thus be omitted. We thus need to add
to the parent's farfield expansion the fields from the parent cell's neighbors'
children (which we will call PNCs), with the understanding that the other
immediate neighbors of the cell in question are also excluded. Figure 3-2
shows a two-dimensional projection in which the immediate neighbors (excluded
from the farfield) and the PNCs (included in the farfield) of a level 3 cell
are shown. All of the shaded cells are with the parent of the level 3 cell or
the parent's neighbor cells. The remaining unshaded level 2 cells are within
the farfield of the parent cell, and are included in its contribution to the
level 3 cell's farfield.
Figure 3-2. Neighbors and PNCs for a level 3 cell.
To start this process, note that all cells at level 1 (the first set of eight
child cells) are immediate neighbors of each other, and thus there is no
farfield contribution from the parent or the PNCs of these cells.
By induction, we can continue generating Taylor expansions all the way to the
leaf cells, at which point every cell has a Taylor expansion representing the
farfield within the cell.
We define the following useful terms:
The contributions from each PNC to a cell's Taylor expansion coefficients are
multipole expansion coefficients are those of the PNC.
The contributions from the parent cell's Taylor series coefficients to one of
its child cells' coefficients are:
here and the
coefficients on the right are from the parent cell.
Once the cell Taylor expansions have been determined, computing the interaction
energy and force due to the farfield at each atom is then simply a matter of
evaluating the Taylor series at the atom's position (given by equation 4 above)
and calculating the interaction of the field and the atomic charge.
Since the farfield changes much more slowly than the nearfield, we have found
it feasible to perform the farfield calculation at intervals instead of every
timestep. The centers and coefficients of the Taylor series expansions
representing the farfield are kept constant during the interval.
The remaining nearfield interactions between an atom and the other atoms in the
same cell and in neighboring cells are computed explicitly, using the
appropriate charge-charge interaction equations (Coulomb or van der Waals).
Next / Previous / Table of Contents