K.-T. Lim Thesis - Chapter 3

Chapter 3. Cell Multipole Method

3.1. Introduction

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) [1] 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 [2], 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:

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 system.
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 calculated explicitly.

3.2. Octree Decomposition

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.

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.

[IMAGE]
Figure 3-1. Two-dimensional representation of octree decomposition.
Cells are labelled with their level.

3.3. Multipole Expansions

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:

$q_cell = \sum_{atoms}{q_atom}; \mu_{cell,a} = \sum_{atoms}{e R_a q_atom}; Q_{cell,ab} = \sum_{atoms}{e (e+2) R_a R_b q_atom - e R^2 q_atom \delta_ab}$ (1)

where $q, \mu_a,$ and Q_ab are the charge, dipole, and quadrupole moments, respectively; $R_a = R_{atom,a} - R_{center,a}$ ; R_atom is the position of the atom; R_center is the position of the center of the expansion; the power-law for the energy is $E = k R^{-e}$ ; and a,b = x, y, z .

3.4. Choice of Expansion Centers

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.

$R_centroid = 1/n_{atoms-in-cell} \sum {R_atom}$ (2)

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 coefficients.

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.

3.5. Combination of Multipoles

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 determined.

$q_parent = \sum_children {q_child}; \mu_{parent,a} = \sum_children {\mu_{child,a} + e R_a q_child}; Q_{parent,ab} = \sum_children {Q_{child,ab} + (e+2) R_a R_b q_child - (2 R.\mu + e R^2 q_child) \delta_ab}$ (3)

where $R = R_{center,parent} - R_{center,child}$ .

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 cell.

3.6. Taylor Series Expansions

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 previous step.

We define the Taylor series expansion of the field at the position of an atom to be

$V(R) = V_0 + \sum_a {V_a R_a} + \sum_ab {V_ab R_a R_b} + ...$ (4)

where R = R_atom - R_center and V_0, V_a, and V_ab 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.

3.7. Computation of Taylor Series Coefficients

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.

[IMAGE]
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:

$R = R_{center,PNC} - R_{center,cell}; \mu = \sum_a {\mu_a R_a}; Q = \sum_ab {Q_ab R_a R_b}$ (5)

The contributions from each PNC to a cell's Taylor expansion coefficients are then:

$V_0 = R^-e q - R^{-e-2} \mu + 1/2 R{-e-4} Q$ (6)

$V_a = eq R^{-e-2} R_a - R^{-e-2} [(e+2) R^-2 R_a \mu - \mu_a] + R^{-e-4} [1/2 (e+4) R^-2 R_a Q - \sum_b {Q_ab R_b}]$ (7)

$V_aa = 1/2 eq R^{-e-2} [(e+2) R^-2 R_a^2 - 1] + R^{-e-4} (e+2) \mu_a R_a - 1/2 (e+2) \mu R^{-e-4} [(e+4) R^-2 R_a^2 - 1] + 1/2 R^{-e-4} Q_aa + 1/4 (e+4) R^{-e-6} \sum_b {Q_ab R_a R_b}$ (8)

$V_ab = e(e+2)q R^{-e-4} R_a R_b + (e+2) R^{-e-4} [\mu_a R_b + \mu_b R_a - (e+4) R^-2 \mu R_a R_b] + R^{-e-4} Q_ab - (e+4) R^{-e-6} \sum_g {Q_ag R_b R_g + Q_gb R_a R_g} + 1/2 (e+4)(e+6) Q R^{-e-8} R_a R_b$ (9)

where the $q, \mu, and Q$ 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:

$V_{child,0} = V_0 + \sum_a {V_a R_a} + \sum_ab {V_ab R_a R_b}; V_{child,a} = V_a + V_aa R_a + \sum_b {V_ab R_b}; V_{child,ab} = V_ab$ (10)

where $R = R_{center,child} - R_{center,parent}$ here and the coefficients on the right are from the parent cell.

3.8. Farfield Evaluation and Nearfield Computation

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

Kian-Tat Lim