The Macroscopic multipole method [36], MMM, employs fast multipole techniques to calculate
electrostatic forces in a system of finite periodic unit cells arranged in a lattice of
cells in
2D or 
in 3D in 
operations.
The method is based on the observation, also realized by Schmidt and Lee [44],
that once the multipole expansion of the unit cell is computed,
translates of this unit cell (i.e. periodic images) are independent of the position of the cell
and hence have the same multipole coefficients as the unit cell.
The algorithm presented here utilizes the Fast Multipole Algorithm mechanisms
[26,9], however it introduces new error criteria to determine the
number of macroscopic images it needs to maintain an appropriate error bound.
In a point charge system, given p (the number of multipole terms) and k (lattice size),
the macroscopic multipole method commences by dividing the simulation cell recursively into
smaller sub-cells as in the FMA method.
The multipole expansion of each sub-cell at the finest level of refinement is calculated
and subcells are grouped into bigger structures up to the unit cell in what is
known as the upward pass. At this point, macroscopic multipole procedures are called to compute
the multipole expansion 
for cell 
, a multiple of the unit cell 
, for 
,
by using multipole to multipole transformations.
The following step starts at the highest level (i=k) and converts the multipole expansion 
to a local expansion about the unit cell center only if the macroscopic region at this level is
``well separated'' from the central unit cell, otherwise, the region is subdivided recursively into
9 cells (2D) or 27 cells (3D) and this step is repeated again. At the deepest level of recursion,
the forces are evaluated directly between cells that are not well separated. Once all macroscopic
cells are dealt with, the downward pass of the FMA can proceed.
In practical simulations, 
and p=8,16 are adequate for average and high accuracy
simulations. This algorithm is also highly efficient incurring only an overhead of about 
over FMA simulation of a single unit cell. The results of this method were within 3-4
significant figures of the Ewald sum results for p=8 and k=4. However, the method is capable of achieving
higher accuracy by increasing p at the expense of longer execution time.
In addition, this method can efficiently handle non-cubic systems, i.e. 
lattice
(
), which allows the study of surfaces, i.e. systems that are finite in one of three dimensions.