Optimizing the regularization in size-consistent second-order Brillouin-Wigner perturbation theory

Formally, MP2 has many desirable properties. For example, it is free of delocalization errors, unlike the widely popular density functional theory (DFT).^1–3 In contrast to both DFT and the direct random phase approximation (RPA),⁴ there is no self-correlation error. Consistent with Pople’s high standards for an approximate model chemistry⁵ (at the heart of which are formal principles which are practically useful in chemical predictions), MP2 is size-consistent, size-extensive, and orbital invariant.⁶ A model is size-consistent if the total energy of a supersystem comprised of noninteracting subsystems is the same as the sum of the energies of the isolated subsystems; this is an essential property when studying phenomena such as bond breaking. Second, a method is size-extensive if the total correlation energy in a linear chain of atoms grows linearly with number of electrons, which is essential for reaching the thermodynamic limit. Third, a method that yields the same correlation energy despite arbitrary orbital rotations in the occupied (or virtual) subspace is considered to be orbital invariant – a property that enables transformations to chemically-relevant bases such as the natural orbital or localized orbital representations.

MP2 routinely outperforms Hartree–Fock (HF) theory across myriad test sets with respect to experimental or near-exact numerical reference values.⁷ The accuracy of MP2 can be very high in the case of closed-shell and small organic molecules, and can exceed the accuracy of popular DFT functionals for important chemical properties such as reaction barrier heights (which are sensitive to delocalization errors).⁸ Indeed, MP2 is the most popular wavefunction component to be incorporated into double-hybrid density functionals, with promising results in many chemically-relevant situations.^9–16 

However, over the years many shortcomings of MP2 have been found. It is well-known that perturbation theory in general is not suitable for multi-reference states, in which higher order (connected) excitations are required for a qualitatively correct description of the wavefunction. In addition, MP2 (and even higher orders of perturbation theory) can fail in certain cases where the reference determinant is severely spin-contaminated.^17–22 In strongly correlated cases, this is a fatal issue; however, in weakly correlated systems, where the spin-symmetry breaking is artificial (i.e., due to deficiencies in the model chemistry’s treatment of dynamic correlation), MP2 with a restricted open-shell (RO) reference determinant can at times remedy this situation.²¹ Approximate Brueckner orbital approaches,^23,24 e.g., orbital optimized (MP2) methods,^25–32 can also clean up spin-symmetry breaking at the level of the Hartree–Fock orbitals.

Despite such notable improvements over conventional MP2, energy-gap dependent protocols for MP2 regularization lack the desired level of transferability required to be widely used in a black-box fashion. For example, with regularization parameters optimized for NC and transition metal systems, the accuracy for main-group thermochemistry (TC) and electronic response properties is notably deteriorated – at times these regularized MP2 approaches are worse than conventional MP2 by factors of 2 or 3.⁷ Similarly, κ-MP2 demonstrated very promising improvements in accuracy relative to MP2 for nuclear magnetic resonance (NMR) chemical shifts only when element-specific κ values were employed.³⁸ The prospect of developing a more transferable approach to regularized second-order perturbation theory, which preserves high accuracy for NC and TM datasets, is the primary motivation for the present work.

Brillouin-Wigner perturbation theory (BWPT)^39–42 is an alternative to the Rayleigh-Schrödinger approach (the latter gives rise to MP2). We recently proposed a size-consistent variant of second-order BWPT, which naturally regularizes the t-amplitudes by shifting the occupied orbital energies in the denominator to lower values, thus increasing the effective orbital energy gaps and damping artificially overestimated amplitudes.⁴³ The single free parameter in our BW-s2 approach was determined such that the dissociation limit of a system with two electrons in two orbitals (e.g., the H₂ molecule in a minimal basis set) is exact. Importantly, this model, which we refer to as BW-s2, is size-consistent, size-extensive, and orbital invariant. While BW-s2 was found to be less accurate than (optimally parameterized) κ-MP2 in cases where exceptionally strong regularization was required, for a wide variety of main group TC its performance is superior to κ-MP2 and conventional MP2. In this work we aim to explore the landscape of the free parameter, which we will call α, by investigating many different data sets representative of NC, large-gap TMTC, main-group TC, barrier heights, and molecular dipoles and polarizabilities.

In our original set of benchmarks,⁴³ we found that BW-s2 consistently outperforms MP2 across myriad chemical problems, which is very encouraging. However, it was evident that specific, optimal choices of κ in κ-MP2 could significantly outperform BW-s2 in problems where strong regularization was required (such as transition metal thermochemistry). How much improvement is possible if we lift the restriction of α = 1, and instead view α as a parameter that controls regularization strength? That is the question that we will investigate here.

