Analysis of the errors in explicitly correlated electronic structure theory

Abstract

The explicitly correlated second order Møller–Plesset (MP2-R12) methods perform well in reproducing the last detail of the correlation cusp, allowing higher accuracy than can be accessed through conventional means. Nevertheless in basis sets that are practical for calculations on larger systems (i.e., around triple- or perhaps quadruple-zeta) MP2-R12 fails to bridge the divide between conventional MP2 and the MP2 basis set limit. In this contribution we analyse the sources of error in MP2-R12 calculations in such basis sets. We conclude that the main source of error is the choice of the correlation factor r₁₂. Sources of error that must be avoided for accurate quantum chemistry include the neglect of exchange commutators and the extended Brillouin condition. The generalized Brillouin condition is found not to lead to significant errors.

---

1. Introduction

It is troubling that the success of the Hylleraas calculation on helium^1–3 has proven so difficult to carry over to more complicated cases. Hylleraas solved the problem of the slow convergence of his configuration interaction expansions for helium⁴ by including terms that depended explicitly on the distance between the two electrons.¹ For larger atoms or molecules such methods become rapidly impractical because integrals over the coordinates of many electrons always appear, and these are both extremely hard to evaluate individually, and extremely numerous.

Despite countless innovations in the field of explicitly correlated electronic structure theory^5–8 there is really only one class of methods that can truly claim to have escaped the restriction to tiny molecules. These methods—now collectively labelled ‘R12’—were introduced in 1985 by Kutzelnigg⁹ and their success pivots on the ingenious use of the resolution of the identity (RI) to avoid the difficult many-electron integrals.⁹

The RI is based on 1-electron projection operators of the form

P̂′ = |p′〉〈p′| ≈ 1̂ (1)

where {p′} is a suitably chosen orthonormal basis, possibly the molecular orbital (MO) basis,⁹ but probably somewhat larger.^10,11 Here and throughout summation over repeated dummy indices is assumed. Such RI operators can be used to approximate many-electron integrals by sums of products of 2-electron integrals: for example,

〈ijm|r⁻¹₁₂r₂₃|mlk〉 ≈ 〈ijm|r⁻¹₁₂P̂^′₂r₂₃|mlk〉

This concept can be applied at a variety of levels of electronic structure theory, leading to the R12 class of methods.^12–19

At the level of second order Møller–Plesset perturbation theory (MP2) the first order pair function is written

|u_ij〉 = t^ij_ab|ab〉 + t^ij_klQ̂₁₂f₁₂|kl〉 (3)

where the correlation factor in the original methods was f₁₂ = r₁₂. The operator Q̂₁₂ ensures the strong orthogonality of the pair function to the occupied space, and various forms can be used to achieve this.^10,11,20,21 These include

Q̂^P₁₂ = (1 − P̂₁)(1 − P̂₂)

where Ô, V̂ and P̂ project, respectively onto the occupied, virtual and full MO spaces.

The last two of these operators, Q̂^O₁₂ and Q̂₁₂, lead to formally equivalent first order ansätze; the first is distinct. In the particular case where the RI is performed in the MO basis all of these operators lead to the same working equations.¹¹ Since details can be found elsewhere^11,13,20 we will not examine the equations that arise in each case; we will however identify those approximations that are commonly made in deriving the working equations.

First it is convenient to assume that the Fock equations have been solved exactly for the occupied space: f̂|i〉 = ε_i|i〉. This approximation is referred to in the literature as the generalized Brillouin condition (GBC),¹³ and allows the simplification of certain expressions by deleting commutators of the form [f̂₁ + f̂₂, Q̂^O₁₂]. If either Q̂^P₁₂ or Q̂₁₂ are used, the analogous commutators can only be eliminated under the more stringent assumption of the extended Brillouin condition (EBC),¹³ namely that the Fock equations hold exactly for all MOs: f̂|p〉 = ε_p|p〉. In this case the EBC has the additional benefit that the coupling between the conventional and R12 parts of the first order wavefunction vanishes.

