Improving MP2 bandgaps with low-scaling approximations to EOM-CCSD Free

Second-order Møller–Plesset perturbation theory (MP2) is the simplest ab initio treatment of dynamical electron correlation. Its low cost makes it especially attractive for large systems, including periodic solids. Although periodic MP2 has been found to perform reasonably well for the description of ground-state properties,^1–12 its performance is less satisfactory for charged excitation energies and bandgaps.^5,13 For example, in Ref. 5, MP2 was applied to 13 semiconductors and insulators and exhibited average errors of 0.5% for lattice constants, 4.1% for bulk moduli, and 0.23 eV for cohesive energies but predicted negative bandgaps for materials that are known to be semiconducting, such as silicon and silicon carbide. This unsatisfactory performance was attributed to the lack of screening in finite-order perturbation theory. Indeed, the GW approximation^14–16 and equation-of-motion coupled-cluster theory^17–21 describe excitation energies with infinite-order perturbation theory and predict accurate bandgaps of semiconductors,^12,16,22 albeit with a computational cost that is higher than that of MP2.

Here, we study the performance of a second-order approximation to equation-of-motion coupled-cluster theory with single and double excitations (EOM-CCSD), first presented in Refs. 23 and 24. Despite making sequential second-order approximations, the method will be seen to be equivalent to the use of a self-energy containing all second-order diagrams and a few third-order diagrams.

Consider the Møller–Plesset partitioning of the many-body Hamiltonian, leading to the Hartree–Fock (HF) orbitals ϕ_p(r) with orbital energies ɛ_p; as usual, we denote the orbitals occupied in the HF determinant by i, j, k, l, those unoccupied by a, b, c, d, and general orbitals by p, q, r, s. The self-energy evaluated to second-order in perturbation theory is

where the antisymmetrized two-electron integrals are defined by $Vpqrs=⟨pq|rs⟩−⟨pq|sr⟩$⁠, with

where x is a combined space and spin variable. Unlike the GW approximation, the MP2 self-energy has exact second-order exchange and is therefore free of self-screening error.

An alternative theory can be obtained by a second-order approximation to EOM-CCSD, leading to a method originally referred to as EOM-MBPT(2)²³ or EOM-CCSD(2).²⁵ In this method, the ground-state CCSD amplitudes are approximated by their MP2 values, avoiding the expensive iterative solution of the CCSD amplitude equations. In this work, we consider the additional approximation of partitioning the EOM Hamiltonian into single and double excitation spaces and perturbatively treating the latter. Under this approximation, the large double excitation block of the similarity-transformed Hamiltonian is a diagonal matrix of orbital energy differences. This method has been referred to as DSO-GF²³ and P-EOM-MBPT(2);^{21,24,26–28} because we always use a Hartree–Fock reference, we will refer to the method as P-EOM-MP2.

Unlike typical Green’s function techniques, the EOM approach yields ionization potentials (IPs) and electron affinities (EAs) from separate eigenvalue calculations. In practice, these eigenvalues are found iteratively using the Davidson algorithm. As shown by Nooijen and Snijders,²³ the P-EOM-MP2 IPs can be equivalently obtained from the self-consistent eigenvalues of a matrix with elements $εiδij+ΣijEOM(ω)$⁠, where

and likewise for the EAs. The above matrix is clearly similar to the MP2 self-energy matrix (1), except for three differences. The first difference is the neglected coupling between the particle and hole spaces. Within the common diagonal approximation to the self-energy, this coupling is irrelevant and we have numerically confirmed that it is a negligible difference in this work. The second difference is the perturbative replacement of ω = ɛ_j in one of the two terms. When this replacement is done in the MP2 self-energy, we find that it makes the results slightly worse and is thus not responsible for the improvement to be shown in the P-EOM-MP2 bandgaps (vide infra).

The third and most important difference is the presence of the intermediate

where the antisymmetrization operator is P₋(kl)A_kl = A_kl − A_lk and

Viewed in terms of the similarity-transformed Hamiltonian $H̄=e−THeT$⁠, the first term in Eq. (3) reflects a renormalization of the one-body interactions due to ground-state correlation and the presence of W in the second term reflects a renormalization of the two-body interactions, i.e., it is a screened Coulomb interaction. Alternatively, the intermediate W can be understood as the inclusion of a few third-order self-energy diagrams, with a perturbative evaluation of the frequency denominator. These conclusions concerning the importance of the various differences between MP2 and P-EOM-MP2 are the same as those observed for the IPs of molecules.²³ In Fig. 1, we show the time-ordered self-energy diagrams included in the MP2 and P-EOM-MP2 Green’s functions beyond first-order (i.e., beyond Hartree–Fock).

