Faraday Symp. Chem.SOC.,1984 19 165-173 Numerical Perturbation Calculations for Diatomic Molecules BY EDWARDA. MCCULLOUGH AND KENTW. RICHMAN JR,* JOHN MORRISON? Department of Chemistry and Biochemistry UMC 03 Utah State University Logan Utah 84322 U.S.A. Received 7th August 1984 A numerical method is described for diatomic-molecule perturbation calculations based on a natural orbital expansion of the pair function. Results are presented for several diatomic molecules and are compared with results obtained with finite basis sets. A discussion of basis-set errors and the advantages and disadvantages of the numerical method is given. For atoms it is common to carry out many-body perturbation calculations entirely numerically. 9 For molecules however nearly all perturbation calculations presently use basis-set expansion^.^ The inaccuracies introduced in this way can be difficult to estimate since the convergence of such expansions is often quite erratic.In this article we report on numerical second-order perturbation calculations for several diatomic molecules. Our purposes are twofold. First we wish to test the feasibility of obtaining numerical solutions of the pair equation fpr diatomic molecules. Secondly we would like to examine some of the best published basis set calculations with the idea of determining the levels of basis-set error in such calculations generally. THEORY In the diagrammatic perturbation expansion of the energy of a closed-shell system2 the second-order direct diagram D,,has the value We consistently shall use a b and c for occupied (core) orbitals r s and t for unoccupied (virtual) orbitals and i and j for unspecified orbitals.Here these orbitals are chosen to be eigenfunctions of the closed-shell Fock operator h. Eqn (1) may be evaluated readily using the first-order pair function lyd(ab-+ rs) = C I rs) (rs 1 r;,l I ab)/(&,+ &t -E -cS). (2) rs An inhomogeneous equation for 'yd is obtained by operating on eqn (2) from the left with E 4-&b-h(1)-h(2). The result is t Present address Computer Science Department Box 1679 Station B Vanderbilt University Nashville Tennessee 37235 U.S.A. 165 NUMERICAL PERTURBATION CALCULATIONS FOR DIATOMICS Similarly the pair function for the exchange diagram satisfies the equation [E +E~ -h(1) -h(2)] W,(ab -+ rs) = -Z I rs) (rs 1 r; I ba).(4) rs If we add eqn (3) and (4) we obtain a determinant Itya ybl,on the right. It is more convenient however to use a symmetry-adapted linear combination of determinants since the spin dependence can then be factored out and spin eliminated from the pair equation. We shall use a bracket notation to denote a symmetry adapted spatial function. For example the functions for i # j are where the upper (lower) sign refers to the singlet (triplet) case. With this notation the symmetry-adapted pair function yUb,satisfies the equation where it is understood that all quantities now are functions only of spatial variables. For atoms eqn (6) can be solved numerically with high efficiency and accuracy.* While the same technique could be applied to molecules the appearance of double partial wave sums would make the calculations quite time consuming.For this reason we shall here pursue another approach. We begin by noting that eqn (6) can be obtained by requiring that the functional be stationary with respect to any first-order variation Styab consistent with strong orthogonality to the core orbitals. The stationary value of the functional cab is the symmetry -adap ted pair energy. One can automatically satisfy strong orthogonality by expanding yabin terms of the virtual orbitals. This approach is often used for basis-set calculations where in the course of obtaining the occupied orbitals one obtains as a by-product an approximate set of virtual orbitals.It is much less attractive for numerical calculations because each orbital (virtual as well as occupied) must be found independently. Also if h is the closed-shell Fock operator all of the virtual orbitals usually lie in the continuum. An alternative approach is to use a natural orbital5 expansion of yUbas suggested by Kutzelnigg,g and make cab stationary with respect to both variations in the expansion coefficients and the forms of the pair natural orbitals (PNO). This leads to a set of equations that are formally quite similar to the multiconfiguration Hartree-Fock equations. The natural orbital expansion is of the form k For a singlet pair in which orbitals a and b have the same spatial symmetry it is possible to choose uk and vk to be identical.This leads to certain simplifications which we shall not consider more fully. Henceforth we shall regard uk and uk as distinct. Substituting eqn (8) into eqn (7) gives where the ck are assumed to be real. E. A. McCULLOUGH JR J. MORRISON AND K. W. RICHMAN The stationary condition with respect to ck gives The stationary condition with respect to say uk consistent with orthonormality of the PNO and strong orthogonality to the core leads to the equation where the d are Lagrange multipliers and Picr)( 1) = vi(1) Idz w*(2) ry,lvi(2) for an arbitrary function cr). An analogous equation can be derived for vk. In addition to variations of individual PNO rotations between pairs of PNO must be allowed. For example consider a rotation of two u uk -+ uk+8uland ul-+ ul-8uk.Expanding &,b in powers of 8 to first order we obtain so that the stationary condition with respect to rotations implies that cf Llk = ct Lkl. If this condition is not satisfied automatically the PNO must be rotated until it is. The optimum rotation angle can be obtained by expanding &abto higher order in 8 and setting deab/d8 equal to zero. One must also consider of course rotations between pairs (uk v,) and (vk vl). The PNO equations also can be derived by substituting eqn (8) directly into eqn (6). This is legitimate only if the PNO expansions could be a solution of eqn (6) i.e. only if the expansion is infinite. Taking the scalar product of the result from the left with say vk(2) we obtain ck(Eu+Eb-h)uk-x Cl((vk~hlvl)ulk(vkIhl ul>vl) 1 = c <'k 1 s> vr (rsI rTtI va vbfvb va>.(14) rs The right-hand-side of eqn (14) can be simplified using the completeness of eigen- functions of h,which allows us to express each sum over excited states Cvt(r)v/F(r') as a delta function minus a sum over core states. After some rearrangement eqn (14) becomes +ck(ca+cb)uk+~ (vcvkIr~~l~avb*vbva) vc. (15) c By comparing eqn (15) with eqn (11) we may obtain expressions for the Lagrange multipliers. In particular we see that ck dlk = -Cl (Vk I hI 01) I # k. (16) It is then shown easily that ci Alk = cf dkl. The Lagrange multipliers also may be evaluated by taking the scalar product of eqn (1 1) from the left with each PNO. Equating these expressions with those obtained NUMERICAL PERTURBATION CALCULATIONS FOR DIATOMICS from eqn (1 5) we can derive consistency conditions for an exact PNO solution of the pair equation such as -C1<UkIhlU1) =Ck<UIIhIUk)+(UIUkIY~~IWaWbfWbWa) Ifk.(17) To the extent that the PNO expansion is exact these consistency conditions are satisfied and cabis automatically stationary with respect to rotations. This does not follow for a truncated expansion however where the effect of rotations must be considered and departures from the consistency conditions can be expected. Our first problem of course is how to solve eqn (I 1) for molecules. For diatomic molecules it is possible to obtain numerical solutions using a partial wave expansion of the PNO in prolate spheroidal coordinates.For instance we write up(<,q,4) =2Xklm(<) qm(v,6) (18) where the qmare spherical harmonics and the Xklm are numerical functions. The symmetry-adapted PNO for diatomic molecules can be labelled with m quantum numbers e.g. upup’. The advantages and disadvantages of this approach as well as computational details have been reported previously for the multiconfiguration Hartree-Fock problem. The principal advantage of eqn (18) is that it converges smoothly which allows us to determine PNO with controlled numerical accuracy. The principal disadvantage is that it is only applicable to diatomic molecules. RESULTS The formal similarity between the pair natural orbital and multiconfiguration Hartree-Fock energy functionals enabled us to use a very slightly modified version of our numerical MCSCF code to solve the PNO equations.All computations were carried out in double precision on a VAX 11/780. Our first example is atomic beryllium which we treated as a heteronuclear diatomic molecule with zero charge on one nucleus and the internuclear distance set arbitrarily to 0.5a0.This is an important test case since very accurate pair energies are available for Be. Table 1 shows the results. The pair energies Eab(mm’),obtained from expansions containing N(mm’) PNO configurations are compared with the essentially exact cab(mm’)calculated from the work of Malinowski et al. (MPJ).8 The MPJ results are equivalent to infinite PNO expansions for each (rnrn’). The results in table 1 are indicative of the general level of accuracy that may be expected from application of the PNO approach to molecules.One sees that fairly short expansions can produce Eab(mm’)correct to d 5 x hartree,* and that the initial convergence rate of the PNO expansion is quite high as csS2(a2)demonstrates. The second-order energy E2(6),obtained by summing our degeneracy-weighted Eab(mm’),is -0.0715 hartree. The MPJ values of E2(6)and E2(w)are -0.0743 and -0.0763 hartree respectively The error in our E2(6)due to truncating the PNO expansions is of the same order as that due to neglect of (rnrn’) >2. Table 2 gives our results for H and LiH. Here only the total pair energies &,b for (mm’)d 2 are shown. For the energetically most important pairs we used two sets of PNO expansion lengths to investigate the convergence with respect to N(mm’).Our best E2(6)for LiH is -0,0696 hartree which may be compared directly with * 1 hartree =4.359814 x lo-’* J. E. A. McCULLOUGH JR J. MORRISON AND K. W. RICHMAN Table 1. Numerical PNO and accurate pair energies (hartree) for Be accurate ab 4 mm’ N(mm’) PNO Eab(mm’) Eab(mm’) 9 -0.0206 -0.0209 6 -0.0164 -0.0168 3 -0.0015 -0.0018 4 -0.0101 -0.0108 6 -0.0163 -0.0168 3 -0.00 17 -0.0020 8 -0.0017 -0.0017 6 -0.001 3 -0.0013 3 -0.0001 -0.000 1 1 -0.0002 -0.0002 1 -0.0004 -0.0005 Table 2. Numerical PNO pair energies (hartree) for H (R= 1.4a,) and LiH (R = 3.015a0) ab -+ mm’ N(mm’) &ab N(mm’) ‘ab H2 10; -+ Cr2,712,s2 9 5 3 -0.0325 13,9 5 -0.0331 LiH 102 -+ 02,712,62 9 6 3 -0.0375 13 9 5 -0.0381 202 48,712 d2 9 6 3 -0.0286 13 9 5 -0.0290 l02o(lC) + u2,n2 4 3 -0.0013 1a2a(3q jog‘ 71.n’ 1 1 -0.0004 the value of -0.0654 hartree obtained by Bartlett and Silverg (BS) using a large Slater basis set.Their basis did not contain any functions beyond 6 but they did sum over all virtual orbitals. Thus their calculation is equivalent to an (mm’)d 2 PNO calculation within their finite basis and the difference between their E2(6)and ours is a measure of the basis-set error in their calculation. We now turn to FH. This molecule is considerably more complex than LiH because of the greater number of electrons and the presence of several spatially proximate pairs. For these reasons FH provides a more realistic and more stringent test of high-accuracy methods.Systematic variation of all PNO expansion lengths was impractical in the case of FH owing to the large number of pairs. For the intrashell ozpairs we chose initial expansion lengths based on our experience with smaller systems. The occupied 7t shell was a new feature; however we were aided by the similarity between correlation in the 71 shells of FH and Ne. We computed several 7t2 pair energies for Ne and compared them with accurate values.l0 We then used these comparisons as a guide in choosing initial N(mm’) for FH 7t2 pairs. One advantage of the PNO approach is that one can estimate roughly the maximum energy that could be gained by adding another PNO configuration from the energy contribution of the least important configuration already present.The relatively greater complexity of FH was revealed clearly in the intershell pair calculations. The iterative solution of the PNO equations requires initial guesses for the PNO and these were often problematic for intershell pairs While the FH occupied NUMERICAL PERTURBATION CALCULATIONS FOR DIATOMICS Table 3. Numerical PNO and basis-set pair energies (hartree) for FH (R= 1.7328~~) ab PNO &,b basis set &,b 202 -0.01 12 -0.0100 2034'Z) 203u(~Z) 302 -0.0156 -0.0028 -0.0273 -0.0138 -0.0026 -0.0251 20ln('I-I) -0.0138 -0.0124 20 174311) 3~17c(~II) 301~(~ll) 17cz (average singlet) -0.0027 -0.0136 -0.0088 -0.0199 -0.0027 -0.0130 -0.0088 -0.0189 1q 3 ~ ) -0.0090 -0.0089 orbitals strongly resemble slightly distorted fluorine orbitals attempts to utilize initial guesses based on atomic symmetry were only moderately successful.For intershell pairs involving 20 or 30 especially the PNO frequently displayed a tendency to move away from the fluorine by large distances in either direction. Plots of these PNO exhibited few vestiges of atomic symmetry. More variation of expansion lengths was carried out for these pairs. In addition to problems arising from choices of initial guesses and expansion lengths convergence of the iterative solution procedure was poorer in general for all aa pairs of FH than for smaller systems. Possibly this was due to the more numerous strong orthogonality constraints. Some pair energies for FH are shown in table 3 where they are compared with the large Slater basis results of BS.Pairs involving la are excluded as the BS basis set contained very few high-exponent functions and could not be expected to represent correlation involving lo very well. Our E2(6)is -0.331 hartree. Based on our experience and limited testing we estimate that this is not less than 95% of the exact E2(6).The BS E2(6)is -0.306 hartree. Almost one third of the difference between their result and ours is due to the la2 pair. This probably is not a serious error since few molecular properties should be sensitive to 1 a2correlation. Discussion of the remaining pairs will be deferred until the next section. For FH we also tested the sensitivity of the pair energies to the accuracy of the occupied orbitals.We chose for this test orbitals computed in a moderate sized Slater basis by Nesbet.ll His wavefunction gives an energy of -100.057 hartree compared with the exact Hartree-Fock energy of -100.071 hartree.12 This basis is of slightly better quality than typically might be used for large molecules although much inferior to the BS basis sets. We computed for Nesbet's wavefunction eight of the energetically most important Eab(mm'),excluding any involving 1a.The largest deviation from our results calculated with the accurate numerical wavefunction was 2% and most differences were smaller than 1%. Basis-set error in the occupied orbitals is unlikely to be important in high-quality molecular-perturbation calculations.DISCUSSION Before embarking on a detailed discussion of the numerical PNO method vis-a-vis basis-set methods we must identify the essential differences between the two approaches. E. A. McCULLOUGH JR J. MORRISON AND K. W. RICHMAN Table 4. Numerical PNO optimized Slater orbital (OSO) and optimized gaussian geminal (OGG) second-order energies (hartree) for H and LiH molecule OSO E2(S) PNO E2(S) OGG E2(cn) H2 -0.0331 -0.0331 -0.0342 LiH -0.0678 -0.0696 -0.0722 Errors in the numerical PNO method are of three types :(a)error in individual PNO (b)error introduced by truncating the PNO expansion for each (mm’) and (c) error due to truncating the sums over (mm’).Any expansion of the pair function in a set of orthonormal orbitals can be reduced to PNO form by a suitable unitary transformation.This includes the conventional sum-over-virtual orbitals basis-set method where both the PNO expansion lengths and the maximum (mm’) are determined by the basis set. The only essential difference then between the numerical PNO method and most basis-set methods is in the representation of the PNO and the only practical advantage of solving the PNO equations numerically is that type (a)errors can be made negligible at least for the relatively short expansions we used. Recently two less conventional basis-set methods have been developed. The method of Adamowicz and Bartlett13 also is restricted to diatomic molecules. In their initial application they employed an expansion of each pair function in simple Slater-orbital products; however they were able to use long expansions and perform complete exponent optimizations.Their results are the best obtainable from a given number of Slater orbitals with fixed (nlm)quantum numbers. The method of Monkhorst and utilizes explicity correlated gaussian geminal expansions of the pair functions while avoiding many of the difficulties normally associated with this approach. Long germinal expansions with extensive optimization are a feature of their work. Unlike almost all other methods this method does not involve truncated (mm’) sums. As yet both of these methods have been applied only to small systems. A comparison with our numerical PNO method is given in table 4. The optimized Slater orbital E2(6)for H is in excellent agreement with ours.Recent work indicates that both are within 1% of the exact E2(6).15The situation for LiH is less satisfactory. Although Adamowicz and Bartlett used 26 completely optimized Slater functionsperpair they only reduced the error in the BS E2(6)by ca. 60%. Note that the PNO representation gives the shortest orbital expansion of a two-electron f~nction.~ Therefore one could not expect to improve on our E2(6)for LiH with any basis containing fewer than 27 functions since that is the number of PNO in our longest expansion. The convergence exhibited by the PNO expansions in tables 1 and 2 indicates that much of the difference between our E2(6)and the optimized geminal E2(co)truly is type (c) error. The type (a)-(b) error in the optimized Slater orbital results for LiH is of the same order as the type (c)error.Type (a)-(b) errors are even more serious in the BS calculations. We emphasize here that we have not singled out the BS calculations for special criticism. On the contrary their work is of very high quality. Basis sets for molecular perturbation calculations are rarely larger and frequently much smaller than those used by BS. Thus it is disconcerting to find that the type (a)-(b) error in their E2(6)for LiH is the most important error in their calculation by a factor of 1.6. Even more disconcerting is the behaviour of some of the FH pairs involving the bonding 30 orbital. Actually BS employed several basis sets. In passing from a 22 NUMERICAL PERTURBATION CALCULATIONS FOR DIATOMICS function basis to 32 functions (their largest) the changes in the 2a3a(%) and 3a2 pair energies were -0.0017 and -0.0020 hartree respectively.As table 3 shows the residual errors for these same pairs are at least -0.0018 and -0.0022 hartree so only about half of the error has been eliminated with 10 additional basis functions. In general triplet correlation is more accurately described than singlet correlation in the BS calculation on FH even when the former is large (e.g. for the 3aln and In2pairs). 10% errors in singlet pair energies are common. The total correlation energy of FH has been estimated to be -0.381 hartree.16 The value for Ne is almost exactly the same and in Ne E2(m)amounts to > 99% of the total.l* Our calculations indicate that E2(6)for FH is not less than -0.331 hartree but probably not more than -0.35 hartree based on our 5% estimate for type (b) errors.Type (c) errors thus should be between -0.03 and -0.05 hartree. The type (a)-(b)error in the BS E2((6)is at least -0.025 hartree and may be as much as -0.04 hartree. Many people have noted that an error in the correlation energy is important only if it affects some molecular property of interest. With respect to this point we would only make two observations. First empirical evidence certainly suggests that some molecular properties can be computed reliably with basis sets much smaller than those used by BS.3 Secondly the difference between E2(m) and the total correlation energy for many light molecules probably does not exceed lo% whereas even the BS basis set for FH leads to second-order energy errors of the order of lo% and the error in the first-order wavefunction itself may be much larger.Extensive error cancellation in higher order seems to be the only way these two observations can be reconciled. Basis-set improvements are straightforward in principle. Type (a)-(b) errors call for more or better functions of a given m.Here the numerical PNO method might make a special contribution since it is possible to project the PNO onto a basis” and determine by hindsight an optimum basis for correlation. Experience gained with a few representative examples might lead to a general prescription for choosing such basis sets. Type (c) errors can be reduced by going to higher m,but the asymptotic convergence rate of the (mm’) sums probably is slow.The evidence for this comes from the slow convergence rate of the analogous (If‘) sums in atoms where extrapolation formulae are often employed to get better estimates of the asymptotic limits.1° Similar formulae might be useful in the diatomic-molecule case. We conclude with some general remarks about the PNO approach in general and the numerical PNO method in particular. The rapid initial convergence of the PNO expansions for fixed (mm’) makes this approach quite attractive if moderate accuracy is the goal. It can certainly be implemented with basis-set methods as well as numerically. The iterative convergence problems we encountered might be overcome with guaranteed convergence algorithms of the type developed for MCSCF calculations.The PNO approach can be extended readily to higher orders. The first-order pair function, eqn (6) suffices for the calculation of third-order energies. This is made more difficult by the non-orthogonality of PNO for different pairs; however the number of PNO is relatively small. An equation for the second-order pair function can be derived for which the first-order pair function is used to calculate a modified right-hand side of eqn (6) in a manner analogous to that developed for atoms.18 This leads to an iterative scheme that could be carried to arbitrarily high order. If one examines the derivation of eqn (1 5) one finds that the right-hand side of eqn (6) is responsible for the terms involving the K operators as well as the terms forcing orthogonality to the occupied orbitals.These terms would have to be re-valuated at each order of perturbation theory. E. A. McCULLOUGH JR,J. MORRISON AND K. W. RICHMAN The PNO approach has certain defects as well. The higher PNO oscillate rapidly which may lead to some difficulty with type (a)errors if very long expansions are used. More seriously the number of PNO giving an energy contribution above a fixed threshold increases as the threshold decreases so that rapid initial convergence is offset by slower asymptotic convergence. Overall truncation error smaller than ca. 5% in E2(m)would probably be difficult to achieve in practical calculations and there is no reason to expect smaller errors at higher orders.Thus iterating to very high order as has been done successfully in atomic coupled-cluster calculation^,^^ seems questionable with the PNO approach. These last remarks of course apply equally well to basis-set methods. For the above reasons we do not view a numerical PNO approach as the optimum method for numerical perturbation calculations even on diatomic molecules. The method is probably most useful in moderate accuracy applications at fairly low orders of perturbation theory. Still there are a number of important problems which could be addressed at this level including basis-set incompleteness questions of the type considered in this paper. We thank the US. National Science Foundation and the Swedish Research Council for partial support of this work.H. Kelly Adv. Chem. Phys. 1969 14 129. I. Lindgren and J. Morrison Atomic Many-body Theory (Springer-Verlag Berlin 1982). R. J. Bartlett Annu. Rev. Phys. Chem. 1981 32 359. J. Morrison J. Phys. B 1973 6 2205 and references therein; A. M. Mirtensson J. Phys. B 1979 12 3995. P-0. Lowdin and H. Shull Phys. Rev. 1956 101 1730. W. Kutzelnigg in Methods of Electronic Structure Theory ed. H. F. Schaefer I11 (Plenum Press New York 1977) chap. 5. E. A. McCullough Jr J. Phys. Chem. 1982 86 2178. P. Malinowski M. Polasik and K. Jankowski J. Phys. B 1979 12 2965. R. J. Bartlett and D. M. Silver J. Chem. Phys. 1975 62 3258. lo K. Jankowski and P. Malinowski Phys. Rev. A 1980 21 45. l1 R. K. Nesbet J. Chem. Phys.1962 36 1518 specifically tables XI1 and XIV. l2 P. A. Christiansen and E. A. McCullough Jr J. Chem. Phys. 1977 67 1877. l3 L. Adamowicz and R. J. Bartlett In?. J. Quantum Chem. to be published. l4 K. Szalewicz B. Jeziorski H. J. Monkhorst and J. G. Zabolitzky J. Chem. Phys. 1983,78 1420; 79 5543. l5 L. Adamowicz and R. J. Bartlett personal communication. C. F. Bender and E. R. Davidson Phys. Rev. 1969 183 23. l7 L. Adamowicz and E. A. McCullough Jr Int. J. Quantum Chem. 1983 24 19. S. Garpman I. Lindgren J. Lindgren and J. Morrison 2.Phys. A 1976,276,167; A. M. Mirtensson Ph.D. Thesis (Chalmers University of Technology Gothenburg 1978). I. Lindgren and S. Salomonson Phys. Scr. 1980,21,335; J. Morrison and S. Salomonson Phys. Scr. 1980 21 343.