Dynamical correlation energy of metals in large basis sets from downfolding and composite approaches Free

Ground-state electronic properties of metallic solids have traditionally been computed using density functional theory (DFT),^1–4 which is partially justified by the fact that many popular functionals are parameterized by numerically exact results on the uniform electron gas (UEG).^5–8 In recent years, interest has grown around the application of ab initio wavefunction-based electronic structure techniques for condensed-phase systems^9–16 since they do not suffer from uncontrolled errors inherent to the DFT exchange–correlation functional.^6,17,18 Promising methods in this direction include the random-phase approximation (RPA)^16,19,20 and coupled-cluster theory.^21–26 Importantly, both these methods preclude the well-known divergences of finite-order perturbation theories, such as second-order Møller–Plesset perturbation theory (MP2) via an infinite-order resummation.^27–31 

Although coupled-cluster theory has been successfully applied to an increasing number of atomistic semiconductors and insulators,^{14,15,32–37} its applicability to metals has been primarily focused around the UEG, also known as jellium.^29,38–45 Despite their reasonable accuracy, these calculations have demonstrated the typical slow convergence of the correlation energy as a function of the number of virtual (unoccupied) orbitals included.^11,40,46,47 This slow convergence is especially problematic because of the high cost of coupled-cluster calculations with large basis sets. For example, coupled-cluster theory with single and double excitations (CCSD) has a computational cost that scales as O(N²M⁴), where N and M are the number of electrons and basis functions, respectively. To date, the results near the complete basis set (CBS) limit have been primarily computed via the extrapolation of results obtained with a finite increasing number of basis functions,^24,40,45 although explicitly correlated⁴⁸ and transcorrelated^49,50 methods provide promising alternative approaches.

Composite methods (sometimes called focal point methods) are a simple, alternative class of approaches for recovering dynamical correlation within large basis sets.^51–53 A common composite scheme combines the results of high-level and low-level theories using three calculations. For example, using CCSD as the high-level theory and MP2 as a low-level theory, the CCSD correlation energy in a large basis is approximated as

where M_act < M is the number of “active” basis functions. References 54 and 55 provide a similar but more sophisticated CCSD/MP2 composite method based on an analysis of the basis set convergence of various diagrammatic contributions to the correlation energy. While such CCSD/MP2 composite approaches have been applied successfully to a number of semiconductors and insulators,^36,54,56,57 their applicability to metals is questionable because of the failures of MP2 theory. One goal of this work is to test the composite CCSD/MP2 approach for three-dimensional metals.

A more theoretically satisfying approach would be to perform a single calculation where low-energy excitations near the Fermi surface are treated with CCSD and are coupled to high-energy excitations treated with MP2. This particular approach, which is similar to tailored CC⁵⁸ and the broader class of active-space CC methods,^59–63 has variously been called CC/PT,⁶⁴ CCSD-MP2,⁶⁵ and multilevel CC.^66,67 Two of us (Lange and Berkelbach) recently tested this method for a few simple atomistic semiconductors and insulators,⁵⁶ and here, we aim to assess its performance for three-dimensional metals, where the differences between CCSD and MP2 are more striking. Since the effects of the frozen high-energy MP2 amplitudes are folded down onto the low-energy CCSD amplitudes (see below), we refer to this method as a “downfolding” approach. In principle, this downfolding CCSD/MP2 method should provide a distinct advantage over the conceptually simpler composite approach, as downfolding does not include the MP2 treatment of low-energy excitations that are responsible for divergence in the thermodynamic limit (TDL). After providing the theoretical details of these two methods, we compare their performance for the UEG at a fixed number of electrons and in the TDL. Before concluding, we also examine the straightforward use of the direct RPA in place of MP2.

Here, we briefly review the theory underlying the downfolding and composite approaches. The N occupied spin-orbitals are indexed by i, j, k, and l; the $M−N$ virtual orbitals by a, b, c, and d; and the M general orbitals by p, q, r, and s. The MP2 and coupled-cluster with double excitations (CCD) correlation energies are given by

where $tijab$ are the amplitudes of the double excitation operator $T2=14∑ijabtijabaa†ab†ajai$ and ⟨pq‖rs⟩ are antisymmetrized two-electron repulsion integrals; contributions from single excitations vanish because the UEG has no capacity for orbital relaxation by symmetry. At the lowest order in perturbation theory,