The computational and storage costs of P-EOM-MP2 are significantly less than those of EOM-CCSD. Assuming N_o occupied orbitals and N_v virtual orbitals with N_v > N_o, the ground-state CCSD calculation has an iterative $O(No2Nv4)$ cost, whereas the ground-state MP2 calculation has a non-iterative $O(No2Nv2)$ cost [or O(N⁵) when including the one-time integral transformation]. The cost of excited-state calculations by both methods is determined by the non-iterative cost of forming matrix elements of the similarity transformed Hamiltonian $H̄$ and the iterative cost of contracting this matrix with a trial vector R. Both methods have the same non-iterative $O(No3Nv3)$ cost associated with the formation of intermediates such as the one shown in Eq. (4), but the cost of this step can be dominated by that of the iterative steps. In EOM-CCSD for EAs, the most expensive contraction is due to the doubles–doubles block,

with $O(NoNv4)$ scaling. In P-EOM-MP2 for EAs, the diagonal approximation to the doubles–doubles block shifts the most expensive contraction to the off-diagonal blocks that couple single and double excitations, e.g.,

with $O(NoNv3)$ scaling. If the one-time O(N⁶) intermediate construction is too expensive, then it can be avoided by reordering the contraction in Eq. (7) to yield iterative $O(No2Nv3)$ scaling. Regarding storage, EOM-CCSD requires four-virtual integrals with $O(Nv4)$ storage, whereas P-EOM-MBPT2 only requires three-virtual integrals with $O(NoNv3)$ storage.

We have applied the above two theories to the calculation of the minimum bandgaps of 11 simple, three-dimensional semiconductors and insulators. Seven have a diamond/zinc-blende crystal structure, Si, SiC, GaP, BP, GaN, C, and BN; two have a rock-salt crystal structure, MgO and LiF; and two have a face-centered cubic crystal structure, Ar and Ne. Calculations were done with periodic boundary conditions using the PySCF software package.^29,30 All calculations were done without pseudopotentials using the all-electron cc-pVTZ basis set except for Ne and Ar, which used the aug-cc-pVTZ basis set. Calculations using larger basis sets (not shown) suggest that our results are close to the basis set limit, consistent with analogous results obtained with the GW approximation.³¹ Two-electron integrals were calculated by periodic Gaussian density fitting³² using JKFIT auxiliary basis sets.³³ 

For charged excitation energies, finite-size effects are large.^12,34 Here, we have included one Madelung constant correction to the occupied orbital energies and another to the final IPs. The former correction has no impact on wavefunction-based theories, such as EOM-CCSD, but does have an impact on finite-order perturbation theories (similar to the differing behaviors of ground-state CCSD and MP2); the latter correction is familiar from periodic calculations of charged systems and can be given a many-body interpretation on the basis of the excited-state structure factor.³⁴ We have performed calculations with N_k = 2³–5³ k-points sampled uniformly in the Brillouin zone. Bandgaps were then extrapolated to the thermodynamic limit assuming an $O(Nk−1/3)$ finite-size error. Other treatments of finite-size effects are possible, but all are expected to exhibit finite-size errors with the same scaling. The rock-salt and face-centered cubic materials have direct bandgaps at the Γ point, and the diamond/zinc-blende materials have indirect bandgaps, which were determined by performing separate IP and EA calculations with k-point meshes shifted to include the valence band maximum and conduction band minimum, respectively.

In Fig. 2, we compare the minimum bandgaps obtained by MP2, P-EOM-MP2, and the GW approximation to experimental values at 300 K. Because the calculations do not account for vibrational effects, we have adjusted the experimental values according to calculated electron–phonon renormalizations from the literature^39,40 based on the Allen–Heine–Cardona framework.^42–44 We only include the zero-point renormalization for all materials except LiF, which has a sizable thermal contribution to the renormalization at 300 K;⁴¹ for the other materials, this latter contribution is relatively small. Lattice expansion is already accounted for because our lattice constants are experimental 300 K values. Precise numbers and crystal geometries are given in Table I. We note that experimental bandgaps and calculated electron–phonon renormalizations vary throughout the literature.

