The “Bubblepole” (BUPO) Method for Linear-Scaling Coulomb Matrix Construction with or without Density Fitting

In this work, we describe the development of a new algorithm for the computation of Coulomb-type matrices using the well-known resolution of the identity (RI) or density fitting (DF) approximation. The method is linear-scaling with respect to system size and computationally highly efficient. For small molecules, it performs almost as well as the Split-RI-J algorithm (which might be the most efficient RI-J implementation to date), while outperforming it for larger systems with about 300 or more atoms. The method achieves linear scaling through multipole approximations and a hierarchical treatment of multipoles. However, unlike in the fast multipole method (FMM), the algorithm does not use a hierarchical boxing algorithm. Rather, close-lying objects like auxiliary basis shells and basis set shell pairs are grouped together in spheres that enclose the set of objects completely, which includes a new definition of the shell-pair extent that defines a real-space radius outside of which a given shell pair can be safely assumed to be negligible. We refer to these spheres as "bubbles" and therefore refer to the algorithm as the "Bubblepole" (BUPO) algorithm, with the acronym being RI-BUPO-J. The bubbles are constructed in a way to contain a nearly constant number of objects such that a very even workload arises. The hierarchical bubble structure adapts itself to the molecular topology and geometry. For any target object (shell pair or auxiliary shell), one might envision that the bubbles "carve" out what might be referred to as a "far-field surface". Using the default settings determined in this work, we demonstrate that the algorithm reaches submicro-Eh and even nano-Eh accuracy in the total Coulomb energy for systems as large as 700 atoms and 7000 basis functions. The largest calculations performed (the crambin protein solvated by 500 explicit water molecules in a triple-ζ basis) featured more than 2000 atoms and more than 33,000 basis functions.