Regularized second-order correlation methods for extended systems

Second-order Møller–Plesset perturbation theory (MP2) holds a unique place in the hierarchy of wavefunction-based electronic structure methods.^1–4 Historically, it was the first correlated method in quantum chemistry that was size-extensive, invariant to unitary orbital rotations, and computationally affordable. Until the development of modern density functional theory (DFT), it was, therefore, the workhorse method of molecular quantum chemistry.⁵ 

Even today, MP2 still plays an important role, e.g., in its spin-component-scaled variants or as a part of double-hybrid DFT methods.^6–8 Despite this popularity, the limitations of MP2-like methods are also well-known. In particular, they fail spectacularly for strongly correlated or metallic systems.^9–11 Furthermore, a strong overestimation of dispersion interactions is observed for large polarizable systems due to the absence of higher-order screening effects.^12,13

With the significant recent interest in implementing and applying wavefunction methods to solids, MP2 is again at the center of attention as it represents the natural first step in such endeavors.^14–19 However, the known problems for large and metallic systems clearly limit the usefulness of MP2 in this context. Indeed, canonical MP2 by construction must fail for metals, since the vanishing bandgaps lead to diverging contributions to the correlation energy. Interestingly, coupled-cluster (CC) methods do not share this problem, despite being closely related to MP2.^10,20 This improved behavior can be interpreted as a renormalization of the bandgap due to the inclusion of screening effects. This makes CC methods highly attractive for condensed matter applications.^17,21 Unfortunately, this advantage comes at a significantly increased computational cost.

In light of these issues, there has been significant interest in obtaining more robust MP2-like methods. In particular, different forms of regularization have been proposed as a way of empirically imitating higher-order screening effects at the MP2 level.^22–25 This concept has proven very fruitful for molecular applications, but to the best of our knowledge, it has so far not been applied to extended systems.

In this paper, we investigate the performance of regularized MP2-like methods for some prototypical periodic systems in which each constitute known issues for canonical MP2. In particular, we consider the homogeneous electron gas (HEG) (a metal prototype), the one-dimensional Hubbard model (a strongly correlated system), and the interaction of graphene with a water molecule (a challenging dispersion-driven interaction). We also propose a new, non-empirical regularization method based on the second-order Brillouin–Wigner (BW) perturbation theory.

In the following, we use indices i, j, k, l for occupied and a, b, c, d for unoccupied spin-orbitals χ, respectively. Antisymmetrized two-electron repulsion integrals in this basis are denoted as ⟨ij‖ab⟩. Using this notation, the MP2 correlation energy can be written as

where the denominator $Δijab=εa+εb−εi−εj$ is computed from the corresponding orbital energies ɛ.

Clearly, this energy must diverge if any $Δijab$ becomes zero, i.e., for metallic systems. In principle, a small constant δ could be added to the denominator to avoid this.²³ Such a level-shift prevents the division by zero, thus regularizing the MP2 correlation energy expression. This raises the question of how large the regularizer δ should be, however. Unfortunately, it has been found that there is no simple answer to this question: no single value of δ can both restore Coulson–Fisher points for single bond breaking and retain good thermochemical performance for weakly correlated systems.²³ In other words, the δ-regularization approach is not flexible enough to fix the divergence of MP2 while retaining its merits.

To address this, Lee and Head-Gordon explored several more sophisticated regularization schemes, in which each contribution to the MP2 energy is individually regularized according to the denominator $Δijab$⁠.²⁴ This allows attenuating the offending terms without affecting the well-behaved ones. The most successful of these schemes is κ-regularization, which uses the expression

For orbital-optimized MP2, this form of regularization (using κ = 1.4 $EH−1$⁠) was found to yield accurate molecular thermochemistry while also restoring the Coulson–Fisher points in single-bond dissociation curves (the absence of the latter being a well-known failure of canonical MP2).²⁴ 

In this work, we also explore a different regularization approach. This is based on Brillouin–Wigner (BW) perturbation theory, where the second-order correlation energy reads²⁶ 

This equation must be solved iteratively, as $EcBW2$ appears on both the left- and right-hand sides.

BW2 has some formal advantages over MP2. It yields the exact correlation energy for two-level systems and avoids the divergence of the energy for vanishing $Δijab$⁠.^26,27 There was therefore considerable interest in the Brillouin–Wigner series in the early days of quantum chemistry, which persists in the context of multi-reference perturbation theory.^28,29 However, BW2 also has a considerable downside, namely, that it is not size-extensive.³⁰ This is obviously problematic for applications to extended systems. Nonetheless, the idea of using a correlation energy-dependent regularization is attractive, since the diverging correlation energies (e.g., beyond Coulson–Fisher points in bond-breaking) are a clear signal that regularization is needed.

