- Email: [email protected]
Incipient manifestation of the shell structure of atoms within the WDA model for the exchange and kinetic energy density functionals
Chemical Physics ELSEVIER
Chemical Physics 196 (1995) 455-463
Incipient manifestation of the shell structure of atoms within the WDA model for the exchange and kinetic energy density functionals M.D. Glossman a, L.C. Balbfis b, J.A. Alonso b a Programa QUINOR, Facultad de Ciencias Exactas, Universidad Nacional de La Plata, Calle 47 esquina 115, Casilla de Correo 962, 1900 La Plata, Argentina b Departamento de Fisica Te6rica, Facuhad de Ciencias, Universidad de Valladolid, E 47011, Valladolid, Spain Received 6 December 1994; in final form 27 March 1995
Abstract
The radial electron density obtained for all the atoms of the main groups of the Periodic Table through the solution of the Euler equation associated with the nonlocal weighted density approximation (WDA) for the exchange and kinetic energy density functionals shows an incipient shell structure which is absent in other calculations using kinetic energy functionals based on the electronic density. The WDA radial density reveals two local maxima and the position of the first maximum correlates with the position of the maximum for the ls orbital in the Hartree-Fock approximation. The cusp condition at the nucleus is fulfilled accurately. Also we study the density-based electron localization function (DELF) as a complementary procedure for the visualization of shells.
1. Introduction One of the challenges in the field of electronic structure calculations is the development of interpretative procedures for the analysis of wavefunctions and charge densities in terms of well-known concepts such as atomic charges, valence, electronegativities of atoms and functional groups, the character of a bond, etc. One of the most significant attempts to meet this challenge involves partitioning the electron density of an atom [1-3] or molecule [4-7]. There has been a continuing interest in exploring possible relationships between the shell structure of atoms and their electronic density distribution [1,2,8-15]. In this respect, attention has focused upon the radial density function, D(r)= 4axrZp(r),
which goes through a series of maxima and minima with increasing radial distance from the nucleus [1,2,8-11,16]. Nevertheless, it has been found [9,11,17] that the radial densities obtained from Hartree-Fock (HF) [18] wavefunctions do not exhibit the correct shell structure for heavy atoms. This can be seen from the absence of certain minima or maxima at distances where they should formally be present [1,2,9,11]. On the other hand, a lot of effort has been put in studying energy functionals fully expressed in terms of the electron density [16,19,20]. The ThomasFermi model and its extensions [16,19,20] have been criticized on the ground that they lead to radial electronic densities where the expected shell structure is absent. There have been some attempts to
0301-0104/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0 3 0 1 - 0 1 04( 9 5 ) 0 0 1 2 9 - 8
456
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
build an energy density functional whose Euler equation resulting from the minimization of the energy with respect to the density should allow to obtain a radial electronic density showing shell structure. Of course, the shell structure of atoms has been introduced in those models using ad hoc procedures, like imposing a piecewise exponential decay on the density [10] or using variational densities constructed as a superposition of hydrogenlike one-electron wavefunctions [21-23]. However, the problem of finding such functionals presenting the proper shell behavior without introducing any constraint remains unsolved. We have recently studied for atoms and clusters of atoms [24-28] the solutions of the Euler equation associated with the weighted density approximation (WDA) model. The WDA was proposed several years ago by Gunnarsson et al. [29] and by Alonso and Girifalco [30] as an improvement over the local density approximation (LDA) for exchange and correlation effects, and it was extended to the kinetic energy by Alonso and Girifalco [30]. Besides the highly improved results of the physical and chemical properties of the studied systems in comparison with those obtained using gradient-corrected ThomasFermi models, it has been noticed [28] that the radial densities obtained through the WDA appear to show incipient indications of shell structure, which, as mentioned above, is absent in gradient expansion results, In this work, we explore in detail this particular feature of the WDA radial electron density of atoms and we find that there are two local maxima at positions close to those for the first two maxima of the HF density and that the relative magnitude of the two maxima is also in qualitative agreement with the HF trend. The position of the first maximum correlates with the maximum of the electronic density for the ls orbital in the HF approximation. Correspondingly, the electronic density can be partitioned into a corelike component describing a charge close to two electrons and a second component which is a relatively slowly varying function of the radial distance in a way resembling that proposed by Politzer and Parr [1] and by Tal and Bader [6]. Additional tests have been performed by studying the density-based electron localization function (DELF) formulated by Becke and Edgecombe [31] and its nonlocal WDA counterpart [32].
2. Weighted density approximation (WDA) A brief description of the WDA model for the kinetic and exchange-correlation energy density functionals is presented here. A full description of the model has been given in our previous paper [28]. The exact exchange-correlation energy [29] of a system of electrons can be written as (Hartree atomic units will be used through the paper) Exc[ p] =
½fp(r) d r f
Vxc( r ) =
8Exc [ P ] ~ p( r ~
p(r')G(r' r') dr'.
(1) I r - r'l This expression can be interpreted as the Coulomb interaction between the charge density p(r) and the exchange-correlation hole charge pxc(r, r') = p(r')G(r, r') which surrounds an electron at r. In this formula G(r, r') is the pair-correlation function of the system. The functional derivative of Exc[ p] with respect to the density is called the exchangecorrelation potential
(2)
The exact pair-correlation function fulfills the symmetry condition G(r, r') = G(r', r). Since the exact G(r, r') is generally unknown, it is necessary to resort to some approximation to it. In the local density approximation (LDA) the product p(r')G(r, r') is replaced by p(r)Gh[ I r -- r' I; p(r)], where Gh[ ] r-- r' I; p(r)] is the pair-correlation function in a homogeneous system with "constant" density equal to the "local" density p(r). If we restrict attention to the exchange-only case (that is, we neglect Coulomb correlation), then Gh[ [ r -- r' 1; p] is exactly known in analytical form 9 [ sin y - y cos y ] 3, Gh( Y; P) = - 2 / y3
)
r' y= Ir[(3'rr2p) 1/3, (3) which leads to the well-known Dirac expression for the LDA exchange energy [33]. On the other hand, in the nonlocal WDA, the correct factor p(r') is preserved in the exchange-correlation hole charge Pxc(r, r') = p(r')G(r, r') and G(r, r') is approximated by GWDA(r, r ' ) = G h [ [ r - - r ' l ; t3(r)], where G h is again the pair-correlation function of a homo-
M.D. Glossmanet al./ ChemicalPhysics196 (1995)455-463 geneous system, but now characterised by an effective "weighted" density /5(r), evaluated at each point r using the sum rule for the exchange hole
We then consider the following energy functional E[ p] = TWDA[ p] +EWDA[ p]
l f fp(r)p(r')
charge:
fp(r')Gh[ I r - r ' l ;
IS(r)] d r ' = - 1 .
+
(4)
becomes: P]
2 J p ( r ) dr x
[ p(r') G r ht lr-r'l; J
dr', (5)
drdr'
Ir-r'[
+fo(r)v,(r)
With the WDA ansatz for G, the exchange energy
E WDA[
457
dr,
(8)
where the last two terms represent the classical electron-electron interaction and the interaction of the electrons with the nucleus of the atom, represented by the potential Vi. The Euler equation for this energy functional becomes 1 V p . Vp P'= VT(r) + 8 p2
1 V2p +Ck [ 2/3 4 7 L
and the corresponding VxwoA is obtained from Eq. (2). A consequence is that the long range behaviour of VxwoA for a neutral atom is vxWDA~ --1/2r. This is an important improvement over the LDA exchange potential, which shows an exponential decay. One can also establish a WDA approximation for the kinetic energy. The key is a relation between the one and two particle density matrices in HartreeFock theory (see Ref. [30] for details). Although this relation is not universal, it is exact for some cases of practical interest, like atoms with closed electronic shells. Using that relation and the WDA approximation for the pair-correlation function G (of course at the exchange-only, or Hartree-Fock, level) one arrives at the following approximate expression for the kinetic energy [30]
TWDA[p]=ftwDAtr;p]dr,
(6)
1 (Vp) 2
1 2 ~V p,
-~/3 815(r')
(7)
8 p with C k = (3/10)(3'rr2) 2/3. The first term is a nonlocal extension of the local Thomas-Fermi term (C k p5/3) and the second is the original Weizs~icker quantum correction [34]. The integral of the last term does not contribute to the total energy in cases of practical interest.
]
+ 3
~p(r) dr'I,
(9)
where VT(r ) = Vi(r ) + Ve(r ) -4- Vx(r ) .
(10)
VT is the total potential, sum of the nuclear, electronic (Hartree)and exchange parts. Defining a function g ( r ) = 3 f ( r ) + ~2 p ( r ) - 2 / 3 /-
× jp(r')[9(r')-l/3_~/5(r') p(r) dr',
(11)
where f ( r ) = [ ~(r)/p(r)] 2/3, the Euler equation (9) takes the form 1 V p . Vp p2
/Z = V T ( r ) + 8
1 __VzP+ g ( r ) 4 p
× 5C 5 k pZ/3,
where the kinetic energy density is given by tWDA[r; p] = Ck p~2/3 +
2
(12)
that closely resembles the form of the Euler equation of the Thomas-Fermi-Weizs~icker functional [20, 35,36]. By putting ~b(r) = pl/2(r) and 5 Vaf(r) = VT(r) + g(r) × 3Ck p2/a(r), (13) Eq. (12) can be viewed as a Schr6dinger-like equation [ -2V 1 2 + Vaf(r)] ~O(r ) = pAb( r ) ,
(14)
458
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463 180
. . . .
i
. . . .
i
. . . .
i
. . . .
~ ,
160
160 i
i
.
.
.
.
=
.
.
.
.
i
.
.
.
.
i
'
'
140 Pb
140 ~.
Rn
120
120
~.
~ loo .--g" loo v)
"-¢::
80
Xe
60
~
n"
Ba 80
10
Kr
60
n,"
Rb
40
40
20
20
0
. . . .
oo
r
o.s
. . . .
i
. . . .
1.0 rl12(a.u.)
I
1.5
. . . .
i
'
Ca
0
2.0
oo
o.5
1.0 rla(a.u.)
1.s
Fig. 1. Radial density (4"rrr 2 p ( r ) ) of noble gas atoms obtained by the WDA model,
Fig. 2. Radial density (4w r 2 p ( r ) ) of Ca, Rb, Ba and Pb atoms obtained by the WDA model.
which we have solved by the conventional self-consistent method for only one fictitious s orbital normalized to Z electrons,
stead of the two maxima of the W D A density. The W D A density of Kr is further compared in Fig. 4 to that obtained from a K o h n - S h a m density functional calculation [16]. This K o h n - S h a m calculation evidently treats the single-particle kinetic energy exactly. On the other hand, exchange (and correlation) effects were included in the LDA. The K o h n - S h a m density of Fig. 4 and the HF density of Fig. 3 are very close, so the conclusions from Figs. 3 and 4 concerning the W D A density are the same.
3. Results and discussion The calculated W D A radial densities ( 4 7 r r 2 p ( r ) ) of the noble gas atoms have been plotted in Fig. 1 and those corresponding to Ca, Rb, Ba and Pb in Fig. 2. An analysis of similar plots for the radial densities of all the atoms in the main groups of the Periodic Table reveals that there are two local maxima at positions close to those for the first two maxima of the HF density. This is true for all but the lighter atoms, where there is a shoulder instead of a well developed second maximum. The evolution of the relative height of the two maxima of the W D A density is also in qualitative agreement with the HF trend. A direct comparison of densities is now presented. In Fig. 3 we represent the radial density of Kr obtained by the WDA, a gradient corrected (T O + T2) functional and the HF approaches. The gradient corrected functional includes the first two terms, TO + 7"2, in an expansion of the kinetic energy in the gradients of the density for a system of slowly vaying density, and it uses the L D A for exchange. It can be seen that the radial density from the gradient expansion functional shows only one maximum in-
6o . . . . . . . . . . . . . . . . . . . . . . . . . so
l ! l..'~"
m 40 = v~ ~ 30 '= "~ ~ 20 a: lo 0
•~1
, ....
w~A
I~
....
GRADIENT FUNCTIONAL
A :.:vi , ' ~/~ I ~. \ \ '...,~ "
.........
oo
0.4
, ....
~
, ..............
o.a
1.2 1.6 2.0 2.4 r1'2(a.u.) Fig. 3. Radial density ( 4 ~ r r Z p ( r ) ) of Kr, obtained from the WDA, gradient expansion (To + T2) and HF methods.
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
To interpret our results it is useful to recall first that the von Weizs~icker term alone, that is the second term on the right-hand side of Eq. (7), without any other contribution gives the exact kinetic energy of a two-electron system at the HF level. In addition, the value o f / 5 ( r ) for He, obtained from the sum rule (4), equals zero for all points r. This leads to G W O A ( r l 2 ) = -- 1 / 2 (see Eq. (3)), and then to 1
EWDA(He) = -- 2
(- ss
p(r)p(r')
dr
I t - r'l
dr'
)
5
, ....
, ....
,,
.'. !A:: WDA
4o
~30 -o 20 t~ ~r ~o
0
oa ~.2 ~B 20 rlJ2(a.u.) Fig. 4. Radial density (47rr2p(r)) of Kr, obtained by the WDA and Kohn-Sham methods, o.o
0.4
.4 ~"
(15)
This argument suggests that the von Weizs~cker term in the W D A functional describes the kinetic
, ....
o
'
focus on another consequence of the result /5(r) = 0 for He. The term in the kinetic energy containing ( / 9 ( r ) ) 2/3 is in this case zero, and T T M for He reduces to the pure yon Weizs~icker term; that is, T wDA is exact for He (or any two-electron ion) at the HF level. For an atom with more than two electrons, ~b(r) :# 0 in the region of the ls core shell, but we can still expect that the von Weizs~icker term dominates the kinetic energy density in this region,
.........
.6
E
which is the exact exchange energy of the He atom. In this particular case E w o A cancels the spurious self-interaction of each electron with itself contained in the classical Coulomb energy of Eq. (8). Although this is an important result that helps to explain the s u c c e s s of the W D A for exchange [37], w e n o w
so
459
0
.2
.~
rmax
.6
(W[)A)a u
Fig. 5. Position of the maximum of the radial density of the ls shell calculated by Desclaux [38] using the HF method, versus the position of the first maximum of the WDA radial density. Each circle represents a different atom. The data set includes atoms
from the main groups of the Periodic Table.
energy of the ls core shell, and that this term plus the nonlocal term f C k p~)2/3 d r account for the kinetic energy of the remaining electrons. According to this interpretation the first oscillation in the radial densities plotted in Figs. 1 and 2 (or in similar plots for other atoms) is to be associated with the ls shell. To substantiate this interpretation in Fig. 5 we represent for several atoms the position of the maximum of the radial density of the ls shell given by the HF calculations of Desclaus [38] against the position of the first maximum of the W D A radial density for the same atoms. The data set includes mainly atoms from the main groups of the Periodic Table, that is, excluding transition metal atoms. The points in this plot fit very well to a straight line rmax(lS)HF = rmaxWDA_ 0.0201 with an overall shift of only 0.0201 a.u. Other evidence for the ls orbital character of the first maximum of the W D A density is provided by ploting (r)ls(HF), the average value of r in the ls orbital obtained by HF calculations, against the position of the first maximum of the W D A radial density for the same atoms. For a hydrogenic ion the relation 3 ( r ) is = ~a, with a = a o / Z , is fulfilled, a 0 being the Bohr radius and Z the atomic number. On the other hand, the maximum in the radial density for the same
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
460 1.0
1.0 0
.8
.8
.6
t~ .6
^ .4 V
F .4
.2
2
0
0 .2
.4
.
.
.
.
.
.6
.
.
.
.
.2
rmax (WDA) au
.4
.6
rmax (HF) au
Fig. 6. Mean value of r in the ls shell, obtained from HF calculations [39], versus the position of the first maximum of the WDA radial density. The data set includes atoms from the main groups of the Periodic Table.
Fig. 7. Mean value of r in the ls shell versus the maximum of the radial density of that shell. All data are from HF calculations [38]. The data set includes atoms from the main groups of the Periodic Table.
hydrogenic ls orbital occurs at rmax(lS)= a, leading to the relation ( r ) i s = 3rmax(lS). In Fig. 6 we represent the values of ( r ) l s given by HF calculations [39] against the position of the first maximum of the W D A radial density for the same atoms. The points in this plot are fitted by a straight line ( r ) l s ( H F ) = WDA 1.57rma x -- 0.0282, SO the slope is very close to 3 / 2 as expected for hydrogenic ls orbitals. The small shift is practically the same as is the linear fit in Fig. 5. The deviation of the slope, 1.57, with respect to the value 1.5 is due to screening effects. The ls shell is not exactly hydrogenic due to the presence of the other electrons. This deviation is also observed in a plot similar to that of Fig. 6, but produced with pure HF data. Such a plot is given in Fig. 7. In this case the points are fitted by the relation ( r ) l s ( H F ) = HF 1.61rma x 0.00513. The deviation of the slope from 1.5 is very similar to that in Fig. 6. Evidently, the fit by a linear relation is better when all the data has been obtained from HF calculations. A plot similar to that of Fig. 6, but now using the K o h n - S h a m results for ( r ) l s , instead of the HF values, is given in Fig. 8. The relation that fits the points is
It has been shown some time ago [40] that the spherically averaged atomic charge density p(r) can be represented in terms of decreasing exponentials, and that there is a change in the rates of the exponentials that takes place at certain points associated with the minima in the radial density. Politzer and Parr [1] have considered the outermost minimum of the radial density to establish a boundary between the core and valence regions in atoms. It is obvious that if the
-
(r)ls(KS)
1,
.
.
.
.
.
.
.
.
10
O
09 05 .~ 07 06
-
=
1
WDA
.6393rmax
-- 0 . 0 3 2 4 9 , SO there is v e r y
little difference in using HF or KS data for the plot. A final comment concerning Figs. 6, 7 and 8 is that the slope of the fitted lines becomes closer to 1.5 if He is not included,
o5 ~ 0, o3 0~ 01 00 00
o o
01
0.2
o
03
04
05
06
07
r,u=,(WDA) (a.u.)
Fig. 8. Mean value of r in the ls shell, obtained from Kohn-Sham calculations, versus the position of the first maximum of the WDA
radial density. The data set includes atoms from the main groups of the Periodic Table.
M.D. Glossmanet al./ ChemicalPhysics 196 (1995) 455-463 electron density can be represented as a sum of decreasing exponentials, then in a logarithmic plot of the density against the radial distance r, the points where there is a change of slope will reflect the minima in the radial density. In Fig. 9 we have plotted the HF and W D A densities for xenon in a logarithmic scale against r in a short interval close to the nucleus (r~< 0.5 a.u.). The W D A density, represented by a dotted line, reflects reasonably well the behavior of the HF density, represented by a solid line. Of course, discrepancies exist because the HF density shows more structure than the W D A density. The change in slope occurring in both curves at r - - 0 . 0 5 a.u. reflects the first minimum in the radial electron density. The change is more pronounced for the HF case (the distinction between the lS shell and the rest is better defined). The HF density shows other changes of slope at r = 0.17 a.u. and r---0.45 a.u., which reflect the completion of Ne-like and Ar-like cores, respectively. Those changes are not reflected in the W D A density which is unable to reproduce the full shell structure of the atom. But the accuracy of the W D A density in a region of the atom close to the nucleus is undeniable. Let us now study in more detail the behavior of the density in this region. The slope of In p(r) for r very close
13 ......... , ......... , ................... , ........ 12 11 lo HF 9 WDA 8 ~ ~' -= 6 5 4 3 2 1 o0
ol
02
03
04
05
r (a.u.) Fig. 9. Electronic density of Xe in a logarithmic scale, obtained from the WDA and HF approaches. Only the region close to the nucleus is shown.
461
100 80
~ 60 -~ g3 40
20 . 20
.
.
.
40
60
80
10o
atomic number Z
Fig. 10. Test of the cusp condition at the nucleus, a is the exponentin a fit of the WDA density near the nucleus to a functionp(r)= poe -"r.
to the nucleus should be - 2 Z according to the Kato's cusp theorem [41]. In this region the main contribution to pWDA(r) can be described by an Is type orbital, as shown in the preceding section. By fitting the W D A density near the origin to an exponential p(r) = po e-~r, we have obtained values of a very close to 2Z, that is, the cusp condition is well satisfied. In Fig. 10 we plot the values of a / 2 against Z for all the atoms that we have studied. The points in this plot fit a straight line a / 2 = 0.955Z + 0.0809 showing the accurate fulfillment of the cusp theorem. Notice the alternance below and above the straight line in going from a main group to the next group of the Periodic Table. Recently, the electron localization function (ELF) proposed by Becke and Edgecombe [31] has been used for visualization of electronic shells. Within this context, Gadre et al. [32] have formulated a densitybased electron localization function (DELF) via the nonlocal weighted density approximation [30]. Thus, this formulation appears to be adequate to study the shell structure within the WDA. The density-based electron localization function (DELF) can be explicitly written [32] 1 DELF = (1 + ~ ( r ) / p (
r))2/3'
(16)
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
462
1.o .. ........ , ......... , ......... , ......... , ......... it:\
~:~ li~ \ -,: li ', \ 08 71~ '\ \ tu' !~ \\ \ o 07 I', \ \ ~\ ,"'. ~\~ . '~'~'~:~'-":~'-" ~ 0.9
,
\
0.6
-
-
Ne
--
Ar Kr - - - xe ......
~
,
~
~
~ ~ ~ "--~.,,,,~
~.-'~...~...~,
05 ............................................... 0.0 0.2 0.4 0.6 o.a 1.0 r (a.u.) Fig. 11. Density-based electron localization function (DELF) of Ne, Ar, Kr and Xe atoms, obtained by the W D A approach (see
Eq. (16)of the text),
where p(r) is the WDA electron density and iS(r) is the effective weighted density obtained through Eq. (4). In Fig. 11 we have represented DELF against r for the neon, argon, krypton and xenon atoms. The kind of information obtained from this plot is consistent with the information obtained from the radial density (Fig. 1) or from the plot of In p (Fig. 9). Notwithstanding, it is noticeable that the transition from one region of the atom to the other is sharper in the DELF representation and consequently the points where this transition takes place can be located in an easier way. To conclude this Section it is worth to analyze the reason why the separation between the ls core and the rest of the electrons is not achieved by the usual extended Thomas-Fermi functionals, either the Thomas-Fermi-von Weizs~icker functional (in which the lf[(Vp)2/p]dr term enters with its full ~1 coefficient) or the functionals arising from a series expansion in the gradients of the density (in those, the Weizs~icker term enters with a coefficient reduced by a factor ~). In all those functionals the ThomasFermi term enters as a fundamental one, providing a large contribution to the energy in the region of the ls shell. On the other hand, in the WDA, the fact that tS(r)<< p(r) in that region leads to a drastic reduction of the Thomas-Fermi-like term compared
to the von Weizs~icker contribution, which is itself large in this region. Results qualitatively similar to the ones in this paper have been obtained recently by Wang and Teter [42]. The work of these authors fits within the spirit of the WDA, as formulated for classical liquids [43,44], although the differences with respect to our work are substantial. Wang and Teter arrive to a kinetic energy functional that also contains the full Weizs~icker term, plus a nonlocal term from which a Thomas-Fermi-like contribution with a reduced coefficient is subtracted. Furthermore, the functional contains a free parameter which affects not only the numerical constants in front of some terms, but also the powers of the electron density. When applied to a very simple model of closed shell atoms (in which both exchange and the classical electron-electron Hartree t e r m w e r e neglected) shells qualitatively similar to those we have reported in Fig. 1 were obtained. It is gratifying that similar results have been found with two different kinetic energy functionals (and two very different treatments of the other energies). The agreement in the main results has to be ascribed to the basic WDA spirit used in both works.
4. Conclusions Several tests have been performed in order to study if the radial density obtained by solving the Euler equation associated with the WDA model for free atoms displays some kind of shell structure. A direct inspection of the radial electronic density shows that, except for light atoms, there are two maxima whose position and relative height vary following the same trend as the first two maxima of the HF radial density. For light atoms only one maximum and a pronounced shoulder at the position of the second are found. The same behavior has been observed when using other different criteria like plots of the density in a logarithmic scale and plots of the density-based electron localization function (DELF). Arguments have been given in favor of the interpretation that the first maximum is associated to the Is shell, and additional tests have been carried out relating the position of the first maximum in the
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
463
W D A radial density with the expectation value o f r
[15] M. Kakkar and K.D. Sen, Chem. Phys. Letters 226 (1994)
for I s orbitals in H F calculations. A l t h o u g h the W D A functional presented in this
241. [16] R.G. Parr and W. Yang, Density functional theory of atoms and molecules (Oxford Univ. Press, New York, 1989). [17] K.D. Sen, M. Slamet and V. Sahni, Chem. Phys. Letters 205 (1993) 313. [18] E. Clementi and C. Roetti, At. Data Nucl. Data Tables 14 (1974) 177. [19] E.S. Kryachko and E.V. Ludefia, Density functional theory of many-electron systems (Kluwer, Dordrecht, 1990). [20] R.M. Dreizler and E.K.U. Gross, Density functional theory. An approach to the quantum many-body problem (Springer, Berlin, 1990). [21] P. Csavinszky, Phys. Rev. A 24 (1981) 1215. [22] P. Csavinszky, Intern. J. Quantum Chem., Quantum Chem. Symp. 23 (1989) 591. [23] P. Csavinszky, Phys. Rev. A 39 (1989) 1519. [24] M.D. Glossman, A. Rubio, L.C. Balb~is and J.A. Alonso, Intern. J. Quantum Chem., Quantum Chem. Symp. 26 (1992) 347. [25] M.D. Glossman, A. Rubio, L.C. Balbfis and J.A. Alonso, New J. Chem. 16 (1992) 115. [26] M.D. Glossman, A. Rubio, LC. Balbfis, J.A. Alonso and LI. Serra, Intern. J. Quantum Chem. 45 (1993) 333. [27] M.D. Glossman, A. Rubio, L.C. Balbfis and J.A. Alonso, Phys. Rev. A 47 (1993) 1804. [28] M.D. Glossman, L.C. Balbfis and J.A. Alonso, Intern. J. Quantum Chem. 49 (1994) 171. [29] O. Gunnarsson, M. Jonson and B.I. Lundqvist, Phys. Rev. B 20 (1979) 3136. [30] J.A. Alonso and L.A. Girifalco, Phys. Rev. B 17 (1978) 3735. [31] A.D. Becke and K.E. Edgecombe, J. Chem. Phys. 92 (1990) 5397. [32] S.R. Gadre, S.A. Kulkarni and R.K. Pathak, J. Chem. Phys. 98 (1993) 3574. [33] P.A.M. Dirac, Proc. Cambridge Phil. Soc. 26 (1930) 376. [34] C.F. Weizs~icker, Z. Physik 96 (1935) 431. [35] W. Yang, Phys. Rev. A 34 (1986) 4575. [36] B.M. Deb and S.K. Ghosh, Intern. J. Quantum Chem. 23 (1983) 1. [37] O. Gritsenko, N. Cordero, A. Rubio, L.C. Balbfis and J. A. Alonso, Phys. Rev. A 48 (1993) 4197. [38] J.P. Desclaux, At. Data Nucl. Data Tables 26 (1973) 311. [39] C. Froese-Fischer, The Hartree-Fock method for atoms (Wiley, New York, 1977); Computer Phys. Commun. 14 (1978) 145. [40] H. Weinstein, P. Politzer and S. Srebenik, Theoret. Chim. Acta 38 (1975) 159. [41] T. Kato, Commun. Pure Appl. Math. 10 (1957) 151. [42] L.W. Wang and M.P. Teter, Phys. Rev. B 45 (1992) 13196. [43] P. Tarazona, Phys. Rev. A 31 (1985) 2672. [44] W.A. Curtin and N.W. Ashcroft, Phys. Rev. A 32 (1985) 2909.
paper is only able to separate the l s shell f r o m the rest o f the electrons of the atom, our results give s o m e c o n f i d e n c e that i m p r o v e d W D A - t y p e functionals, containing the v o n Weizs~icker term as a fundamental ingredient m a y be able to reproduce the full shell structure o f the atoms in the future. If this goal is achieved, it will constitute a significant step in density functional theory.
Acknowledgment
This w o r k has b e e n supported by D G I C Y T (Grant P B 9 2 - 0 6 4 5 - C 0 3 - 0 1 ) and Junta de Castilla y Le6n. M.D.G. is m e m b e r of the R e s e a r c h Career o f the C o n s e j o N a c i o n a l de Investigaciones Cientlficas y T6cnicas ( C O N I C E T ) , R e p u b l i c o f Argentina, f r o m w h i c h support is gratefully a c k n o w l e d g e d .
References [1] P. Politzer and R.G. Parr, J. Chem. Phys. 64 (1976) 4634. [2] R.J. Boyd, J. Chem. Phys. 66 (1977) 356. [3] R.F.W. Bader and M.E. Stephens, J. Am. Chem. Soc. 97 (1975) 7391. [4] R.F.W. Bader, P.M. Beddall and P.E. Cade, J. Am. Chem. Soc. 93 (1971) 3095. [5] R.F.W. Bader, Accounts Chem. Res. 8 (1975) 34. [6] Y. Tal and R.F.W. Bader, Intern. J. Quantum Chem. S12 (1978) 153. [7] P. Politzer, J. Reuther and G.T. Kasten, J. Chem. Phys. 67 (1977) 2385. [8] J.T. Waber and D.T. Cromer, J. Chem. Phys. 42 (1965) 4116. [9] R.J. Boyd, J. Phys. B 9 (1976) L69. [10] W.P. Wang and R.G. Parr, Phys. Rev. A 16 (1977) 891. [11] A.M. Simas, R.P. Sagar, A.C.T. Ku and V.H. Smith Jr., Can. J. Chem. 66 (1988) 1923. [12] R.P. Sagar, A.C.T. Ku, V.H. Smith Jr. and A.M. Simas, J. Chem. Phys. 88 (1988) 4367. [13] Z. Shi and R.J. Boyd, J. Chem. Phys. 88 (1988) 4375. [14] P. Politzer, J.S. Murray, M.E. Grice, T. Brinck and S. Ranganathan, J. Chem. Phys. 95 (1991) 6699.
Chemical Physics 196 (1995) 455-463
Incipient manifestation of the shell structure of atoms within the WDA model for the exchange and kinetic energy density functionals M.D. Glossman a, L.C. Balbfis b, J.A. Alonso b a Programa QUINOR, Facultad de Ciencias Exactas, Universidad Nacional de La Plata, Calle 47 esquina 115, Casilla de Correo 962, 1900 La Plata, Argentina b Departamento de Fisica Te6rica, Facuhad de Ciencias, Universidad de Valladolid, E 47011, Valladolid, Spain Received 6 December 1994; in final form 27 March 1995
Abstract
The radial electron density obtained for all the atoms of the main groups of the Periodic Table through the solution of the Euler equation associated with the nonlocal weighted density approximation (WDA) for the exchange and kinetic energy density functionals shows an incipient shell structure which is absent in other calculations using kinetic energy functionals based on the electronic density. The WDA radial density reveals two local maxima and the position of the first maximum correlates with the position of the maximum for the ls orbital in the Hartree-Fock approximation. The cusp condition at the nucleus is fulfilled accurately. Also we study the density-based electron localization function (DELF) as a complementary procedure for the visualization of shells.
1. Introduction One of the challenges in the field of electronic structure calculations is the development of interpretative procedures for the analysis of wavefunctions and charge densities in terms of well-known concepts such as atomic charges, valence, electronegativities of atoms and functional groups, the character of a bond, etc. One of the most significant attempts to meet this challenge involves partitioning the electron density of an atom [1-3] or molecule [4-7]. There has been a continuing interest in exploring possible relationships between the shell structure of atoms and their electronic density distribution [1,2,8-15]. In this respect, attention has focused upon the radial density function, D(r)= 4axrZp(r),
which goes through a series of maxima and minima with increasing radial distance from the nucleus [1,2,8-11,16]. Nevertheless, it has been found [9,11,17] that the radial densities obtained from Hartree-Fock (HF) [18] wavefunctions do not exhibit the correct shell structure for heavy atoms. This can be seen from the absence of certain minima or maxima at distances where they should formally be present [1,2,9,11]. On the other hand, a lot of effort has been put in studying energy functionals fully expressed in terms of the electron density [16,19,20]. The ThomasFermi model and its extensions [16,19,20] have been criticized on the ground that they lead to radial electronic densities where the expected shell structure is absent. There have been some attempts to
0301-0104/95/$09.50 © 1995 Elsevier Science B.V. All rights reserved SSDI 0 3 0 1 - 0 1 04( 9 5 ) 0 0 1 2 9 - 8
456
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
build an energy density functional whose Euler equation resulting from the minimization of the energy with respect to the density should allow to obtain a radial electronic density showing shell structure. Of course, the shell structure of atoms has been introduced in those models using ad hoc procedures, like imposing a piecewise exponential decay on the density [10] or using variational densities constructed as a superposition of hydrogenlike one-electron wavefunctions [21-23]. However, the problem of finding such functionals presenting the proper shell behavior without introducing any constraint remains unsolved. We have recently studied for atoms and clusters of atoms [24-28] the solutions of the Euler equation associated with the weighted density approximation (WDA) model. The WDA was proposed several years ago by Gunnarsson et al. [29] and by Alonso and Girifalco [30] as an improvement over the local density approximation (LDA) for exchange and correlation effects, and it was extended to the kinetic energy by Alonso and Girifalco [30]. Besides the highly improved results of the physical and chemical properties of the studied systems in comparison with those obtained using gradient-corrected ThomasFermi models, it has been noticed [28] that the radial densities obtained through the WDA appear to show incipient indications of shell structure, which, as mentioned above, is absent in gradient expansion results, In this work, we explore in detail this particular feature of the WDA radial electron density of atoms and we find that there are two local maxima at positions close to those for the first two maxima of the HF density and that the relative magnitude of the two maxima is also in qualitative agreement with the HF trend. The position of the first maximum correlates with the maximum of the electronic density for the ls orbital in the HF approximation. Correspondingly, the electronic density can be partitioned into a corelike component describing a charge close to two electrons and a second component which is a relatively slowly varying function of the radial distance in a way resembling that proposed by Politzer and Parr [1] and by Tal and Bader [6]. Additional tests have been performed by studying the density-based electron localization function (DELF) formulated by Becke and Edgecombe [31] and its nonlocal WDA counterpart [32].
2. Weighted density approximation (WDA) A brief description of the WDA model for the kinetic and exchange-correlation energy density functionals is presented here. A full description of the model has been given in our previous paper [28]. The exact exchange-correlation energy [29] of a system of electrons can be written as (Hartree atomic units will be used through the paper) Exc[ p] =
½fp(r) d r f
Vxc( r ) =
8Exc [ P ] ~ p( r ~
p(r')G(r' r') dr'.
(1) I r - r'l This expression can be interpreted as the Coulomb interaction between the charge density p(r) and the exchange-correlation hole charge pxc(r, r') = p(r')G(r, r') which surrounds an electron at r. In this formula G(r, r') is the pair-correlation function of the system. The functional derivative of Exc[ p] with respect to the density is called the exchangecorrelation potential
(2)
The exact pair-correlation function fulfills the symmetry condition G(r, r') = G(r', r). Since the exact G(r, r') is generally unknown, it is necessary to resort to some approximation to it. In the local density approximation (LDA) the product p(r')G(r, r') is replaced by p(r)Gh[ I r -- r' I; p(r)], where Gh[ ] r-- r' I; p(r)] is the pair-correlation function in a homogeneous system with "constant" density equal to the "local" density p(r). If we restrict attention to the exchange-only case (that is, we neglect Coulomb correlation), then Gh[ [ r -- r' 1; p] is exactly known in analytical form 9 [ sin y - y cos y ] 3, Gh( Y; P) = - 2 / y3
)
r' y= Ir[(3'rr2p) 1/3, (3) which leads to the well-known Dirac expression for the LDA exchange energy [33]. On the other hand, in the nonlocal WDA, the correct factor p(r') is preserved in the exchange-correlation hole charge Pxc(r, r') = p(r')G(r, r') and G(r, r') is approximated by GWDA(r, r ' ) = G h [ [ r - - r ' l ; t3(r)], where G h is again the pair-correlation function of a homo-
M.D. Glossmanet al./ ChemicalPhysics196 (1995)455-463 geneous system, but now characterised by an effective "weighted" density /5(r), evaluated at each point r using the sum rule for the exchange hole
We then consider the following energy functional E[ p] = TWDA[ p] +EWDA[ p]
l f fp(r)p(r')
charge:
fp(r')Gh[ I r - r ' l ;
IS(r)] d r ' = - 1 .
+
(4)
becomes: P]
2 J p ( r ) dr x
[ p(r') G r ht lr-r'l; J
dr', (5)
drdr'
Ir-r'[
+fo(r)v,(r)
With the WDA ansatz for G, the exchange energy
E WDA[
457
dr,
(8)
where the last two terms represent the classical electron-electron interaction and the interaction of the electrons with the nucleus of the atom, represented by the potential Vi. The Euler equation for this energy functional becomes 1 V p . Vp P'= VT(r) + 8 p2
1 V2p +Ck [ 2/3 4 7 L
and the corresponding VxwoA is obtained from Eq. (2). A consequence is that the long range behaviour of VxwoA for a neutral atom is vxWDA~ --1/2r. This is an important improvement over the LDA exchange potential, which shows an exponential decay. One can also establish a WDA approximation for the kinetic energy. The key is a relation between the one and two particle density matrices in HartreeFock theory (see Ref. [30] for details). Although this relation is not universal, it is exact for some cases of practical interest, like atoms with closed electronic shells. Using that relation and the WDA approximation for the pair-correlation function G (of course at the exchange-only, or Hartree-Fock, level) one arrives at the following approximate expression for the kinetic energy [30]
TWDA[p]=ftwDAtr;p]dr,
(6)
1 (Vp) 2
1 2 ~V p,
-~/3 815(r')
(7)
8 p with C k = (3/10)(3'rr2) 2/3. The first term is a nonlocal extension of the local Thomas-Fermi term (C k p5/3) and the second is the original Weizs~icker quantum correction [34]. The integral of the last term does not contribute to the total energy in cases of practical interest.
]
+ 3
~p(r) dr'I,
(9)
where VT(r ) = Vi(r ) + Ve(r ) -4- Vx(r ) .
(10)
VT is the total potential, sum of the nuclear, electronic (Hartree)and exchange parts. Defining a function g ( r ) = 3 f ( r ) + ~2 p ( r ) - 2 / 3 /-
× jp(r')[9(r')-l/3_~/5(r') p(r) dr',
(11)
where f ( r ) = [ ~(r)/p(r)] 2/3, the Euler equation (9) takes the form 1 V p . Vp p2
/Z = V T ( r ) + 8
1 __VzP+ g ( r ) 4 p
× 5C 5 k pZ/3,
where the kinetic energy density is given by tWDA[r; p] = Ck p~2/3 +
2
(12)
that closely resembles the form of the Euler equation of the Thomas-Fermi-Weizs~icker functional [20, 35,36]. By putting ~b(r) = pl/2(r) and 5 Vaf(r) = VT(r) + g(r) × 3Ck p2/a(r), (13) Eq. (12) can be viewed as a Schr6dinger-like equation [ -2V 1 2 + Vaf(r)] ~O(r ) = pAb( r ) ,
(14)
458
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463 180
. . . .
i
. . . .
i
. . . .
i
. . . .
~ ,
160
160 i
i
.
.
.
.
=
.
.
.
.
i
.
.
.
.
i
'
'
140 Pb
140 ~.
Rn
120
120
~.
~ loo .--g" loo v)
"-¢::
80
Xe
60
~
n"
Ba 80
10
Kr
60
n,"
Rb
40
40
20
20
0
. . . .
oo
r
o.s
. . . .
i
. . . .
1.0 rl12(a.u.)
I
1.5
. . . .
i
'
Ca
0
2.0
oo
o.5
1.0 rla(a.u.)
1.s
Fig. 1. Radial density (4"rrr 2 p ( r ) ) of noble gas atoms obtained by the WDA model,
Fig. 2. Radial density (4w r 2 p ( r ) ) of Ca, Rb, Ba and Pb atoms obtained by the WDA model.
which we have solved by the conventional self-consistent method for only one fictitious s orbital normalized to Z electrons,
stead of the two maxima of the W D A density. The W D A density of Kr is further compared in Fig. 4 to that obtained from a K o h n - S h a m density functional calculation [16]. This K o h n - S h a m calculation evidently treats the single-particle kinetic energy exactly. On the other hand, exchange (and correlation) effects were included in the LDA. The K o h n - S h a m density of Fig. 4 and the HF density of Fig. 3 are very close, so the conclusions from Figs. 3 and 4 concerning the W D A density are the same.
3. Results and discussion The calculated W D A radial densities ( 4 7 r r 2 p ( r ) ) of the noble gas atoms have been plotted in Fig. 1 and those corresponding to Ca, Rb, Ba and Pb in Fig. 2. An analysis of similar plots for the radial densities of all the atoms in the main groups of the Periodic Table reveals that there are two local maxima at positions close to those for the first two maxima of the HF density. This is true for all but the lighter atoms, where there is a shoulder instead of a well developed second maximum. The evolution of the relative height of the two maxima of the W D A density is also in qualitative agreement with the HF trend. A direct comparison of densities is now presented. In Fig. 3 we represent the radial density of Kr obtained by the WDA, a gradient corrected (T O + T2) functional and the HF approaches. The gradient corrected functional includes the first two terms, TO + 7"2, in an expansion of the kinetic energy in the gradients of the density for a system of slowly vaying density, and it uses the L D A for exchange. It can be seen that the radial density from the gradient expansion functional shows only one maximum in-
6o . . . . . . . . . . . . . . . . . . . . . . . . . so
l ! l..'~"
m 40 = v~ ~ 30 '= "~ ~ 20 a: lo 0
•~1
, ....
w~A
I~
....
GRADIENT FUNCTIONAL
A :.:vi , ' ~/~ I ~. \ \ '...,~ "
.........
oo
0.4
, ....
~
, ..............
o.a
1.2 1.6 2.0 2.4 r1'2(a.u.) Fig. 3. Radial density ( 4 ~ r r Z p ( r ) ) of Kr, obtained from the WDA, gradient expansion (To + T2) and HF methods.
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
To interpret our results it is useful to recall first that the von Weizs~icker term alone, that is the second term on the right-hand side of Eq. (7), without any other contribution gives the exact kinetic energy of a two-electron system at the HF level. In addition, the value o f / 5 ( r ) for He, obtained from the sum rule (4), equals zero for all points r. This leads to G W O A ( r l 2 ) = -- 1 / 2 (see Eq. (3)), and then to 1
EWDA(He) = -- 2
(- ss
p(r)p(r')
dr
I t - r'l
dr'
)
5
, ....
, ....
,,
.'. !A:: WDA
4o
~30 -o 20 t~ ~r ~o
0
oa ~.2 ~B 20 rlJ2(a.u.) Fig. 4. Radial density (47rr2p(r)) of Kr, obtained by the WDA and Kohn-Sham methods, o.o
0.4
.4 ~"
(15)
This argument suggests that the von Weizs~cker term in the W D A functional describes the kinetic
, ....
o
'
focus on another consequence of the result /5(r) = 0 for He. The term in the kinetic energy containing ( / 9 ( r ) ) 2/3 is in this case zero, and T T M for He reduces to the pure yon Weizs~icker term; that is, T wDA is exact for He (or any two-electron ion) at the HF level. For an atom with more than two electrons, ~b(r) :# 0 in the region of the ls core shell, but we can still expect that the von Weizs~icker term dominates the kinetic energy density in this region,
.........
.6
E
which is the exact exchange energy of the He atom. In this particular case E w o A cancels the spurious self-interaction of each electron with itself contained in the classical Coulomb energy of Eq. (8). Although this is an important result that helps to explain the s u c c e s s of the W D A for exchange [37], w e n o w
so
459
0
.2
.~
rmax
.6
(W[)A)a u
Fig. 5. Position of the maximum of the radial density of the ls shell calculated by Desclaux [38] using the HF method, versus the position of the first maximum of the WDA radial density. Each circle represents a different atom. The data set includes atoms
from the main groups of the Periodic Table.
energy of the ls core shell, and that this term plus the nonlocal term f C k p~)2/3 d r account for the kinetic energy of the remaining electrons. According to this interpretation the first oscillation in the radial densities plotted in Figs. 1 and 2 (or in similar plots for other atoms) is to be associated with the ls shell. To substantiate this interpretation in Fig. 5 we represent for several atoms the position of the maximum of the radial density of the ls shell given by the HF calculations of Desclaus [38] against the position of the first maximum of the W D A radial density for the same atoms. The data set includes mainly atoms from the main groups of the Periodic Table, that is, excluding transition metal atoms. The points in this plot fit very well to a straight line rmax(lS)HF = rmaxWDA_ 0.0201 with an overall shift of only 0.0201 a.u. Other evidence for the ls orbital character of the first maximum of the W D A density is provided by ploting (r)ls(HF), the average value of r in the ls orbital obtained by HF calculations, against the position of the first maximum of the W D A radial density for the same atoms. For a hydrogenic ion the relation 3 ( r ) is = ~a, with a = a o / Z , is fulfilled, a 0 being the Bohr radius and Z the atomic number. On the other hand, the maximum in the radial density for the same
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
460 1.0
1.0 0
.8
.8
.6
t~ .6
^ .4 V
F .4
.2
2
0
0 .2
.4
.
.
.
.
.
.6
.
.
.
.
.2
rmax (WDA) au
.4
.6
rmax (HF) au
Fig. 6. Mean value of r in the ls shell, obtained from HF calculations [39], versus the position of the first maximum of the WDA radial density. The data set includes atoms from the main groups of the Periodic Table.
Fig. 7. Mean value of r in the ls shell versus the maximum of the radial density of that shell. All data are from HF calculations [38]. The data set includes atoms from the main groups of the Periodic Table.
hydrogenic ls orbital occurs at rmax(lS)= a, leading to the relation ( r ) i s = 3rmax(lS). In Fig. 6 we represent the values of ( r ) l s given by HF calculations [39] against the position of the first maximum of the W D A radial density for the same atoms. The points in this plot are fitted by a straight line ( r ) l s ( H F ) = WDA 1.57rma x -- 0.0282, SO the slope is very close to 3 / 2 as expected for hydrogenic ls orbitals. The small shift is practically the same as is the linear fit in Fig. 5. The deviation of the slope, 1.57, with respect to the value 1.5 is due to screening effects. The ls shell is not exactly hydrogenic due to the presence of the other electrons. This deviation is also observed in a plot similar to that of Fig. 6, but produced with pure HF data. Such a plot is given in Fig. 7. In this case the points are fitted by the relation ( r ) l s ( H F ) = HF 1.61rma x 0.00513. The deviation of the slope from 1.5 is very similar to that in Fig. 6. Evidently, the fit by a linear relation is better when all the data has been obtained from HF calculations. A plot similar to that of Fig. 6, but now using the K o h n - S h a m results for ( r ) l s , instead of the HF values, is given in Fig. 8. The relation that fits the points is
It has been shown some time ago [40] that the spherically averaged atomic charge density p(r) can be represented in terms of decreasing exponentials, and that there is a change in the rates of the exponentials that takes place at certain points associated with the minima in the radial density. Politzer and Parr [1] have considered the outermost minimum of the radial density to establish a boundary between the core and valence regions in atoms. It is obvious that if the
-
(r)ls(KS)
1,
.
.
.
.
.
.
.
.
10
O
09 05 .~ 07 06
-
=
1
WDA
.6393rmax
-- 0 . 0 3 2 4 9 , SO there is v e r y
little difference in using HF or KS data for the plot. A final comment concerning Figs. 6, 7 and 8 is that the slope of the fitted lines becomes closer to 1.5 if He is not included,
o5 ~ 0, o3 0~ 01 00 00
o o
01
0.2
o
03
04
05
06
07
r,u=,(WDA) (a.u.)
Fig. 8. Mean value of r in the ls shell, obtained from Kohn-Sham calculations, versus the position of the first maximum of the WDA
radial density. The data set includes atoms from the main groups of the Periodic Table.
M.D. Glossmanet al./ ChemicalPhysics 196 (1995) 455-463 electron density can be represented as a sum of decreasing exponentials, then in a logarithmic plot of the density against the radial distance r, the points where there is a change of slope will reflect the minima in the radial density. In Fig. 9 we have plotted the HF and W D A densities for xenon in a logarithmic scale against r in a short interval close to the nucleus (r~< 0.5 a.u.). The W D A density, represented by a dotted line, reflects reasonably well the behavior of the HF density, represented by a solid line. Of course, discrepancies exist because the HF density shows more structure than the W D A density. The change in slope occurring in both curves at r - - 0 . 0 5 a.u. reflects the first minimum in the radial electron density. The change is more pronounced for the HF case (the distinction between the lS shell and the rest is better defined). The HF density shows other changes of slope at r = 0.17 a.u. and r---0.45 a.u., which reflect the completion of Ne-like and Ar-like cores, respectively. Those changes are not reflected in the W D A density which is unable to reproduce the full shell structure of the atom. But the accuracy of the W D A density in a region of the atom close to the nucleus is undeniable. Let us now study in more detail the behavior of the density in this region. The slope of In p(r) for r very close
13 ......... , ......... , ................... , ........ 12 11 lo HF 9 WDA 8 ~ ~' -= 6 5 4 3 2 1 o0
ol
02
03
04
05
r (a.u.) Fig. 9. Electronic density of Xe in a logarithmic scale, obtained from the WDA and HF approaches. Only the region close to the nucleus is shown.
461
100 80
~ 60 -~ g3 40
20 . 20
.
.
.
40
60
80
10o
atomic number Z
Fig. 10. Test of the cusp condition at the nucleus, a is the exponentin a fit of the WDA density near the nucleus to a functionp(r)= poe -"r.
to the nucleus should be - 2 Z according to the Kato's cusp theorem [41]. In this region the main contribution to pWDA(r) can be described by an Is type orbital, as shown in the preceding section. By fitting the W D A density near the origin to an exponential p(r) = po e-~r, we have obtained values of a very close to 2Z, that is, the cusp condition is well satisfied. In Fig. 10 we plot the values of a / 2 against Z for all the atoms that we have studied. The points in this plot fit a straight line a / 2 = 0.955Z + 0.0809 showing the accurate fulfillment of the cusp theorem. Notice the alternance below and above the straight line in going from a main group to the next group of the Periodic Table. Recently, the electron localization function (ELF) proposed by Becke and Edgecombe [31] has been used for visualization of electronic shells. Within this context, Gadre et al. [32] have formulated a densitybased electron localization function (DELF) via the nonlocal weighted density approximation [30]. Thus, this formulation appears to be adequate to study the shell structure within the WDA. The density-based electron localization function (DELF) can be explicitly written [32] 1 DELF = (1 + ~ ( r ) / p (
r))2/3'
(16)
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
462
1.o .. ........ , ......... , ......... , ......... , ......... it:\
~:~ li~ \ -,: li ', \ 08 71~ '\ \ tu' !~ \\ \ o 07 I', \ \ ~\ ,"'. ~\~ . '~'~'~:~'-":~'-" ~ 0.9
,
\
0.6
-
-
Ne
--
Ar Kr - - - xe ......
~
,
~
~
~ ~ ~ "--~.,,,,~
~.-'~...~...~,
05 ............................................... 0.0 0.2 0.4 0.6 o.a 1.0 r (a.u.) Fig. 11. Density-based electron localization function (DELF) of Ne, Ar, Kr and Xe atoms, obtained by the W D A approach (see
Eq. (16)of the text),
where p(r) is the WDA electron density and iS(r) is the effective weighted density obtained through Eq. (4). In Fig. 11 we have represented DELF against r for the neon, argon, krypton and xenon atoms. The kind of information obtained from this plot is consistent with the information obtained from the radial density (Fig. 1) or from the plot of In p (Fig. 9). Notwithstanding, it is noticeable that the transition from one region of the atom to the other is sharper in the DELF representation and consequently the points where this transition takes place can be located in an easier way. To conclude this Section it is worth to analyze the reason why the separation between the ls core and the rest of the electrons is not achieved by the usual extended Thomas-Fermi functionals, either the Thomas-Fermi-von Weizs~icker functional (in which the lf[(Vp)2/p]dr term enters with its full ~1 coefficient) or the functionals arising from a series expansion in the gradients of the density (in those, the Weizs~icker term enters with a coefficient reduced by a factor ~). In all those functionals the ThomasFermi term enters as a fundamental one, providing a large contribution to the energy in the region of the ls shell. On the other hand, in the WDA, the fact that tS(r)<< p(r) in that region leads to a drastic reduction of the Thomas-Fermi-like term compared
to the von Weizs~icker contribution, which is itself large in this region. Results qualitatively similar to the ones in this paper have been obtained recently by Wang and Teter [42]. The work of these authors fits within the spirit of the WDA, as formulated for classical liquids [43,44], although the differences with respect to our work are substantial. Wang and Teter arrive to a kinetic energy functional that also contains the full Weizs~icker term, plus a nonlocal term from which a Thomas-Fermi-like contribution with a reduced coefficient is subtracted. Furthermore, the functional contains a free parameter which affects not only the numerical constants in front of some terms, but also the powers of the electron density. When applied to a very simple model of closed shell atoms (in which both exchange and the classical electron-electron Hartree t e r m w e r e neglected) shells qualitatively similar to those we have reported in Fig. 1 were obtained. It is gratifying that similar results have been found with two different kinetic energy functionals (and two very different treatments of the other energies). The agreement in the main results has to be ascribed to the basic WDA spirit used in both works.
4. Conclusions Several tests have been performed in order to study if the radial density obtained by solving the Euler equation associated with the WDA model for free atoms displays some kind of shell structure. A direct inspection of the radial electronic density shows that, except for light atoms, there are two maxima whose position and relative height vary following the same trend as the first two maxima of the HF radial density. For light atoms only one maximum and a pronounced shoulder at the position of the second are found. The same behavior has been observed when using other different criteria like plots of the density in a logarithmic scale and plots of the density-based electron localization function (DELF). Arguments have been given in favor of the interpretation that the first maximum is associated to the Is shell, and additional tests have been carried out relating the position of the first maximum in the
M.D. Glossman et al. / Chemical Physics 196 (1995) 455-463
463
W D A radial density with the expectation value o f r
[15] M. Kakkar and K.D. Sen, Chem. Phys. Letters 226 (1994)
for I s orbitals in H F calculations. A l t h o u g h the W D A functional presented in this
241. [16] R.G. Parr and W. Yang, Density functional theory of atoms and molecules (Oxford Univ. Press, New York, 1989). [17] K.D. Sen, M. Slamet and V. Sahni, Chem. Phys. Letters 205 (1993) 313. [18] E. Clementi and C. Roetti, At. Data Nucl. Data Tables 14 (1974) 177. [19] E.S. Kryachko and E.V. Ludefia, Density functional theory of many-electron systems (Kluwer, Dordrecht, 1990). [20] R.M. Dreizler and E.K.U. Gross, Density functional theory. An approach to the quantum many-body problem (Springer, Berlin, 1990). [21] P. Csavinszky, Phys. Rev. A 24 (1981) 1215. [22] P. Csavinszky, Intern. J. Quantum Chem., Quantum Chem. Symp. 23 (1989) 591. [23] P. Csavinszky, Phys. Rev. A 39 (1989) 1519. [24] M.D. Glossman, A. Rubio, L.C. Balb~is and J.A. Alonso, Intern. J. Quantum Chem., Quantum Chem. Symp. 26 (1992) 347. [25] M.D. Glossman, A. Rubio, L.C. Balbfis and J.A. Alonso, New J. Chem. 16 (1992) 115. [26] M.D. Glossman, A. Rubio, LC. Balbfis, J.A. Alonso and LI. Serra, Intern. J. Quantum Chem. 45 (1993) 333. [27] M.D. Glossman, A. Rubio, L.C. Balbfis and J.A. Alonso, Phys. Rev. A 47 (1993) 1804. [28] M.D. Glossman, L.C. Balbfis and J.A. Alonso, Intern. J. Quantum Chem. 49 (1994) 171. [29] O. Gunnarsson, M. Jonson and B.I. Lundqvist, Phys. Rev. B 20 (1979) 3136. [30] J.A. Alonso and L.A. Girifalco, Phys. Rev. B 17 (1978) 3735. [31] A.D. Becke and K.E. Edgecombe, J. Chem. Phys. 92 (1990) 5397. [32] S.R. Gadre, S.A. Kulkarni and R.K. Pathak, J. Chem. Phys. 98 (1993) 3574. [33] P.A.M. Dirac, Proc. Cambridge Phil. Soc. 26 (1930) 376. [34] C.F. Weizs~icker, Z. Physik 96 (1935) 431. [35] W. Yang, Phys. Rev. A 34 (1986) 4575. [36] B.M. Deb and S.K. Ghosh, Intern. J. Quantum Chem. 23 (1983) 1. [37] O. Gritsenko, N. Cordero, A. Rubio, L.C. Balbfis and J. A. Alonso, Phys. Rev. A 48 (1993) 4197. [38] J.P. Desclaux, At. Data Nucl. Data Tables 26 (1973) 311. [39] C. Froese-Fischer, The Hartree-Fock method for atoms (Wiley, New York, 1977); Computer Phys. Commun. 14 (1978) 145. [40] H. Weinstein, P. Politzer and S. Srebenik, Theoret. Chim. Acta 38 (1975) 159. [41] T. Kato, Commun. Pure Appl. Math. 10 (1957) 151. [42] L.W. Wang and M.P. Teter, Phys. Rev. B 45 (1992) 13196. [43] P. Tarazona, Phys. Rev. A 31 (1985) 2672. [44] W.A. Curtin and N.W. Ashcroft, Phys. Rev. A 32 (1985) 2909.
paper is only able to separate the l s shell f r o m the rest o f the electrons of the atom, our results give s o m e c o n f i d e n c e that i m p r o v e d W D A - t y p e functionals, containing the v o n Weizs~icker term as a fundamental ingredient m a y be able to reproduce the full shell structure o f the atoms in the future. If this goal is achieved, it will constitute a significant step in density functional theory.
Acknowledgment
This w o r k has b e e n supported by D G I C Y T (Grant P B 9 2 - 0 6 4 5 - C 0 3 - 0 1 ) and Junta de Castilla y Le6n. M.D.G. is m e m b e r of the R e s e a r c h Career o f the C o n s e j o N a c i o n a l de Investigaciones Cientlficas y T6cnicas ( C O N I C E T ) , R e p u b l i c o f Argentina, f r o m w h i c h support is gratefully a c k n o w l e d g e d .
References [1] P. Politzer and R.G. Parr, J. Chem. Phys. 64 (1976) 4634. [2] R.J. Boyd, J. Chem. Phys. 66 (1977) 356. [3] R.F.W. Bader and M.E. Stephens, J. Am. Chem. Soc. 97 (1975) 7391. [4] R.F.W. Bader, P.M. Beddall and P.E. Cade, J. Am. Chem. Soc. 93 (1971) 3095. [5] R.F.W. Bader, Accounts Chem. Res. 8 (1975) 34. [6] Y. Tal and R.F.W. Bader, Intern. J. Quantum Chem. S12 (1978) 153. [7] P. Politzer, J. Reuther and G.T. Kasten, J. Chem. Phys. 67 (1977) 2385. [8] J.T. Waber and D.T. Cromer, J. Chem. Phys. 42 (1965) 4116. [9] R.J. Boyd, J. Phys. B 9 (1976) L69. [10] W.P. Wang and R.G. Parr, Phys. Rev. A 16 (1977) 891. [11] A.M. Simas, R.P. Sagar, A.C.T. Ku and V.H. Smith Jr., Can. J. Chem. 66 (1988) 1923. [12] R.P. Sagar, A.C.T. Ku, V.H. Smith Jr. and A.M. Simas, J. Chem. Phys. 88 (1988) 4367. [13] Z. Shi and R.J. Boyd, J. Chem. Phys. 88 (1988) 4375. [14] P. Politzer, J.S. Murray, M.E. Grice, T. Brinck and S. Ranganathan, J. Chem. Phys. 95 (1991) 6699.