and Eq. (2) gives the MP2 correlation energy. The high density of states at the Fermi surface, the long-ranged nature of the Coulomb potential, and the dimensionality are together responsible for the divergence of MP2 in the TDL. We note that a vanishing energy denominator alone is not responsible for the divergence. For example, the two-dimensional electron gas has a well-defined MP2 energy but a divergent MP3 energy.⁷ By contrast, the CCD amplitudes solve a system of nonlinear equations

where H is the electronic Hamiltonian. A standard approach for reaching the CBS limit is to perform a series of calculations with increasing M and use a M⁻¹ extrapolation.

In both the composite and downfolding approaches, we partition the orbitals into a set of M_act active orbitals, composed of all occupied orbitals and low-energy virtual orbitals, and a set of $M−Mact$ frozen (inactive) orbitals, composed of high-energy virtual orbitals. In principle, occupied orbitals can also be partitioned, but typically they do not significantly contribute to the computational cost. In the composite CCD/MP2 approach, the correlation energy is calculated according to Eq. (1). Importantly, for metals, the low-energy active space double excitations are treated by CCD and not by MP2, so we expect the method to be well-behaved in the thermodynamic limit.

In the downfolding CCD/MP2 approach, the double excitation operator T₂ is partitioned into internal excitations fully contained within the active space and external excitations that involve at least one frozen orbital, $T2=T2(int)+T2(ext)$⁠. Fixing the $T2(ext)$ amplitudes to their MP2 values via Eq. (3), the downfolding method involves first solving Eq. (4) for only the internal amplitudes and then evaluating the correlation energy expression Eq. (2) using both the internal and external amplitudes. Compared to the O(N²M⁴) cost of full CCD, the composite approach has $O(N2Mact4)+O(N2M2)$ cost and the downfolding approach has $O(N2Mact2M2)$ cost, which can provide significant savings, depending on the practical value of the ratio M_act/M.

Let us now provide more insight into the “downfolding” perspective. Note that because the internal and external excitation operators commute, the defining energy and amplitude equations of the downfolding approach can also be written as

where $H̄=e−T2(ext)HeT2(ext)$ (with fixed $T2(ext)$ as detailed above) and $EMP2(ext)=14∑ijab(ext)tijab⟨ij‖ab⟩$ is the MP2 correlation energy due to external excitations. These resemble ordinary CCD energy and amplitude equations within the active space only, except that the bare Hamiltonian H is replaced by an effective Hamiltonian $H̄$ that is similarity-transformed by the external excitation amplitudes. The effective Hamiltonian within the active space can be expressed as

where {⋯} indicates the normal ordering of the operators. This effective Hamiltonian can be seen to contain effective one- and two-body interactions that are frequency independent,²¹ in contrast to other downfolding approaches, such as the constrained random-phase approximation.^68,69 For example, the all-occupied two-body interaction becomes

where the primed summation indicates that one or both a and b are inactive virtual orbitals. The frequency independence can be understood because our observable is the total energy rather than a spectral function. To summarize, an approach that solves the internal CCD amplitude equations in the presence of frozen external amplitudes is equivalent to a CCD calculation in an active space of orbitals using an effective (downfolded) Hamiltonian that is similarity-transformed by the external excitation operator.

It is straightforward to show that the composite approach, normally understood as a three-step procedure as shown in Eq. (1), is equivalent to Eq. (5) but where the effective Hamiltonian $H̄$ is replaced by the bare Hamiltonian H in the amplitude equation (5b). From this perspective, the performance differences between the downfolding and composite approaches are attributable to the screening of the integrals in the effective Hamiltonian when determining the internal amplitudes. Nevertheless, we reiterate that the composite CCD/MP2 approach is expected to perform well because it replaces the problematic MP2 treatment of low-energy, internal double excitations with a well-behaved CCD treatment.

We study the UEG as the simplest model of metals. A brief review of the UEG model in finite cells with finite plane-wave basis sets is given in the  Appendix, and we refer the reader to the literature for more details.^31,40,70 To illustrate the performance of the composite and downfolding methods, we focus on the Wigner–Seitz radius r_s = 4 (corresponding to the approximate valence electron density of metallic sodium), where CCD has been found to recover about 85% of the correlation energy.^40,43 We use a twisted boundary condition by performing calculations at the Baldereschi point,⁷¹ which has been shown to provide smoother convergence to the TDL.^72,73

