A density-fitting implementation of the density-based basis-set correction method

This work reports an efficient density-fitting implementation of the density-based basis-set correction (DBBSC) method in the MOLPRO software. This method consists in correcting the energy calculated by a wave-function method with a given basis set by an adapted basis-set correction density functional incorporating the short-range electron correlation effects missing in the basis set, resulting in an accelerated convergence to the complete-basis-set limit. Different basis-set correction density-functional approximations are explored and the complementary-auxiliary-basis-set single-excitation correction is added. The method is tested on a benchmark set of reaction energies at the second-order Møller – Plesset (MP2) level and a comparison with the explicitly correlated MP2-F12 method is provided. The results show that the DBBSC method greatly accelerates the basis convergence of MP2 reaction energies, without reaching the accuracy of the MP2-F12 method but with a lower computational cost.


| INTRODUCTION
One of the main goals of quantum chemistry is the accurate prediction of molecular properties, which requires to tackle the electron correlation problem.For this, there are two main families of computational electronic-structure methods: wave-function theory (WFT) 1 which targets the complicated many-electron wave function, and densityfunctional theory (DFT) 2 which uses the simpler one-electron density.
While DFT has become the workhorse of quantum chemistry thanks to its appealing balance between computational cost and accuracy, the lack of a systematic scheme to improve the quality of densityfunctional approximations has renewed the interest in the development of WFT methods in the last few decades.
A serious limitation of WFT methods is their slow convergence of the correlation energy with the size of the one-electron basis set.This slow convergence originates from the short-range singularity of the Coulomb electron-electron repulsion which induces a derivative discontinuity in the exact eigenstate wave functions, known as the electron-electron cusp condition. 3There are two main approaches for dealing with this problem.The first approach consists in extrapolating the results to the complete-basis-set (CBS) limit by using increasingly large basis sets. 4,5e second approach consists in using explicitly correlated R12 or F12 methods which incorporate in the wave function a correlation factor reproducing the electron-electron cusp (see, e.g., References 6-16).
An alternative approach to accelerate basis-set convergence was recently proposed, which we will refer as the density-based basis-set correction (DBBSC) method. 17It consists in correcting the energy calculated by a WFT method with a given basis set by an adapted basis-set correction density functional incorporating the short-range electron correlation effects missing in the basis set, resulting in an accelerated convergence to the CBS limit.In practice, this basis-set correction density functional is constructed from range-separated DFT 18 by defining a basis-dependent local range-separation parameter which provides a local measure of the incompleteness of the basis set.This DBBSC method was validated for configuration-interaction and coupled-cluster calculations of atomization energies, [19][20][21] excitation energies, 22 dissociation energy curves, 23 and dipole moments. 24,25It was also extended to GW calculations 26 and to linear-response theory, 27 and some mathematical aspects of the method were studied in detail on a one-dimensional model system. 28 this work, we report an efficient implementation of the DBBSC method in the MOLPRO software [29][30][31] in which density fitting 32 is used to alleviate the computational bottleneck of the method, namely the calculation of the local range-separation parameter.This allows us to use the DBBSC method on larger molecular systems than what was previously possible.We thus apply the DBBSC method for correcting the basis-set errors in the molecular reaction energies of the FH51 benchmark set 33,34 at the second-order Møller-Plesset (MP2) level.
We also test different basis-set correction density-functional approximations, as well as the addition of a single-excitation correction for one-electron basis-set errors.Finally, we compare the performance of the DBBSC method with the explicitly correlated MP2-F12 method. 11e paper is organized as follows.In Section 2, we explain the theory of the present implementation of the DBBSC method.
Section 3 provides computational details for the calculations on the FH51 benchmark set.In Section 4, we give and discuss our results.
Finally, Section 5 contains our conclusions.

| THEORY
For simplicity, we give the equations for closed-shell states and we assume real-valued HF spatial orbitals fφ p g.

| The DBBSC method at the MP2 level
Given the MP2 total energy E B MP2 in a basis set B, we apply the nonself-consistent basis-set correction 17,19 as where E B ½n B HF is the basis-correction density functional evaluated at the active HF density n B HF (i.e., excluding core orbitals in case of frozen-core calculations).In order not to affect the CBS limit, this functional E B ½n must be such that it vanishes when the basis set B is complete.Moreover, provided a good enough approximation is used for E B ½n, the basis-set corrected MP2 energy, referred to as "MP2+DFT," is expected to converge faster to the MP2 CBS limit than the uncorrected MP2 energy.