In this work, we benchmark the performance of various values of α against a variety of data sets in an effort to tune the accuracy of BW-s2 [henceforth, the empirical variant will be referred to as BW-s2(α)]. Notably, the particular value of α does not influence the size-consistency of the method, but it may be a determining factor in the overall quality of the results. We will assess the transferability of the α parameter across various chemical problems, and attempt to make a recommendation for a broadly applicable α value.

The results for all benchmark sets apart from electronic properties are shown in Table I and are plotted individually as a function of α for each data set in Figs. S6–S9. These data include NC for sets of small dimers such as A24,⁵⁴ S22,⁵⁵ S66,⁵⁶ and the non-I-containing subset of X40 (hereafter referred to as X31),⁵⁷ along with the large π-stacked dimers of L7.⁵⁸ TC is assessed on H-atom transfer (HTBH38) and non-H-atom transfer (NHTBH38) sets,^59,60 along with the more comprehensive single-reference subset of W4-11.⁶¹ As compared to our original work, we extend our coverage of TMTC with reaction energies from MOR39³⁶ (a subset of MOR41 with triple-ζ reference values),⁴³ MC09,⁶² and a set of 13 Au, Pt, and Ir reaction energies that we call AuPtIr13.⁶³ Finally, we also include the ACONFL set of relative alkane conformational isomer energies.⁶⁴ 

The MP2 results always improve on those obtained from HF, and gap-regularized κ-MP2 improves further on these results in all benchmark sets, apart from TC where the results degrade by up to 3.2 kcal/mol. A similar trend emerges for BW-s2(α), with noticeable improvements over MP2 for NC, TMTC, and ACONFL data sets, but the results for barrier heights degrade by only about half as much as κ-MP2 (0.03–1.9 kcal/mol less accurate for modest parameters in the range 1 ≤ α ≤ 4). Furthermore, BW-s2(α) performs roughly 1 kcal/mol better than MP2 on the W4-11 benchmark set regardless of the particular value of α, whereas κ-MP2 performs slightly (0.7 kcal/mol) worse. The improvements in NC, TMTC, and ACONFL sets with minimal degradation in the results for TC suggest that the BW-s2(α) α-parameter is more transferable than the κ in κ-MP2.^7,43

Regarding the transferability argument, it is instructive to consider electronic properties such as dipole moments and polarizabilities that are shown in Table II. Whereas κ-MP2 doubles the errors relative to MP2 for both dipoles and polarizabilities, BW-s2(α) exhibits an exceptional flatness in the errors as a function of α, peaking at 4% for the most severe α = 8.0 where the errors nonetheless remain lower than κ-MP2. Reference 7 reports κ-MP2 errors (1.6 ≥ κ ≥ 1) for dipoles that span the range 4.7%–7.5% while polarizability errors span 4.2%–5.9%. Not only are these close to the largest errors that we report for BW-s2(α), but BW-s2(α = 1) actually improves the results for dipole moments relative to MP2, whereas κ-MP2 errors monotonically increase as κ decreases. These results for electronic properties suggest that the α parameter in BW-s2(α) is indeed much more transferable between classes of chemical problem than gap-dependent regularizers.

Some particularly interesting highlights are that BW-s2(α) can reduce errors relative to MP2 in the L7 data set from 9.5 to 1.3 kcal/mol. TMTC data can also be improved by a factor of 2–3 relative to the MP2 results, reducing errors from 14 kcal/mol to 4–6 kcal/mol for MOR39 and MCO9 sets if moderate to large α parameters are applied. Finally, errors in alkane conformational energies can be reduced from $∼1$ kcal/mol with MP2 to just 0.1 kcal/mol with BW-s2(α), achieving something close to chemical accuracy. Of course, excellent performance for particular kinds of chemical problem does not suggest a “universal” α value, and there is likely no α parameter that is entirely satisfactory in all chemical contexts. However, we make the recommendation of α = 4 based on the analysis presented alongside Fig. 1. Taking a closer look, the BW-s2(α = 4) error statistics suggest considerable improvements relative to MP2 for NC, main-group TC (W4-11), TMTC, and ACONFL sets, while minimal damage is done to the results for H-atom/non-H-atom transfer barrier heights and electronic properties.