In Fig. 1, we show the basis set convergence of the correlation energy for a finite UEG with N = 332 electrons. The uncorrected MP2 and CCD correlation energies exhibit their typical slow convergence with increasing basis set size and show asymptotic behavior where the basis set error decays as M⁻¹. Extrapolation to the CBS limit yields E_c/N = −0.0401 E_h for MP2 and E_c/N = −0.0262 E_h for CCD. At the largest finite basis shown, M = 1502, the results exhibit a significant basis set error of about 0.01 E_h for both methods, highlighting the challenge of recovering dynamical correlation in metals with large basis sets. Importantly, we emphasize that the MP2 correlation energy does not diverge for any finite system but only upon extrapolation to the TDL (see below).

Recall that the CCD/MP2 composite and downfolding approaches involve both a “target” number of orbitals M and an active number of orbitals M_act. In Fig. 1, we show the results obtained for M = 1502 as M_act is varied. By construction, both methods yield the target MP2 correlation energy when there are no active virtual orbitals (M_act = N) and the target CCD correlation energy when all orbitals are active (M_act = M). We observe that both methods converge smoothly to the target CCD result and that the downfolding approach exhibits a faster convergence due to its coupling between the internal and external excitation spaces. We also see a similar behavior for other numbers of electrons and densities (not shown), indicating that neither finite-size effects nor the specific metallic density changes the overall picture.

To better quantify the rate of convergence, in Fig. 2, we plot the absolute deviation of the correlation energy from the “target” CCD result obtained with M = 1502. The error is plotted as a function of the difference between the inverse number of active orbitals and the inverse number of total orbitals and analyzed in terms of the power law $|ΔEc|∝Mact−1−M−1α$⁠. We compare the convergence of traditional CCD, the composite approach, and the downfolding approach. For plain CCD, we see the linear convergence of the correlation energy, with α ≈ 1, over a large range of basis set sizes. The composite method exhibits an early, rapid convergence reaching a maximum scaling of around α ≈ 2 before slowing to the same α ≈ 1 convergence as M_act approaches M. The rapid convergence for small M_act is responsible for absolute errors that are about one order of magnitude better than those obtained by simple truncation. In fact, the plain CCD result does not obtain mE_h accuracy until essentially all orbitals are correlated, whereas the composite result achieves this accuracy when only 50% of the virtual orbitals are correlated in the expensive CCD calculation; this results in a speedup of a factor of 16 compared to the full CCD calculation. Finally, the downfolding result exhibits rapid but non-monotonic convergence, making it difficult to extract a power law. Before slightly overshooting the “target,” the power law exponent reaches α ≈ 3 or better, a significant improvement over the composite CCD/MP2 and standard CCD approaches. This fast rate of convergence provides mE_h accuracy when about one third of the virtual orbitals are correlated, giving a speedup of a factor of 9.

The good performance of the composite approach indicates that MP2, while an inapplicable theory for three-dimensional metals, is safe to use for basis set corrections. As discussed above, the reason for this applicability can be understood by considering the MP2 correction that is applied in Eq. (1). This correction is a difference between two MP2 correlation energies, both of which correlate orbitals near the Fermi surface. These two MP2 energies are separately divergent in the TDL, but their difference is not; moreover, this difference is precisely $EMP2(ext)$ defined previously.

The reliable scaling of the MP2 basis set error at large M suggests that the MP2 CBS limit can be obtained by extrapolation for any given number of electrons. In contrast, the asymptotic scaling regime for the CCD correlation energy cannot always be reached. Thus, for any calculation performed with a given M, we propose to add the MP2 correlation energy difference $δ(2)M=EMP2∞−EMP2M$⁠, where $EMP2∞$ is obtained by M⁻¹ extrapolation. Having obtained an estimate of the CBS limit for a given number of electrons, the finite-size extrapolation to the TDL can be done separately. We expect this scheme to be not only more reliable than extrapolating CCD on its own but also less costly since it involves more MP2 and fewer CCD calculations.