| Local range-separation parameter
The dependence on the basis set of the basis-correction density functional E B ½n comes from the local range-separation parameter μ B ðrÞ.It is defined as 17,19 μ B ðrÞ ¼ where W B ðrÞ is the on-top value of the effective interaction localized with the HF wave function In Equation (3), n B 2,HF ðrÞ is the HF on-top pair density with the active HF density n B HF ðrÞ ¼ 2 P act i φ i ðrÞ 2 , and f B HF ðrÞ has the expression where p and q run over all (occupied + virtual) HF spatial orbitals, i and j run over active HF spatial orbitals, and ðφ p φ i jφ q φ j Þ are the twoelectron Coulomb integrals in chemists' notation.We recall that by active orbitals we mean occupied orbitals without the frozen-core orbitals, in case of frozen-core calculations.
The local range-separation parameter μ B ðrÞ provides a local measure of the incompleteness of the basis set.A straightforward calculation of f B HF ðrÞ in Equation ( 5) requires to first calculating the molecular-orbital two-electron integrals ðφ p φ i jφ q φ j Þ with a dominant scaling of OðN act N 4 all Þ, and then performing the sums at each grid point which scales as OðN 2 act N 2 all N grid Þ, where N act is the number of active orbitals, N all is the total number of orbitals in the basis, and N grid is the number of spatial grid points.This is the computational bottleneck of the basis-set correction calculation.This scaling can be reduced by density fitting. 32,35Introducing an auxiliary fitting basis set fχ A g, the orbital product is approximated as where d pi A are the Coulomb-fitting coefficients with and Orthonormalizing the auxiliary fitting basis functions with respect to the metric J, we can approximate the two-electron integrals as and the quantity f B HF ðrÞ in Equation ( 5) as Thus, with density fitting, there is no need to build explicitly the two-electron integrals anymore and the calculation of f B HF ðrÞ in Equation ( 12) now scales as OðN act N all N fit N grid Þ where N fit is the number of auxiliary fitting basis functions.In practice, the same auxiliary fitting basis sets optimized for density fitting in MP2 can be used here.

| Approximate basis-correction density functional
We approximate the basis-correction density functional with the local form 19 E B ½n ≈ ð e sr c,md ðnðrÞ, rnðrÞ, μ B ðrÞÞdr, ð13Þ where e sr c,md ðn, rn, μÞ is the complementary multideterminant shortrange correlation functional energy density 19,36 e sr c,md ðn, rn, μ B Þ ¼ e c ðn, rnÞ ÞÞ=3 and n 2 ðnÞ is a model of the on-top pair density.In Equation ( 14), e c ðn, rnÞ is a standard Kohn-Sham correlation functional energy density.As in previous works, the default choice is the PBE correlation functional. 37In this work, we also test using the LDA, 38 LYP, 39 TPSS, 40 and SCAN 41 correlation functionals.
Note that the TPSS and SCAN functionals are meta-GGA functionals, that is, they depend also on the non-interacting positive kinetic energy density τðrÞ ¼ ð1=2Þ P act i jrφðrÞj 2 , and thus constitute a slight extension of Equations ( 13) and ( 14).
The default choice 19 for n 2 ðnÞ is to use the on-top pair density of the uniform-electron gas (UEG) where the on-top pair-distribution function g 0 ðnÞ is parametrized in Eq. ( 46) of Reference 42.In this work, we also explore two other on-top pair-density models.[45] n CS 2 ðnÞ ¼ where and where q is an empirical parameter.The second one is the Hollett-Pegoretti (HP) model 46 n HP 2 ðnÞ ¼ where We may choose the value of the parameter q, for example, by impos- Hylleraas-type wave function, [47][48][49] we find q ¼ 1:88 for the CS model and q ¼ 2:05 for the HP model.When these on-top pair-density models are used with the PBE correlation functional in Equation ( 14), we call the resulting basis-set correction functionals PBE-CS and PBE-HP, respectively.