We, therefore, propose a simple modification of BW2, which restores size-extensivity (termed xBW2 in the following):

with the number of electrons N_e. In other words, the correlation energy in the denominator is replaced by the correlation energy per electron. In this way, a constant regularizer is obtained in the limit of infinite systems so that the size-extensivity of MP2 is recovered. This is shown numerically for He chains in Fig. 1. Note that, while we introduce xBW2 here as a form of regularized MP2, BW is a perturbation theory in its own right. For more details on the relation of MP2, BW2, and xBW2 from the perspective of formal perturbation theory, see  Appendix B.

Furthermore, $EcNe≪Δijab$ for weakly correlated systems, so that MP2 is approximately recovered in this limit. This is important, as we do not want to destroy the good performance of MP2 in such cases, e.g., for main group reaction energies. To illustrate this numerically, the correlation between MP2 and xBW2 reaction energies is shown in Fig. 2 for 17 isomerization and 98 bimolecular reactions involving 66 small main-group molecules from the W4-11-GEOM database³¹ (full data in the supplementary material).³² This shows that xBW2 and MP2 are nearly identical in this case with 8 and 9 meV mean absolute deviation for the isomerization and bimolecular reaction energies, respectively.

To the best of our knowledge, this method has not been proposed before in the literature. It is closely related to several other approaches, however. On one hand, the trick of using the correlation energy per electron instead of the total correlation energy can also be used to derive the averaged coupled-pair functional (ACPF) from the configuration interaction approach.^33,34 In this sense, xBW2 can be seen as a second-order approximation to ACPF. On the other hand, a similar second-order expression was also derived by Zhang, Rinke, and Scheffler from the Bethe–Goldstone equation (BGE2).³⁵ In this method, pair correlation energies are used for regularization instead of the correlation energy per atom. Consequently, each electron-pair ij has a different regularizer. In this sense, BGE2 is related to the coupled electron-pair approximation (CEPA)³⁶ and xBW2 can be considered an “averaged” BGE2. This averaging has the advantage that (unlike BGE2 and κ-MP2) xBW2 is invariant to unitary orbital transformations, although the BGE2 regularization is arguably more flexible. More generally, the iterative nature of the xBW2 equation resembles coupled-cluster theory, where the coupling between individual amplitudes is replaced by an average coupling term that is the same for all amplitudes.

As this discussion shows, a series of regularized MP2 methods have been proposed in the literature, using constant and dynamic regularization terms. In the following, we will focus on two dynamic schemes, namely, κ-MP2 (as a prototypical semi-empirical regularization based on $Δijab$⁠) and xBW2 (representing a non-empirical, energy-dependent regularization scheme).

As a first model system, we consider the homogeneous electron gas (HEG), which plays a central role in understanding the properties of simple metals and is of essential importance to the foundations of density functional theory. The divergence of second-order perturbation theory for the infinite HEG has long been proven analytically.^9,37,38 More recently, finite periodic electron gases have emerged as important numerical benchmark systems for many-body theories.^10,20,39,40 Importantly, this showed that MP2 also diverges for these systems, even though they display significant energy gaps. Since the divergence of canonical MP2 for the HEG is thus well established, it is of particular interest to understand how the regularized methods behave in comparison.¹⁰ 

In Fig. 3, the per electron correlation energies of κ-MP2 and xBW2 are plotted against the corresponding MP2 values. Here, each point corresponds to an electron gas with N_e = 14–1598 electrons in a cubic supercell, at a density corresponding to the Wigner–Seitz radius of 1 a.u. The calculations are performed with a finite plane-wave basis, as detailed in Ref. 20. For comparison, we also include data for the “mosaic” coupled-cluster doubles method (mCCD) and random phase approximation (RPA) with screened second-order exchange (RPA+SOSEX) from that reference.

RPA+SOSEX represents a canonical example of a convergent theory in this context. While the per electron correlation energy continues to increase with the size of the supercell for MP2, RPA+SOSEX quickly approaches a nearly constant value. For our purposes, the mCCD method is also an interesting benchmark. Here, the amplitude equations of CCD are reduced to the pure driver term (also present in MP2) and the “mosaic” diagrams, which can be interpreted as a renormalization of the single-particle Hamiltonian.¹⁰ From this perspective, mCCD, thus, resembles an ab initio analog to the regularized MP2 methods studied herein. It can be seen that mCCD also deviates strongly from MP2 for larger supercells, indicating an apparent convergence (or at least slower divergence).

