Communication: Quasi-robust local density fitting

The mean-field or correlated treatment of the electronic Schrödinger equation using the LCAO (linear combination of atomic orbital) approach necessitates the evaluation of four-centre two-electron repulsion integrals

where χ_μ are atomic orbitals, ϕ_p = C_μpχ_μ are molecular orbital, and C_μp are the molecular orbital coefficients. The density fitting approximation was introduced to reduce the computational effort required to build the integrals (rs$|$pq) through a reduced rank formulation.^1,2 The central idea is to fit each charge density χ_μ(r₁)χ_ν(r₁) with an auxiliary basis set of atom-centred density fitting functions ϕ_P(r₁) so that

The coefficients of the approximate charge density $|μν̃$⁠) are determined by minimising the error in the Coulomb energy

This enforces that the error in the integral is bilinear in the fitting error ϵ_J = (ϵ_κλ$|$ϵ_μν) and the approximation is termed robust.³ The optimal coefficients are computed by solving

This formulation has the unsatisfactory feature that although the AOs $|$μν) belong to at most two atoms, the auxiliary functions $|$P) used to fit that charge density are drawn from all of the atoms of the system and the coefficients decay in the worst case only as the reciprocal of the distance between the charge centre of $|$μν) and the fitting function $|$P).⁴ 

Local density fitting formulations circumvent this extraneous long-range behaviour either (a) by replacing Eq. (5), which is a least squares fit to the density using the Coulomb metric, by a least squares fit with a short-ranged metric such as the overlap,^5–7 Yukawa-Coulomb⁴ or attenuated⁸ or Gaussian damped⁹ Coulomb metrics, or (b) by using a domain-based fit that excludes contributions from distant fitting functions a priori.^10–14 In both cases, the error ϵ_J in the integral Eq. (4) is then linearly dependent on the fitting error and, when using standard libraries of atom-centred density fitting functions, the magnitude of linear error makes Eq. (4) impractical for general application.^5,15–18 Local density fitting therefore either follows a pair-domain scheme, with standard fitting in the unified domain of (κλ$|$ and $|$μν), or uses Dunlap’s robust formula for the integral³ that removes the linear error by adding the linear terms

In all cases, the dependence on integrals (Q$|$μν) with (Q$|$ distant from $|$μν) remains, appearing either in the fitting formula or in the integral expression.

Consider again Eq. (6), which determines the fitting coefficients. The functions $|$P) fit the charge density $|$μν) and the functions (Q$|$ act as test functions where the Coulomb potential from $|μν̃$⁠) should be accurate or, more precisely, yield accurate Coulomb energies $(Q|μν̃)$⁠. If we insert an approximate resolution of the identity (RI) $1≈∑PR|PSPR−1R|$⁠, where S_PR is the overlap metric between auxiliary basis functions, we obtain

Here (Rμν) is a three-centre one-electron overlap integral. The decay of the coefficients $CμνP$ follows $SPR−1(Rμν)$ and is therefore exponential in the distance between the charge centre $|$μν) and the function $|$P). The test functions (Q$|$⁠, however, span the whole molecule. Equation (8) suggests $CμνP=SPR−1(Rμν)$⁠, which is understood to be impractical as a direct means for determining $CμνP$ because the incompleteness errors in the approximate RI are substantial, leading to large linear errors in the Coulomb energy.⁶ For functions $|$μ) and $|$ν) located on different atoms, the charge centre $|$μν) lies in the vacuum between the atoms, whereas the fitting functions are all atom-centred and therefore often not well suited to fitting $|$μν). As a consequence, the approximate charge density $|μν̃$⁠) resulting from Eq. (6) acquires unphysical contributions from $|$P) far away from $|$μν) that cancel the linear errors and yield the correct Coulomb energy $(Q|μν̃)=(Q|μν)$⁠. The long-range behaviour of standard density fitting is thus seen to be a consequence of the requirement that the linear error in the integral vanishes entirely when the RI is incomplete. To summarise, the auxiliary functions in Eq. (5) have a dual role: they are used both to approximate a charge distribution $|$μν) and to test the accuracy of the approximation $|μν̃$⁠). It is the second role that is responsible for the non-locality inherent to standard density fitting.

This paper communicates a new fully separable local density fitting formulation where a charge density $|$μν) is fit only using functions $|$P) near by and where it is not necessary to compute the Coulomb interaction between $|$μν) and distant functions (Q$|$⁠. The method is derived by decoupling the two roles of the auxiliary basis functions and by permitting linear error that is made small by the fitting procedure. It will be shown to be stable and almost as accurate as standard robust density fitting without local approximations in practical calculations using standard density fitting basis sets.

The target integrals are computed as

For a given set of μν, a set of fitting functions P_μν ∈ {P} local to $|$μν) are selected according to the criteria

with S_PP = 1. A set of test functions X_μν ∈ {P} are selected according to the criteria

where $f(A,B)=αβα+β|A−B|2$ for two functions with centres A, B and (smallest primitive) exponents α, β. The coefficients $CμνP$ are determined so that the Coulomb potential is accurate at the test points, by minimising the errors in the Coulomb energies between the charge density $|μν̃$⁠) and the test densities $|$X),