In Fig. 3(a), we plot the MP2 correlation energy as a function of the inverse number of electrons for finite UEG systems containing N = 90–2392 electrons. For each system size, we performed MP2 calculations at different basis set sizes (gray squares), and the top four gray curves connect systems with a similar ratio of M/N. We then performed M⁻¹ extrapolations at each particle number using data from M/N ≈ 3.0, 3.5, 4.0, and 4.5 to obtain the CBS limit at each system size (black stars). Upon approaching the TDL, the MP2 correlation energy diverges, as seen most easily in our largest calculations for the smaller values of M/N. Despite the divergence of its components, the MP2 CBS correction δ⁽²⁾, plotted explicitly in Fig. 3(a) for M/N ≈ 4 (pink thin diamonds), does not diverge and is thus safe for use in metallic systems.

To confirm that the success of the composite CCD/MP2 method depends on approaching the TDL, in Fig. 3(b), we plot results for CCD at finite M/N (green circles) and composite CCD/MP2 (orange crosses). For the composite CCD/MP2 results, we applied the MP2 CBS correction $δ(2)M$ to the corresponding CCD results at each finite value of M/N. For comparison, we also show CCD results where the CBS limit is estimated by simple M⁻¹ extrapolation (dark green stars). In contrast to MP2, CCD is well defined in the TDL; however, at each system size, convergence to the CBS limit is slow, as shown by the green datasets. Simple M⁻¹ extrapolation gives a very large and potentially unreliable CCD CBS correction. In contrast, we observe that composite CCD/MP2 has much faster convergence to the CBS limit, as shown by the orange datasets. Performing an N⁻¹ extrapolation of the CBS-corrected composite CCD/MP2 with M/N ≈ 4.0 gives us a CBS and TDL extrapolated correlation energy of E_c = −0.0293 E_h. Using simple M⁻¹ and N⁻¹ extrapolations of the CCD results yields a slightly smaller correlation energy of E_c = −0.0273 E_h, giving an approximate error bar of 2 mE_h on our extrapolated results. For comparison, the exact value⁷⁴ is E_c = −0.0318 E_h, indicating that CCD recovers 85%–90% of the correlation energy, in general agreement with past CCD results.⁴⁰ The same finite-size extrapolation of the CBS-corrected M/N ≈ 3.0 data gives a correlation energy that differs by only 0.7 mE_h, indicating the excellent convergence of the composite method.

Before concluding, we recognize that a variety of other low-level theories can be combined with CCD, in both the composite and downfolding manner. Of particular interest is the direct RPA (dRPA), which is more appropriate for three-dimensional metals than MP2. The dRPA amplitudes are the solution of the CCD equations where only selected terms are retained and only direct (non-antisymmetrized) electron repulsion integrals are used.^70,75 Based on its good performance and lower cost, we focus only on the composite method.

The improved accuracy of dRPA over MP2 might make its higher cost worthwhile if the composite method can be converged with a smaller active space. In Fig. 4, we compare how the choice of low-level theory affects the performance of the composite method for the same N = 332 UEG from Fig. 1. In order to test the transferability of our conclusions, we show results at two densities, r_s = 4 and r_s = 1. At r_s = 4, the dRPA result closely tracks the CCD result for all the values of M, which unsurprisingly yields a composite method with excellent performance. At r_s = 1, the behavior is qualitatively similar, although a non-negligible deviation between the dRPA and CCD correlation energy occurs for large M values.

We now aim to better quantify and compare the performance of the composite CCD/MP2 and CCD/dRPA approaches, especially with the overall computational cost in mind. Moreover, we would like to assess the dependence of the results on the large, “target” basis set size. In Fig. 5, we plot the percentage of virtual orbitals in the CCD active space required to achieve an error of 1 mE_h as a function of the target basis set size M/N. We note that the success of a composite method depends on the parallelism of the correlation energies of the high- and low-level theories as a function of M⁻¹. When the target basis set is small, the MP2 and CCD correlation energies exhibit very different behaviors, and thus, the composite CCD/MP2 method is quite poor. Almost all the virtual orbitals need to be included in the CCD active space to achieve an error below 1 mE_h. However, as the target M increases, the required percentage drops monotonically due to the near-perfect parallelism of the MP2 and CCD correlation energies. By studying smaller system sizes (not shown), we have checked that this behavior continues out to larger values of M/N.

