- Email: [email protected]
Exploration of Brueckner orbital trial wave functions in diffusion Monte Carlo calculations
Chemical Physics Letters 644 (2016) 117–120
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Exploration of Brueckner orbital trial wave functions in diffusion Monte Carlo calculations Michael J. Deible 1 , Kenneth D. Jordan ∗ Department of Chemistry, University of Pittsburgh, Pittsburgh, PA 15260, United States
a r t i c l e
i n f o
Article history: Received 29 July 2015 In final form 29 November 2015 Available online 9 December 2015
a b s t r a c t This study explores the utility of Brueckner orbitals as trial wave functions for diffusion Monte Carlo (DMC) calculations. Comparison is made with Hartree–Fock (HF) and density functional theory (DFT) orbitals allowing for an assessment of how well the three sets of orbitals describe the nodal surfaces. For the neutral test systems, PBE0 orbitals or Brueckner orbitals give DMC energies that are appreciably lower than those obtained using Hartree–Fock orbitals. For a CO2 − anion test case, a significantly lower DMC energy is obtained when using Brueckner orbitals rather than DFT orbitals as the trial function. © 2015 Elsevier B.V. All rights reserved.
1. Introduction The diffusion Monte Carlo (DMC) [1,2] method is capable of giving exact electronic energies within the fixed-node approximation [3,4]. A trial wave function, generally taken to be a Slater determinant of orbitals from density-functional theory (DFT) or Hartree–Fock (HF) calculations, is used to enforce the fixed-node approximation. This condition is required to prevent the ground state wave function from collapsing to a Bosonic state. If the trial wave function were to exactly describe the nodal surface for exchange of electrons, the DMC method, if used in all-electron calculations and run for sufficiently large number of moves and corrected for time step bias, would give the exact non-relativistic ground state energy [5]. In general, trial functions employing a Slater determinant of DFT orbitals give lower total energies than do trial functions employing a Slater determinant of Hartree–Fock orbitals, which implies that the use of DFT orbitals provides a better description of the nodal surface [6,7]. In the present study, we explore the use of trial functions comprised of a Slater determinant of Brueckner orbitals (BO) [8,9]. One might expect that such trial wave functions would give a nodal surfaces superior to a Slater determinant of DFT orbitals, as the Slater determinant of BOs is that with the maximum overlap with the exact wave function for the basis set employed [10]. There has been considerable discussion in the literature concerning the similarity and differences between HF, DFT, and
∗ Corresponding author. E-mail address: [email protected] (K.D. Jordan). 1 Present Address: Department of Civil and Mechanical Engineering, United States Military Academy, West Point, NY 10996, United States. http://dx.doi.org/10.1016/j.cplett.2015.11.053 0009-2614/© 2015 Elsevier B.V. All rights reserved.
Brueckner orbitals. Scuseria [11] has shown that the DFT equations can be derived through approximations made to the Brueckner equations. Lindgren and Salmonson [12] have argued that DFT and Brueckner orbitals are closely related. Heßelmann and Jansen have shown that Brueckner orbitals offer an improvement over Hartree–Fock orbitals when calculating first-order intermolecular interaction energies [13], whereas DFT orbitals may or may not give an improved description of the first-order Coulomb and exchange energies, with the performance depending on the functional used to generate the orbitals [14]. Jankowsi et al. [15,16] calculated the distance between orbital subspaces, and based on this measure, concluded that DFT orbitals can differ appreciably from both HF and Brueckner orbitals. An alternative criterion for testing the quality of various types of orbitals is to determine how well they describe the nodal surface for exchange of electrons which can be evaluated by assessing their performance when used as trial functions for DMC calculations. In this study, we investigate the performance of trial functions using HF, DFT, and Brueckner orbitals in DMC calculations on a series of diatomic molecules and on a bent CO2 − ion to determine whether the use of Brueckner orbitals leads to lower DMC energies than obtained using trial functions comprised of DFT orbitals. 2. Methodology The diatomic molecules studied include BeO [17], N2 [18], O2 [18], F2 [18], and CN [18], with the equilibrium geometries being taken from the respective references. To obtain atomization energies, calculations were also carried out on the atoms in their ground electronic states. For these test systems, both pseudopotential and all-electron calculations were carried out. In the all-electron calculations the orbitals were expanded in terms of the cc-pVTZ basis
118
M.J. Deible, K.D. Jordan / Chemical Physics Letters 644 (2016) 117–120
set [19]. In the calculations of the dimer atomization energies using pseudopotentials, the Trail-Needs AREP pseduopotentials [20,21] and basis sets of Xu et al. [22] were used for all atoms with the exception of Be. For N, F, and C, the spd portions of the triplezeta basis set were used and for oxygen the sp portion of the quintuple basis set was combined with the d functions from the triple-zeta basis set. For Be, the pseudopotential and valence triplezeta basis set of Burkatzki, Filippi, and Dolg (BFD) were used [23]. The calculations on CO2 and its anion made use of the Trail-Needs pseudopotentials [20,21]. For both C and O 8s8p3d contracted Gaussian basis sets formed by adding to the 6s6p portion of the quintuple-zeta basis set of Xu et al. [22] two diffuse s and two diffuse p functions with exponents three and nine times smaller than those of the most diffuse primitive functions in the 6s6p set as well as d functions with exponents of 0.5586, 0.2271, and 0.1024 and 1.2192, 0.4904, and 0.2053 for carbon and oxygen atoms, respectively. The hybrid B3LYP [24,25] functional was used for the DFT calculations on the all-electron systems, and both the B3LYP and PBE0 [26] functionals were used for the pseudopotential calculations of the diatomics. The Brueckner orbitals were obtained from coupled-cluster calculations in which the orbitals were rotated so as to eliminate single excitations to all orders [27,28]. For one system, N2 , we also considered a trial function based on PBE [29] orbitals and also employed the more flexible cc-pVQZ-g basis set for the all-electron calculations [19]. In addition to the diatomic test cases described above, we also considered CO2 and CO2 − with CO bond lengths of 1.215 A˚ and an OCO angle of 147◦ . This geometry was chosen because earlier studies have shown electron correlation effects are important for the vertical electron detachment energy (VDE) which increases from 0.18 to 0.32 eV in going from the Hartree–Fock to CCSD(T) method [30]. DFT methods, on the other hand, drastically overbind the excess electron at this geometry, with, for example, the B3LYP value of the VDE being 0.80 eV. All open-shell systems were described using spinunrestricted orbitals. The trial functions were generated using the Gaussian 09 code [31], and the diffusion Monte Carlo calculations were carried out using the CASINO code [32]. The correction scheme of Ma et al. [33] was used in the DMC calculations to account for electron-nuclear cusps in the all-electron calculations. The DMC calculations were preceded by variational Monte Carlo (VMC) calculations to optimize (via energy minimization) the parameters in the Jastrow factors [34] which include explicit electron-nuclear (e-n), electron–electron (e–e), and e-e-n terms. The all-electron DMC calculations were carried out for time steps of 0.001, 0.003, 0.005, and 0.007 a.u., and the resulting energies were extrapolated to zero time step by use of quadratic fits. The calculations were carried out with 40 000 walkers and for sufficient number of steps so as to reduce the statistical errors in the
extrapolated DMC energies of the diatomics and their atomization energies to under 0.3 and 0.5 kcal/mol, respectively. The errors in the atomization energies are defined as the differences between the experimental values, corrected for vibrational zero-point energy (ZPE), and the corresponding DMC results. For N2 , O2 , F2 , and CN, the experimental atomization energy values and zero-point energies are taken from reference [35]. For BeO, the experimental atomization energy and ZPE are taken from cccbdb.nist.gov. The parameters of the Jastrow factors for the diatomics with pseudopotntials and for the CO2 test system were optimized via variance minimization. Time steps of 0.005, 0.0075, and 0.01 a.u. were used in the DMC calculations on the diatomic species with pseudopotentials, and time steps of 0.0125, 0.003, and 0.005 a.u. were used for the CO2 test system. Linear fits were used to extrapolate to zero time step. The T-move procedure was used in the pseudopotential calculations to correct for the localization error [36].
3. Results The total energies from the all-electron DMC calculations on the diatomic species and associated atoms are summarized in Table 1. For the atoms, DMC calculations with the Hartree–Fock, B3LYP, and Brueckner orbitals give total energies that essentially agree to within statistical error. This is consistent with earlier studies [7] that found that for atoms DMC energies were largely insensitive to whether DFT or HF orbitals are used for the trial functions. The situation is quite different for the diatomics: going from HF orbitals to B3LYP orbitals leads to energy decreases ranging from 1.63 kcal/mol for BeO to 13.8 kcal/mol for CN. The DMC energies obtained using B3LYP and Brueckner orbitals agree to within one standard deviation for N2 and BeO, but for O2 , F2 , and CN, significantly lower DMC energies are obtained when using B3LYP orbitals in the trial function. It should be noted that the use of spin-unrestricted calculations introduces spin contamination in the wave function. This has a larger effect on the Hartree–Fock orbitals than it does for the B3LYP and Bruckner orbitals. The largest spin contamination is found for the CN molecule, which has S2 expectation values of 1.158, 0.758, and 0.764 for the HF, B3LYP, and Brueckner orbital wave functions, respectively, compared to the exact value of 0.750. While using RHF rather than UHF orbitals can result in lower DMC energies for open-shell systems [7,37] it would not change the qualitative result that lower DMC energies result when using B3LYP and Brueckner orbitals than when using HF orbitals (whether RHF or UHF), and that use of B3LYP orbitals generally give a lower DMC energy, thus superior nodal surface, than use of Brueckner orbitals. Our DMC calculations using Hartree–Fock
Table 1 Total energies from DMC calculationsa using Hartree–Fock, B3LYP, and Bruckner orbitals. Species
Beryllium Carbon Nitrogen Oxygen Fluorine N2 O2 F2 BeO CN a b c
Energyb (a.u.) Hartree–Fock
B3LYP
Brueckner
Brueckner (cc)c
−14.6575(1) −37.8296(2) −54.5759(2) −75.0512(2) −99.7161(3) −109.5007(6) −150.2808(5) −199.4740(5) −89.8823(4) −92.6668(5)
−14.6572(1) −37.8301(2) −54.5765(2) −75.0518(2) −99.7169(3) −109.5047(4) −150.2873(5) −199.4850(2) −89.8849(3) −92.6888(5)
−14.6575(1) −37.8296(2) −54.5757(2) −75.0518(3) −99.7161(3) −109.5042(4) −150.2856(5) −199.4812(5) −89.8845(3) −92.6876(5)
– – – – −99.7163(2) – – −199.4816(3) – –
Results extrapolated to dt = 0 as described in the text. Statistical errors (one standard deviation) are given in parentheses. Results obtained using Brueckner orbitals obtained from coupled cluster calculations including core correlation.
M.J. Deible, K.D. Jordan / Chemical Physics Letters 644 (2016) 117–120 Table 2 Errors in the atomization energies calculated using the DMC method with single determinant trial functions and different orbital choices. Atomization energy errora (kcal/mol)
Species
N2 O2 F2 BeO CN a
Table 3 Total energies from DMC pseudopotential calculationsa using Hartree–Fock, B3LYP, PBE0, and Bruckner orbitals. Species
Hartree–Fock
B3LYP
Brueckner
−9.52 (45) −8.26 (40) −11.97 (49) 3.47 (29) −17.08 (36)
−7.76 (35) −4.94 (40) −6.07 (40) 4.91 (23) −3.97 (36)
−7.07 (35) −6.00 (49) −7.45 (49) 4.47 (27) −3.90 (36)
N2 O2 F2 CN BeO a
Statistical errors (one standard deviation) are given in parentheses.
b c
Atomization energy error (kcal/mol)
and B3LYP orbitals for the atoms give slightly higher DMC energies than obtained by Per et al. [7] using the same trial wave functions. However, for the diatomics, Per et al’s DMC energies are lower than our energies by up to 2.4 kcal/mol. These differences are due to the greater flexibility of the atomic basis sets used to represent the orbitals in the study of Per et al. Indeed, for N2 , the DMC energy obtained using a Slater determinant of B3LYP orbitals as the trial function is 1.7 kcal/mol lower in energy with the cc-pVQZ-g basis set than with the cc-pVTZ basis set. Table 2 and Fig. 1 report the errors in the calculated atomization energies as previously defined. Consistent with the results discussed above for the total energies, the errors in the atomization energies tend to be smaller when calculated using either DFT or Brueckner orbitals than HF orbitals, with the errors obtained using the DFT orbitals being somewhat smaller. The exception is BeO, for which the DMC calculations using a Hartree–Fock determinant as the trial function give the smallest error in the atomization energy. As seen from Table 2, of the molecules considered, this is the only one for which the DMC calculations overestimate the atomization energy. The reason for this is well understood, namely, that the Be atom has significant multiconfigurational character, which is important for describing its nodal surface in all-electron DMC calculations [38]. Thus, using a single determinant trial function results in too high a value of the DMC energy of the Be atom and an overestimation of the atomization energy of BeO. In analyzing these results, it is important to note that the Brueckner orbitals were obtained from calculations using frozen 1s cores. As a check on whether relaxing this constraint significantly impacts the DMC energies in the case of F2 , we also carried out DMC calculations using Brueckner orbitals generated by correlating all electrons and using the cc-pCVTZ basis set [19,39] which includes functions for correlating the core. The error in the DMC value of the atomization energy of F2 obtained using as the trial function a Slater determinant of Brueckner orbitals generated correlating all electrons is essentially identical to that using Brueckner orbitals generated in the frozen-core approximation. As noted in the Introduction, a Slater determinant of Brueckner orbitals is the Slater determinant wave function that has the 10 5
HF
B3LYP
BD
0 -5 -10 -15 -20
N2
O2
F2
BeO
CN
Figure 1. Errors in the dimer atomization energies for different trial functions. HF and BD denote Hartree–Fock and Brueckner doubles, respectively.
119
Energyb (a.u.) Hartree–Fock
B3LYP
PBE0
BDc
−19.8688 (2) −31.8668 (3) −48.2722 (3) −15.4266 (2) −17.0119 (3)
−19.8715 (2) −31.8707 (3) −48.2768 (3) −15.4476 (2) −17.0142 (3)
−19.8723 (3) −31.8715 (3) −48.2781 (3) −15.4480 (2) −17.0146 (3)
−19.8721 (3) −31.8717 (3) −48.2788 (3) −15.4478 (2) −17.0149 (3)
Results extrapolated to dt = 0 as described in the text. Statistical errors (one standard deviation) are given in parentheses. BD denotes Brueckner doubles.
maximum overlap with the exact wave function for the basis set employed. To illustrate the extent of the change of the overlap when going from Hartree–Fock to Brueckner orbitals, in the case of BeO, we carried out full CI calculations using Hartree Fock and Brueckner orbitals represented in terms of the cc-pVDZ basis set [19]. The coefficient of the leading Slater determinant in the expansion goes from 0.925 to 0.942 when changing from Hartree–Fock to Brueckner orbitals. For all diatomic molecules considered, the single-determinant based DMC calculations are unable to achieve chemical accuracy (±1 kcal/mol) in the atomization energies regardless of the orbitals used [7,40,41]. The failure of the DMC calculations to accurately reproduce the atomization energies is a result of the fixed node approximation. Part of the error is due to the limitations of the basis set used to represent the orbitals. For example, in the case of N2 , the error in the atomization energy at a DMC time step of 0.001 a.u. is reduced from −7.83 to −6.28 kcal/mol in going from the cc-pVTZ to the cc-pVQZ-g basis set for representing the B3LYP orbitals in the single determinant trial function. However, most of the error in the description of the nodal surfaces is due to the use of single determinant trial functions. Indeed, it has been found that significantly improved values of the atomization energies are obtained by use of full-valence CASSCF trial functions, particularly when optimizing the CI coefficients, orbitals, and Jastrow factors simultaneously [37]. Chemical accuracy can also be achieved by using large trial functions from CI calculations employing natural orbitals [42]. We now consider the results for BeO, CN, O2 , F2 , and N2 obtained from the calculations employing pseudopotentials. The total energies obtained from DMC calculations using various types of orbitals in the single determinant trial functions are summarized in Table 3. From this table it is seen that for the N2 , O2 , and F2 molecules the DMC calculations using trial functions employing Brueckner orbitals give energies lower than those employing B3LYP orbitals, but statistically similar energies to those obtained using PBE0 orbitals. For CN and BeO the total DMC energies calculated using the three sets of orbitals agree to within statistical uncertainty. In contrast to the all-electron calculations, when pseudopotentials are employed to eliminate the 1s core orbitals, lower DMC energies are obtained from trial functions using Brueckner orbitals rather than using B3LYP orbitals. Table 4 summarizes the results of the DMC calculations on CO2 and CO2 − . It is seen from the table that for neutral CO2 , the DMC calculations using trial functions of DFT orbitals (B3LYP, PBE0, or PBE), give comparable or slightly lower energies than obtained using a trial function represented in terms of Brueckner orbitals. In contrast, for the anion, the DMC calculations using a trial function comprised of Brueckner orbitals give an appreciably lower energy than the calculations using trial functions represented in terms of DFT orbitals, with the energy difference being much more pronounced for the case of PBE than for the B3LYP or PBE0 orbitals.
120
M.J. Deible, K.D. Jordan / Chemical Physics Letters 644 (2016) 117–120
Table 4 DMC energies (a.u.) of CO2 and CO2 − at an OCO angle of 147◦ .a , b Orbitals
Anion
Neutral
Energy difference (eV)
B3LYP PBE0 PBE BDc HFc
−37.6811 (2) −37.6809 (3) −37.6781 (3) −37.6825 (2) −37.6745 (2)
−37.6691 (2) −37.6696 (2) −37.6691 (2) −37.6688 (2) −37.6608 (2)
−0.327 (8) −0.31 (1) −0.24 (1) −0.378 (8) −0.373 (8)
a b c
˚ The CO bond lengths are set to 1.215 A. Statistical errors (one standard deviation) are given in parentheses. HF and BD denote Hartree–Fock and Brueckner doubles, respectively.
4. Conclusions In summary, we have demonstrated that DMC calculations using trial functions with a single Slater determinant of Brueckner orbitals gives atomization energies of a set of test diatomics much closer to experiment than DMC calculations employing Slater determinants of HF orbitals. When pseudopotentials are employed the DMC energies obtained using Brueckner orbitals are essentially identical to those employing PBE0 orbitals and give lower energies than those using B3LYP orbitals. It is well known that DFT methods do not properly describe certain types of anions, and for these there can be a significant advantage to using Brueckner rather than DFT orbitals in the trial function. For example, for a bent CO2 − test system, the DMC calculations using as a trial function a Slater determinant of Brueckner orbitals give a lower energy than when employing DFT orbitals thereby establishing that a Slater determinant of Brueckner orbitals better describes the nodal surface in this case.
[4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
[18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]
Acknowledgement This research was supported by grant CHE136234 from the National Science Foundation. Computational resources were provided by the Center of Simulation and Modeling at the University of Pittsburgh, with hardware purchased under the NSF MRI award 1229064. References [1] W.M. Foulkes, L. Mitas, R.J. Needs, G. Rajagopal, Rev. Mod. Phys. 73 (2001) 33–83. [2] B.M. Austin, D.Y. Zubarev, W.A. Lester, Chem. Rev. 112 (2012) 263–288. [3] J.B. Anderson, J. Chem. Phys. 65 (1976) 4121–4127.
[32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]
J.B. Anderson, J. Chem. Phys. 63 (1975) 1499–1503. D. Ceperley, J. Stat. Phys. 63 (1991) 1237–1267. J. Ma, D. Alfè, A. Michaelides, E. Wang, J. Chem. Phys. 130 (2009) 154303. M.C. Per, K.A. Walker, S.P. Russo, J. Chem. Theory Comput. 8 (2012) 2255–2259. K.A. Brueckner, Phys. Rev. 96 (1954) 508–516. R.K. Nesbet, Phys. Rev. 109 (1958) 1632–1638. S. Larsson, V.H. Smith Jr., Phys. Rev. 178 (1969) 137–152. G. Scuseria, Int. J. Quant. Chem. 55 (1995) 165–171. I. Lindgren, S. Salomonson, Int. J. Quantum Chem. 90 (2002) 294–308. A. Heßelmann, G. Jansen, J. Chem. Phys. 112 (2000) 6949–6952. A. Heßelmann, G. Jansen, Chem. Phys. Lett. 357 (2002) 464–470. K. Jankowski, K. Nowakowski, J. Wasilewski, Chem. Phys. Lett. 389 (2004) 393–399. K. Jankowski, K. Nowakowski, J. Wasilewski, Struct. Chem. 15 (2004) 437–445. K.P. Huber, G. Herzberg, G. Van Nostrand, Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules, Reinhold Company, New York, NY, 1979, pp. 732. D. Feller, K.A. Peterson, D.A. Dixon, J. Chem. Phys. 129 (2008) 204105. T.H. Dunning, J. Chem. Phys. 90 (1989) 1007–1023. J.R. Trail, R.J. Needs, J. Chem. Phys. 122 (2005) 174109. J.R. Trail, R.J. Needs, J. Chem. Phys. 122 (2005) 014112. J. Xu, M.J. Deible, K.A. Peterson, K.D. Jordan, J. Chem. Theory Comput. 9 (2013) 2170–2178. M. Burkatzki, C. Filippi, M. Dolg, J. Chem. Phys. 126 (2007) 234105. A.D. Becke, J. Chem. Phys. 98 (1993) 5648–5652. C. Lee, W. Yang, R. Parr, Phys. Rev. B 37 (1988) 785–789. C. Adamo, V. Barone, J. Chem. Phys. 110 (1999) 6158–6170. C.E. Dykstra, Chem. Phys. Lett. 45 (1977) 466–469. N. Handy, J. Pople, M. Head-Gordon, K. Raghavachari, G.W. Trucks, Chem. Phys. 164 (1989) 185–192. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. V. Voora, T. Sommerfeld, K.D. Jordan, unpublished results. M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J.A. Montgomery Jr., J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, J.J. Dannenberg, S. Dapprich, A.D. Daniels, O. Farkas, J.B. Foresman, J.V. Ortiz, J. Cioslowski, D.J. Fox, Gaussian 09, Revision B.01, Gaussian, Inc., Wallingford, CT, 2010. R.J. Needs, M.D. Towler, N.D. Drummond, P. López Ríos, J. Phys. Condens. Matter 22 (2010) 023201. A. Ma, M.D. Towler, N.D. Drummond, R.J. Needs, J. Chem. Phys. 122 (2005) 224322. N.D. Drummond, M.D. Towler, R.J. Needs, Phys. Rev. B 70 (2004) 235119. D. Feller, K.A. Peterson, J. Chem. Phys. 110 (1999) 8384–8396. M. Casula, Phys. Rev. B 74 (2006) 1–4. F.R. Petruzielo, J. Toulouse, C.J. Umrigar, J. Chem. Phys. 136 (2012) 124116. D. Bressanini, Phys. Rev. B 86 (2012) 115120. D.E. Woon, T.H. Dunning, J. Chem. Phys. 103 (1995) 4572–4585. N. Nemec, M.D. Towler, R.J. Needs, J. Chem. Phys. 132 (2010) 034111. J.C. Grossman, J. Chem. Phys. 117 (2002) 1434–1440. M.A. Morales, J. McMinis, B.K. Clark, J. Kim, G.E. Scuseria, J. Chem. Theory Comput. 8 (2012) 2181–2188.
Contents lists available at ScienceDirect
Chemical Physics Letters journal homepage: www.elsevier.com/locate/cplett
Exploration of Brueckner orbital trial wave functions in diffusion Monte Carlo calculations Michael J. Deible 1 , Kenneth D. Jordan ∗ Department of Chemistry, University of Pittsburgh, Pittsburgh, PA 15260, United States
a r t i c l e
i n f o
Article history: Received 29 July 2015 In final form 29 November 2015 Available online 9 December 2015
a b s t r a c t This study explores the utility of Brueckner orbitals as trial wave functions for diffusion Monte Carlo (DMC) calculations. Comparison is made with Hartree–Fock (HF) and density functional theory (DFT) orbitals allowing for an assessment of how well the three sets of orbitals describe the nodal surfaces. For the neutral test systems, PBE0 orbitals or Brueckner orbitals give DMC energies that are appreciably lower than those obtained using Hartree–Fock orbitals. For a CO2 − anion test case, a significantly lower DMC energy is obtained when using Brueckner orbitals rather than DFT orbitals as the trial function. © 2015 Elsevier B.V. All rights reserved.
1. Introduction The diffusion Monte Carlo (DMC) [1,2] method is capable of giving exact electronic energies within the fixed-node approximation [3,4]. A trial wave function, generally taken to be a Slater determinant of orbitals from density-functional theory (DFT) or Hartree–Fock (HF) calculations, is used to enforce the fixed-node approximation. This condition is required to prevent the ground state wave function from collapsing to a Bosonic state. If the trial wave function were to exactly describe the nodal surface for exchange of electrons, the DMC method, if used in all-electron calculations and run for sufficiently large number of moves and corrected for time step bias, would give the exact non-relativistic ground state energy [5]. In general, trial functions employing a Slater determinant of DFT orbitals give lower total energies than do trial functions employing a Slater determinant of Hartree–Fock orbitals, which implies that the use of DFT orbitals provides a better description of the nodal surface [6,7]. In the present study, we explore the use of trial functions comprised of a Slater determinant of Brueckner orbitals (BO) [8,9]. One might expect that such trial wave functions would give a nodal surfaces superior to a Slater determinant of DFT orbitals, as the Slater determinant of BOs is that with the maximum overlap with the exact wave function for the basis set employed [10]. There has been considerable discussion in the literature concerning the similarity and differences between HF, DFT, and
∗ Corresponding author. E-mail address: [email protected] (K.D. Jordan). 1 Present Address: Department of Civil and Mechanical Engineering, United States Military Academy, West Point, NY 10996, United States. http://dx.doi.org/10.1016/j.cplett.2015.11.053 0009-2614/© 2015 Elsevier B.V. All rights reserved.
Brueckner orbitals. Scuseria [11] has shown that the DFT equations can be derived through approximations made to the Brueckner equations. Lindgren and Salmonson [12] have argued that DFT and Brueckner orbitals are closely related. Heßelmann and Jansen have shown that Brueckner orbitals offer an improvement over Hartree–Fock orbitals when calculating first-order intermolecular interaction energies [13], whereas DFT orbitals may or may not give an improved description of the first-order Coulomb and exchange energies, with the performance depending on the functional used to generate the orbitals [14]. Jankowsi et al. [15,16] calculated the distance between orbital subspaces, and based on this measure, concluded that DFT orbitals can differ appreciably from both HF and Brueckner orbitals. An alternative criterion for testing the quality of various types of orbitals is to determine how well they describe the nodal surface for exchange of electrons which can be evaluated by assessing their performance when used as trial functions for DMC calculations. In this study, we investigate the performance of trial functions using HF, DFT, and Brueckner orbitals in DMC calculations on a series of diatomic molecules and on a bent CO2 − ion to determine whether the use of Brueckner orbitals leads to lower DMC energies than obtained using trial functions comprised of DFT orbitals. 2. Methodology The diatomic molecules studied include BeO [17], N2 [18], O2 [18], F2 [18], and CN [18], with the equilibrium geometries being taken from the respective references. To obtain atomization energies, calculations were also carried out on the atoms in their ground electronic states. For these test systems, both pseudopotential and all-electron calculations were carried out. In the all-electron calculations the orbitals were expanded in terms of the cc-pVTZ basis
118
M.J. Deible, K.D. Jordan / Chemical Physics Letters 644 (2016) 117–120
set [19]. In the calculations of the dimer atomization energies using pseudopotentials, the Trail-Needs AREP pseduopotentials [20,21] and basis sets of Xu et al. [22] were used for all atoms with the exception of Be. For N, F, and C, the spd portions of the triplezeta basis set were used and for oxygen the sp portion of the quintuple basis set was combined with the d functions from the triple-zeta basis set. For Be, the pseudopotential and valence triplezeta basis set of Burkatzki, Filippi, and Dolg (BFD) were used [23]. The calculations on CO2 and its anion made use of the Trail-Needs pseudopotentials [20,21]. For both C and O 8s8p3d contracted Gaussian basis sets formed by adding to the 6s6p portion of the quintuple-zeta basis set of Xu et al. [22] two diffuse s and two diffuse p functions with exponents three and nine times smaller than those of the most diffuse primitive functions in the 6s6p set as well as d functions with exponents of 0.5586, 0.2271, and 0.1024 and 1.2192, 0.4904, and 0.2053 for carbon and oxygen atoms, respectively. The hybrid B3LYP [24,25] functional was used for the DFT calculations on the all-electron systems, and both the B3LYP and PBE0 [26] functionals were used for the pseudopotential calculations of the diatomics. The Brueckner orbitals were obtained from coupled-cluster calculations in which the orbitals were rotated so as to eliminate single excitations to all orders [27,28]. For one system, N2 , we also considered a trial function based on PBE [29] orbitals and also employed the more flexible cc-pVQZ-g basis set for the all-electron calculations [19]. In addition to the diatomic test cases described above, we also considered CO2 and CO2 − with CO bond lengths of 1.215 A˚ and an OCO angle of 147◦ . This geometry was chosen because earlier studies have shown electron correlation effects are important for the vertical electron detachment energy (VDE) which increases from 0.18 to 0.32 eV in going from the Hartree–Fock to CCSD(T) method [30]. DFT methods, on the other hand, drastically overbind the excess electron at this geometry, with, for example, the B3LYP value of the VDE being 0.80 eV. All open-shell systems were described using spinunrestricted orbitals. The trial functions were generated using the Gaussian 09 code [31], and the diffusion Monte Carlo calculations were carried out using the CASINO code [32]. The correction scheme of Ma et al. [33] was used in the DMC calculations to account for electron-nuclear cusps in the all-electron calculations. The DMC calculations were preceded by variational Monte Carlo (VMC) calculations to optimize (via energy minimization) the parameters in the Jastrow factors [34] which include explicit electron-nuclear (e-n), electron–electron (e–e), and e-e-n terms. The all-electron DMC calculations were carried out for time steps of 0.001, 0.003, 0.005, and 0.007 a.u., and the resulting energies were extrapolated to zero time step by use of quadratic fits. The calculations were carried out with 40 000 walkers and for sufficient number of steps so as to reduce the statistical errors in the
extrapolated DMC energies of the diatomics and their atomization energies to under 0.3 and 0.5 kcal/mol, respectively. The errors in the atomization energies are defined as the differences between the experimental values, corrected for vibrational zero-point energy (ZPE), and the corresponding DMC results. For N2 , O2 , F2 , and CN, the experimental atomization energy values and zero-point energies are taken from reference [35]. For BeO, the experimental atomization energy and ZPE are taken from cccbdb.nist.gov. The parameters of the Jastrow factors for the diatomics with pseudopotntials and for the CO2 test system were optimized via variance minimization. Time steps of 0.005, 0.0075, and 0.01 a.u. were used in the DMC calculations on the diatomic species with pseudopotentials, and time steps of 0.0125, 0.003, and 0.005 a.u. were used for the CO2 test system. Linear fits were used to extrapolate to zero time step. The T-move procedure was used in the pseudopotential calculations to correct for the localization error [36].
3. Results The total energies from the all-electron DMC calculations on the diatomic species and associated atoms are summarized in Table 1. For the atoms, DMC calculations with the Hartree–Fock, B3LYP, and Brueckner orbitals give total energies that essentially agree to within statistical error. This is consistent with earlier studies [7] that found that for atoms DMC energies were largely insensitive to whether DFT or HF orbitals are used for the trial functions. The situation is quite different for the diatomics: going from HF orbitals to B3LYP orbitals leads to energy decreases ranging from 1.63 kcal/mol for BeO to 13.8 kcal/mol for CN. The DMC energies obtained using B3LYP and Brueckner orbitals agree to within one standard deviation for N2 and BeO, but for O2 , F2 , and CN, significantly lower DMC energies are obtained when using B3LYP orbitals in the trial function. It should be noted that the use of spin-unrestricted calculations introduces spin contamination in the wave function. This has a larger effect on the Hartree–Fock orbitals than it does for the B3LYP and Bruckner orbitals. The largest spin contamination is found for the CN molecule, which has S2 expectation values of 1.158, 0.758, and 0.764 for the HF, B3LYP, and Brueckner orbital wave functions, respectively, compared to the exact value of 0.750. While using RHF rather than UHF orbitals can result in lower DMC energies for open-shell systems [7,37] it would not change the qualitative result that lower DMC energies result when using B3LYP and Brueckner orbitals than when using HF orbitals (whether RHF or UHF), and that use of B3LYP orbitals generally give a lower DMC energy, thus superior nodal surface, than use of Brueckner orbitals. Our DMC calculations using Hartree–Fock
Table 1 Total energies from DMC calculationsa using Hartree–Fock, B3LYP, and Bruckner orbitals. Species
Beryllium Carbon Nitrogen Oxygen Fluorine N2 O2 F2 BeO CN a b c
Energyb (a.u.) Hartree–Fock
B3LYP
Brueckner
Brueckner (cc)c
−14.6575(1) −37.8296(2) −54.5759(2) −75.0512(2) −99.7161(3) −109.5007(6) −150.2808(5) −199.4740(5) −89.8823(4) −92.6668(5)
−14.6572(1) −37.8301(2) −54.5765(2) −75.0518(2) −99.7169(3) −109.5047(4) −150.2873(5) −199.4850(2) −89.8849(3) −92.6888(5)
−14.6575(1) −37.8296(2) −54.5757(2) −75.0518(3) −99.7161(3) −109.5042(4) −150.2856(5) −199.4812(5) −89.8845(3) −92.6876(5)
– – – – −99.7163(2) – – −199.4816(3) – –
Results extrapolated to dt = 0 as described in the text. Statistical errors (one standard deviation) are given in parentheses. Results obtained using Brueckner orbitals obtained from coupled cluster calculations including core correlation.
M.J. Deible, K.D. Jordan / Chemical Physics Letters 644 (2016) 117–120 Table 2 Errors in the atomization energies calculated using the DMC method with single determinant trial functions and different orbital choices. Atomization energy errora (kcal/mol)
Species
N2 O2 F2 BeO CN a
Table 3 Total energies from DMC pseudopotential calculationsa using Hartree–Fock, B3LYP, PBE0, and Bruckner orbitals. Species
Hartree–Fock
B3LYP
Brueckner
−9.52 (45) −8.26 (40) −11.97 (49) 3.47 (29) −17.08 (36)
−7.76 (35) −4.94 (40) −6.07 (40) 4.91 (23) −3.97 (36)
−7.07 (35) −6.00 (49) −7.45 (49) 4.47 (27) −3.90 (36)
N2 O2 F2 CN BeO a
Statistical errors (one standard deviation) are given in parentheses.
b c
Atomization energy error (kcal/mol)
and B3LYP orbitals for the atoms give slightly higher DMC energies than obtained by Per et al. [7] using the same trial wave functions. However, for the diatomics, Per et al’s DMC energies are lower than our energies by up to 2.4 kcal/mol. These differences are due to the greater flexibility of the atomic basis sets used to represent the orbitals in the study of Per et al. Indeed, for N2 , the DMC energy obtained using a Slater determinant of B3LYP orbitals as the trial function is 1.7 kcal/mol lower in energy with the cc-pVQZ-g basis set than with the cc-pVTZ basis set. Table 2 and Fig. 1 report the errors in the calculated atomization energies as previously defined. Consistent with the results discussed above for the total energies, the errors in the atomization energies tend to be smaller when calculated using either DFT or Brueckner orbitals than HF orbitals, with the errors obtained using the DFT orbitals being somewhat smaller. The exception is BeO, for which the DMC calculations using a Hartree–Fock determinant as the trial function give the smallest error in the atomization energy. As seen from Table 2, of the molecules considered, this is the only one for which the DMC calculations overestimate the atomization energy. The reason for this is well understood, namely, that the Be atom has significant multiconfigurational character, which is important for describing its nodal surface in all-electron DMC calculations [38]. Thus, using a single determinant trial function results in too high a value of the DMC energy of the Be atom and an overestimation of the atomization energy of BeO. In analyzing these results, it is important to note that the Brueckner orbitals were obtained from calculations using frozen 1s cores. As a check on whether relaxing this constraint significantly impacts the DMC energies in the case of F2 , we also carried out DMC calculations using Brueckner orbitals generated by correlating all electrons and using the cc-pCVTZ basis set [19,39] which includes functions for correlating the core. The error in the DMC value of the atomization energy of F2 obtained using as the trial function a Slater determinant of Brueckner orbitals generated correlating all electrons is essentially identical to that using Brueckner orbitals generated in the frozen-core approximation. As noted in the Introduction, a Slater determinant of Brueckner orbitals is the Slater determinant wave function that has the 10 5
HF
B3LYP
BD
0 -5 -10 -15 -20
N2
O2
F2
BeO
CN
Figure 1. Errors in the dimer atomization energies for different trial functions. HF and BD denote Hartree–Fock and Brueckner doubles, respectively.
119
Energyb (a.u.) Hartree–Fock
B3LYP
PBE0
BDc
−19.8688 (2) −31.8668 (3) −48.2722 (3) −15.4266 (2) −17.0119 (3)
−19.8715 (2) −31.8707 (3) −48.2768 (3) −15.4476 (2) −17.0142 (3)
−19.8723 (3) −31.8715 (3) −48.2781 (3) −15.4480 (2) −17.0146 (3)
−19.8721 (3) −31.8717 (3) −48.2788 (3) −15.4478 (2) −17.0149 (3)
Results extrapolated to dt = 0 as described in the text. Statistical errors (one standard deviation) are given in parentheses. BD denotes Brueckner doubles.
maximum overlap with the exact wave function for the basis set employed. To illustrate the extent of the change of the overlap when going from Hartree–Fock to Brueckner orbitals, in the case of BeO, we carried out full CI calculations using Hartree Fock and Brueckner orbitals represented in terms of the cc-pVDZ basis set [19]. The coefficient of the leading Slater determinant in the expansion goes from 0.925 to 0.942 when changing from Hartree–Fock to Brueckner orbitals. For all diatomic molecules considered, the single-determinant based DMC calculations are unable to achieve chemical accuracy (±1 kcal/mol) in the atomization energies regardless of the orbitals used [7,40,41]. The failure of the DMC calculations to accurately reproduce the atomization energies is a result of the fixed node approximation. Part of the error is due to the limitations of the basis set used to represent the orbitals. For example, in the case of N2 , the error in the atomization energy at a DMC time step of 0.001 a.u. is reduced from −7.83 to −6.28 kcal/mol in going from the cc-pVTZ to the cc-pVQZ-g basis set for representing the B3LYP orbitals in the single determinant trial function. However, most of the error in the description of the nodal surfaces is due to the use of single determinant trial functions. Indeed, it has been found that significantly improved values of the atomization energies are obtained by use of full-valence CASSCF trial functions, particularly when optimizing the CI coefficients, orbitals, and Jastrow factors simultaneously [37]. Chemical accuracy can also be achieved by using large trial functions from CI calculations employing natural orbitals [42]. We now consider the results for BeO, CN, O2 , F2 , and N2 obtained from the calculations employing pseudopotentials. The total energies obtained from DMC calculations using various types of orbitals in the single determinant trial functions are summarized in Table 3. From this table it is seen that for the N2 , O2 , and F2 molecules the DMC calculations using trial functions employing Brueckner orbitals give energies lower than those employing B3LYP orbitals, but statistically similar energies to those obtained using PBE0 orbitals. For CN and BeO the total DMC energies calculated using the three sets of orbitals agree to within statistical uncertainty. In contrast to the all-electron calculations, when pseudopotentials are employed to eliminate the 1s core orbitals, lower DMC energies are obtained from trial functions using Brueckner orbitals rather than using B3LYP orbitals. Table 4 summarizes the results of the DMC calculations on CO2 and CO2 − . It is seen from the table that for neutral CO2 , the DMC calculations using trial functions of DFT orbitals (B3LYP, PBE0, or PBE), give comparable or slightly lower energies than obtained using a trial function represented in terms of Brueckner orbitals. In contrast, for the anion, the DMC calculations using a trial function comprised of Brueckner orbitals give an appreciably lower energy than the calculations using trial functions represented in terms of DFT orbitals, with the energy difference being much more pronounced for the case of PBE than for the B3LYP or PBE0 orbitals.
120
M.J. Deible, K.D. Jordan / Chemical Physics Letters 644 (2016) 117–120
Table 4 DMC energies (a.u.) of CO2 and CO2 − at an OCO angle of 147◦ .a , b Orbitals
Anion
Neutral
Energy difference (eV)
B3LYP PBE0 PBE BDc HFc
−37.6811 (2) −37.6809 (3) −37.6781 (3) −37.6825 (2) −37.6745 (2)
−37.6691 (2) −37.6696 (2) −37.6691 (2) −37.6688 (2) −37.6608 (2)
−0.327 (8) −0.31 (1) −0.24 (1) −0.378 (8) −0.373 (8)
a b c
˚ The CO bond lengths are set to 1.215 A. Statistical errors (one standard deviation) are given in parentheses. HF and BD denote Hartree–Fock and Brueckner doubles, respectively.
4. Conclusions In summary, we have demonstrated that DMC calculations using trial functions with a single Slater determinant of Brueckner orbitals gives atomization energies of a set of test diatomics much closer to experiment than DMC calculations employing Slater determinants of HF orbitals. When pseudopotentials are employed the DMC energies obtained using Brueckner orbitals are essentially identical to those employing PBE0 orbitals and give lower energies than those using B3LYP orbitals. It is well known that DFT methods do not properly describe certain types of anions, and for these there can be a significant advantage to using Brueckner rather than DFT orbitals in the trial function. For example, for a bent CO2 − test system, the DMC calculations using as a trial function a Slater determinant of Brueckner orbitals give a lower energy than when employing DFT orbitals thereby establishing that a Slater determinant of Brueckner orbitals better describes the nodal surface in this case.
[4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
[18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]
Acknowledgement This research was supported by grant CHE136234 from the National Science Foundation. Computational resources were provided by the Center of Simulation and Modeling at the University of Pittsburgh, with hardware purchased under the NSF MRI award 1229064. References [1] W.M. Foulkes, L. Mitas, R.J. Needs, G. Rajagopal, Rev. Mod. Phys. 73 (2001) 33–83. [2] B.M. Austin, D.Y. Zubarev, W.A. Lester, Chem. Rev. 112 (2012) 263–288. [3] J.B. Anderson, J. Chem. Phys. 65 (1976) 4121–4127.
[32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]
J.B. Anderson, J. Chem. Phys. 63 (1975) 1499–1503. D. Ceperley, J. Stat. Phys. 63 (1991) 1237–1267. J. Ma, D. Alfè, A. Michaelides, E. Wang, J. Chem. Phys. 130 (2009) 154303. M.C. Per, K.A. Walker, S.P. Russo, J. Chem. Theory Comput. 8 (2012) 2255–2259. K.A. Brueckner, Phys. Rev. 96 (1954) 508–516. R.K. Nesbet, Phys. Rev. 109 (1958) 1632–1638. S. Larsson, V.H. Smith Jr., Phys. Rev. 178 (1969) 137–152. G. Scuseria, Int. J. Quant. Chem. 55 (1995) 165–171. I. Lindgren, S. Salomonson, Int. J. Quantum Chem. 90 (2002) 294–308. A. Heßelmann, G. Jansen, J. Chem. Phys. 112 (2000) 6949–6952. A. Heßelmann, G. Jansen, Chem. Phys. Lett. 357 (2002) 464–470. K. Jankowski, K. Nowakowski, J. Wasilewski, Chem. Phys. Lett. 389 (2004) 393–399. K. Jankowski, K. Nowakowski, J. Wasilewski, Struct. Chem. 15 (2004) 437–445. K.P. Huber, G. Herzberg, G. Van Nostrand, Molecular Spectra and Molecular Structure IV. Constants of Diatomic Molecules, Reinhold Company, New York, NY, 1979, pp. 732. D. Feller, K.A. Peterson, D.A. Dixon, J. Chem. Phys. 129 (2008) 204105. T.H. Dunning, J. Chem. Phys. 90 (1989) 1007–1023. J.R. Trail, R.J. Needs, J. Chem. Phys. 122 (2005) 174109. J.R. Trail, R.J. Needs, J. Chem. Phys. 122 (2005) 014112. J. Xu, M.J. Deible, K.A. Peterson, K.D. Jordan, J. Chem. Theory Comput. 9 (2013) 2170–2178. M. Burkatzki, C. Filippi, M. Dolg, J. Chem. Phys. 126 (2007) 234105. A.D. Becke, J. Chem. Phys. 98 (1993) 5648–5652. C. Lee, W. Yang, R. Parr, Phys. Rev. B 37 (1988) 785–789. C. Adamo, V. Barone, J. Chem. Phys. 110 (1999) 6158–6170. C.E. Dykstra, Chem. Phys. Lett. 45 (1977) 466–469. N. Handy, J. Pople, M. Head-Gordon, K. Raghavachari, G.W. Trucks, Chem. Phys. 164 (1989) 185–192. J.P. Perdew, K. Burke, M. Ernzerhof, Phys. Rev. Lett. 77 (1996) 3865–3868. V. Voora, T. Sommerfeld, K.D. Jordan, unpublished results. M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J.A. Montgomery Jr., J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, J.J. Dannenberg, S. Dapprich, A.D. Daniels, O. Farkas, J.B. Foresman, J.V. Ortiz, J. Cioslowski, D.J. Fox, Gaussian 09, Revision B.01, Gaussian, Inc., Wallingford, CT, 2010. R.J. Needs, M.D. Towler, N.D. Drummond, P. López Ríos, J. Phys. Condens. Matter 22 (2010) 023201. A. Ma, M.D. Towler, N.D. Drummond, R.J. Needs, J. Chem. Phys. 122 (2005) 224322. N.D. Drummond, M.D. Towler, R.J. Needs, Phys. Rev. B 70 (2004) 235119. D. Feller, K.A. Peterson, J. Chem. Phys. 110 (1999) 8384–8396. M. Casula, Phys. Rev. B 74 (2006) 1–4. F.R. Petruzielo, J. Toulouse, C.J. Umrigar, J. Chem. Phys. 136 (2012) 124116. D. Bressanini, Phys. Rev. B 86 (2012) 115120. D.E. Woon, T.H. Dunning, J. Chem. Phys. 103 (1995) 4572–4585. N. Nemec, M.D. Towler, R.J. Needs, J. Chem. Phys. 132 (2010) 034111. J.C. Grossman, J. Chem. Phys. 117 (2002) 1434–1440. M.A. Morales, J. McMinis, B.K. Clark, J. Kim, G.E. Scuseria, J. Chem. Theory Comput. 8 (2012) 2181–2188.