The number of test functions X_μν is larger than the number of fitting functions P_μν and the appropriate solution to this system of equations is achieved through QR factorisation of the rectangular matrix (P_μν$|$X_μν).

The criteria for selecting the domain P_μν of fitting functions follows directly from the RI in Eq. (8). The criteria for selecting the test functions X_μν is based on Gauss’s theorem: If the Coulomb potential is correct at functions around $|μν̃$⁠) that do not overlap with $|μν̃$⁠), then it will also be correct at functions further away. In this work, the test functions have been drawn from the density fitting set, but any set could in principle be used. Using delta functions, for example, is equivalent to probing the Coulomb potential on a grid. The current choice of test set is preferable since it retains a direct link to robust density fitting. This can be seen by noticing that Eq. (12) reduces to Eq. (6) when {X_μν} = {P_μν} = {P}.

The above local density fitting scheme has several favourable features: In the limit of vanishing T and infinite R, standard robust density fitting is recovered; the coefficients $CμνP$ are in no way dependent on κλ and the formulation is in this sense fully separable and straightforward to parallelise; the coefficients $CμνP$ can be computed with costs that scale linearly with system size; the only components with quadratic scaling over the lengthscale determined by the slow 1/r decay are the integrals (P$|$Q), which only need to be computed once at the start of a calculation.

The numerical results reported in Sec. III will demonstrate that contrary to the predominant view in the literature, it is possible to use non-robust local formulas and obtain an accuracy comparable to that of standard robust density fitting for Coulomb and exchange contributions to the Fock matrix and for the integrals required for the MP2 energy and that the naturally local formulation of density fitting presented above achieves this. The implications for efficient algorithms for computing Coulomb, exchange, and MP2 integrals are discussed in Sec. IV.

The accuracy of the quasi-robust local density fitting approximation for Coulomb, exchange, and MP2 integrals was first assessed by comparing the resulting Coulomb (J), exchange (K), and MP2 correlation energies to those computed exactly and to those computed using standard robust density fitting. First, a Hartree–Fock calculation was performed using the def2-TZVPP basis set¹⁹ and the RI-JK approximation as implemented in the Turbomole program package.²⁰ This involves standard robust density fitting for both the Coulomb and exchange contributions to the energy and Fock matrix. Weigend’s RI-JK fitting basis sets²¹ and RI-MP2 sets²² were used. Using the converged orbitals and the same basis sets, the J, K, and MP2 correlation energies were computed using quasi-robust local density fitting (QR-J, QR-K, and QR-MP2). The energies were also computed using standard robust density fitting^20,23 and analytically computed 4c-2e electron repulsion integrals. All calculations were performed with a locally modified version of the Turbomole program.²⁴ 

Figure 1 displays the density fitting error in the QR-J, QR-K, and QR-MP2 energies as a function of the overlap threshold T for the alkane chain C₂₄H₅₀ in an idealised linear arrangement. Overlap screening was performed for sets of μν belonging to the same shell pair and a conservative threshold of $R=40a02$ was chosen to select the test functions. Figure 2 displays the corresponding distribution of distances between the fitting functions P_μν and the centre of the charge density $|$μν), summed over all shell pairs μν, and reveals the spatial locality of the fitting set as a function of T.

Each of the J, K, and MP2 energies converges to the standard RI values as the threshold T is reduced. The K and MP2 energies converge rapidly and with a threshold of 1 × 10⁻⁵ or 1 × 10⁻⁶, and local density fitting is as accurate as robust density fitting. The QR-J energy converges more slowly, requiring a tighter threshold of 1 × 10⁻⁷ or 1 × 10⁻⁸. Since QR-DF is non-robust, the QR-K energies are not an upper bound, which can result in a variational collapse if the errors are too large.²⁵ The threshold of 1 × 10⁻⁶ was tight enough to avoid variational collapse in all cases. Linear alkanes represent a worst case scenario for local density fitting⁴ due to the large surface area, which corresponds to regions at the edge of the molecule where the atom-centred RI basis is poorly saturated. Examining the corresponding Fig. 2 reveals that not only are distant fitting functions not required for an accurate fit to the potential but also more than half of the functions nearby can also be safely discarded on the basis of poor overlap with the target charge density.

Table I lists the signed fitting for the J, K, and MP2 energies using the QR and standard RI methods for a sequence of alkanes each in an idealised linear geometry, for a cluster of 21 water molecules taken from the Cambridge Cluster Database,²⁶ for Cu$(C3H6O)4+$⁠, the largest Cu⁺ complex from the CuCNDe33 dataset,²⁷ for the methotrexate molecule (C₂₀N₈O₅H₂₂) and the 4T-MeOH cluster model for the adsorption of methanol on the H-ZSM-5 zeolite, taken from Ref. 28. The def2-TZVPP basis set combinations was used in all cases and fitting errors are reported in micro-Hartree per occupied orbital for J and K and per valence orbital for MP2. Conservative thresholds T = 1 × 10⁻⁶ for exchange and MP2 correlation and T = 1 × 10⁻⁸ for the Coulomb energy were used. The MP2 values for C₃₂H₆₆ and methotrexate are missing because the calculation without density fitting required more memory than was available on the computers used for these calculations. The accuracy of the local density fitting method matches that of standard density fitting for the progression of alkane chains with increasing length, and the density fitting error per occupied orbital is constant for both methods. For the more globular structures, the accuracy of the local fitting method is equally impressive.

