Regularized Second-Order Correlation Methods for Extended Systems

While many-body wavefunction theory has long been established as a powerful framework for highly accurate molecular quantum chemistry, these methods have only fairly recently been applied to extended systems in a significant scale. This is due to the high computational cost of such calculations, requiring efficient implementations and ample computing resources. To further aggravate this, second-order Møller-Plesset perturbation theory (MP2) (the most cost effective wavefuntion method) is known to diverge or fail for some prototypical condensed matter systems like the homogeneous electron gas (HEG). In this paper, we explore how the issues of MP2 for metallic and strongly correlated systems can be ameliorated through regularization. To this end, two regularized second-order methods (including a new, size-extensive Brillioun-Wigner approach) are applied to the HEG, the one-dimensional Hubbard model and the graphene-water interaction energy. We find that regularization consistently leads to improvements over the MP2 baseline and that different regularizers are appropriate for metallic and strongly correlated systems, respectively.


I. INTRODUCTION
Second-order Møller-Plesset perturbation theory (MP2) holds a unique place in the hierarchy of wavefunctionbased electronic structure methods. 1-4 Historically, it was the first correlated method in quantum chemistry that was sizeextensive, invariant to unitary orbital rotations and computationally affordable. Until the development of modern density functional theory (DFT) it was therefore the workhorse method in molecular quantum chemistry. 5 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][7][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][10][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 endeavours. [14][15][16][17][18][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 vanishing band-gaps 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 behaviour can be interpreted as a renormalization of the band-gap 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][23][24][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, which each constitute known issues for canonical MP2. In particular, we consider the homogeneous electron gas (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 second-order Brillouin-Wigner perturbation theory.

II. 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 i j||ab . Using this notation, the MP2 correlation energy can be written as: where the denominator ∆ ab i j = ε a + ε b − ε i − ε j is computed from the corresponding orbital energies ε.
Clearly, this energy must diverge if any ∆ ab i j becomes zero, i.e. for metallic systems. In principle a small constant δ could be added to the denominator to avoid this. 23 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. 23 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 ∆ ab i j . 24 This allows attenuating the offending terms without affecting the well behaved ones. The most successful of these schemes is κ-regularization, which uses the expression: (2) For orbital-optimized MP2, this form of regularization (using κ=1.4 E −1 H ) 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 know failure of canonical MP2). 24 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: 26 This equation must be solved iteratively, as E BW2 c appears on both the left and right-hand sides.
This method 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 ∆ ab i j . 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. 30 This is obviously problematic for applications to extended systems. Nonetheless, the idea of using a correlation energy dependent regularization is attractive, since 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 sizeextensivity of MP2 is recovered. This is shown numerically for He chains in Fig. 1.
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. 31,32 In this sense, xBW2 can be seen as a secondorder 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). 33 FIG. 1. Correlation energies per atom for evenly spaced (d=3 Å) Helium chains of increasing length, using the cc-pVDZ basis. A slope larger than zero indicates a non-size-extensive method.
In this method, pair correlation energies are used for regularization instead of the correlation energy per atom. Consequently, each electron-pair i j has a different regularizer. In this sense, BGE2 is related to the coupled electron pair approximation (CEPA) 34 and xBW2 can be considered an 'averaged' BGE2. This averaging has the advantage that (unlike BGE2) xBW2 is invariant to unitary orbital transformations, although the BGE2 regularization is arguably more flexible. More generally, the interative 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 ∆ ab i j ) and xBW2 (representing a non-empirical, energy-dependent regularization scheme).

III. RESULTS
Homogeneous Electron Gas: 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,35,36 More recently, finite periodic electron gases have emerged as important numerical benchmark systems for many-body theories. 10,20,37,38 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. 10 In Fig. 2, 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 atomic unit. 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 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. 10 From this perspective, mCCD thus resembles an ab initio analogue 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 behaviour relative to MP2, somewhat resembling the mCCD behaviour though with a different rate of convergence. This is quite remarkable given that the regularization parameter κ=1.4 E −1 H 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 finitesized HEG supercells, even though these do not display a vanishing band gap. Therefore, the simplistic notion that fixing MP2 in the ∆ ab i j → 0 limit is sufficient to describe metals is not correct. Instead, higher-order screening effects are essen-tial 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. One-Dimensional Hubbard Model: 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 39 , Kanamori 40 and Gutzwiller 41 and in chemistry by Pariser, Parr and Pople [42][43][44] . These model Hamiltonians cover the behaviour of a wide range of correlated systems such as high T c superconductivity 45 , magnetism 40,41 and the Mott metal-insulator transition 46,47 .
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-neighbour 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 45 , as well as approximate many-body methods such as the random phase approximation (RPA) 11 , truncated CC methods 11,48,49 , dynamical mean-field theory (DMFT) 47 and Quantum Monte Carlo (QMC) methods 45,49 . To evaluate the performance of correlated many-body methods, the one-dimensional half-filled model is especially instructive, as it can be solved exactly. This solution shows that the one-dimensional half-filled Hubbard model is insulating in the strongly correlated limit. 46 In the following, we therefore examine the one-dimensional spinless periodic six-site Hubbard model at half-filling. Figure 3 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 excitations (VCCD) 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. 50 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 behaviour 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. 4. 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 repro-  duced 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. Graphene-Water Interaction: 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. [51][52][53][54][55][56][57] As a case in point, MP2 significantly overestimates this interaction, as is commonly observed for non-covalent interactions involving large, polarizable systems. 12 Consequently, only computationally demanding many-body treatments at the CCSD(T) or QMC level offer predictive accuracy on this system. Table I shows the corresponding interaction energies, calculated at various levels of theory for a 4 × 4 graphene supercell (see Appendix A for computational details). Here, CCSD(T) provides an accurate benchmark. As expected, canonical MP2 significantly overbinds (by 18 meV corresponding to 22.5 %), due to the absence of screening effects. Perhaps surprisingly, the CCSD method in turn underestimates the interaction energy by 25 meV, underscoring the challenging nature of this system. 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 an 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 reference [12], in the context of the RPA method.

IV. CONCLUSIONS
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. 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). Fitting the regularization parameter to a representative set of condensed matter systems would certainly lead to even better performance. Indeed, even a systemspecific 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. 58 In practice, Kohn-Sham or self-consistently optimized orbitals are often found to be more accurate for molecular systems. 23,24,59,60 We may expect the same for extended systems. This will be explored in future work. ACKNOWLEDGMENTS We wish to acknowledge the support of the International Max Planck Research School for Elementary Processes in Physical Chemistry.

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

V. APPENDIX A: COMPUTATIONAL DETAILS
Helium Chains: xBW2 and BW2 were implemented in NumPy using pySCF to generate molecular integrals. 61 The cc-pVDZ basis was used for the Helium chain calculations. 62 Homogeneous Electron Gas: HEG calculations are performed with a modified version of the UEGCCD code. 10 Hubbard Model: Calculations on the Hubbard Models were performed with custom NumPy implementations of (regularized) MP2 and VCCD. FCI calculations were performed with pySCF. 61 Water-Graphene Interaction: Water physisorption is considered in the 2-leg geometry as discussed in Ref [51]. 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, with a vacuum gap of 14.79 Åto ensure that the monolayer does not interact with its periodic image. HF orbitals are expanded in a plane-wave basis within the PAW framework using a cutoff energy of 500 eV, as implemented in the VASP code. [63][64][65][66] For the correlated calculations we use frozen natural orbitals (FNOs) computed at the MP2 level. 67,68 . The number of virtual orbitals is truncated by selecting the N v FNOs with the largest occupation number.
Recently, Irmler et al. 69 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. 70