Consistent with Ref. 5, we find that MP2 systematically underestimates the bandgap and predicts negative bandgaps for Si, GaP, BP, and SiC (our MP2 bandgaps are similar to those of Ref. 5, but some differ by as much as 0.5 eV, which we attribute to differences in the treatment of core electrons, basis set effects, and finite-size effects). Remarkably, the P-EOM-MP2 bandgaps are a significant improvement and show good agreement for all materials. The mean signed error (MSE) is +0.28 eV, and the mean unsigned error (MUE) is 0.40 eV. The largest signed error is for Si (+0.96 eV), which has the smallest gap of all materials considered.

In Fig. 2 and Table I, we also compare results to those calculated by the G₀W₀ approximation with a Perdew–Burke–Ernzerhof (PBE) reference. For all materials except LiF, we compare to all-electron, full-frequency calculations by Zhu and Chan,³¹ which were performed with PySCF using identical treatments of core electrons and identical Gaussian basis sets. The result for LiF is from Ref. 35. Remarkably, the P-EOM-MP2 and GW approximation perform similarly well, despite their underlying physical differences. Roughly speaking, the GW approximation performs better for materials with the smallest gaps, while P-EOM-MP2 performs better for those with the largest gaps. The largest errors for the GW approximation are for the large-gap insulators, whose bandgaps are underestimated by about 1 eV or more, which we attribute to the use of a PBE starting point and the absence of second-order exchange. For the GW approximation, the MSE is −0.62 eV and the MUE is 0.64 eV.

We performed additional calculations to estimate the effect of the various diagrams included in the P-EOM-MP2 self-energy shown in Fig. 1. The HF bandgap is always much too large, and the direct second-order diagrams yield a large negative correction, which is largely responsible for the performance exhibited by MP2. Consistent with Ref. 5, we find that the second-order exchange diagrams make a small contribution (≲0.2 eV) for small-gap materials but a larger contribution (≳0.5 eV) for large-gap materials. This can be attributed to the more localized nature of the electronic states in large-gap insulators. The diagram in Fig. 1(c) makes a large contribution (1 eV or more) for all materials and is most responsible for the significant improvement of P-EOM-MP2 over MP2. The final diagram in Fig. 1(d), a vertex correction beyond the GW approximation, typically raises the gap by about 0.2 eV.

In Fig. 3, we show the convergence of the bandgap toward the thermodynamic limit for four of the materials considered here. As mentioned earlier, the finite-size error is large and must be removed by extrapolation. Interestingly, although MP2 and P-EOM-MP2 give similar bandgaps for small k-point meshes, they exhibit very different convergence to the thermodynamic limit. This difference is largest for materials with small bandgaps. We attribute this behavior to our use of Madelung constant corrections in the HF orbital energies. These corrections cause the HF gap to converge to the thermodynamic limit from above such that the systems with smaller k-point meshes are more weakly correlated and the importance of third-order diagrams in the self-energy is diminished. On approach to the thermodynamic limit, the system becomes more strongly correlated and the results of the two methods deviate.

In conclusion, we have shown that the P-EOM-MP2 approach, a second-order approximation to EOM-CCSD, yields solid-state bandgaps that are a significant improvement over those predicted by MP2. The success of P-EOM-MP2 contradicts the conventional wisdom that infinite-order screening is necessary for quantitative accuracy in bandgap prediction. Rather, P-EOM-MP2 represents an affordable balance of low-order screening and exchange, yielding semiquantitative accuracy for materials with a wide range of bandgaps. By starting from Hartree–Fock theory and including antisymmetrization in all interaction vertices, the method is completely ab initio and free of self-interaction and self-screening errors. We note that P-EOM-MP2 is very similar to CC2,⁴⁵ and we therefore expect similar performance from the latter, which also includes some amount of orbital relaxation. Although P-EOM-MP2 has been found to perform well for three-dimensional materials, it will be interesting to apply it to low-dimensional semiconductors, where screening is more complicated. The present results should be compared to those produced with the more expensive EOM-MP2 (without the diagonal approximation) and EOM-CCSD methods. We have recently done this comparison for neutral excitations calculated with EE-EOM-CCSD,⁴⁶ and work is currently in progress for band structures and bandgaps.

This work was supported, in part, by the National Science Foundation under Grant No. CHE-1848369. M.F.L. was supported, in part, by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1746045. 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) Contract No. C090171, both awarded on April 15, 2010. We also acknowledge computing resources from the Flatiron Institute. The Flatiron Institute is a division of the Simons Foundation.