| CABS single-excitation correction
For small basis sets B, the HF energy can have a substantial basis-set error.This HF basis-set error is not corrected by the approximate basis-set correction functionals in Section 2.3 since they only correct for missing short-range correlation.The HF basis-set error can however be easily corrected by using the complementary auxiliary basis set (CABS) (see, e.g., References 10,12,50-54) used in explicitly correlated R12/F12 methods.In this approach, a large orthonormal basis set is formed by the occupied+virtual HF orbitals obtained in the normal basis set B and an additional set of virtual orbitals obtained from the CABS.The HF energy correction due to the addition of the CABS is estimated by second-order perturbation theory, leading to the expression, in a closed-shell formalism, 12,50,51 where i runs over active HF orbitals and α runs over all virtual orbitals (obtained in the normal basis set B and from the CABS).In Equation ( 21), f α i are Fock matrix elements and t i α are single-excitation coefficients found by solving the first-order perturbation equations The correction is often referred to as the CABS single-excitation correction.Note that a similar correction is used in the dual basis-set approach proposed by Wolinski and Pulay 55 for improving HF energies and by Liang and Head-Gordon 56 for Kohn-Sham DFT energies.
The total basis-set corrected MP2 energy is thus and will be referred to as "MP2+CABS+DFT".For comparison, we will also present MP2 results only corrected by the CABS singleexcitation correction, which will referred to as "MP2+CABS."

| COMPUTATIONAL DETAILS
0][31] We have performed tests on the FH51 benchmark set.The FH51 set 33,34 is a set of 51 reaction energies for various organic molecules.It is included in the GMTKN55 database. 57e FH51 set contains a large variety of molecules of different sizes (from 2 to 29 atoms).It is thus suitable to test the DBBSC method over systems of different sizes.Moreover, many systems are large enough so that the present density-fitting implementation has a large impact on the performance of the method.As regards the basis set B, we use the aug-cc-pVnZ basis sets 58 for first-row atoms and the aug-cc-pV(n+d)Z basis sets 59 for second-row atoms, which we jointly abbreviate as avnz, for n ¼ 2 (d), 3 (t), 4 (q), and 5.
We perform canonical-orbital density-fitting HF 60 and densityfitting MP2 32 calculations with the frozen-core approximation.We calculate the basis-set correction with different functionals evaluated at the active HF density, and including the CABS single-excitation correction. 12,50,51The basis-set correction is consistently calculated in the frozen-core approximation, corresponding to using only active orbitals in Equation ( 5) and in the HF density used in Equation (1).For the n ¼ 2, 3, for comparison, we also perform canonical-orbital density-fitting MP2-F12 (in the default 3C(F) variant with a Slatertype geminal exponent of 1 bohr À1 ) 11 calculations, implicitly including the CABS single-excitation correction.The MP2/CBS reference values are estimated as the MP2-F12 values with the av5z basis set.
For a given basis set B, the density-fitting basis sets used are the corresponding B/JKFIT and B/MP2FIT basis sets of Weigend et al. 61,62 (and their extensions 51 ) for the HF and MP2 calculations, respectively.The B/JKFIT basis set is also used as CABS for the CABS single-excitation correction.We have checked the density-fitting errors and found them to be insignificant.For large systems, densityfitting calculations of the basis-set correction can be more than an order of magnitude faster than non-density-fitting calculations (computation times can be found in the Appendix S1).

| RESULTS AND DISCUSSION
As a first test, we compare in Figure 1  In Table 1, we report the mean absolute errors (MAEs) on the reaction energies of the FH51 set with respect to MP2/CBS obtained with the methods already discussed, as well as with additional basisset correction functionals, namely LDA, LYP, TPSS, SCAN, PBE-CS (q ¼ 1:88), and PBE-HP (q ¼ 2:05).For the methods already discussed, the mean errors are consistent with the observations made previously.
For the avdz basis set, we go from a MAE of 2.08 kcal/mol for uncorrected MP2 to a MAE of 0.62 kcal/mol for MP2+CABS+PBE and a MAE of 0.33 kcal/mol for MP2-F12.For the avtz basis set, we go from a MAE of kcal/mol for uncorrected MP2 to a MAE of 0.21  With the other basis-set correction functionals tested, the MAEs are very similar, except for the LYP correlation functional which gives much larger basis errors.We have also tested optimizing the parameter q in the CS and HP on-top pair density-density models in Equation (18) and the parameter c in front of the on-top pair density in Equation ( 14), but we did not obtain significant improvements.Thus, if we set aside LYP, we find a rather small sensitivity of method to the underlying correlation functional for calculating reaction energies.
In the Appendix S1, we report additional statistical indicators which confirm this conclusion.
Finally, as regards the computational cost of the DFT-based basis-set correction in comparison with MP2-F12, we consistently observe, for all basis sets, that MP2+CABS+PBE is approximately 10 times faster than MP2-F12 in the default 3C variant.However, we note that MP2-F12 can be made faster using the 3*A approximation 11 without losing much accuracy in most cases, and MP2+CABS+PBE is only approximately three to four times faster than this cheaper MP2-F12 variant.Of course, the relative gains in computational cost would be much less for more expensive wave-function methods such as CCSD(T).