Developing the formalism one also find terms of the form [f̂₁ + f̂₂, f₁₂]. Clearly the commutators involving the local potential terms vanish, and one has the identity

[f̂₁ + f̂₂, f₁₂] = [t̂₁ + t̂₂ − K̂₁ − K̂₂, f₁₂].

Given the approximately local form of exchange (K̂ ≈ cρ^1/3) one can make the approximation [K̂₁ + K̂₂, f₁₂] = 0 to obtain

[f̂₁ + f̂₂, f₁₂] ≈ [t̂₁ + t̂₂, f₁₂]. (5)

Further arguments for and against the deletion of terms of the form [K̂₁, f₁₂] can be found around eqn. (5.14) of ref. 13. Those methods that include the exchange commutators are labelled B and those that neglect them are labelled A.¹³

These three approximations (GBC, EBC, [K̂₁, f₁₂] = 0) are present in various combinations in different flavours of MP2-R12 theory, as illustrated in Table 1. The methods that use the projector Q̂^P₁₂ (ansatz1) are less accurate and will not be considered here. Of the remaining methods there is an increasing computational cost as one moves down the table. The highest scaling process in MP2-R12/2*A′ scales like [scr O, script letter O](o⁸), but this typically remains modest compared with other O(N⁵) steps. However in MP2-R12/2A′, an additional step that scales O(o⁶v²) (eqn. (79) of ref. 11) is needed and this step quite quickly dominates the computational work. The original method MP2-R12/2B¹¹ adds to this problem by requiring 4-index integrals with two indices drawn from the (generally quite large) RI basis set. The latter problem has recently been circumvented.²²

Table 1 Various flavours of MP2-R12 theory, and the approximations made in them

Method	GBC	EBC	[K̂1, f12] = 0	SOP^a
a Strong orthogonality projector.
R12/1A′	✓	✓	✓	Q̂ ^P12
R12/2*A′	✓	✓	✓	Q̂ ^O12
R12/2A′	✓	✗	✓	Q̂ ^O12
R12/2B	✓	✗	✗	Q̂ ^O12

---

The R12 methods perform very well in describing the last detail of the correlation cusp, efficiently making up the difference between large-basis conventional methods and the basis set limit. For more extended systems, though, large atomic orbital (AO) basis sets become impractical, and one is restricted (at the moment) to basis sets of triple-zeta or at most quadruple-zeta size. Unfortunately the performance in these smaller basis sets is much less satisfactory. For neon in cc-pVTZ the MP2 and MP2-R12/2*A′ correlation energies are −264 mE_h and −292 mE_h, respectively, compared to the basis set limit −320 mE_h. It is dissappointing that MP2-R12 methods only make up about half of the energy gap to the basis set limit, and this is quite typical.

In the preceding text two distinct sources for this 30 mE_h error have been identified: (1) the RI for many-electron integrals; (2) the additional approximations GBC, EBC, [K̂₁, f₁₂] = 0; and to these we can add a third possibility, namely (3) the quality of the ansatz for the first order wavefunction.

In this paper we assess the relative impact of these three sources of error, by scrutinizing the existing literature and by performing new benchmark calculations, and produce clear recommendations for superior explicitly correlated MP2 methods. It is anticipated that the central conclusions for the MP2 level of theory will be applicable to higher levels of theory such as CCSD-R12^16,18 and CCSD(R12).¹⁹

Our benchmarks are based on valence correlation energies of the test set of 20 small first-row molecules from the textbook of Helgaker et al.²³ Most calculations have been performed with our DF-MP2-R12 and DF-MP2-F12 codes^21,24,25 using density fitting to compute all 4-index quantities. In such cases we have used the large cc-pV5Z density fitting basis sets of Hättig.²⁶ The errors arising from density fitting are much smaller than the other errors we are considering, as shown elsewhere.²⁴ Computations without the assumption of the EBC and GBC were carried out with the developmental version of the Massively Parallel Quantum Chemistry Program (MPQC).²⁷

2. The accuracy of 3-electron integrals

