NATURAL ORBITAL BASEDENERGY CALCULATION FOR HeH
AND
LiH
2675
A Natural Orbital Based Energy Calculation for Helium Hydride and
Lithium Hydride
by Charles F. Bender and Ernest R. Davidson University of Washington, Seattle, Washington 98106
(Received January 86,1966)
An iterative procedure for simultaneous calculation of the natural orbitals and energy has been developed. This has been applied to the determination of the wave function for the ground electronic state of HeH and LiH near an internuclear separation of 3.0 bohrs. The energies obtained were -3.38280 hartrees for HeH (experimental approximately -3.3849) and -8.0606 hartrees for LiH (experimental -8.0705).
Introduction The LCAO-MO-SCF procedurk developed by Roothaanl has now been applied to a great many molecules. I t is clear that the wave function produced by this method is a good zeroth-order approximation to the true wave function. Attempts to improve this approximation by configuration interaction have, however, been singularly unsuccessful. Ebbing12 for instance, carried out a calculation on LiH using 53 configurations and obtained only 65% of the correlation energy. This lack of success has led to a resurgence of interest in the nonorthogonal or valence bond formalism. Recently, Harris and Taylor3 have published a four-configuration wave function which is nearly as good as Ebbing's, and Browne and Matsen4 have published a 28-configuration function which gave 85% of the correlation energy. These latter calculations have the serious disadvantage that the resulting functions are difficult to visualize and (because of nonorthogonality) difficult to extend to large numbers of electrons. In a landmark paper, Lowdid introduced the concept of natural spin orbitals and the equations for computing them. These orbitals, which diagonalize the first-order density matrix, form an orthonormal set. The occupation number of each orbital gives a direct measure of its importance in the wave function. The natural orbital expansion for two-electron systems has been much s t ~ d i e d . ~ -It~ is now clear that only four configurations built from natural orbitals are required to obtain 90% of the correlation energy of Hz.6
Further, the first of these is essentially the SCF wave function. Judging from Sinanotlu'slo arguments concerning the types of correlation terms that are important, one might expect to nearly duplicate the result of Browne and Matsen for LiH with only seven configurations. One would be the SCF function, three would be inner-shell correlation terms, and three would be outer-shell terms. This prediction has been partially borne out in a recent paper by Ebbing and Henderson" in which a 41-configuration function for LiH built from a-type basis functions was analyzed into a five-configuration function built from natural orbitals with a loss of only 0.001 in the energy.12 The energy they obtained was consistent with extrapolations from two-electron results. From the two-electron results, it should also be expected that the remaining 10% of the correlation energy will be exceedingly difficult to obtain with any orbital basis. (1) C. C. J. Roothaan, Rev. Mod. Phys., 23, 69 (1951). (2) D. D. Ebbing, J. Chem. Phys., 36, 1361 (1962). (3) F. E. Harris and H. 9. Taylor, Physica,'dO, 105 (1964). (4) J. C. Browne and F. A. Matsen, Phys. Rev., 135, A1227 (1964). (5) P. 0. Lowdin, ibid., 97, 1474 (1955). (6) E. R. Davidson and L. L. Jones, J . Chem. Phys., 37, 1966 (1962). (7) E. R. Davidson, ibid., 39, 875 (1963). (8) 5. Hagstrom and H. Shull, Rev. Mod. Phys., 35, 624 (1963). (9) P. 0. Lowdin and H. Shull, Phys. Rev., 101, 1730 (1956). (10) 0. Sinanoglu, J . Chem. Phys., 36, 706 (1961). (11) D. D. Ebbing and R. C. Henderson, ibid., 42, 2225 (1965). (12) All energies are in hartree units and distances are in bohr radii unless otherwise specified.
Volume 70, Number 8 August 1966
2676
The original paper by Lowdin gave two ways for finding the natural orbitals. The first was to solve a complicated set of integrodifferential equations. Although this has now been done for two electrons,13 it does not seem feasible for many electrons. The second was to diagonalize the density matrix. This has the obvious difficulty that the density matrix, and hence the wave function, must be already known. Nevertheless, this latter rnethod has been extensively used for analyzing the wave functions of He, Hz, Be, and LiHI4 into conceptually simpler forms. I n this paper a new procedure for obtaining the natural orbitals is used which makes it possible to use them for calculating the wave function without actually solving the integrodiff erential equations for the orbitals. This procedure consists simply of the following steps: (a) make an initial guess to the natural orbitals in terms of some basis set; (b) construct an approximate wave function from a reasonable number of configurations formed from these orbitals; (c) compute the natural orbitals for this wave function; and (d) using tfiese orbitals as a guess, repeat steps b and c until the wave function and orbitals converge. Results from the use of this procedure are presented in this paper for HeH and LiH. In the preliminary attempts a t this calculation many difficulties were encountered, but these have (for the most part) been solved. The most encouraging aspect of these calculations was that, for reasonable initial guesses and choices of configurations, the iteration procedure was strongly convergent and the energy of each iteration improved until the limit was approached.
Details of the Method In step c, the calculation of the natural orbitals, the spinless density matrix was diagonalized. This avoided the necessity of using different sets of molecular orbitals for different spin. In principle, the transformation to natural spin orbitals could be made, but this would give twice as many natural orbitals as space orbitals and slow down the transformation of integrals to the new basis. Because of the conceptual advantages of keeping spin and space coordinates separated, it was not thought worthwhile to introduce the extra complication of spin orbitals. This choice also avoided a great deal of difficulty in projection of spin eigenfunctions from Slater determinants based on these orbitals. Another difficulty related to step c was that it is not generally possible to form as many natural orbitals as there are nonorthogonal basis functions because of linear dependency problems. The transformation of electron repulsion integrals to an orthogonal basis is The JouTnal of Physical Chemistry
CHARLES F. BENDER AND ERNEST R. DAVIDSON
particularly sensitive to lack of accuracy in the coefficients. The transformed orbitals which were retained in the energy calculation were all orthogonal and normalized to within an error of in their overlap matrix elements. The main difficulty in the selection of appropriate configurations for step b was, of course, the enormous number of possible configurations. I n a typical calculation there are billions of configurations which can be formed with the correct symmetry. For an unfortunate choice in step a all of these might be equally important, but for a good initial guess only a few will really contribute to the wave function. The selection of important configurations was done by scanning all single and double excitations of the dominant term for important contributions. Using Sinanoglu's suggestions,l0 the coefficients for triple and higher excitations were estimated from the coefficients of the singles and doubles. Then the 50 configurations which seemed most important were used in a variational calculation to get the wave function. ,As convergence was approached, successive iterations always had the first 35 to 45 configurations the same (except, of course, for the small changes in the orbitals themselves). Convergence in the expectation values of operators other than the energy followed the expected pattern. The one-electron operators r S - l , n,-l, 2,2*,and X 2 were given correctly by the dominant configuration and were relatively insensitive to correlation corrections to the wave function. Any error in the expectation value of these operators is due mainly to errors in the dominant configuration. The average value of the kinetic energy and electron repulsion are greatly affected by correlation corrections. I n general, the dominant configuration gave an average electron repulsion which was high by about twice the correlation energy and a kinetic energy which was low by about the correlation energy. Thus, in a very real sense, the effect of the correlation terms is to build up a partial wave expansion of the cusp in the wave function a t all points where two electrons of opposite spin are in contact. The first few terms in such an expansion represent long-range correlation (over dimensions comparable to the size of the orbital) and are usually called by such names as in-out, angular, left-right, etc. The higher terms are of shorter and shorter wave length and are exceedingly difficult to represent accurately with any of the usual basis sets. Probably no more than three or four orthogonal orbitals (13) W. Kutzelnigg, Theoret. Chim. Acta (Berlin), 1, 343 (1963). (14) G. P. Barnett, J. Linderberg, and H. Shull, J . Chem. Phys., 43, S80 (1965).
NATURAL ORBITALBASEDENERQY CALCULATION FOR HeH AND LiH
2677
Table I Part 1 Michels and Harrisa Basis for HeH at R = 3.0 3r s(lr)
ZC
n
j
m
a
B
0 0
0 0 0
3.227 1.836 1 ,447
-3.352 -1.669 1.590
0
ISH
- 0.4345 -0.7190 2.4498
- 15.2874
-2.0386 - 0.0289 1.476 2.610 0.825 1.9914
-1.609 3.614 1.925 1 ,0002
1.428 3.417 2.190 0.0084
- 3.7979
8.3862 0.8649
Part 2 Configuration
1u22u 3u22a1 lu23u1 lu'3u2 1 u'2a2 2u23a1 lu12u13u1 a
Dipole moment, C1
0.995597 0.003843 0.000002 0.000142 0.000001 0.000351 0,000065
E
D.
(Z9
- 3.34300 - 3.35676
0.399 He+H-
2.948
-3.35835
0.402 He+H-
2.947
(X*
+ Y*)
( T)
(Vne)
(Vee)
1.192
3.4178 3.4336
-9.0770 -9.0781
1.6496 1.6211
1.196
3.4352
-9.0778
1.6177
See ref 16.
of the same symmetry can be localized in a small region of space without complete loss of all significant figures in an energy calculation. It is precisely because the higher natural orbitals should be localized in the same region of space as the occupied orbitals that previous SCF-CI calculations have failed to represent the correlation correction to the ground state. The virtual SCF orbitals represent excited orbitals and are localized at great distances from the occupied orbitals. Thus, the virtual orbitals are essentially zero where the error in the SCF ground-state function is largest. Methods other than solving the SCF equations are available for constructing an initial set of orthonormal orbitals. Many of these were tried, but they generally led to unacceptably slow convergence in the CI calculation. The Lowdin transf~rmation'~ (to the orthogonal set with maximum similarity of the nonorthogonal basis) was particularly bad, probably because of the large number of very similar basis functions used in this calculation. The most satisfactory simple scheme was a judicious use of the Gram-Schmidt procedure. After a few trial calculations it became clear that the only really satisfactory solution was to solve the Roothaan SCF equations for the occupied orbitals. In the final calculation on LiH the initial guesses to the
correlation orbitals were chosen to diagonalize the SCF exchange integral operator and to be orthogonal to each other and the occupied orbitals. The correlation orbitals corresponding to very small exchange integrals with the charge density were indeed unimportant (as expected from second-order perturbation theory). The correlation orbitals corresponding to fairly large exchange integrals were important in about the same order as the square of their exchange integrals. Orbitals chosen by the above procedure changed very little during the natural orbital calculation. The SCF orbitals underwent a slight unitary transformation among themselves but remained essentially unchanged. Adjacent correlation orbitals in the ordered list mixed slightly (less than 10% in most cases). This procedure for the initial orthogonalization is the best of those considered and was satisfactory for the molecules treated here. The final answer is not too sensitive to the initial choice of nonorthogonal basis provided it is reasonable and fairly complete. We initially solved the sequence of two-electron atomic ions to get a ls2 core function expressed in elliptical basis functions. As split shells greatly increase the time required for a calculation, we (15) P. 0.Lowdin, J . Chem. Phya., 18, 365 (1950).
Volume 70, Number 8 Aupuet 1966
CIIARLES F. BENDER AND ERNEST R. DAVIDSON
2678
E
a
The J O U Tof~Phllgical Chemistry
oooooooooo
74-44
E
H N m N
NATURAL ORBITALBASEDENERGY CALCULATION FOR HeH AND LiH
2679
Volums 70,Number 8 Avgust 1066
2680
CHARLES F. BENDERAND ERNESTR. DAVIDSON
Table 111 Part 1 Natural Orbitals of LiH a t R = 3.015
n
j
m
U
0 0 1 1 0 2 0 0 1 1 2 0 1 1 1 1 0
0 1 0 1 2 0 0 1 0 1 0 2 0 1 0 1 1
0
4.709 4.709 4.709 4.709 4.709 4.709 1.800 1.800 1.800 1.800 1.800 1.800 1.280 1.280 1.000 1,000 1,000
0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
1.J
2.J
3r
B
ISL~
ISH
S(Z0)
4.709 4.709 4.709 4.709 4.709 4.709 -1.500 -1.500 -1.500 -1.500 -1.500 -1.500 0 0 1.000 1.ooo 1.000
11.9202 -0.7298 0.7799 19.6953 9.0733 10.6898 0.0640 0.0200 0.0194 0,0327 0.0083 0.0054 - 0.0037 - 0.0022 0.0016 -0.0099 0.0404
-0.8514 -2.1152 1.8818 3.0686 0.4259 -2.9503 2.9048 0.7987 0.2720 0.6050 0,4720 0.8299 0.5354 0,1734 0.0382 0.0178 -0.5135
-1.5122 2.4402 0.3023 1.9941
1.1819 3.0238 3.5924 1,9498
(Z) (Z2) (XZ n
+ Y?
n
j
m
0 0 1 1 0 2 0 0 1 1 2 0 1 1 1 1 0
0 1 0 1 2 0 0 1 0 1 0 2 0 1 0 1 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
4.709 4.709 4.709 4.709 4.709 4.709 -1.500 -1.500 -1.500 -1.500 -1.500 -1.500 0 0 1.000 1.000 1.000
n
142.014 190.320 107.783 118.743 -7.258 14.035 9.320 43.768 16.643 -35.009 11.866 -4.345 -1.490 0.303 2.536
-
+ Y*) j
m
a
B
0 0 1 0 1 0
0 1 0 0 0 0
1 1 1 1 1 1
6.757 6.757 6.757 0.980 0.980 2.200
6.757 6.757 6.757 0.980 0.980 -2.200
The Journal of Phyaical Chemistry
1s P(2.J)
0.5396 10.6597 11.9824 0.0927 -0,0191 - 5.5295
-
4.9936 0.6897 2.2090 - 4.9445 1.8640 3.1863 8.6171 3.8829 4.8138 - 2.0088 2.4372 0.4136 -0.2472 - 1.2481
21.7554 15.5027 - 47.3274 21.7173 18.2419 12,8923 -11.9401 0.1656 2.7419 -9.0184 -1.7652 0.9628 2.9039 0.8977 - 0.4524 0.4812 - 1.7429
-13.2911 - 70.5978 42.0115 - 116.2663 72.0269 66.7547 - 5.5772 - 4.8946 - 1.8042 -6.2069 - 1.1050 -3.0160 2.6755 0.7679 -0.6479 0,9836 - 2.6697
131.4339 78.3317 - 140.7811 - 101.5543 44,6850 75.9089 - 6.6466 -9.8221 1.8225 0,7942 -0.1481 3.8667 0,3806 0.2779 -0.1229 0.2401 -0.6716
0.8385 2.9641 5.1946 0.0227
1 3641 4.7634 1.8989 0.0063
- 1.2988
- 1.3774
3.0717 0.8892 0.0022
2.4284 0.3054 0.0010
6.5406
- 7.4827 - 8.5567
80
126.455 2.124 12.528 58.344 -20.716 65.154 29.552 17.669 16.587 0.192 -15.092 0.673 -0.137 5.413
- 161.725
5r s(lu)
I
-
-
0.7345 2.9586 1.8880 5 x 10-5 2s P(1.J)
13.4846
- 96.0328 158.9882 0.2358 -0.0467 5.3472
-
-
I
9u (1u)
1o.J (10)
-7.469
10.304 310.043 200.932 103.032 122.395 111.077 9.010 11.647 -22.355 - 65.133 -16.481 -37.334 22.663 2.770 -2.985 0.175 -0.207
(20)
- 44 037 - 78.739
168,268
1.9065 8.5021 3 ' 5754 0.0001
n
-
- 151.978
(2)
(Z? (XZ
-
7r (2.J)
B
U
4.709 4.709 4.709 4.709 4.709 4.709 1.800 1.800 1.800 1.800 1.800 1.800 1.280 1.280 1.000 1,000 1.000
-
40 P(20)
-93.663 67.619 -748.424 572.789 252,425 -27.269 -1.718 7.771 45,430 -21.433 -24.608 22.311 -7.211 -3.759 4.209 -6.405
-
-0.3915 3.0367 1.9890 2 x 10-6 3s p(2r)
-21.1594 20.1752 -0.3985 0.3480 - 0.0372 65.4621
-
-
-
110
12u (2U)
-905.775 905.645 228.501 408.548 - 478,582 -28.548 - 61.749 25.364 30.592 4.82% 24.425 12.913 4,505 2.967 -6.427 15.825
-
- 1.1462
- 1.5130
7.6955 8.3230 1 x 10-6
6.9277 3.8443 1 x 10-6
-61.7286 92,7669 -0.4317 0.8127 -0.1509 -9.1068
-
(1u)
- 174.440
4s d@o)
6u
P(1U)
179.435
- 142,797 -157.018 -637.936 465.543 325,104 -21.551 - 104.056 -4.647 39.171 18.462 - 49.338 -9.943 14.310 0.692 0.299 -7.294 0.6043 4.1163 3.4712 2 x 10-6 67 d(1a)
5r
P(l,J)
1326.183 -291.301 903.525 0.173 -0.038 25.558
-
387.527 -1171.037 501.531 1.318 -0.263 10.983
-
NATURAL ORBITALBASEDENERGY CALCULATION FOR HeH AND LiH
2681
Table I11 (Continued) Part 1 (Continued) a
2.200 2.200 2.000 2.000
1* P(W
27? P(10)
-4.5354 - 1.7087 - 2.3970 -0.5911
- 5.3495
1.0620 2.3608 3.9159 0.0207
- 1.5018
B
-2.200 -2.200 0 0
a
B
9.924 9.924 3.000 3.000 9.924
9.924 9.924 -3.000 -3.000 9.924
16 d(2u)
-1.4415 3.5530 - 9.9959 2.4812 0.4192 0.0022 26
d(1r)
- 189,266 215.296 19.521 0.747 0 1.4968 2.7494 3.0478 0.0002
1919.743
- 3733.547 -0.237 -0.370 0
-1.5814 2.5361 0.2206 3 x 10-6
4*
3r
5* P(1U)
Pm)
d(2d
54.8084 -7.5758 -35.1808 26.3310
- 38.4799 - 14.2725
1.8974 5.4285 3.4380 0.0002
1.3701 3.1893 2.9198 0.0002
-
15.2227 26.0206
24.328 3.527 -24.440 31.941
-438.615 488.081 47.760 83.660 0
11734.39
- 1.2782
2.8384 0.5360 0.0001
2.1436 0.9905 5 x 10-6
714.33
-,17776.17
-0.90 -1.66 0
-1.15 -1.86 14196.94
- 1.2613
- 1.6856
1.6792 0.4427
2.9367 0.4237 4 x 10-6
6
-4.909 5.468 -31.549
56 (1-J)
-,13541.94
0.9874 2.2336 3.3711 3 x 10-6
- 15.023
- 1.5776
46
38 d(2u)
6r d(1u)
x
Part 2 Energy of LiH a t R = 3.0150 Config
1u22a2 1u23u2
2u22n'iZr' 1U 2 1 7 r ' G ' 2~15~2 2~16~2 1u24u2 2U24d5U' 2cr26a16rr' 2~25~1%' 2 0226'5' 1U240'5U' la24a1G' 2U23U15U' lU216W
lu23x1~' 2u24u'6u' 2u24ua 1U'2U1502 2u29u2 ld2u'27r'2x' 1u'2u'4u16u' 1u27u2 3u22a1g1 2U W ' 4 6 ' lu12u~4u'5u' 2U2(lS'G' 4- i227r') 1u28u2 17r'G'27r'2?r' 2u256'56'
0
E
( T)
0.971786 0.011405 0.001015 0.010384 0.000844 0.000483 0.00289 1 0.000220 0.000024 0.000036 0.000014 0.0001 03 0.000087 0.000058 0.000106 0.000106 0.000032 0.000039 0.000023 0.000008 0.000019 0.000032 0.000049 0.000012 0.000003 0.000022 0.000035 0.000010 0.00001 1 0.000002
-7.98711 - 8.00227 -8.01614 - 8.02889 - 8.03911 - 8.04511 - 8.04937 - 8.05140 - 8.05232 - 8.05322 - 8.05392 - 8.05453 - 8.05503 - 8.05556 -8.05598 - 8.05633 - 8.05662 - 8.05687 - 8.05716 - 8.05751 - 8.05777 - 8.05796 - 8.05811 -8.05828 - 8.05842 - 8.05856 - 8.05889 - 8.05900 -8.05915 - 8.05927
7.9885 8.0197 8.0285 8.0331 8.0444 8.0475 8.0521 8.0551 8.0558 8.0563 8.0568 8.0556 8.0556 8.0563 8.0565 8.0565 8.0567 8.0570 8.0572 8.0574 8.0574 8.0576 8.0576 8.0581 8.0582 8.0584 8.0586 8.0586 8.0588 8.0589
(Vd
- 20.4717 - 20.4783 - 20.4734 - 20.4667 - 20.4681 - 20.4658 - 20.4683 - 20.4695 - 20.4695 - 20.4694 - 20.4694 - 20.4675 - 20.4675 - 20.4678 - 20.4678 - 20.4678 - 20.4678 - 20.4679 - 20.4678 - 20.4679 - 20.4676 - 20.4676 20.4676 - 20.4676 -20.4676 -20.4677 - 20.4676 -20.4676 - 20.4674 20.4674
-
-
(Vd
3.4911 3.4613 3.4337 3.4097 3.3896 3.3781 3,3718 3.3680 3.3664 3.3649 3.3636 3.3623 3.3618 3.3609 3.3604 3.3599 3.3595 3.3590 3.3584 3.3579 3.3573 3.3570 3.3568 3.3562 3.3560 3.3557 3.3552 3.3551 3.3545 3.3543
Volume 70,Number 8 August 1966
CHARLES F. BENDER AND ERNEST R. DAVIDSON
2682
Table I11 (Continued)
Part 2 (Continued) C* 0.000011 0.000010 0.000027 0.000015 0.000006 0.000009 0.000000 0.000013 0.000020 0.000006 0.000010 0.000005 0.000001 O.oooO06 0.000001
E
( T)
(vne)
( Vee)
- 8.05940
8.0593 8.0596 8.0594 8.0595 8.0596 8.0595 8.0595 8.0595 8.0589 8.0591 8.0592 8.0593 8.0593 8.0593 8.0590
- 20.4676 - 20.4675 -20.4670 -20.4670 - 20.4670 - 20.4668 -20.4668 -20.4668 -20.4659 - 20.4650 -20.4659 - 20.4658
3.3538 3.3534 3.3529 3.3527 3.3525 3.3523 3.3523 3.3522 3.3518 3.3515 3.3512 3.3510 3.3510 3.3508 3.3508
- 8.05952 - 8.05967 -8.05978
- 8.05989
- 8.06001 - 8.06001 -8,06008 - 8.06018 - 8.06026 - 8.06043 - 8.06050
-8.06052 -8.06059 -8,06062
tried to use only one effective charge for each symmetry type in the inner shell. With 6u-, 3a-,and 26type orbitals, we were able to calculate the energies of all ions for 2 from 1 to 11 within 0.0050 hartree. The only interesting facet of these results was that basis functions of the form r2 exp(-kr) are essential if the shell is not split. With this term the split-shell answer is obtained, while with only exp(-kr) and r exp(-kr) the answer is no better than with exp( - kr) alone. In carrying out the HeH and LiH calculations, we have relied on the numerous past calculations for the choice of valence shell orbital parameters. As we were primarily interested in developing techniques for CI calculations in terms of molecular orbitals, not much time was spent on optimizing parameters or variation of internuclear distance. Probably our basis sets for these calculations are needlessly large and could be simplified at no great loss in accuracy.
Results for HeH The ground-state energy of HeH was calculated at R 3.00, primarily as a check on our program. In the first calculation, the three a-basis functions reported by Michels and Harris16 were used. Since there were only seven possible spatial distributions (we define a spatial distribution as a particular selection of spacetype orbitals with all possible arrangements of spin), all configurations were included. The resulting wave function was then analyzed into natural orbitals and the energy recomputed. Table I shows the result of this calculation. I t will be noticed that two natural configurations are capable of reproducing the answer Michel and Harris obtained with one split-shell configuration. The full seven-term result is not appreci=
The Journal of Physical Chemistry
-20.4658 -20.4657 - 20.4654
ably better since the other configurations are not very appropriate. In all tables the basis functions are listed as (n,j , m, a,P ) = 2 ~ - ' / a2 - ' / 2 x n7)j(A2 - ,)IrnI/Z x (1 - T2)"n"2exp(-aolX ~7 im+)
+
+ +
where X = (r, T b ) R - ' , 7 = (r, - Tb)&', I$ is the cylindrical angle of rotation around the internuclear axis, r, and rb are distances from the two nuclei, and R is the internuclear separation. The coefficients for transforming the original basis to the natural basis are given in the first part of the table along with the orbital designation (la, 2a, . . .) and the orbital description s(1a) designates an s-type orbital used for local in-out correlation in the la shell, p (la) designates an orbital used for local left-right or angular correlation. In some cases the average values of 2, Z2,and X 2 Y 2are also given for each orbital (2 is along the internuclear axis and measured from the center of the molecule). The number n for each orbital is the occupation number normalized so that the sum of the occupation numbers is the number of electrons. In a very crude sense, n is the number of electrons "occupying" the orbital. The second part of the table contains the list of spatial distributions and the sum of the squares of all configurations involving that distribution when the full list is used in a variational calculation. I n the other columns the expectation values of other operators are given for a variational wave function using all configurations down to and including the one on the same line as the particular entry. The operators T,V,,, and V,, refer to the
+
(16) H. H. Michels and F. E. Harris, J. Chem.Phys., 39,1464 (1963).
NATURAL ORBITALBASEDENERGY CALCULATION FOR HeH
total kinetic energy, the total nuclear-electronic potential, and the total electron-electron repulsion potential, respectively. The other expectation values N
are defined by 0 =
N-'C 4-1
Oi except that
p
refers to
the total dipole moment in debye units. The repulsive effect in HeH for the Michels and Harris function comes primarily from the increase in kinetic energy of the 1s orbital on the H due to its orthogonality to the virtually unperturbed 1s orbital on He. The orthogonality requirement results in a net localization of the orbital away from the He nucleus. This can also be seen from the dipole moment which arises, not from charge transfer, but just from polarization of the 1s H atom orbital. Table I1 gives the result of a similar calculation using an extended basis set. In this calculation the GramSchmidt procedure was used for the initial orthogonalization and several iterations were required for convergence. The final result is by far the best result reported on this molecule. The energy (in hartrees) obtained (-3.38280) may be compared with results by Michels and Harris16 (- 3.35670), Taylor and Harris" (- 3.3618), or Mason1* (-3.3295). This calculation sets the upper bound on the repulsive energy ( E - E,) at R = 3 to be 0.0209 in comparison with the experimental result 0.0188 of Amdur and Mason.19 The formula of Amdur and Mason for the energy gives a slope, dE/dR, of -0.0207 a t this point, while the virial forE ) gives -0.0204. Scans mula dE/dR = -R-'(T of other configurations not included in the best 50 indicate that the energy limit of this basis set is about -3.3843. If the "calculated energy a t R = m " is used instead of the exact energy, Michels and Harris predict a repulsive energy of 0.01904, Taylor and Harris get 0.0302, and Mason gets 0.0182. The core basis used in this calculation gave an He atom energy in error by 0.0029 so the corrected repulsion from the present calculation is 0.0180. This latter argument is not completely valid, of course. Larger discrepancy results, for instance, if the basis limit energy is used to give a corrected repulsion of 0.0166. Part of the difficulty may be due to the fact that comparison of the one configuration energy a t R = 3 with the SCF results for He gives an estimated repulsion of only 0.0158. Thus, if the experimental result is correct, the correlation energy in the separated atoms is larger (by 0.003) than in the molecule. I t should be realized that the list of configurations contains many terms which are unimportant. The last 20 terms in this table contribute only 0.001 to the energy altogether and could be neglected for most purposes. These latter terms are certainly not neces-
+
AND
LiH
2683
sarily representative of the natural expansion of the true wave function. The first 11 terms were persistent and probably belong in the true function. Of these 11, only two do not arise from the product of an inner-shell geminal with an outer-shell orbital. These two, which had similar coefficients in Table I, can best be thought of as arising from a configuration of the type [la2 X(3a p 2 ~ > ~ ] (+rla)',where3a+p2aisnotorthog2u onal to 2a but is more like the free He atom correlation orbital. Similarly, 2a 71 a is more like a free H atom orbital than is 2a alone. Thus, these two intershell correlation terms arise from antisymmetrizing a simple product of He and H wave functions and then expressing the result in a set of orthogonal orbitals.
+
+
+
Results for LiH The results for our most extensive calculation on LiH a t R = 3.015 are presented in Table 111. The first two natural orbitals agree well with those given by Ebbing and Henderson." The energy of the first configuration is slightly lower than Ebbing's estimate2 of the SCF limit at R = 3.000. The first natural orbital is nearly the Li+ SCF orbital. The second resembles an H- SCF orbital polarized toward the Li nucleus and orthogonalized to the Li 1s orbital. The 5u, 6a, and 2n natural orbitals have occupation numbers in agreement with Kutzelnigg's results13 for Li+. These orbitals were localized near the Li nucleus and were of the shape expected for local in-out, leftright, and angular correlation, respectively. The energy contributions, as judged from configurations 5u22u2,6a22a2,and 2s22a2,were nearly the same as Kutzelnigg obtained with corresponding terms for Lif except that 5a22u2gave only 0.0102 to the energy and Kutzelnigg obtained 0.0137 from the first in-out correlation orbital for Lif. The energy contributions, from all configurations which appeared to be associated with inner-shell correlation, were 0.202 for 2-type terms, 0.0160 for II-type, and 0.0010 for A-type. The results of WeissZ0led to estimates of 0.0243 for 2, 0.0167 for II, and 0.0012 for A with the residual of 0.001 assigned to higher angular correlation. The energy of Li+ calculated with 6a, 3n, and 26 basis functions was in error by 0.0035. These results indicate that (a) the ls2 Li core in LiH is nearly the same as in Li+ and (b) the major deficiency in the present core function is in the inadequate description of in-out correlation by the 5a orbital. (17) H. S.Taylor and F. E. Harris, Mol. Phys., 7, 287 (1964). (18) E. A. Mason, J. Chem. Phys., 25, 626 (1956). (19) I. Amdur and E. A. Mason, ibid.,25, 630 (1956). (20) A. Weiss, Phya. Rev.. 122, 1826 (1961).
Volume 70,Number 8 August 1966
CHARLES F. BENDER AND ERNEST R. DAVIDSON
2684
A comparable analysis of Ebbing's wave function" shows that his in-out correlation orbital was better as it gave an energy lowering of 0.0128. RIost of the error in his core function seems to arise from lack of enough a basis functions and use of basis functions which were too diffuse for angular correlation. The outer-shell results are more difficult to understand. The total correlation energy of the molecule minus the known correlation energy of the Li+ ion indicates an outer-shell correlation energy equal to that of H-. However, the partitioning of this energy into a, r, and 6 parts does not follow the pattern of H-. Our calculation gave 0.0204 from a-type terms, 0.0138 from a-type, and 0.0005 from &type with about 0.005 more unaccounted for and presumed to be a. The results of WeisszOfor H- led to estimates of 0.0307 for a, 0.0084 for a, and 0.0003 for 6. However, Kutzelnigg13 obtained 0.092 from one ir-type natural orbital (assuming r is two-thirds of p) and we obtained 0.0127. Apparently the s function energy limit of H- is significantly lower than the s contribution to the s p limit. I n any case, the outer-shell angular correlation in LiH is much higher than previous estimates have indicated it should be. An analysis of Ebbing's outer-shell function shows that his in-out correlation orbital gave a smaller contribution to the energy than did ours and was somewhat different in appearance. Our 3a orbital resembles a 2s hydrogen atom orbital orthogonalized to the 1s Li orbital. It has a spherical nodal surface surrounding the H atom a t a radius of 1.9 bohrs. The convergence of the wave function is even slower than the results of Barnett, Linderberg, and Shull14 would indicate. There were numerous configurations with coefficients in excess of 0.001 (we have included 44 and omitted a t least that many more). The selection of Configurations to be included was necessarily arbitrary. It also was complicated by the fact that the occupation numbers of the correlation orbitals tend to be proportional to the reciprocal of the square of the effective charge while their energy contribution is relatively independent of charge. For example, a configuration like 5a22a2would be expected to have a coefficient about one-third as large as 1a23azbut would make almost the same energy contribution. I n order to pick terms which gave most rapid convergence of the energy, configurations were chosen on the basis of ( H n n - E)Cnz rather than sirnply Cn2. This choice deliberately omits contributions to the wave function which do not affect the energy and, hence, probably slows the convergence of the exceDted values of oDerators unrelated to the energy. It is possible that even more rapid energy convergence (and correspondingly slower convergence of
+
The Journal of Physical Chemistry
other properties) could be obtained if the matrix diagonalized to find the natural orbitals were chosen to be an energy density of some type, instead of the charge density. The one-electron properties of the molecule were given quite well by the SCF function. The one-configuration result for the dipole moment was 5.946 D., while with 45 spatial distributions, it was 5.965. The dipole moment seems to have converged to within kO.01 although the experimental dipole momentz1 (extrapolated to Re) is about 5.83. The nuclear-electron attraction, Vne,changes little in going from 1 to 45 configurations. It is doubtful, however, if the 45-configuration result is more than one significant figure better than the SCF result. The virial quantity, T E, varies from 0.0016 to 0.0001 with no apparent trend toward convergence. However, the ratio T / E is never wrong by more than 0.02'%. The energy itself converges more rapidly than it did in previous results by others. Thus, our seven-term result, -8.04937, is better than all previous calculations with the exception of Browne and i\latsenJ4 who obtained -8.0561. We required only 16 configurations to equal the result Browne and hLatsen obtained with 28 configurat,ions built directly from nonorthogonal basis functions. With a total of 45 spatial distributions we obtained -8.0606 out of an apparent energy
+
Table IV : Properties of LiH This calm
-E Re De, ev pet D. (Xi Y2)
+
(22)
(2) (TLi) 0-H) (TLi-') (TH-') (rLi-*P1(Cos h i ) ) (rH-*pl(cOs OH)) ( ~ L ~ - ~ P ~h (i c) )o s (rH-3p2(cos OH)) qLi qH
8.0606 3.0147 0.0825 5.9650 1.9563 2.7315 -0.1670 1.9573 2.5252 1.5187 0,5601 0.0287 - 0.0828 0.0137 0.0216 -0.0364 0.0460
Exptl
8.0705 3.0148" 0.0924b 5 . 83c
a G. Herzberg, "Spectra of Diatomic Molecules," 2nd ed, D. Van Nostrand Co., Inc., Princeton, N. J., 1950, Table 39. R. Valasco, Can. J . Phys., 35, 1024 (1957). See ref 21.
'
(21) L. Wharton, L. Gold, and W. Klemperer, J. Chem. Phys., 37, 2149 (1962).
2685
limit for our hasis set of -8.064. The energy of the first five configurations involving u-basis functions only was -8.02418 which was better than Ebbing's 41configuration u-only result." The energy limit from our u basis was about -8.032 and the true u limit would seem to be around -8.038. Acknowledgment. The authors gratefully acknowl-
edge the financial support of the National Science Foundation which made these calculations possible. We are also indebted to the University of Washington Research Computing Center for a generous grant of computing time on their IBM 7094 system and to Mr. S. Rothenberg for the use of several of his programs and many valuable discussions.
NOTES
Acetylene Production in the Radiolysis of Methane
by R. W. Hummel Wantage Research Laboratory ( A E R E ) , U.K . Atomic Energy Authority, Wantage. Berkshire, England (Received November 22, 1966)
In a recent note' in this journal, some experiments were described which were designed to measure any acetylene produced by the radiolysis of methane. A number of parameters were varied, but acetylene could not be detected. Emphasis was placed on the possibility of a dose rate dependence, but no acetylene was found even with 10-Mev protons "where the high LET has an effect similar to a very high local dose rate." Thus our previous report,2 giving initial G values for acetylene which decreased from 0.25 to 0.01 with increasing dose when using electrons from a linear accelerator at a pulse dose rate of about 9 X loz3 ev min-l g-l, could not be confirmed. The present note describes a few experiments designed to check our previous work. In the first experiment, purified CHI in a 200-ml Pyrex bulb at about atmospheric pressure and 20" was irradiated -6th @Coy rays a t 1.7 X lo1?ev min-I g-l; total dose 4.6 X 10% ev g-l. Acetylene could not be detected using a 6.5-ml sample of the irradiated gas injected into a 2-m column of 3% diethylhexyl sebacate on silica gel (Perkin-Elmer column S) at 40". A flame ionization detector was used, and the sensitivity was such that 1.4 X 10-4 mole yo of acetylene in a sample taken for analysis gave a peak area of 0.4 in.2, about six times greater than the minimum
detectable amount. Elution times in minutes on the column were as follows: C2H6, 2.7; C2H4,4.0; C3H8, 6.5; C2H2, 10.8; ~-CBHIO, 15; n-C4Hlo and CI&, 18. These analytical conditions were the same as those used in the previous work.2 At the dose given in this experiment, G(C2H2)= 0.03 would have been expected on the basis of the earlier work2 at the higher dose rate. The upper limit on G(C2H2) set by this experiment is 0.002. In the second experiment, 1% of propylene was added to the methane. Otherwise, conditions and methods were as in the first experiment. Samples were taken for analysis at three doses up to 2.3 X lom ev 8-1 and peaks were obtained at the elution time corresponding to that of acetylene. A good linear yielddose curve passing through the origin gave G(C2H2)= 0.045 f 0.005. This value is about one-third of that previously obtained2 in the presence of propylene using 4-Mev electrons at high pulse dose rates. In experiments 3 and 4,eight Pyrex bulbs of 2Wml capacity, containing pure CH, at about atmospheric pressure, were irradiated at 25" with 4-Mev electrons from a linear accelerator at a pulse dose rate of 3 X lOZ3 ev min-' g-I (3 X lozoev min-1 g-I during a duty cycle). Using liquid nitrogen, the hydrocarbon products were separated from methane and hydrogen, and the combined products from the eight bulbs passed through the gas chromatograph. The peak eluting at the retention time corresponding to acetylene was trapped at -196" and, after pumping away the helium carrier gas, analyzed on an MS-10 mass spectrometer (Associated Electrical Industries, Ltd., Manchester, England). After correcting for ad(1) R. H. Johnsen, J. Phys. Chem., 69, 3218 (19651. (2) R. W. Hummel, Discussions Faraday Sac., 36, 75 (1963).
Volume 70,Number 8 August 1966