While there is no universal parameter, BW-s2(α) stands out from gap-dependent regularizers like κ-MP2 (and the similarly-performing σ-MP2 and σ²-MP2 methods)⁷ in the sense that it is clearly more transferable across different chemical problems. This may be due to the fact that BW-s2(α) defines a valid second order BW correction for each α. As a consequence it incorporates the full set of t amplitudes in the regularizer, whereas gap-dependent schemes rely only on the orbital energy gaps. The self-consistent nature of BW-s2(α) may also act to further refine the orbital energy gap, introducing a feedback loop that fine-tunes the resultant amplitudes.

As a final test for the robustness of our parameterization, we consider a secondary free parameter, β, that directly modulates the amount of BW-s2(α) correlation energy such that, E = E_HF + βE_BW-s2(α). The results in Sec. S1 show that the optimal β parameter generally hovers in the range 0.9 ≤ β ≤ 1.1. Furthermore, when α nears its optimal value, β → 1.0 with the exception of non-H-atom barrier heights in NHTBH38 (Fig. S2) where β = 1.1 when α = 1. A β > 1 implies systematic under-correlation, and points to an optimal α for NHTBH38 that is less than 1. In stark contrast to this, the landscape of the parameter space for TMTC in Fig. S3 features an optimal β = 0.7 at low α = 1, which increases to β = 1 only when α → 8. This implies a significant over-correlation for transition-metal systems that is tempered only by larger α parameters.

The NC, W4-11, and ACONFL data sets in Figs. S1, S2, and S4, respectively, show a relatively flat slope defined by the line tracing $minα,βError(α,β)$⁠. For these sets, the optimal β is very close to 1 across α parameters, suggesting that BW-s2(α) offers a balanced description of correlation for NC, main-group TC, and conformational isomers. Overall, the relatively low slopes across the parameter space and the proximity of β to 1 across various α both speak to the transferability of the BW-s2(α) approach. Thus, moving forward we suggest the single parameter BW-s2(α = 4) approach for general chemical applications.

All calculations were performed in a development version of Q-Chem v6.0.2.⁶⁶ All calculations (aside from evaluations of electronic properties) feature SCF convergence thresholds that were set to 10⁻⁸ root-mean-square error. The correlation energy was considered to be converged at a change of 10⁻⁸ Ha between iterations for all calculations except for those of the L7 dataset, where this was relaxed to 10⁻⁵ Ha. Relevant derivatives with respect to electric fields for properties such as dipoles and polarizabilities were evaluated via finite difference. Because finite difference results are especially sensitive to numerical errors, the SCF convergence and correlation energy thresholds were set to 10⁻¹¹. To achieve complete basis set limit extrapolations for NC, TC, and ACONFL we follow the protocol in Ref. 67, which has been verified to perform well with the heavy-aug-cc-pVDZ/heavy-aug-cc-pVTZ basis sets used for L7.⁶⁸ For electronic response properties, we use the same extrapolation method reported in Ref. 7.

Since W depends on the t amplitudes, which themselves depend on the modulation of the energy gap supplied by the W matrix, the BW-s2 equations must be solved self-consistently. We begin each BW-s2 calculation with canonical Hartree–Fock orbitals and an MP2 guess for the initial t amplitudes, though we note the possibility of obtaining a strictly non-divergent initial guess by means of Davidson’s repartitioning of the one-electron Fock operator.⁶⁹ To accelerate these calculations, our implementation uses the resolution-of-the-identity (RI) approximation for the two-electron integrals,^70,71 resulting in a formal scaling of $m×O(N5)$⁠, where m is the number of iterations (typically between 4 and 6) and N is the number of basis functions. Due to computational limitations, the I functions in the auxiliary RI basis sets were removed for transition-metal calculations in the MOR39, MCO9, and AuIrPt13 data sets.

Additional figures pertaining to scans over α, and the secondary β parameter in BW-s2(α). Detailed data are provided for NC, main-group TC, TMTC, ACONFL, and electronic response properties.

This work was supported by the Director, Office of Science, Office of Basic Energy Sciences, of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. K.C.-F. acknowledges support from the National Institute Of General Medical Sciences of the National Institutes of Health under Award No. F32GM149165. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Martin Head-Gordon is a part-owner of Q-Chem, which is the software platform used to perform the developments and calculations described in this work.

Kevin Carter-Fenk: Data curation (lead); Investigation (lead); Methodology (equal); Software (lead); Writing – original draft (equal); Writing – review & editing (equal). James Shee: Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Martin Head-Gordon: Conceptualization (lead); Funding acquisition (lead); Supervision (lead); Writing – review & editing (equal).

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