Until 2002 only two methods were widely used for the evaluation of 3-electron integrals: exact computation^5,6,20,28 and RI in the atomic orbital (AO) basis.^9,13 The former method is practical for molecules only when both the atomic orbitals and the correlation factor (or geminal function) are expanded in Gaussians (see for example refs. 8 and 29). The latter uses the identity operator in the AO basis to resolve (approximately) the 3-electron integrals into sums of products of 2-electron integrals, as in eqn. (2).

Over the past few years, though, there have been several developments. First amongst these in chronology as well as importance is the use of RI approximations in an auxiliary basis set¹¹ (ABS). This allows one to study the effect of the accuracy of the integrals independently, and, from a more pragmatic point of view, allows one to converge the accuracy of the many electron integrals whilst retaining a reasonably modest basis set for the molecular orbitals.

The idea has been extended to formulations that need only RIs in the orthogonal complement of the AO basis.¹⁰ Density fitting (DF) can also be used as an alternative to RIs for the 3-electron integrals³⁰ although this does not appear to offer significant advantages; however combined RI–DF approaches offer enhanced efficiency²⁴ and accuracy.³¹ Finally one can use numerical quadrature for the many electron integrals as shown by Boys and Handy⁷ and recently rediscovered by Ten-no.³²

Although we do not intend to benchmark all of these methods here, we will compare calculations using different ABS expansions to estimate the magnitude of the error that can be expected from the use of the conventional RI in an auxiliary basis. Valence MP2-R12/2*A′ correlation energies are plotted in Fig. 1 for AO basis set cc-pVTZ using uncontracted cc-pVnZ sets for the RI basis, with n = T, Q, 5. To highlight the differences between the RI basis sets, the data are presented relative to the MP2/cc-pVTZ correlation energies. Analogous data are presented in Fig. 2 for the cc-pVQZ AO basis. It can be seen that the differences between the two largest RI basis sets are at most of the order of a few mE_h. This is using the most straightforward ABS RI implementation¹¹ rather than more refined models;^10,31 and we have also used basis sets that were not specifically optimized for this role.

Fig. 1 MP2-R12 corrections using the cc-pVTZ orbital basis, and using uncontracted cc-pVnZ basis sets for the RI with n = T, Q, 5.

Fig. 2 MP2-R12 corrections using the cc-pVQZ orbital basis, and using uncontracted cc-pVnZ basis sets for the RI with n = Q, 5.

Concerning the evaluation of 3-electron integrals we observe that: (1) the errors using the straightforward ABS method in standard basis sets leads are small; (2) the errors can be further reduced by more advanced methods;^10,31 (3) the computational cost rises only linearly with RI basis size.²²

It is therefore both technically straightforward and computationally feasible to reduce the RI errors way below the target accuracies required for standard chemical applications. Having said that, even with the currently available extensions to the RI approach, it will probably not be feasible to treat heavy elements. Elements with occupied f-orbitals require RI basis sets saturated to [small script l] = 6 with the RI-DF technique³¹ and to [small script l] = 9 with conventional RI approaches.

3. Effect of the GBC, EBC and [K̂₁, f₁₂] = 0

As discussed above, several approximations beyond the RI are made in MP2-R12 theory. If the RI is not responsible for our 30 mE_h error, perhaps these approximations are to blame.

The effect of deleting exchange commutators is not negligible. The mean difference in valence correlation energies between the MP2-R12/2A′ and MP2-R12/2B methods for eight small first row systems is 5.5 mE_h and 2.4 mE_h for cc-pVTZ and cc-pVQZ, respectively.¹¹ From this we can derive two conclusions: (1) the discrepancy is certainly large enough to be of concern; (2) it is not the principal source of the 30 mE_h short-fall in correlation energy for the neon atom.