Turning to the composite CCD/dRPA results, we see a very rapid decrease in the percentage of virtual orbitals required in the CCD active space. Beyond M/N ≈ 2, less than 5% of the total number of virtual orbitals needs to be treated by CCD. Despite the increased cost of dRPA over MP2, this reduction in the number of orbitals to be treated by the more expensive CCD more than compensates and the dRPA-based correction is preferred. However, this improved performance is unfortunately not universal. The deviation between dRPA and CCD seen at large M values in Fig. 4(b) (at r_s = 1) spoils the composite method. When this deviation occurs, around M/N ≈ 3.5, the percentage of virtual orbitals required in the active space jumps to about 50% and continues to increase due to the growing deviation between the two methods. The results at smaller N (not shown) confirm that this behavior continues to larger values of M/N. Therefore, at least for this density and at large basis set sizes, the composite CCD/MP2 method is both cheaper and more accurate than the composite CCD/dRPA method.

We have described and analyzed two approaches for eliminating basis set error in the CCD correlation energy of metals using the simple UEG model. Our results indicate that these methods allow for aggressive freezing of virtual orbitals or approximation of external amplitudes, leading to significant reductions in computational cost. Although the downfolding CCD/MP2 approach is slightly more accurate, we find that the simpler composite CCD/MP2 approach is surprisingly effective because divergent contributions near the Fermi surface do not contribute to the basis set correction. For typical basis set sizes with M/N ≳ 4, the CCD/MP2 approach exhibits errors less than 1 mE_h when about half of the virtual orbitals are included in the more expensive CCD calculation. Given the scaling of CCD, this yields a factor of ∼16 speedup or better. In most cases, the composite CCD/dRPA approach was shown to perform even better, typically requiring that only a small percentage of the virtual orbitals need to be included in the CCD calculation. However, we observed a failure of the method at large basis set sizes for one of the two densities studied. These findings tentatively suggest that the composite CCD/dRPA approach is less expensive (for a desired accuracy) but less robust than the CCD/MP2 approach.

In this work, we have addressed post-Hartree–Fock basis set errors in the canonical orbital basis, but the methods we presented could also be straightforwardly applied to a basis of localized orbitals.^{18,57,76–79} Localized orbital basis sets mix orbitals near and far from the Fermi surface, so it will be interesting to test how the composite and downfolding approaches perform when using these localized orbitals for metals.

Future work will focus on applying these techniques to atomistic metals using natural orbitals,⁸⁰ to the excited-state properties of metals,^41,44 and to higher-level theories of correlation.⁴³ For example, we imagine that a composite CCSDT/CCSD or CCSDT/CCSD(T) approach would provide quantitative accuracy for metals while precluding the failure³¹ of the otherwise successful treatment of perturbative triple excitations.

Note added in proof. After this work was submitted, Refs. 81 and 82 were posted to arXiv and provide complementary approaches to basis set corrections for periodic CCSD, especially as applied to metals.

We thank Verena Neufeld and Sandeep Sharma for their comments on this manuscript. This work was supported, in part, by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1644869 (M.F.L.), the Department of Defense through the National Defense Science and Engineering Graduate (NDESG) Fellowship Program (J.M.C.), and the National Science Foundation under Grant No. CHE-1848369 (T.C.B.). We acknowledge computing resources from Columbia University’s Shared Research Computing Facility project, which is supported by NIH Research Facility Improvement Grant No. 1G20RR030893-01, and associated funds from the New York State Empire State Development, Division of Science Technology and Innovation (NYSTAR), under Contract No. C090171, both awarded April 15, 2010. The Flatiron Institute is a division of the Simons Foundation.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Working in a single-particle basis of plane waves with momenta k = (2π/L)(n_x, n_y, n_z) and cell of volume L³ with periodic boundary conditions, the UEG Hamiltonian is given by

where the primed summation requires k₁ + k₂ = k₃ + k₄. The two-electron repulsion integrals are given by

where the Ewald potential is

where v_M = 2.837 297 479/L is the Madelung constant of the cell.⁸³ The N-electron reference determinant has the lowest-energy N/2 plane wave orbitals doubly occupied, and the HF orbital energies are given by $ε(k)=k2/2−∑k′≤kFv(k−k′)$⁠, where k_F is the Fermi momentum. We restrict our calculations to closed-shell configurations, which allows only certain “magic numbers” of electrons and orbitals. At the Baldereschi point, the first few magic numbers are 2, 8, 14, 22, 34, 40, and 52.