Volume 157.
number 4
CHEMICAL
PHYSICS LETTERS
KINETICALLY BALANCED GEOMETRIC FOR RELATIVISTIC MANY-ELEtX’RON
Ajaya K. MOHANTY
12 May 1989
GAUSSIAN BASIS SET CALCULATIONS ATOMS WITH FINITE NUCLEAR SIZE
and E. CLEMENTI
IBM Corporation, ScScrentific l+zgvwerrng Computations, Data Systems Division, Department 48B/MS 428. NetghborhoodRoad, Kmgston, NY 12401. USA Received
3 1January
1989: in final form 23 February
I989
We have performed DiraccFock atomic structure calculations on many-electron atoms using kinetically balanced, geometric Gaussian basis sets. Our results show that “variational collapse” or bound failure does not occur even at 2= 86. Furthermore, with a finite-size nuclear model, we obtain results close to the relativistic Hartree-Fock limit, with a smooth convergence from above. This opens up the possibility of optimisation of exponents in Gaussian basis set calculations within the Dirac-Fock approach for many-electron atoms.
1. Introduction Recently there has been a growing interest in solving the Dirac equation for many-electron atoms using basis set expansion methods. Following Roothaan’s [ I] proposal of basis set expansion in nonrelativistic quantum mechanics, Kim [ 21 developed a scheme for relativistic Hat-tree-Fock-Roothaan calculations for closed-shell atoms and Kagawa [ 3 ] extended these calculations to open-shell atoms. Despite carlier difticultics associated with the basis set methods #I, it is now generally agreed that, with proper boundary conditions imposed on the basis set, Dirac-Fock (DF) self-consistent field (SCF) calculations can be performed on many-electron atoms using techniques similar to non-relativistic methods. Crucial to the success of the basis set approach is the imposition of kinetic balance [ 5-7 ] between the large and the small component basis. Grant et al. [ 8-111 used kinetically balanced Slater-type orbitals (STOs) to perform calculations on H-, He-, Be-, and Ne-like atoms from low to moderately high atomic charge, and obtained results comparable in accuracy (one part in lo*) to numerical methods [12,13]. Recently the same authors have also reported [ 111 cal-
“’ For a review, see ref.
348
[ 4 1.
culations for mercury (Z= 80) and other high-Z atoms. Although STOs have the correct functional form to describe the relativistic wavefunction of atoms at the origin, they are not particularly suitable for analytic SCF molecular calculations. On the other hand, Gaussian-type orbitals (GTOs) are useful in the evaluation of multi-center integrals in molecules but they do not possess the correct functional behavior at the origin. However, with a finite nuclear size mode1 [ 141 the Gaussian basis sets provide a natural description of the relativistic wavefunction within the nucleus. Aerts et al. [ 15,16 ] used Gaussian basis sets to solve the DF equation for one- and manyelectron systems. Malli [ 171 proposed GTOs with both integer and non-integer powers of r for relativistic DF calculations of atomic properties for a wide range of atomic charges including xenon, in the pointnucleus approximation. Recently, Matsuoka et al. [ 18,191 reported kinetically balanced well-tempered basis set calculations and gave configuration average Dirac-Hartree-Fock (DHF) energies of various atoms, comparable to the numerical limit values of Desclaux [ 20 1. In this Letter we have used kinetically balanced geometric Gaussian basis sets with integer powers of r to calculate DHF total and orbital energies of atoms up to atomic charge 2~86.
0 009-2614/89/$ ( North-Holland
In
ref. [ 171 Malli
03.50 @ Elsevier Science Publishers Physics Publishing Division )
B.V.
Volume 157, number 4
CHEMICAL
used Gaussians with non-integer powers of r to calculate DHF total and orbital energies for Xe in the point-nucleus approximation, and the results compared favorably with the numerical values obtained by using the program of Desclaux [ 12 1. However, our results show that the use of geometric Gaussian basis sets with integer powers of r together with kinetic balance works well, particularly for a nucleus with finite extent. The same conclusion was drawn by Matsuoka et al. in ref. [ 19 1, although they used well-tempered Gaussian basis sets with kinetic balance to obtain configuration-averaged DHF energies up to xenon. We have obtained DHF results for relativistic closed- and open-shell configurations of atoms up to radon, and by making direct comparison with the numerical DHF results of Grant et al. [ 121, we show that we can achieve comparable accuracy. Furthermore, variational collapse does not occur even when one goes to very large exponents of the order of a few hundred million. With a finite nucleus size one converges faster to the true DHF value from above with fewer basis functions - a convergence pattern consistent with the results of Ishikawa et al.
1141.
PHYSICS LETTERS
number. P3/2,
d312,
I2 May 1989
For example 612,
f5,2,
bj7
the symmetries have
kappa
quantum
s~,~, P~,~, num-
bers - I, 1, -2, 2, - 3, 3, -4, respectively. The radial functions are expanded in terms of the basis sets as (3a) (3b) where the summation index i runs over the number of basis functions N, gki( r) and gpj( r) are Gaussian basis functions belonging to the large and the small component respectively, and CE and C”,; are the corresponding expansion coefficients. For the large component, the GTOs are of the form g~~(r)=Cf;v”‘exp(-cyir2),
(4)
where nK takes integer values 1, 2, 3, 4, for s~/~; fTjz symmetries respecPI/Z; ~3/2; G/2> d,,2; G/z, tively, and C’k is the normalisation factor. We have constrained our small component radial basis according to the kinetic balance condition [ 571
2. Choice of basis set In the central field approximation the solution to the Dirac-Fock Hamiltonian is given by the oneelectron orbitals [ 121 (atomic units (au) where e = m = fi = 1 are used throughout )
(1) Here, in referring to the individual orbitals, indices K and m are used to characterise the angular state. The quantum number K is given by Ic=-2(j-l)lj+1/2),
(21
where I is the orbital quantum number and j=1+ l/2 or I- l/2 is the total angular momentum quantum number. In the non-relativistic case the 1 quantum number is used to describe the symmetry of various orbitals. Similarly, symmetry species in the relativistic case are characterised by the K quantum
where Cc is the normalisation factor for the small component. Hence, the same set of exponents has been used for the large as well as the small component. The exponents ay, are generated according to ~i=cy,pf-‘,
i=l,
2, .... N,
(6)
where (Ye and /J are parameters determined by the basis set size N and the particular atomic system being considered. For xenon (Z= 54) N was taken to be of order 30. These geometric basis sets were originally introduced by Clementi et al. [ 2 I ] for non-relativistic atomic and molecular computations. In the present work we have employed similar geometric Gaussian basis sets with the kinetic balance and proper boundary conditions [ 8- 111 to calculate total DHF energy and orbital energies of atoms up to atomic charge 2~86. The derivation of the SCF equations for relativistic atoms in the closed-shell case are given in refs. [2,8]. For the open-shell case we use the same equation as in ref. [ 3 1. The expressions for all the required two-electron radial inte349
CHEMICAL
Volume 157, number 4
grals can easily be obtained in terms of incomplete beta functions [ 171. We have, however, adopted a previously developed algorithm [ 221. All calculations were performed in the finite nucleus model with a uniform charge distribution of radiusRz2.2677X IOW5A”‘, where A is the atomic mass number. The speed of light is taken to be 137.0373 au. We have also obtained a similar set of results using the point nucleus approximation.
3. Numerical results Calculations were performed on hydrogen-like systems up to 2=92, and for several open- and closed-shell atoms up to Z= 86 using kinetically balanced GTGs with integer powers of r. The results for a few representative atoms are presented in tables l4. The Breit interaction is not included in our calculation. For atoms up to Z=54 it was possible to find a
basis set which gave energies within a few millihartree of the numerical values. However, for radon (2~86) the DHF total energy obtained was within 0.029 hartree. The total DHF energy was always found to approach the true value from above. The relativistic virial theorem [ 23 ] CT>=-(0
He Li Be B C N 0 F Ne Na MS Al Si P S Cl Ar
(energies are in au)
State ‘r
-es1
-_t c)
Is2
2.86181335 7.4335332 14.5758919 24.5366169 37.6574196 54.3 169626 74.X393234 99.50230 128.69194
2.86 181328 7.4335329 14.5758909 24.5366132 37.6574128 54.3169528 74.8393043 99.50226 128.69 193 162.07809 199.93505 242.33112 289.44985 341.48891 398.60868 460.93977 528.68373
2&,2,,,2 2s2 2Pl,,,,,,,, 2p:,* 2P31:3,w* zp:,,,, 2p:.,,2,3,2 2P$a
381:1,2>1.2 3z.2 3Pl,,,,,,,, 3p:,, 3!&,2,3,L 3p::2,2 3P::,,2w 3P&*
162.07810 199.93508 242.33114 289.44987 34 1.48896 398.60875 460.93987 528.68384
(7)
was tested in the point nucleus case and was satisfied to at least one part in 10 ‘. In table 1 we present the atomic total DHF energies (in au) for He-Ar and compare them with the numerical DHF energies computed by using the program of Grant et al. [ 121. The results indicate that even with Gaussian basis sets, one can produce energies and orbital eigenvalues close to the numerical limits. In table 2 we give examples of orbital energies and the total DHF energy for xenon, which differ from the numerical values within a few millihartree. In table 3 we present our basis set for xenon and also show the convergence pattern for xenon in the
Table I Orbital and total DHF energies for He-Ar with the finite nucleus approximation Atom
12 May 19X9
PHYSICS LETTERS
-c d)
14.575878
128.69163 199.93463
528.68254
a) For open-shell cases, we adopt a notation similar to that of Grant et al. [ 121, i.e. nl;,,,,, where n denotes the principal quantum number. 1labels the orbital quantum number, 4 denotes the occupation number, v denotes the seniority number, j denotes the resultant quantum number of the open shell and x is the total quantum number of the state. The closed-shell structures are denoted by the last shell and its occupation number. b, Numerical DHF value [ 121 obtained with the average level (AL) option (Breit interaction not included). ‘) This work. ‘) Some selected closed-shell energies obtained by Matsuoka and Huzinaga basis sets.
350
[ 191by using non-relativistic optimised well-tempered Gaussian
Volume 157, number 4
CHEMICAL
Table 2 Orbital and total DHF energies of ground state Xe with the finite nucleus approximation (energies are in au)
Orbital lSl/Z 2s 112 3s 112 4s 1/J 5s l/2
2P,,, 3~~2 4P,,2 5P,(, 2PW 3P,/2
4PW 5P,,2 3d3/2 4d 311 3d,i, 4d5,z DHF total
-_Ea)
_tb)
1277.25880 202.46510 43.01043 8.42988
1277.25885
1.01012 189.67954 37.65989 6.45245 0.49256 177.70469
1.00999 189.67965 37.65992 6.45229 0.49240 177.70466
202.465 18
43.01049 8.42970
35.32536 5.98283 0.4398 1 26.023452 2.71114 25.53723 2.63382 -7446.9002
35.32568 5.98259 0.43966 26.023450 2.71115 25.53720 2.63359 -7446.8961
a) Same as footnote b of table 1. (The DHF energy given by Desclaux [20] is - 7446.9006 au.) b’ This work. (The DHF energy given by Matsuoka and Huzinaga [ 191 is -7446.8987 au.) Table 3 Convergence pattern on the DHF energy for the ground state of xenon in the tinitc nucleus size approximation. au
Energies are in
IvPI/Zb,
-E,,,”
25 28 28 28 28 30 31 31
20 22 24 27
7446.862525 7446.893170 7446.895598 7446.896125 7446.896 I34 7446.896 136 7446.896 136 7446.896 136
=) N I,,L is the number of geometric basis large component of s,,2 symmetry. b, NW is the number of geometric basis large component of pi,) symmetry. cl In this calculation we have kept the ND,,2= 20, A$,,, = 15, N,,,, = 15. The atedbyusingcu,=0.061248andP=2.13.
12 May 1989
Table 4 Orbital and total DHF energies ofground state Rn with the finite nucleus approximation (energies are in au) Orbital
-e a)
“_Eb)
lSI,, 2s I/Z 3s I/Z 4s,,* 5s I/2 6s l/2 2P,,, 3P,,, 4P,,, 5PI,, 6P,,, 2P3,, 3P,,, 4P,,, SP,,, 6~ J/Z 3d,,, 4d3,, 5d,,z 3ds,, 4d,,, 5ds,> 4fV2 4f7,* DHF total
3641.153 668.803 166.831 41.313 8.409 1.071 642.328 154.894 36.019 6.409 0.5404 541.103 131.731 30.121 5.176 0.3840 112.5699 2 I .548 2.190 107.760 20.439 2.017 9.194 8.929 -23601.974
3641.154 668.806 166.836 41.315 8.411 1.070 642.330 154.896 36.021 6.411 0.5407 541.105 131.734 30.123 5.178 0.3843 112.5673 21.550 2.191 107.762 20.441 2.018 9.197 8.931 -23601.945
a) Same as footnote b of table Desclaux [20] is -23601.978
I. (The DHF energy given by au.)
‘) This work.
N 6111”
28 28 29 30
PHYSICS LETTERS
functions
used for the
functions
used for the
following exponents
as constant: were gener-
finite nucleus case. Our results for xenon confirm the conclusion drawn by Ishikawa et al. [ 141 from their results for H-like mercury that convergence is much
faster
in the finite
nucleus
case than
the point-nu-
cleus calculation. From table 3 we see that by the time we have reached N,,,, =28 (number of basis functions for sl/2 symmetry) and N,,,, = 28 (number of basis functions for P,,~ symmetry), the DHF energy has converged to a limiting total energy of -7446.8961 au as against the numerical value of - 7446.9002 au. Even if we go beyond N,,,,, =28 and N zP,,z ~28 the DHF energy does not change to any appreciable degree. We consider that the kinetic balance condition given by (5), along with proper boundary conditions imposed on the basis set, contributes to the variational stability, as has been demonstrated by previous authors [ 8- 11,19 1. This type of convergence pattern seems to be a simple generalisation of the result of Ishikawa et al. [ 141 to the many-electron case. Table 4 shows the DHF energy and the orbital 351
Volume 157, number 4
CHEMICAL
PHYSICS LETTERS
energies for radon (the heaviest atom reported here). It can be seen that the orbital energies of both the inner and outer shells are obtained to an accuracy comparable with the numerical method; the maximum discrepancy is 0.005 au for the 3~,,~ orbital. We also predict the energy difference between the 4f,,, and 4f,,, orbitals to be within 0.001 au of the numerical DHF results. In passing, we should mention that the convergence pattern for radon is similar to that presented in table 3 for xenon in the finite nucleus approximation although we do not present the results explicitly. To summarise our results we can state that kinetic balance works well with GTOs, particularly for a nucleus with finite extent. The accuracy of our results up to and including radon is comparable to that of numerical methods and “variational collapse” does not occur even when one uses exponents as large as several hundred million. With a finite nucleus size there is faster convergence to the true DHF value from above with fewer exponents. This type of convergence pattern seems to be consistent with the results of Ishikawa et al. [ 141.
4. Conclusion In conclusion, we have demonstrated that with the use of kinetically balanced geometric Gaussian basis sets, Dirac-Fock SCF calculations in many-electron atoms can yield results comparable in accuracy to that of numerical DHF methods, particularly in the finite nucleus model. Our results show that “variational collapse” does not occur if we use a matched basis set satisfying kinetic balance and proper boundary conditions for bound state orbitals - a point which has been emphasised by Grant et al. [S111. Our calculations support the view that the finite-size nucleus approximation facilitates the use of GTOs in Dirac-Fock calculations for heavy atoms: the DHF total energy converges more rapidly to the true value from above and one avoids the necessity for a very large basis set, However, the results presented in this paper are by no means “final”. More work is being done to develop compact Gaussian basis sets within the Dirac-Fock approach [ 24 ] and to include 352
Breit interactions
in SCF calculations.
12 May 1989
Acknowledgement Interactions with Dr. 1.P. Grant, Professor G. Malli and Dr. Subhas Chakravorty regarding basis set calculations are gratefully acknowledged.
References [ I] C.C.J. Roothaan, Rev. Mod. Phys. 23 ( 1951) 69; 32 (1960) 179. Y.K. Kim, Phys. Rev. 154 (1967) 17; 159 (1967) 190. T. Kagawa, Phys. Rev. A 12 (1975) 2245;22 (1480) 2340. W. Kutzelnigg, Intern. J. Quantum Chem. 25 ( 1984) 107. Y. Ishikawa, R.C. Rinning and K.M. Sando, Chrm. Phys. Letters 101 ( 1983) 11 1. [6] K.G. Dyall, I.P. Grant and S. Wilson, J. Phys. B 17 (1984) 493,120l; L45. [7] R.E. Stanton and S. Havriliak, J. Chem. Phys. 81 (1984) 1910. [S] I.P. Grant, J. Phys. B 19 (1986) 3187. [9] I.P. Grant and H.M. Quiney, Advan. Atom. Mol. Phys. 23 (1988) 37. [IO] H.M. Quiney, I.P. Grant andS. Wilson, J. Phys. B 20 (1987) 1413. [ I1 P.M. Quiney, 1.P. Grant and S. Wilson, On the Relativistic Many-Body Perturbation Theory of Atomic and Molecular Electronic Structure, in: Lecture notes in chemistry, ed. U. Kaldor (Springer, Berlin), to be published. [ 121 I.P. Grant, B.J. McKenzie, P.H. Norrington, D.F. Mayers and N.C. Pyper, Comput. Phys. Commun. 2 1 ( 1980) 207. [ 13 ] J.P. Desclaux, Comput. Phys. Commun. 9 ( 1975) 3 1. [ 141 Y. Ishikawa, R. Baretty and R.C. Binning Jr., Chem. Phys. Letters 121 (1985) 130. [ 151 P.J.C. Aerts and W.C. Nieuwpoort, Chem. Phys. Letters 113 (1985) 16.5. [ 1610. Visser, P.J.C. Aerts, D. Hegarty and W.C. Nieuwpoort, Chem. Phys. Letters 134 ( 1987) 34. [ 171 G. Malli, Chem. Phys. Letters 68 (1979) 529. [IS] 0. Matsuoka, M. Klobukowskl and S. Huzinaga, Chem. Phys. Letters 113 (1985) 395. [I91 0. Matsuoka and S. Huzinaga, Chem. Phys. Letters 140 (1987) 567. [ 201 J.P. Desclaux, At. Nucl. Data Tables 12 ( 1973 ) 3 I 1. [21] E. Clementi and G. Corongiu. Chcm. Phys. Letters 90 (1982) 359; IBM Techn. Rept. POK-11 (1982). [ 22 ] B. ROOS, C. SaIez, A. Veillard and E. Clement], Techn. Rept. RJ518, IBM Research (August 12, 1968). [23] J. Maly and M. Houssonnois, Theoret. Chim. Acta 28 (1973) 363. [24] AK. Mohanty and E. Clementi, IBM Technical Report, Relativistic Atomic SCF Program, to be published. [2] [3] [4] [S]