This error may instead arise from the Brillouin conditions. Fortunately it is now possible to eliminate both the GBC and EBC assumptions from MP2-R12 theory³³ so we are in a position to quantify the impact of these two approximations—details of the theory will appear separately.³³ MP2-R12 energy corrections are shown in Figs. 3 and 4 for the cc-pVTZ and cc-pVQZ basis sets. Data are presented for the case where neither Brillouin condition is assumed (labelled ‘exact’); for the case where only the GBC is assumed; and for the case where both the GBC and EBC are assumed. It can be seen that the impact of the GBC in these calculations is completely negligible. In the cc-pVTZ basis the mean absolute error arising from the GBC is less than 0.1 mE_h and the maximum error is only 0.4 mE_h. In the cc-pVQZ basis set the errors are marginally larger; this arises from the fact that the triple-zeta basis set appears to hit a ‘sweet-spot’, with the GBC errors being significantly smaller than those found with either cc-pVDZ or cc-pVQZ basis sets. This phenomenon need not distract us here, and will be discussed elsewhere.³³

Fig. 3 MP2-R12 corrections to valence correlation energies of the test molecules in mE_h, using the cc-pVTZ AO basis set. The first set of data comes from MP2-R12 calculations without the assumption of either the GBC or the EBC (exact); the second and third rely on the GBC and both the GBC and EBC, respectively. Details of the calculations are given in the text.

Fig. 4 As for Fig. 3 but using the cc-pVQZ AO basis set.

The effect of the EBC can be seen to be significantly larger than that of the GBC, but certainly not large enough to account for the disappointing performance of MP2-R12 methods. For both basis sets the mean absolute error from the combined effect of the GBC and EBC is around 1.5 mE_h. In any case the EBC can readily be avoided by choosing the projection operator Q̂^O₁₂ given in eqn. (4) and following the procedure of Klopper and Samson.¹¹

4. Form of the correlation factor

There is a considerable history relating to the development of optimal correlation factors, but a detailed review of the literature is neither practical nor necessary here. The interested reader is referred to the quantum Monte Carlo literature,³⁴ where use of a variety of Jastrow factors³⁵ is commonplace. There has been considerable recent activity in attempting to use correlation factors other than r₁₂ in MP2-R12-like theories. In 2004 two of us presented results that showed that even a single Gaussian function exp(−αr²₁₂) could perform much better than r₁₂ at reproducing the MP2 basis set limit.³⁶ At the same time Ten-no demonstrated that excellent results could be obtained by using a single Slater function exp(−ζr₁₂) for the correlation factor.^37,38

Quite a lot is now known about the evaluation of the necessary 2-electron integrals that arise with these and other correlation factors. For example, formulae have been published for integrals with the following 2-electron kernels: a simple Gaussian exp(−αr²₁₂);^5,6 the damped factor r₁₂exp(−αr²₁₂);³⁹ the damped Coulomb form r⁻¹₁₂exp(−αr²₁₂);²¹ arbitrary powers and a Gaussian r^k₁₂exp(−αr²₁₂);⁴⁰ and a simple Slater exp(−αr₁₂).³⁸

Two of us have previously published a description of the MP2-F12 method in which an arbitrary correlation factor is chosen in place of r₁₂,²¹ built as a frozen linear combination of Gaussian geminals:

(Ten-no has previously described a transcorrelated method that relies on a correlation factor of the same form.)⁴¹ Here we use nine Gaussians (with exponents in an even tempered sequence with ratio 3 and centre 8) and determine the coefficients c_γ by a least squares fit to a Slater function with unit exponent. We have tried other forms for the correlation factor, and different exponents, but the Slater exp(−r₁₂) does appear to be amazingly effective for valence correlation energies. A more extensive study of the effect of using different correlation factors is in progress in our group and elsewhere.⁴²

Our calculations are performed at the MP2-{R12,F12}/2*A′ levels of theory using cc-pVTZ and cc-pVQZ for the AO basis and uncontracted cc-pV5Z for the RI. The MP2-F12 and MP2-R12 corrections are compared with the difference between conventional MP2 and the MP2 basis set limit in Figs. 5 and 6. The basis set limit numbers were obtained by extrapolating cc-pV5Z and cc-pV6Z correlation energies using the cubic formula of Helgaker et al.⁴³

