Boxing Scheme

There exist a number of O(N) algorithms for calculating the energy due to a sum of pair-wise short-range interactions [197], all of which involve dividing the system into cubical boxes (or some other space-filling shape [198]) and then interacting the contents of each box with only those in neighbouring boxes. As the number of significant interactions per particle, M becomes O(1) when $M \ll N$, the algorithm scales as O(N). The algorithm used here is the most common in the field -- the standard linked-cell method [199]. The additional complexity of alternative methods is unlikely to offer significant benefit, and the CFMM implementation in the Q-CHEM program, which uses the linked-cell method, has already been highly optimized.

The first step is to determine the number of boxes in each Cartesian direction. The sidelength of each box is taken to be the minimum distance over which two point charges, interacting via the CASE operator is negligible, that is, the value of x when
$$\begin{align}\frac{\erfc( \omega x)}{x} = \mathrm{Desired \; integral \; accuracy}, \end{align}$$
which can easily be found (given $\omega$ and the integral accuracy required) using the Newton-Raphson method [200]. The number of boxes in a direction is then the length of the molecule in that direction, divided by the sidelength.

Following the CFMM, each contracted shell-pair is assigned to a box. If the contracted shell-pair contains primitives with centers in different boxes, then the contracted shell-pair is split into two (or more) shell-pairs, one for each box that the original shell-pair centers straddle.

The box size above is fine for point charges, but the CASE integrals deal with continuous distributions, which can be much larger than the box itself. To solve this problem the CFMM notion of a well-separated index (WS) [157], which determines the distance above which a distribution behaves like a point charge (to within the integral accuracy required) was introduced. The Coulombic interaction between two spherical Gaussian charge distributions can be represented (using the notation of the previous section) as
$$\begin{align}I = R^{-1} \erf \left(R/\sqrt{\frac{1}{\zeta}+\frac{1}{\eta}}\right). \end{align}$$
The $\erf$ factor rapidly approaches 1 as R grows, and the two distributions then interact as classical point charges. Thus, the extent of a distribution is defined as
$$\begin{align}r_{\mathrm{ext}} = \sqrt{\frac{2}{\zeta}} \erf^{-1}(1-\epsilon) \end{align}$$
where $\epsilon$ is the desired precision. A Taylor expansion is made for $\erf$, allowing the extent to take the simple form
$$\begin{align}r_{\mathrm{ext}} = \frac{1}{2} \sqrt{\frac{2}{\zeta} \ln(\epsilon)}. \end{align}$$

In molecular calculations the charge is represented as a product of Gaussians, which contain a prefactor U, dependent on the distance AB between the original shells. This can be included in the definition of extent, allowing for shorter interaction distances. The extent is thus given by
$$\begin{align}r_{\mathrm{ext}} = \frac{1}{2} \sqrt{\frac{2}{\zeta} \ln(\epsilon) - AB^{2}}. \end{align}$$
The WS index is then defined as
$$\begin{align}WS = 2 \left\lceil \frac{r_{\mathrm{ext}}}{l} \right\rceil \end{align}$$
where l is the box size.

The boxes are then looped over, and every shell-pair of a box is interacted with the contents of all boxes within the range WS+1 multiples of the box length. This is exactly the same as the near-field part of the CFMM algorithm.