| CONCLUSION
We have reported an efficient density-fitting implementation of the DBBSC method in the MOLPRO software using different basis-set correction density-functional approximations and including the CABS single-excitation correction.We have tested the method on the FH51 benchmark set of reaction energies at the MP2 level and provided a comparison with the explicitly correlated MP2-F12 method.
For the smallest basis sets, the CABS single-excitation correction provides an important correction on reaction energies which is not included in the basis-set correction density-functional approximations.
The basis-set corrected reaction energies are quite insensitive to the choice of the basis-set correction density-functional approximation, with the notable exception of the LYP functional which gives much worse results.This point should be further analyzed in the future.
Overall, the basis-set corrected MP2 reaction energies calculated with a n-zeta basis set are of slightly higher quality than uncorrected MP2 reaction energies calculated with ðn þ 1Þ-zeta quality.However, the explicitly correlated MP2-F12 method is consistently more accurate, with reaction energies calculated with a n-zeta basis set being of slightly lower quality than uncorrected MP2 reaction energies calculated with ðn þ 2Þ-zeta quality.We believe that the DBBSC method is still valuable for accelerating the basis convergence of MP2 due to the fact that it has a lower computational cost than MP2-F12.
Finally, let us mention that the present implementation of the DBBSC method can be applied to any other wave-function methods, such as CCSD(T), with expected similar gains in accuracy.After completion of the present work, we became aware of a very similar independent work that has just been published. 63

n exact 2 ðrÞ
that the integral of the model on-top pair density equals the integral of the exact on-top pair density, with a highly accurate 418-term

- 1 F
the different basis-set correction density-functional approximations for the basis-set convergence of the ground-state MP2 correlation energy of the tetramethylpentane molecule (C 9 H 20 ).We see that all the density-functional approximations lead to a quite similar acceleration of the convergence of MP2 correlation energy toward its CBS limit.At the level of the correlation energy, all the proposed density-functional approximations thus provide a reasonable basis-set correction with quite a substantial acceleration of the basis-set convergence, albeit not as impressive as the one obtained with MP2-F12.The errors on the reaction energies of the FH51 set with respect to MP2/CBS calculated with MP2, MP2+CABS, MP2+CABS+PBE, and MP2-F12 are reported in Figure2.With the avdz basis set, MP2 I G U R E 1 Basis-set convergence of the ground-state MP2 correlation energy of the tetramethylpentane molecule (C 9 H 20 ) with different basis-set correction density-functional approximations (evaluated at the HF density) and with MP2-F12 using avnz basis sets.can have quite large basis errors for some reaction energies, up to about 13 kcal/mol.Obtaining MP2 reaction energies with all basis errors below 1 kcal/mol requires the use of the av5z basis set.The CABS single-excitation correction is crucial to reduce the largest basis errors on MP2 reaction energies obtained with the avdz basis set.Even with larger basis sets, the CABS single-excitation correction still helps to reduce the basis errors for some reaction energies.Adding the PBE-based basis-set correction further reduces the basis errors, albeit not always in a systematic way since there are a few cases where the basis error increases.It is noteworthy that the basis errors of the MP2+CABS+PBE reaction energies are all smaller than 1 kcal/mol with the avtz basis set and larger basis sets.MP2-F12 globally outperforms MP2+CABS+PBE, giving reaction energies with basis errors below about 1 kcal/mol already with the avdz basis set.
Mean absolute errors (in kcal/mol) in reaction energies of the FH51 set with respect to MP2/CBS with avnz basis sets.