κ-MP2 displays a significantly improved behavior relative to MP2, somewhat resembling the mCCD behavior though with different rates of convergence. This is quite remarkable given that the regularization parameter κ = 1.4 $EH−1$ was empirically optimized to small molecule thermochemistry data, which is completely unrelated to the HEG. This could indicate that a somewhat stronger regularization might be adequate for metallic systems but that the functional form of κ-MP2 is, in principle, well suited for this type of system.

In contrast, xBW2 only displays a marginal improvement over MP2. It is, therefore, not trivially the case that regularization allows applying MP2-like methods to metallic systems. Indeed, it is worth reemphasizing that MP2 diverges for finite-sized HEG supercells, even though these do not display a vanishing bandgap. Therefore, the simplistic notion that fixing MP2 in the $Δijab→0$ limit is sufficient to describe metals is not correct. Instead, higher-order screening effects are essential for describing uniform electron densities with or without gaps. These effects are apparently captured by the κ-MP2 regularizer to some extent but not by xBW2.

Next, we investigate the performance of regularized MP2-like methods for the strongly correlated Hubbard model. The Hubbard model and related Hamiltonians were independently proposed in physics by Hubbard,⁴¹ Kanamori,⁴² and Gutzwiller⁴³ and in chemistry by Pariser, Parr, and Pople.^44–46 These model Hamiltonians cover the behavior of a wide range of correlated systems, such as high T_c superconductivity,⁴⁷ magnetism,^42,43 and the Mott metal–insulator transition.^48,49

The Hubbard Hamiltonian reduces the electron interactions in extended systems to a short-range repulsion U ≥ 0 on a lattice of single orbital sites. The nearest-neighbor sites are connected by hopping matrix elements t. Therefore, the physics of the model is governed by the correlation strength U/t. In the case of U ≫ t, the energy penalty U for a double occupancy of sites outweighs the kinetic energy gain t and the system becomes strongly correlated. In this limit, the performance of various many-body methods has been tested. These include exact diagonalization schemes, such as the Lanczos algorithm,⁴⁷ as well as approximate many-body methods, such as the random phase approximation (RPA),¹¹ truncated CC methods,^11,50,51 dynamical mean-field theory (DMFT),⁴⁹ and Quantum Monte Carlo (QMC) methods.^47,51 To evaluate the performance of correlated many-body methods, the one-dimensional half-filled model is especially instructive, as it can be solved exactly, revealing that it is insulating in the strongly correlated limit.⁴⁸ In the following, we, therefore, examine the one-dimensional spinless periodic six-site Hubbard model at half-filling.

Figure 4 depicts the ground-state energy curve of xBW2 as a function of U/t. Here, we compare xBW2 to MP2, Variational Coupled Cluster with Double (VCCD) excitations, and the exact full Configuration Interaction (FCI) method. VCCD is included here as it represents the best possible energy that can rigorously be obtained with an MP2-like wavefunction.⁵² Indeed, VCCD only slightly underestimates the correlation energy in the strongly correlated limit, while MP2 diverges. As with the HEG, this is despite the fact that the energy gap retains its finite value for all correlation strengths. Interestingly, xBW2 displays a massive improvement over MP2, essentially curing the strongly divergent behavior of the latter. The xBW2 curve, in fact, displays excellent quantitative agreement with the reference methods, somewhat fortuitously falling between the VCCD and FCI curves.

An analogous plot is shown for κ-MP2 in Fig. 5. As the Hubbard model Hamiltonian uses an arbitrary energy scale given by the hopping parameter t, there is no meaningful way to translate the empirically optimized value of κ to this system. We, therefore, consider a range of regularization strengths. For κ → ∞, the original diverging MP2 curve is recovered, while for κ → 0, the restricted Hartree–Fock (RHF) curve is obtained, as the correlation energy contribution vanishes. Taking these limits of κ into account, several intermediate values were considered, so that the regularized calculations reproduced the FCI energy at the correlation strengths U/t = 8, 13, 18, respectively. Interestingly, none of these methods outperforms xBW2, despite the fact that they were explicitly fitted to the FCI results. Indeed, all κ-MP2 curves display significant curvature, so that they are either over- or underregularized in some range and ultimately diverge. Notably, RPA completely fails for this system, yielding complex correlation energies beyond correlation strength U/t = 2.¹¹ 

Finally, we test the regularized methods for the interaction of graphene with a single water molecule. This system has been intensively studied as a highly challenging benchmark for many-body methods.^53–59 As a case in point, MP2 significantly overestimates this interaction, as is commonly observed for non-covalent interactions involving large, polarizable systems.¹² Consequently, only computationally demanding many-body treatments such as the coupled cluster method with singles, doubles, and perturbative triples [CCSD(T)] or QMC offer predictive accuracy on this system.