Fig. 5 MP2-R12 and MP2-F12 corrections to valence correlation energies of the test molecules in mE_h using the cc-pVTZ AO basis set. For comparison the difference between MP2/V[5,6]Z and MP2/VTZ correlation energies are also shown. Details of the calculations are given in the text.

Fig. 6 As for Fig. 5 but using the cc-pVQZ AO basis set.

It can be seen that substituting r₁₂ by exp(−r₁₂) recovers almost all of the difference between MP2-R12 and the basis set limit. For cc-pVTZ the MP2-R12 correlation errors are in error by an average of 21 mE_h; this is reduced to just 1.8 mE_h using MP2-F12. In the cc-pVQZ basis, the mean error is reduced from 7.5 mE_h to 0.6 mE_h.

5. Conclusions

MP2-R12 theory fails to deliver accurate MP2 correlation energies in AO basis sets that are currently practical for large molecules. This could arise from several sources of error in MP2-R12 methods: (1) the use of RIs to form 3-electron integrals; (2) the use of other approximations ([K̂₁, f₁₂] = 0, GBC, EBC); (3) the defects inherent in the ansatz.

It has been shown conclusively that the last factor is by far the most significant. For explicitly correlated methods that aspire to chemical accuracy in modest (triple- or quadruple-zeta) basis sets it is therefore essential that r₁₂ is replaced by a more appropriate correlation factor.

The exact form of this factor could be determined in various ways. The most obvious is to construct the correlation factor as a linear combination of suitable functions—e.g. Gaussian geminals—and minimize the Hylleraas functional with respect to the coefficients. This has to be done independently for each pair to maintain extensivity and orbital invariance:

|u_ij〉 = t^ij_ab|ab〉 + t^ij_kl,γQ̂₁₂e^{−γr2₁₂}|kl〉. (7)

However the evidence of Figs. 5 and 6 would suggest that this additional complexity would be unnecessary, since one gets so far by using a single Slater-type correlation factor. A compromise would be to have different Slater exponents for the singlet and triplet pairs, optimized over a suitable test set of molecules.

The errors that arise from the deletion of exchange commutators and from the EBC are too large to ignore in triple- and quadruple-zeta basis sets. It does not appear to be necessary to drop the GBC so that the main structure of Klopper and Samson’s theory¹¹ appears to remain extremely useful. We are therefore currently implementing the MP2-F12/2B method in the hope that it will yield higher accuracy than is currently possible.

Acknowledgements

FRM and AJM are grateful to the EPSRC for funding this research (GR/R93704/01). FRM thanks the Royal Society for funding through a University Research Fellowship scheme and is indebted to Professor W. Klopper and Professor S. Ten-no for stimulating discussions on the explicitly correlated electronic structure theory.

References

