Computer-predicted ionization energy of carbon within 1 cm −1 of the best experiment
aa r X i v : . [ phy s i c s . a t o m - ph ] J un Computer-predicted ionization energy of carbon within cm − of the best experiment Nike Dattani ∗ Harvard-Smithsonian Center for Astrophysics, Atomic and Molecular Physics Division, 02138, Cambridge, MA, USA, andMcMaster University, 606-8103, Hamilton, ON, Canada,
Giovanni LiManni † Max Planck Institute for Solid State Systems, Department of Electronic Structure Theory, Suttgart, Germany.
David Feller ‡ Washington State University, Pullman, Washington 99164-4630, USA, andUniversity of Alabama, Tuscaloosa, Alabama 35487-0336, USA,
Jacek Koput § Adam Mickiewicz University, 61–614 Poznan, Poland. (Dated: 25th June 2020)We show that we can predict the first ionization energy of the carbon atom to within 0.872 cm − of the experimental value. This is an improvement of more than a factor of 6.5 over the precedingbest prediction in [Phys. Rev. A , 022503], and opens the door to achieving sub-cm − accuracyfor ab initio predictions in larger elements of the periodic table. In the last seven years, ionization energies (IEs) havebeen calculated with unprecedented precision for the Liatom [1, 2], Be atom [3] and B atom [4]. Tight varia-tional bounds for non-relativistic ground state energiesassuming a clamped, point-sized nucleus have reached 49digits in units of Hartree for He [5], 16 digits for Li [1],12 digits for Be [6], and 11 digits for B [4] (see Table I).Calculated IEs have been made in agreement with exper-iment to within 10 − cm − for Li, 10 − cm − for Beand 1 cm − for B (see Table II).For the C atom, before this present study, no high-precision calculation had been reported to predict an IEto within ∼ − agreement with experiment, but twonew experimental papers have been published on this IEvery recently [7, 8] as well as a new measurement of theelectron affinity with more than an order of magnitudebetter precision than the best previous experiment [9].Table I shows that the method used for the smaller atomsup to boron has not had success for carbon. With twiceas many variationally optimizable parameters, one fewerdigit was obtained for the B atom than for the Be atom,which also suggests that it would be very difficult to vari-ationally optimize a fully explicitly correlated wavefunc-tion ansatz for atoms and molecules coming from the restof the periodic table.The best known variational bound for the non-relativistic, clamped, point-sized nucleus (NR-CPN)ground state energy for C was calculated in 2015 usingfixed-node diffusion Monte Carlo (FN-DMC) with thenodes of the electronic wavefunction fixed at the loca- ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] tions of a CISD/cc-pV5Z wavefunction, and the statisti-cal uncertainty based on the stochastic fluctuations was ± µE Hartree [10]. However, the IE for C predicted byFN-DMC was in discrepancy with experiment by morethan 40 cm − . In this paper, the approach we use to cal-culate the NR-CPN energy of the ground state of C isFCIQMC (full configuration interaction quantum MonteCarlo) with basis sets as large as aug-cc-pCV8Z. Table Ishows that our NR-CPN energy is at least 76 µE Hartree higher than the variational upper bound obtained fromFN-DMC; but since in our approach, imperfections inthe description of the wavefunction for the neutral atomare almost the same as in the cation, the individual er-rors almost completely cancel when calculating the en-ergy difference. Therefore, with our approach we achieveagreement with experiment that is comparable within anorder of magnitude to what has been seen with the ex-plicitly correlated approach for atoms as big as (but notexceeding) boron.After adding relativistic and quantum electrodynamics(QED) corrections, and corrections to the clamped nu-cleus approximation, we obtained an IE for the groundstate of C which is in only 0.872 cm − disagreement withthe best known experimental estimate. While this isnot as impressive as the method of variationally optimiz-ing parameters in an explicitly correlated wavefunctionansatz has proven to be for Li and Be, the disagreementwith experiment has the same order of magnitude as thelatter approach for B (see Table II). We finally note thatthe approach used in this paper, of calculating FCIQMCon a basis set of non-explicitly correlated orbitals has suc-cessfully treated systems with far more electrons (tran-sition metal atoms [11], diatomics [12], multi-referencepolyatomics such as ozone [13], larger molecules such asbutadiene [14], and even solid state systems [15]), so itis conceivable that the approach used in this paper mayin the near future be able to determine (with fair ac-curacy) the IEs which at present remain experimentallyelusive or poorly known according to NIST’s atomic spec-tra database. These include arsenic (whose experimentalIE has an uncertainty of ± − ), Pm, Pa, Fm, Md,No, Sg, Bh and Hs (whose IEs are only known basedon extrapolations of other experimental data and haveuncertainties between ±
140 cm − and ± − ), Rfand Db (whose IEs are only known from theoretical cal-culations), and Mt, Ds, Rg, Cn, Nh, Fl, Mc, Lv, Ts,and Og (for which no IE is given in NIST’s most recentdatabases). I. METHODOLOGY
We begin with our main result in Table III, whichshows that our computer-predicted ionization energycomes mainly from the NR-CPN Hamiltonian. This en-ergy was calculated in four stages which we describe inthe sub-sections below: (A) We developed larger core-valence (CV) basis sets than previously available for car-bon, (B) we calculated the 1- and 2-electron integrals inthese basis sets, (C) we solved the NR-CPN Schrödingerequation at the FCI level in our finite-sized basis setsof two different sizes, and (D) we extrapolated the finitebasis set results to estimate the energies at the c omplete b asis s et (CBS) limits. Finally, sub-section (E) describeshow we added the corrections due to special relativity,QED, and due to the atom having an unclamped, zero-radius, nucleus. Table III. Summary of our main result. All energies arebetween the fine centre of gravity (fcog) of C (cid:0) P (cid:1) and thefcog of C + (cid:0) P (cid:1) so the experimental spin-orbit lowering of12.6725 cm − (calculated in our Supplemental Material, andbased on measurements reported in [8]) needs to be subtractedfrom all numbers to obtain the C (cid:0) P (cid:1) ← C + (cid:0) P / (cid:1) energy.The experimental uncertainty is a 68% confidence interval,meaning that there is a 32% chance that the true energy isoutside the range spanned by the uncertainty.Hamiltonian Ionization Energy (Calc - Obs)[ cm − ] [ cm − ]NR-CPN 90 863.037X2C -30.023Breit & QED -0.48DBOC -0.235Total (theory) Present 90 832.299Experiment 2017 [8] 90 833.021(9) -0.722Experiment 1998 [25] 90 833.171(15) -0.872Experiment 1966 [26] 90 833.122(100) -0.823Theory 2010 [27] 90 838.75 5.74Theory 2017 [28] 90 840.16 7.15Theory 2015 [10] 90 786.66 -46.35 A. Optimization of ‘tight function’ exponents for theaug-cc-pCV7Z and aug-cc-pCV8Z basis sets
The largest orbital basis sets known for C prior tothis work were the (aug-cc-pV X Z, X =7,8,9) sets usedby Feller in 2016 [29]. These basis sets did not contain‘tight’ exponent functions for capturing the effects of thecorrelation between the core (1 s , s ) electrons and thevalence electrons (2 p ). The largest known basis set forcarbon prior to this work including the CV (core-valence)correction was the aug-cc-pCV6Z [30] set. In this workwe start by optimizing the ‘tight’ exponents for the CVcorrection to Feller’s 2016 aug-cc-pV7Z and aug-cc-pV8Zbasis sets, yielding the first aug-cc-pCV7Z basis set forcarbon, and the first aug-cc-pCV8Z basis set known forany element.The final aug-cc-pCV X Z basis sets have X new tightfunctions of s -type, X − p -type, X − d -type,and so forth, up to the final i -type function for X = 7and the final k -type function for X = 8. The j th ex-ponent corresponding to a function of type L is named γ X,L,j , and is assumed to follow an ‘even-tempered’model: γ X,L,j = α X,L,j β j − X,L,j .In the non-linear optimization procedure to obtain α ,L,j and β ,L,j , the starting values were chosen to bethe α ,L,j and β ,L,j values that were already optimizedin [30]. These were then treated as free parameters tominimize the difference between the frozen core and all-electron CISD energies of the carbon atom with all otherexponent functions fixed. The MOLPRO program [31] wasused to calculate the CISD energies, and the L - BFGS - B program of [32] was used to optimize the free parame-ters. For X = 7, the s -type functions were added first,then once they were optimized they were held fixed whilethe p -type functions were added and optimized. Thenboth the s - and p -type functions were held fixed whilethe d -type functions were added, and so on up to thesingle i -type function. The procedure for X = 8 was thesame, except the procedure continued to k -type func-tions, and the starting values came from the newly op-timized X = 7 case rather than the X = 6 case from[30]. MOLPRO does not support k -functions, so to opti-mize the k -function we calculated the CISD energy atthree points using GAUSSIAN [33] and estimated the valueof α , , yielding the lowest energy by using a quadraticfit.The tight exponents optimized in this work for aug-cc-pCV7Z and aug-cc-pCV8Z are presented in the Supple-mental Material. B. Calculation of 1- and 2-electron integrals including k -and l - functions The calculation of the 1- and 2-electron integrals for(aug)-cc-p(C)V X Z basis sets with X ≥ k - and l - functions, butfor first row elements, k -functions appear in X = 7 ba- Table I. Upper bounds for total non-relativistic electronic energies. VO stands for variational optimization (parameters in awavefunction ansatz are optimized to yield the lowest energy). Hylleraas-Log indicates the use of Hylleraas functions sup-plemented with auxiliary log functions, and ECG( M ) stands for explicitly correlated Gaussian ansatz with M optimizableparameters. Numbers in parentheses are estimated uncertainties within the method used, so for FN-DMC does not includefixed-node error and for FCIQMC does not include basis set error. No numbers were obtained with basis set extrapolations.Total non-relativistic clamped point-nucleus (NR-CPN) energy [Hartree] Method/Ansatz type ReferenceH 1 -0.5 Analytic & Exact 1926 SchrödingerHe 2 -2.903 724 377 034 119 598 311 159 245 194 404 446 696 925 309 838 VO/Hylleraas-Log 2006 Schwartz [5]Li 3 -7.478 060 323 910 134 843 VO/Hylleraas 2017 Wang [1]Be 4 -14.667 356 494 9 VO/ECG(4096) 2013 Puchalski [3]B 5 -24.653 867 537 VO/ECG(8192) 2015 Puchalski [4]C 6 -37.844 48(2) FN-DMC/CISD/cc-pV5Z
FCIQMC/aug-cc-pCV8Z
Present Work -C 6 −37.843 333 VO/ECG(1000) 2013 Bubin [16]Table II. The most precisely calculated ionization energies for the first six atoms, compared to the best known experimentalmeasurements to date. The last column indicates that if aiming for the best precision, an experimental measurement is stillthe best way to obtain the energy for most atoms, but for Be, the energy has been obtained more precisely in silico than inany experiment to date. The value for carbon of 90 832.299 cm − was calculated in the present work.Transition Experiment Theory Calc - Obs | Calc - ObsUncertainty in obs | More precise[cm − ] [cm − ] [cm − ] H H + (cid:0) S (cid:1) ← H (cid:0) S (cid:1)
109 678.771 732(23) [17] 109 678.771 743 07(10) a He He + (cid:0) S (cid:1) ← He (cid:0) S (cid:1)
198 310.666 37(2) [20] 198 310.665 07(1) [21] -0.001 3 65.00 Theory Li Li + (cid:0) S (cid:1) ← Li (cid:0) S (cid:1)
43 487.159 40(18) [22] 43 487.159 7(7) [1] -0.000 3 1.66 Experiment Be Be + (cid:0) S (cid:1) ← Be (cid:0) S (cid:1)
75 192.64(6) [23] 75 192.699(7) [3] 0.059 0.98 Theory B B + (cid:0) S (cid:1) ← B (cid:0) P (cid:1)
66 928.036(22) [24] 66 927.91(21) [4] -0.126 5.73 Experiment C C + (cid:0) P (cid:1) ← C (cid:0) P (cid:1)
90 833.171(15) b [25] 90 832.299 - -0.872 58.13 Experiment a This number is based on the data in [18] although it is not explicitly written anywhere there. Two of the authors of [18] havepresented the number in Table III of [19]. b After the completion of this work two new values in disagreement with each other,90 833 . . sis sets and l -functions appear when X = 8. To calcu-late these integrals, we have used a locally modified ver-sion of MOLCAS [34] in order to support larger basis sets.The 1- and 2-electron integrals for C and C + were eval-uated in the basis of the optimized CASSCF(6,5) andCASSCF(5,5) orbitals respectively, with the five activeorbitals being the 1 s , 2 s , 2 p x , 2 p y and 2 p z of the Catom/ion. This active space is the minimal active spaceincluding all electrons, that is able to provide balancedorbitals for the three degenerate states of the P state ofthe C atom, or the P state of the C + ion. C. Calculation of NR-CPN energies in finite basis setswithout truncating the possible excitation levels (FCIQMC)
A deterministic FCI (full configuration interaction)calculation for the 5e − C + ion in the aug-cc-pCV7Zbasis set would require almost 55 TB of RAM, and forthe neutral atom would require more. Therefore we useFCIQMC for all NR-CPN calculations. The method was introduced in [35], and we use the initiator method firstdescribed in [36], and the semi-stochastic method as de-scribed in [37]. The calculations are performed using thedeveloper version of the software NECI [38].Within a given Hamiltonian (in this case the NR-CPNHamiltonian) and basis set, there are three sources oferror in the FCIQMC energy calculations:1. Trial wavefunction error (∆ E trial ), which ap-proaches zero in the limit where the numberof determinants used in the trial wavefunctionapproaches the number of determinants in theFCIQMC wavefunction;2. Initiator error (∆ E initiator ), which approaches zeroin the limit where the number of walkers N walkers gets sufficiently large; and3. Stochastic error (∆ E stoch ), which for a given num-ber of walkers and trial wavefunction determi-nants is estimated as the square root of the un-biased variance among different estimates E i ofthe energy from their mean ¯ E after different num-bers N of Monte Carlo macro-iterations (deter-mined using the Flyvbjerg-Petersen blocking analy-sis [39]) after the walkers have reached equilibrium:∆ E stochastic ≈ r P Ni =1 ( E i − ¯ E ) N − = O ( / √ N ).Our goal was to obtain all energies to a precision of ± ǫ where ǫ ≤ µE Hartree ≈ . − (within the basissets used). To ensure that ∆ E initiator can be neglected,we used a sufficiently large value of N walkers for every en-ergy calculation, so that the energy difference between us-ing N walkers and N walkers was smaller than 1 µE Hartree .Likewise, to ensure that ∆ E trial can be neglected, weused a sufficiently large number of determinants in thetrial wavefunction for every energy calculation, such that∆ E trial would also be smaller than 1 µE Hartree . We thenran every calculation for enough macro-iterations N suchthat ∆ E stoch was smaller than ∆ E trial and ∆ E initiator .Further details are presented in the Supplemental Mate-rial, including tables which show that all three sources oferror in our final numbers are not larger than we claim. Table IV. Final NR-CPN energies. The break-down of howthese energies were obtained, and the Hartree-Fock energiesthat were used for the extrapolations are available in the Sup-plemental Material. Numbers within parentheses indicate un-certainties in the last digit(s) shown, and their determinationis described in the Supplemental Material. C (cid:0) P (cid:1) C + (cid:0) P (cid:1) P → P [ E Hartree ] [ E Hartree ] [cm − ]aug-cc-pCV7Z -37.844 251 5(05) -37.430 345 1(01) 90 841.955(028)aug-cc-pCV8Z -37.844 355 5(08) -37.430 412 5(05) 90 849.987(054)Eq.(3), n = 3 . n = 4 -37.844 514 2 -37.430 514 3 90 862.471Mean -37.844 521 4 -37.430 519 0 90 863.037 D. Extrapolations to the CBS (complete basis set) limit
We use two different families of formulas to extrapo-late the correlation energies (FCIQMC energies with theHartree-Fock energies subtracted out) from E X − and E X to E CBS : E CBS = E X − AX n , (1) E CBS = E X − A ( X + / ) n . (2)If we set n = 3 in Eq. (1), we recover the formula orig-inally proposed in [40]. If we set n = 4 in Eq. (2), werecover the formula originally proposed in [41]. If wehave values for E X at two different X values, we can eliminate A in both cases, so Eq.(1) leads to Eq.(3) andEq.(2) leads to Eq.(4): E CBS = X n E X − ( X − n E X − X n − ( X − n , (3) E CBS = E X + (2 X − n ( E X − E X − )(1 + 2 X ) n − (2 X − n . (4)As explained on page 5 of [29], extrapolations to theCBS limit using n = 3 in Eq.(1) tend to over-shoot theCBS limit. The value of n = 3 . E CBS obtained from using n = 3 . n = 4 in Eq.(4), for X = 8 were added to the Hartree-Fock energies for X = 8 and are presented in Table IV.The final NR-CPN energy was taken as the mean of bothvalues obtained from extrapolating the correlation energyand adding it to the Hartree-Fock energy with X = 8. E. Estimation of relativistic, QED, finite nuclear mass, andfinite nuclear size corrections
Scalar relativistic corrections were calculated by com-paring the energies using the spin-free version of the 1e − X2C (exact 2-component) Hamiltonian [42], to the ener-gies of the NR-CPN Hamiltonian. The integrals of ourX2C Hamiltonian with ROHF orbitals were done in the
CFOUR program, and were calculated at various levels ofcoupled cluster theory with the
MRCC program [43].Further scalar relativistic effects were included byadding the Breit and QED corrections (including the vac-uum polarization and the self-energy terms that togethercomprise a Lamb-like shift) from the state-averagedDirac-Fock calculations done in [27]. The overall contri-bution from the Breit and QED correctons (combined)to the IE for C was -0.48 cm − .Diagonal Born-Oppenheimer breakdown corrections(DBOC) to the clamped nucleus approximation [44] werecalculated using CFOUR , with
MRCC used for the coupled-cluster part. Our value of -0.235 cm − is about triple thevalue of -0.08 cm − estimated in [27], due to includinghigher levels of correlation.The basis set and correlation convergence of the X2Cand DBOC corrections is shown in the Supplemental Ma-terial, along with our calculation of the finite nuclearsize correction of 0.00543 cm − . The final correctionsthat contributed to our final computer-predicted ioniza-tion energy are presented in Table III. II. CONCLUSION
Table III summarizes the various contributions to ourvalue of the IE, and compares our final value to experi-ment and to three recent theoretical estimates. Our valueis 0.872 cm − smaller than the best experimental value.The best theoretical estimate of the IE before this workwas in [27], and was in disagreement with experiment bymore than a factor of 6.5 more than our present result.We believe that this could have been due to any or allof three things: (1) approximations inherent to the F12approach used for their NR-CPN energy, (2) the pertur-bative nature of their scalar relativistic corrections (i.e.using the mass velocity and Darwin terms, rather thanthe X2C Hamiltonian used in the present work) and (3)the CCSD approximation made in their DBOC correc-tion (as opposed to the CCSDTQ used in the presentwork which we have shown appears to be converged tothe FCI limit). III. ACKNOWLEDGMENTS
We wish to thank Mariusz Puchalski, KrzysztofPachucki, Robert Moszynski, Jacek Komasa, GordonDrake, Michal Lesiuk, Michal Przybytek, and Wim Klop-per for helpful discussions, comments and suggestions.We also thank Alexander Kramida and Kunari Harisof NIST for information about their recent pre-print ona newer experimental ionization energy for carbon (men-tioned in the footnote to Table II), Alexander Kramidafor information about the experimental ionization of hy-drogen used in Table II, and Barry Taylor and PeterMohr of NIST for information about the theoretical ion-ization energy of hydrogen used in Table II.
REFERENCES [1] L. M. Wang, C. Li, Z.-C. Yan, and G. W. F. Drake,Physical Review A , 032504 (2017); G. W. Drake,Private Communication from 4 to 5 June 2018.[2] M. Puchalski, D. Kędziera, and K. Pachucki, Phys. Rev.A , 062509 (2010).[3] M. Puchalski, J. Komasa, and K. Pachucki, PhysicalReview A , 030502 (2013).[4] M. Puchalski, J. Komasa, and K. Pachucki, PhysicalReview A , 062501 (2015).[5] C. Schwartz, ArXiv Mathematical Physics e-prints(2006), math-ph/0605018.[6] M. Puchalski, D. Kedziera, and K. Pachucki, PhysicalReview A , 032503 (2013).[7] W. L. Glab, K. Haris, and A. Kramida, Journal ofPhysics Communications , 055020 (2018).[8] K. Haris and A. Kramida, (2017), arXiv:1704.07474.[9] D. Bresteau, C. Drag, and C. Blondel, Phys. Rev. A ,013414 (2016).[10] Y. Yang, I. Kylänpää, N. M. Tubman, J. T. Krogel,S. Hammes-Schiffer, and D. M. Ceperley, The Journalof Chemical Physics , 124308 (2015).[11] R. E. Thomas, G. H. Booth, and A. Alavi, PhysicalReview Letters , 033001 (2015).[12] D. Cleland, G. H. Booth, C. Overy, and A. Alavi, Journalof Chemical Theory and Computation , 4138 (2012).[13] A. D. Powell, N. S. Dattani, R. F. K. Spada, F. B. C.Machado, H. Lischka, and R. Dawes, The Journal ofChemical Physics , 094306 (2017).[14] C. Daday, S. Smart, G. H. Booth, A. Alavi, and C. Fil-ippi, Journal of Chemical Theory and Computation ,4441 (2012).[15] G. H. Booth, A. Grüneis, G. Kresse, and A. Alavi, Na-ture , 365 (2012). [16] S. Bubin, M. Pavanello, W.-C. Tung, K. L. Sharkey, andL. Adamowicz, Chemical Reviews , 36 (2013).[17] A. Kramida, Atomic Data and Nuclear Data Tables ,586 (2010).[18] U. D. Jentschura, S. Kotochigova, E.-O. Le Bigot, P. J.Mohr, and B. N. Taylor, Phys. Rev. Lett. , 163003(2005).[19] P. J. Mohr, D. B. Newell, and B. N. Taylor, Rev. Mod.Phys. , 035009 (2016).[20] D. Z. Kandula, C. Gohle, T. J. Pinkert, W. Ubachs, andK. S. E. Eikema, Physical Review A , 062512 (2011).[21] K. Pachucki, V. Patkóš, and V. A. Yerokhin, PhysicalReview A , 062510 (2017).[22] B. A. Bushaw, W. Nörtershäuser, G. W. F. Drake, andH.-J. Kluge, Phys. Rev. A , 052503 (2007).[23] R. Beigang, D. Schmidt, and P. J. West, Le Journal dePhysique Colloques , C7 (1983).[24] A. E. Kramida and A. N. Ryabtsev, Physica Scripta ,544 (2007).[25] E. S. Chang and M. Geller, Physica Scripta , 326(1998).[26] L. Johansson, Arkiv Fysik , 201 (1966).[27] W. Klopper, R. A. Bachorz, D. P. Tew, and C. Hättig,Phys. Rev. A , 022503 (2010).[28] D. Feller, The Journal of Chemical Physics , 34103(2017).[29] D. Feller, The Journal of Chemical Physics , 014105(2016).[30] A. K. Wilson, T. van Mourik, and T. H. Dunning,Journal of Molecular Structure: THEOCHEM , 339(1996).[31] H.-J. Werner et al. , 550 (1997).[33] M. J. Frisch et al. , “Gaussian 16 Revision A.03,” (2016),Gaussian Inc. Wallingford, CT.[34] F. Galvan et al. , Journal of Chemical Theoryand Computation , 5925 (2019), pMID: 31509407,https://doi.org/10.1021/acs.jctc.9b00532.[35] G. H. Booth, A. J. W. Thom, and A. Alavi, The Journalof Chemical Physics , 054106 (2009).[36] D. Cleland, G. H. Booth, and A. Alavi, The Journal ofChemical Physics , 041103 (2010).[37] N. S. Blunt, S. D. Smart, J. A. F. Kersten, J. S. Spencer,G. H. Booth, and A. Alavi, The Journal of ChemicalPhysics , 184107 (2015).[38] “Public version is available here: https://github.com/ghb24/NECI_STABLE ,”.[39] H. Flyvbjerg and H. G. Petersen, The Journal of Chem-ical Physics , 461 (1989).[40] W. Kutzelnigg and J. D. Morgan, The Journal of Chem-ical Physics , 4484 (1992).[41] J. M. Martin, Chemical Physics Letters , 669 (1996).[42] K. G. Dyall, The Journal of Chemical Physics , 9618(1997); L. Cheng and J. Gauss, ibid . , 244104 (2011).[43] M. Kallay et al. , “MRCC, a quantum chemical programsuite.”; M. Kallay and P. R. Surjan, The Journal ofChemical Physics , 2945 (2001).[44] W. Kutzelnigg, Molecular Physics , 909 (1997);J. Gauss, A. Tajti, M. Kallay, J. F. Stanton, and P. G.Szalay, The Journal of Chemical Physics125