Table I presents the corresponding interaction energies, calculated at various levels of theory for a 4 × 4 graphene supercell (see  Appendix A for computational details). Note that these calculations are analogous to the ones in Ref. 53. However, since we here focus on benchmarking different electronic structure methods, finite-size effects beyond the 4 × 4 supercells are not taken into account. The methods are thus directly compared for a finite supercell, as in the case of the HEG.

CCSD(T) provides an accurate benchmark for this system. As expected, canonical MP2 significantly overbinds (by 18 meV, corresponding to 22.5%), due to the absence of screening effects. Perhaps surprisingly, coupled cluster with singles and doubles only (CCSD), in turn, underestimates the interaction energy by 25 meV. We also include results for the direct Random Phase Approximation (dRPA) using Hartree–Fock (HF) and Perdew–Burke–Ernzerhof (PBE) reference determinants. These turn out to be quite similar to CCSD in this case, though dRPA@PBE represents a slight correction in the right direction. This underscores the challenging nature of this system, which even infinite order methods like RPA and CCSD struggle with. It should be noted, however, that the RPA results can be significantly improved by including the contribution from GW single excitations.⁵³ 

As in the previous cases, regularization consistently leads to an improvement over the MP2 results. In the case of xBW2, this only amounts to a minor adjustment, however, as the interaction energy is merely lowered by 4 meV. In contrast, the improvement for κ-MP2 is more substantial yielding interaction energy within 7 meV of the CCSD(T) value.

These results indicate that higher-order screening effects can (to some extent) be captured by the κ-MP2 regularizer and less effectively by xBW2. This mirrors the findings for the HEG case discussed above. Indeed, a relation between the HEG and dispersion interactions in large polarizable systems has also been discussed in Ref. 12, in the context of the RPA method.

In this work, we have explored the potential of regularized second-order perturbation theories for predicting the correlation energies of extended systems, using κ-MP2 and the newly proposed xBW2 as representative examples. While neither of these methods is a silver bullet that cures all deficiencies of canonical MP2, regularization consistently leads to an improvement. In particular, κ-MP2 appears well suited to describe screening effects in metallic and polarizable systems. Meanwhile, the xBW2 method provides a remarkably effective remedy for the breakdown of MP2 in the strongly correlated limit of the one-dimensional Hubbard model.

Notably, we have not attempted the empirical adjustment of the regularization to the scrutinized systems (except in the inevitable case of κ-MP2 for the Hubbard model). Instead, the literature reported κ parameter (adjusted to molecular thermochemistry) was used, while xBW2 is a non-empirical method (though it could be converted into an empirical one by scaling the regularization term). However, as recently shown by Head-Gordon and co-workers, a stronger regularization improves the description of dispersion interactions and transition metal thermochemistry with κ-MP2.⁶⁰ Fitting the regularization parameter to a representative set of condensed matter systems would thus certainly lead to even better performance. Indeed, even a system-specific choice of the regularization parameter would be possible, in analogy to the DFT + U method. However, this would require an objective and transferable protocol for determining this value.

Finally, it should be noted that a Hartree–Fock reference determinant has been used throughout. This is the canonical reference for MP2, but not necessarily the optimal choice.⁶¹ In practice, Kohn–Sham or self-consistently optimized orbitals are often found to be more accurate for molecular systems.^23,24,62,63 We may expect the same for extended systems, as can be seen for the dRPA results in Table I. This will be explored in future work.

See the supplementary material for individual MP2 and xBW2 reaction energy data for Fig. 2.

E.K. acknowledges support from the IMPRS for Elementary Processes in Physical Chemistry.

The authors have no conflicts of interest to disclose.

The data that support the findings of this study are available within the article and its supplementary material.

xBW2 and BW2 were implemented in NumPy using pySCF to generate molecular integrals.⁶⁴ The cc-pVDZ basis was used for the helium chain calculations.⁶⁵ 

HEG calculations are performed with a modified version of the UEGCCD code.¹⁰ 

Calculations on the Hubbard models were performed with custom NumPy implementations of (regularized) MP2 and VCCD. FCI calculations were performed with pySCF.⁶⁴ 

Water physisorption is considered in the two-leg geometry as discussed in Ref. 53. The interaction energy E_int at the equilibrium distance d is obtained as