1. E. A. Hylleraas, Z. Phys., 1929, 54, 347.
2. T. Helgaker and W. Klopper, Theor. Chim. Acta, 1999, 103, 180 , for a modern perspective on ref. 1.
3. H. Hettema, Quantum Chemistry: Classic Scientific Papers, in 20th Century Chemistry, World Scientific, Singapore, 2000, vol. 8, for an English translation of ref. 1.
4. E. A. Hylleraas, Z. Phys., 1928, 48, 469.
5. S. F. Boys, Proc. R. Soc. London, Ser. A, 1960, 258, 402.
6. K. Singer, Proc. R. Soc. London, Ser. A, 1960, 258, 412.
7. S. F. Boys and N. C. Handy, Proc. R. Soc. London, Ser. A, 1969, 310, 43.
8. K. Szalewicz, B. Jeziorski, H. J. Monkhorst and J. G. Zabolitzky, Chem. Phys. Lett., 1982, 91, 169.
9. W. Kutzelnigg, Theor. Chim. Acta, 1985, 68, 445.
10. E. F. Valeev, Chem. Phys. Lett., 2004, 395, 190.
11. W. Klopper and C. C. M. Samson, J. Chem. Phys., 2002, 116, 6397.
12. W. Klopper and W. Kutzelnigg, Chem. Phys. Lett., 1987, 134, 17.
13. W. Kutzelnigg and W. Klopper, J. Chem. Phys., 1991, 94, 1985.
14. V. Ternath, W. Kutzelnigg and W. Klopper, J. Chem. Phys., 1991, 94, 2002.
15. W. Klopper and W. Kutzelnigg, J. Chem. Phys., 1991, 94, 2020.
16. J. Noga, W. Kutzelnigg and W. Klopper, Chem. Phys. Lett., 1992, 199, 497.
17. R. J. Gdantiz, Chem. Phys. Lett., 1993, 210, 253.
18. J. Noga and W. Kutzelnigg, J. Chem. Phys., 1994, 101, 7738.
19. H. Fliegl, W. Klopper and C. Hättig, J. Chem. Phys., 2005, 122, 84107.
20. P. Wind, W. Klopper and T. Helgaker, Theor. Chem. Acc., 2002, 107, 173.
21. A. J. May and F. R. Manby, J. Chem. Phys., 2004, 121, 4479.
22. W. Klopper, J. Chem. Phys., 2004, 120, 10890.
23. T. Helgaker, P. Jørgensen and J. Olsen, Molecular Electronic Structure Theory, Wiley, Chichester, 2000.
24. F. R. Manby, J. Chem. Phys., 2003, 119, 4607.
25. R. D. Amos, A. Bernhardsson, A. Berning, P. Celani, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, C. Hampel, G. Hetzer, P. J. Knowles, T. Korona, R. Lindh, A. W. Lloyd, S. J. McNicholas, F. R. Manby, W. Meyer, M. E. Mura, A. Nicklass, P. Palmieri, R. Pitzer, G. Rauhut, M. Schütz, U. Schumann, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson and H.-J. Werner, MOLPRO, a package of ab initio programs designed by H.-J. Werner and P. J. Knowles, Version 2002.9 (development version), 2005.
26. C. Hättig, Phys. Chem. Chem. Phys., 2005, 7, 59.
27. C. L. Janssen, I. B. Nielsen, M. L. Leininger, E. F. Valeev and E. T. Seidl, The Massively Parallel Quantum Chemistry Program (MPQC): Version 2.3 (alpha), Sandia National Laboratories, Livermore, CA, USA, 2005, http://www.mpqc.org/.
28. S. Ten-no, Chem. Phys. Lett., 2000, 330, 175.
29. K. Szalewicz, B. Jeziorski, H. J. Monkhorst and J. G. Zabolitsky, J. Chem. Phys., 1983, 78, 1420.
30. A. J. May and F. R. Manby, 2004, unpublished work.
31. S. Ten-no and F. R. Manby, J. Chem. Phys., 2003, 119, 5358.
32. S. Ten-no, J. Chem. Phys., 2004, 121, 117.
33. E. F. Valeev, 2005, to be submitted.
34. W. M. C. Foulkes, L. Mitas, R. J. Needs and G. Rajagopal, Rev. Mod. Phys., 2001, 73, 33.
35. R. J. Jastrow, Phys. Rev., 1955, 98, 1479.
36. F. R. Manby and A. J. May, Abstr. Pap. Am. Chem. Soc., 2004, 228, 273.
37. S. Ten-no, Abstr. Pap. Am. Chem. Soc., 2004, 228, 275.
38. S. Ten-no, Chem. Phys. Lett., 2004, 398, 56.
39. C. C. M. Samson, W. Klopper and T. Helgaker, Comput. Phys. Commun., 2002, 149, 1.
40. V. Weber and C. Daul, Comput. Phys. Commun., 2004, 158, 1.
41. S. Ten-no, Chem. Phys. Lett., 2000, 330, 169.
42. D. Tew and W. Klopper, J. Chem. Phys., 2005 , submitted.
43. T. Helgaker, W. Klopper, H. Koch and J. Noga, J. Chem. Phys., 1997, 106, 9639.

---