Table II reports the differences between the local density fitted and robust density fitted isomerisation energies for a subset of the ISOL (large isomerisation) test set of molecules,²⁹ which are between 24 and 42 atoms in size and large enough for the local approximations to take effect. The same TZVPP basis set combination and thresholds T, R were used for these calculations as those in Table I, but here the isomerisation energies were computed using the converged self-consistent Hartree–Fock energies, employing density fitting for both J and K. Correspondingly, the RI-MP2 energy was computed using orbitals converged using RI-JK and the QR-MP2 energy was computed using the orbitals converged using the QR-JK method. No variational collapse was observed in any of the HF calculations performed. The values in Table II show that the converged HF energies obtained by local density fitting are as accurate as the energies obtained by standard density fitting, with an RMS error compared to HF isomerisation energies computed using analytic 4c-2e integrals of 0.07 kJ/mol compared to 0.06 kJ/mol for standard fitting.

The QR-MP2 values are less accurate. Since Table I demonstrates that the QR-MP2 energy is accurate when the correct orbitals and eigenvalues are provided, the slightly poorer MP2 energies in Table II most likely result from the deterioration in the virtual orbitals and eigenvalues due to the local fitting at the SCF level. Indeed, the cases where a large MP2 deviation is observed are also the cases where the orbital relaxation from RI-HF to QR-HF is significant. Nonetheless, the maximum local density fitting error in the computed QR-MP2 isomerisation energies is only 2 kJ/mol.

The numerical results show that accurate HF and MP2 energies can be obtained using the quasi-robust local density fitting method. From a formal perspective, this approach is computationally attractive since the coefficients $CμνP$ are zero unless $|$μ), $|$ν), and $|$P) have non-negligible overlap and the number of $CμνP$ is therefore small and grows linearly with system size. Similarly the only (Q$|$μν) integrals that are required are those where (Q$|$ and $|$μν) are connected by overlap criteria, which also results in an early onset of linear scaling with system size for the 3c-2e integral evaluation. The cost of solving the fitting equations for each block of μν is independent of system size but carries an increased prefactor since QR factorisation (or divide and conquer single value decomposition) must be performed in the place of Cholesky decomposition.

In the context of efficient linear-scaling algorithms for computing HF exchange, quasi-robust density fitting is expected to be both more accurate and more efficient than ARI (atomic RI), which although restricts the fitting functions $|$P) to those near to $|$μν) uses the whole set of test functions (Q$|$ and requires the same 3c-2e evaluation as the full RI.^11,12 Since QR-K is as accurate as RI-K, it is more accurate than the COSX pseudospectral method,³⁰ but to properly compare the relative costs efficient QR-K implementations must be designed and tested, for example, one that combines quasi-robust local density fitting with the LinK density screened exchange algorithm^31,32 and uses a low-scaling proxy for Eq. (10).

Replacing standard density fitting with quasi-robust density fitting for the evaluation of the HF Coulomb contribution does not change the scaling since the cubic cost of evaluating the fitting coefficients can be avoided through early contraction with the density matrix.³³ The overall scaling for the Coulomb term is in both cases determined by the number of required (κλ$|$μν), which is quadratic over the lengthscale of non-negligible Coulomb interactions. Any computational saving of using the QR fitting in place of standard fitting would come from replacing the Cholesky decomposition of the full (Q$|$P) matrix, at $O(N3)$ cost with $O(N)$ QR factorisations of small matrices. The size regime where this becomes advantageous is also the regime where it is likely to be more efficient to use the existing efficient linear-scaling methods based on multipole moments for far-field contributions.^34,35

Quasi-robust local density fitting brings distinct advantages in the context of local correlation methods. Most implementations of PAO (projected atomic orbital) and PNO (pair natural orbital) local correlation methods^14,36,37 use pair-domains for local density fitting, which requires integrals 3c-2e integrals (Q_μi$|$μi) and (Q_νj$|$μi) for all pairs of localised occupied orbitals i, j that have significant pair energies. In the quasi-robust formulation, the density fitting can be performed for i independently of j and the Coulomb integral formula does not require the cross terms (Q_νj$|$μi). Therefore, once the fitting coefficients have been computed, constructing the integrals is significantly simpler and cheaper. Moreover, the simplicity of the indexing required for the integral formula makes it significantly more straightforward to cast local correlation methods as a sequence of sparse tensor operations so that the state-of-the-art sparse tensor libraries suitable can be used to exploit modern massively parallel computer architectures.

D.P.T. thanks the Royal Society for financial support. This article is dedicated to the memory of Professor Reinhart Ahlrichs.