where E(d) is the total energy of the system with the water at a distance d = 3.37 Å from the graphene monolayer, and E(d_far) is the energy of the non-interacting system, with the water molecule at a distance d_far = 7.395 Å from the graphene layer. A 4 × 4 graphene sheet is employed, containing 32 carbon atoms. The simulation cell contains a vacuum gap of 14.79 Å, to ensure that the monolayer does not interact with its periodic image. We stress that since we are primarily interested in benchmarking different regularized MP2 methods, finite-size effects beyond the 4 × 4 graphene supercell are neglected. HF orbitals are expanded in a plane-wave basis within the PAW framework using cutoff energy of 500 eV, as implemented in the Vienna Ab initio Simulation Package (VASP).^66–69 For the correlated calculations, we use frozen natural orbitals (FNOs) computed at the MP2 level.^70,71 The number of virtual orbitals is truncated by selecting the N_v FNOs with the largest occupation number.

Recently, Irmler et al.⁷² proposed a simple approximation to correct for the basis set incompleteness error (BSIE) of CCSD. This effectively corresponds to a rescaled pair-specific MP2 term. Herein, we correct for the BSIE of correlated methods using the CBS limit of MP2, such that

where the superscript denotes the total number of virtual orbitals N_v, and E_c is the correlation energy calculated within the different methods. All many-body calculations are performed with the CC4S code.⁷³ RPA calculations employing MP2 FNOs were performed with VASP using N_v = 6000 orbitals. Additional RPA calculations with PBE orbitals were performed using the full plane-wave basis with a 500 eV energy cutoff, and the results were extrapolated to the infinite basis set limit using the internal VASP extrapolation.⁷⁴ 

As indicated in the main manuscript, xBW2 was derived as a size-extensive modification of BW2, i.e., by identifying the root of its non-extensive behavior and renormalizing the offending term. To provide a more general perspective, we here consider how xBW2 can be derived on an equal footing with MP2 and BW2 from formal perturbation theory.^75,76 To this end, we begin by splitting the Hamiltonian $Ĥ$ into a zeroth-order approximation $Ĥ0$ and a perturbation $V̂$⁠,

where $Ĥ0$ is chosen as the sum of Fock operators. For a single-determinant reference wavefunction ϕ₀, we, thus, obtain

and

with $fii=⟨χi|Ĥ0|χi⟩$⁠.

In the following, we use intermediate normalization, so that

for the exact many-body wavefunction Ψ, which implies that

Here, the determinants ϕ_k form the orthogonal complement to ϕ₀ and c_k are the corresponding expansion coefficients. Note that the determinants ϕ_k are typically grouped according to their excitation rank with respect to ϕ₀, i.e., into singles, doubles, triples, etc. We will assume that all ϕ_k are also eigenfunctions of $Ĥ0$ (i.e., that $Ĥ0ϕk$ = E_k,0ϕ_k).

By projecting ϕ₀ onto the Schrödinger equation, we obtain

with the exact ground-state energy E.

Having established these basic properties of $Ĥ0$ and $V̂$⁠, we can proceed to construct a perturbative expansion of Ψ and E. To this end, we define a projection operator as

and rewrite the Schrödinger equation as

where E_LS is an arbitrary level-shift parameter (more on this below).

We can now define the perturbative corrections to the energy and wavefunction as

which are given by the following equations:

From this, the well-known first-order energy expression can be obtained,

Regardless of how E_LS is selected, the first-order energy correction is thus always the difference between the total energy of ϕ₀ (E_HF) and the sum of occupied orbital eigenvalues (E₀).

Introducing the resolvent operator $R̂0$ as

we can express the second-order energy corrections as

Here, the latter two terms vanish (because $R̂0ϕ0=0$⁠), yielding

We can now return to the nature of the arbitrary level-shift parameter E_LS. Different choices for this parameter can be used to obtain different perturbation theories, with different formal and convergence properties. Most prominently, we can set E_LS = E₀ to obtain Møller–Plesset perturbation theory. With this choice (and assuming that ϕ₀ is the canonical Hartree–Fock determinant), we obtain the well-known MP2 correlation energy expression given in Eq. (1).

The other obvious choice is to use E_LS = E, which yields the BW series. Here, additional approximations are required, however, since the exact ground-state energy E is unknown. Common choices are E ≈ E₀ + E⁽¹⁾ + E⁽²⁾ and E ≈ E₀ + E⁽²⁾, the latter of which leads to the BW2 expression in Eq. (3). From this analysis, it becomes clear the xBW2 correlation energy expression corresponds to the choice $ELS=E0+E(2)Ne$⁠.

While this derivation is somewhat post hoc, it opens interesting perspectives, allowing the derivation of equations for general reference determinants and higher-order corrections. This will be explored in future work.