Study of interplanar binding in graphite by extended Thomas-Fermi theory

A model of a graphite crystal is used which consists of a set of parallel slabs of positive charge 0 immersed in an electron sea. Each slab, about 1 A wide, contains the charge of the nucleus and five electrons per carbon atom, homogeneously distributed in the volume of the slabs. The electron density in the region between slabs is calculated from Thomas-Fermi-Dirac theory including corrections for inhomogeneity to the kinetic energy and correlation energy. Also, a calculation is reported with the electron density obtained by a minimization of the Thomas-Fermi-Dirac-Kirzhnits functional. The results are in semiquantitative agreement with empirical data.


I. INTRODUCTION Density-functional
theory has allowed a detailed inter- pretation of the electronic structure of many solids.In general, the study involves the solution of a set of self- consistent one-body Schrodinger equations and gives a quantitative agreement with empirical data.(A thorough review is given in the recent book by Lundqvist and   March.') In view of the success of these methods, it can be questioned whether the simpler, but less accurate, Thomas-Fermi (TF) theory and its extensions are still useful.However, there are two reasons for using TF theory.In the first place, there are systems so complex that a more refined calculation may be too involved.In the second place, theories resting upon a single equation, such as the TF theory and its extensions, may give the general trend for the variation of properties which, to some extent, is lost in more detailed calculations.The paradigmatic example is the TF model of the atom.
One of the most difficult tasks of TF and related theories is the quantitative interpretation of the chemical bond.In fact, although extended TF theory is able to predict the different energy contributions (kinetic, direct electrostatic, exchange, etc.) of a many-electron system with small errors (of a few percent), it is well known that binding comes from a delicate balance between these terms, and errors may accumulate rather than cancel.Actually, examples of a quantitative success of TF theory in the prediction of bond properties are scarce, although many qualitative agreements have been reported.
The first problem that appears in the TF theory of the chemi- cal bond is that pure TF equation does not predict binding between atoms, as was rigorously proved by Teller.
Furthermore, the inclusion of the local density (Dirac) correction for exchange is not enough, and gradient terms must be included.However, gradient corrections should be used carefully as they represent the first terms in a series expansion not convergent in general.On the other hand, no practical alternative has been found other than an exact treatment of the kinetic energy term, which leads again to one-body equations, thus departing from TF theory.
The calculation of the interlayer force in graphite is one of the few examples in which extended TF theory gives an accurate interpretation of the bond properties.Graphite has the unique property of having two entirely different types of interaction binding the structure.In the basal plane, the C atoms are held in a two-dimensional hexago- nal lattice by strong covalent bonds.These bonds are highly directional and account for a binding energy of or- der 5 eV/C-atom.On the other hand, the planes are held to each other by much weaker energies, of order 0.05 eV/C-atom, which are neither covalent nor directional.
Fully ab initio band-structure calculations have been per- formed for graphite ' which have obtained excellent results for the electronic charge density, but these studies did not examine the total energy of the graphite system.
On the other hand, several calculations with extended TF theory have been reported ' and the present calculation represents a new contribution on this line.Graphite has several advantages with respect to other systems.In the first place, there are many reliable experi- mental data on the structural properties of graphite avail- able for direct comparison with theory, like x-ray diffraction, " low-energy electron diffraction, " neutron scatter- ing, ' and hydrostatic pressure.' ' In the second place, graphite and its intercalated compounds are materials with important practical applications.Finally, the bind- ing between layers can be quantitatively studied with an extremely simple "jellium" model, which gives rise to equations that can be solved very easily, even analytically in some cases.This model, which is described in detail in the next section, represents a kind of theoretical laborato- ry for the test of different improvements to TF theory.

31
The purpose of the present calculation is to try to give an answer to a number of questions that arise in the com- parison of two previous Thomas-Fermi (or density- functional) calculations of interplanar binding in graphite (papers by Santos and Villagra and by DiVincenzo, Mele, and Holzwarth' ).Both papers used a density functional of the form Here, the first term is the usual Thomas-Fermi energy Urz[p]= f dr[(3rrR /10m)(3/~) ~p ~3+(8'. )'E~], (2) p(r) being the number of electrons per unit volume and E(r) the total electric field at r.The term Uz is the first gradient correction to the kinetic energy, which, as calcu- lated by Kirzhnits, ' is (3) The last term of (3) gives no contribution for atoms, mole- cules, and solids with pointlike nuclei, so that it is possible to use, instead of (3), Besides the form of the density functional, the methods used in the commented papers differ in the model of graphite and the procedure used to get the density.DiVincenzo et al. used a simple superposition of densi- ties, each one calculated for an isolated layer of carbon atoms in hexagonal array, and no minimization of the functional was made.In contrast, Santos and Villagra, considered a much simpler jellium model of graphite, but calculated the density by minimizing the energy Uzz[p]+ Uz[p] for all electrons (6 per atom), which is equivalent to solve the Thomas-Fermi-Dirac equation.
The results of both calculations are similar, but the agreement with experimental data is better in the paper by Santos and Villagra.For instance, 4% error in the inter- layer distance compared with 15%%uo, ' both by defect; or 30/o error for the compressibility against 180%, ' both by excess.The main conclusion of this comparison seems to be that minimization of the functional is rather impor- tant for the calculation of the density.The conclusion, however, is preliminary due to the difference in the models used in both papers and the fact that Santos and Villagra neglected correlation energy.In the present pa- per we study the influence of both facts by performing new calculations similar to those of Santos and Villagra, but with a more elaborate graphite model (presented in Sec.II) and including correlation energy (Sec.III).In or- der to test the importance of a correct minimization of the functional, we have also made a calculation starting from the Thomas-Fermi-Dirac-Kirzhnits equation, which corresponds to minimizing Uzz+ U~+ U~(Sec.IV).It is worth noting that the initial purpose of the present calcu- lation was trying to improve the early results of Santos and Villagra rather than making a comparison with the work by DiVincenzo et al. and, actually, when that work was published, the present calculations were practically concluded.This explains why we have not used the same functional as DiVincenzo et al. for the correlation energy.

II. MODEL OF GRAPHITE
On a graphite crystal it is possible to consider three kinds of electrons placed at different regions of space: (1) There are electrons ls of the carbon atoms, which are situated inside spheres centered at each nucleus, with a radius of order ao/Z-0. 1 A, ao being the Bohr radius and Z the atomic number (i.e. , Z =6).These spheres are very small in comparison with distances between nuclei (1.42 A).
(2) There are electrons of the o bonds between carbon atoms in the same layer.These electrons lie in a zone, with thickness about 1 A, the center of which is the layer containing the nuclei.
(3) The remaining, 2p"electrons (one per atom) occupy chiefly the region halfway between layers, a zone more than 1.5 A wide.These are the electrons which contribute to the binding between layers and they are the only ones that we must explicitly consider in our calculation.
It is to be expected that the electron density changes quickly in the neighborhood of nuclei, but more slowly far from them.Also, we may assume that the density in the third zone considered above will change mainly in the direction normal to the layers, but it will scarcely change in the directions parallel to the layers.This is due to the fact that the distance between layers (3.5 A) is much larger than the distance between nearby nuclei in a layer (1.4 A).In this way, we arrive at a picture of a graphite crystal as a kind of metal having "ions" in the form of slabs of positive charge (with a homogeneous distribution equivalent to one elementary charge per atom) immersed in a sea of 2p, electrons which glue the "ions." This pic- ture helps us understand the metallic properties of graphite more easily than the usual picture, which sees the graphite crystal as a set of planar macromolecules bound by van der Waals forces.
We will study forces between layers of graphite by con- sidering the variation of the energy per electron of a system of parallel "ionic slabs" when the distance R between slabs changes.Our problem is just to evaluate the energy per electron, E(R), of an electron system placed at a re- gion bounded by two parallel planes at a distance R. The number of electrons per unit area in a plane agrees with the number of atoms per unit area of a graphite layer because we consider one electron per atom contributing to binding.
Then, the surface density of electrons is o =0.1067 a.u.[from now on, we will use hartree units (a.u. ) such that 8=m =e = 1].Hence, the laws of electro- statics plus symmetry arguments fix the electric field in each boundary plane to be 2no.=0.6706 a.u.After that, we shall evaluate the electric potential at each point in the region by means of the Thomas-Fermi-Dirac equation, which also allows for the calculation of the energy of the electron system.Before doing that, we must evaluate the width of each "ionic slab" in order to specify completely the model of graphite.
The calculation of the width of an "ionic slab" rests upon the assumption that the electron density near a nu- cleus does not change very much in going from the free atom to graphite except that we consider atoms with structure ls 2s2p (4 valence electrons) instead of the ground-state structure 1s 2s 2p .Then, the density around the nucleus is p(«)=2 I @i.I'+ I P2.I'+3 I Ap I' where P are the atomic wave functions averaged to spheri- cal symmetry.We take the atomic functions to be AJ kl -A, p p' P&, --ae, g2, --be -+ere beled R, the distance between graphite layers in our model will be R +2z0.In the following sections we will evaluate the energy per atom as a function of the interlayer dis- tance,

THOMAS-FERMI-DIRAC EQUATION
%'e assume that the total energy of a graphite crystal is the sum of additive contributions associated, respectively, with the "ionic slabs" and the regions between slabs.The former are assumed to be a constant, which can be set equal to zero by adjusting the scale of energies.The ener- gy associated to one of the regions between slabs is evaluated by the Thomas-Fermi method.The calculation parallels closely one presented in a paper by Santos and   Villagra where a more crude model of graphite was used which, however, led to a very similar calculational prob- lem.Therefore, we will not repeat the details here, and we give only a sketch of the procedure used and the results obtained.
We start with the Thomas-Fermi-Dirae (TFD) equation where the effective charges X; and the coefficients (all in a.u.These functions are normalized for the volume element r dr (i.e. , without the factor 4n. ).Now, we consider two parallel planes, each one at a distance z0 from the nucleus, in such a way that a net charge of 5 electrons of the atom remains between the planes and 1 electron outside.That is, we have f p(r)d r =5, where p(r) is given as above.
After a straightforward calculation we obtain z0 --1.1826 a.u.Therefore, the width of the "ionic slab" is 2z0 or about 1 A as was estimated at the beginning of this section.As the distance between "ionic slabs" was la- where P is the electric potential and x a coordinate per- pendicular to the layers.The problem is one dimensional due to the assumed translational invariance.That equa- tion is solved with the boundary conditions dX dX =0.6706 a.u., (10) where x=8/2 is a point halfway between layers (there the electric field dgldx is zero by symmetry); x =0 is the boundary of an "ionic slab, " in which an electric field ex- ists as explained in Sec.EI.Equation ( 9) with boundary conditions (10) does not have a solution when R ~8. 1 a.u., a fact related to the finite radius of the atom in the TFD theory.This prevents us from calculating accurately the exfoliation energy of graphite because we would need the energy for R ~oo.We will estimate that quantity by an extrapolation procedure explained below.
Once the function P(x) is obtained we calculate the en- ergy per atom through the TFD functional where o is the surface density of atoms in a layer and p(x) the electron density; which is related to the potential through Poisson's equation d =4mp .dx Note that our method of calculating the densityin par- ticular, the second boundary condition (10)implies that the valence electron density is not permitted to leak into the slab region.
The calculated energy is a function of the distance be- tween graphite layers.That function does not have a minimum, which means that the TFD theory does not predict a binding between layers.A similar result was also obtained in previous calculations.
However, a bind- ing is predicted if we add the correction for inhomogenei- ty to the kinetic energy of the electrons in the form of When this energy.is added to Eq. ( 12) a function is ob- tained which has a minimum at a distance Ao.The prac- tical method to do the calculations closely parallels the one used in Ref. 8. In particular, the integrals (ll) and ( 13) were evaluated numerically for several values of R, which allowed the calculation of the interlayer equilibri- um distance do, and the force constant k by adjusting a suitable second-order polynomial near the minimum.The calculated values are do=Ro+2zo=6 36 a u k:-[d (UgFD+ Ux )/dx ]"g --0.0023 a.u. ,( 14) to be compared with the empirical values 6.34 a.u.(Ref. 13) and 0.0019 a.u. (Ref. 12).As stated before, the calcu- lation of the exfoliation energy of graphite in TFD theory is ambiguous.We might estimate it from the difference between the TFD energy at the minimum and at the larg- est separation between layers which gives a continuous solution to Eq. ( 9).The exfoliation energy estimated this way is 0.00198 a.u. to be compared with the empirical' value 0.0016 a.u.The agreement is fairly good for an al- most ab initio calculation of such a complex structure as a graphite crystal.[The only empirical parameters used have been the o. bond length in graphite and the atomic wave functions of the carbon atoms ( 6).] We shall see later that the agreement is not as good as it seems.
A more correct procedure to calculate the relevant quantities is to do a fitting of the calculated TFD energies at different values of R by means of a suitable function and, then, to obtain the quantities from that function.We have chosen the function and the values obtained are reported in the first row of Table I.The agreement with empirical values is not as good as before.
Correlation energy can be included in the TF theory by adding a suitable term to the energy.We have used the %'igner' expression which, in our case, is

Rf2
Uc(R)= 2cr -' f ap i (1+bp'i ) 'dx ( 16) with a =0.44(4m/3)', b =7.8(4m./3)' Putting the density calculated from Eqs. ( 9) and ( 12) in Eq. ( 16) a numerical integration gives a term which must be added to UrFD(R).As it should be, we obtain a func- tion with a deeper minimum than without correlation en- ergy.A procedure similar to the first one explained above the best fitting of the function U(d)= -Ad +8 exp( -ad)+C to the calculated values between d=4 and 9 a.u.The equi1ibrium interlayer distance dp, force constant k =(d U/dR ), and exfoliation energy per atom are calculated from that function.A11 quanti- ties are given in hartree units (A=e =m =1).UTFD is the energy calculated from Thomas-Fermi-Dirac theory, U& ( UJc ) is the corn- plete (incomplete) inhomogeneity correction to the kinetic energy in the form of Kirzhnits, and Uc the correlation energy.U»DK is the energy calculated from the Thomas-Fermi-Dirac-Kirzhnits theory.The first 7 rows refer to the model of graphite described in Sec.II, rows 8 to 12 refer to the intermediate model described in Sec.III, rom 10 reports results of Ref.Again, a more correct procedure is the fitting of the func- tion (15) to the energies calculated for the different values of R. The parameters obtained this way are given in row 3 of Table I.As we see, the correlation energy shortens the interlayer equilibrium distance by about 3% and con- tributes to binding by about 25%.However, the agreement with empirical values is worse than before, which seems to indicate that the excellent agreement of Eqs. ( 14) is somewhat accidental.
With respect to the inhomogeneity correction a problem appears with the second term of Eqs.
(3) and ( 13).This term does not contribute in the TF theory of atoms and, more generally, when the external charges (nuclei) are considered pointlike.In our model of graphite this is not the case and, in our opinion, this term should be included.Indeed, it is quite important.In order to estimate the contribution of the term we have made also the calcula- tions without it.The results are reported in rows 2 and 4 of Table I, where the inhomogeneity correction (3') to the kinetic energy is labeled Uz .
The adequacy of our model of graphite can be studied through a comparison with the "jellium" model of Ref. 8.Here we consider a graphite crystal as a set of parallel "ionic slabs" of width 1.18 a.u.immersed in an electron sea with 1 electron per carbon atom.The "jellium" model of Ref. 8 can be seen as similar, but with "ionic slabs' having zero width and the electron sea having 6 electrons per carbon atom.The results with both calculations are not too different (see Table I), which gives some justifica- tion for the model proposed.It is also possible to consider an intermediate model in which the electron sea has 4 electrons per atom.That is, we put the nuclei and the 1s electrons in the "ionic slabs" but we treat all valence elec- trons (4 per atom) by the Thomas-Fermi method.The width of the ionic slab can also be calculated easily by putting 2 instead of 5 electrons per atom in each "ionic slab" [see Eq. ( 8)].The result is zo --0.175 a.u.The re- sults of the calculations with this intermediate model have also been made, and are reported in Table I.
For comparison, the results of the calculation by DiVincenzo et al. ' are also included in Table I.It can be seen that, after including correlation, our results are a lit- tle closer to theirs.

IV. CALCULATION FROM THOMAS-FERMI-DIRAC-KIRZH NITS EQUATION
In view of the influence of a correct minimization of the functional, we have also made calculations with a den- sity obtained by a minimization of the functional that in- cludes the Kirzhnits term (4) besides those previously in- cluded, ( 2) and (3).The Euler-I.agrange equation of that variational problem gives which leads to a fourth-order equation when combined with Poisson's Eq. ( 8).As a consequence, we need two boundary conditions, instead of only one, at each boun- dary.We choose, besides the electric field E given by Eq.
(6), the density p at the surface of each "ionic slab" and the gradient of the density halfway between slabs, this coming from symmetry considerations, i.e.
+CX e being an integration constant.This constant cannot be determined directly from (10) and (19) because we do not know d p/dE at the boundaries.Then, we used a con- sistency condition by testing several values of a for each fixed pz&2, until a numerical integration of (20) gives the correct po [Eq. ( 19) for Eo [Eq. ( 10)].This procedure gives numerically the density as a function of the electric field.On the other hand, the distance from the middle plane between slabs is easily obtained from x = (2m ) ' I p(E') 'dE',   (21)   which, at the end, gives numerically p as a function of x.
After that, the numerical evaluation of the energy from the functional is straightforward.
A difficulty appears with the choice between (3) and (3')   for the gradient correction.If the energy is calculated from the functional (1) over the whole volume, including the "ionic slabs, " then the second term of (3) gives no con- tribution, so that (3) and (3') are equivalent (except for a global constant in the energy).However, if the functional is applied only to the region outside the slabs, then the second term of (3) gives an important contribution to the energy.It seems that the first procedure is more correct and, indeed, we checked that it gives results closer to the empirical data.Without including correlation, we ob- tained the data reported on rom 5 of Table I.With corre- lation energy included we got those of row 6.
The density po is calculated by averaging the atomic den- sity (5) on the plane z =zo.In comparison with the calcu- lation of the previous section, we see that an additional constraint is imposed by demanding continuity of the den- sity at the surface of the "ionic slabs." Equations ( 18) and ( 12) give rise to a fourth-order dif- ferential equation, that can be reduced to a third-order one by using the electric field modulus E as variable, in- stead of the coordinate x.This equation can be integrated once and leads to It can be seen that the more complex minimization pro- cedure described in this section does not improve the agreement with empirical data.This seems to be a conse- quence of the smaller freedom allowed by the additional constraint (19).Then, we tested a minimization of the functional without the first condition (19), i.e. , choosing for each interlayer distance the value of po that minimizes the functional.The results appear in row 7 of Table I.
The conclusion of all these calculations is that a minim- ization of the Thomas-Fermi-Dirac-Kirzhnits functional, giving rise to a fourth-order differential equation, is not reliable.An alternative procedure to get the density from a second-order equation that includes, nevertheless, gra- dient corrections is the equation derived recently by Schwinger and generalized by Santos and Leal.' (Another generalization has been given very recently by Englert and Schwinger. ) We hope that these equations will allow more reliable calculations of binding properties within extended TF theory and work in this direction is in progress.In any case, the results of the present calcula- tion encourage the use of TF theory for more complex problems, like the study of intercalated graphite com- pounds, which has already been initiated.
form taken by DiVincenzo et al. ' However, in the jellium model of Santos and Villagra, where the pointlike nuclei were replaced by layers of positive charge, the form (3) was used.The term U",[p] contains ex- change and correlation energy.Santos and Villagra.neglected correlation and used for exchange the Dirac ex- pression U~[p] =( 3e /4)(3/n-)' p (4) while DiVincenzo et al. considered both exchange and correlation by means of the expression proposed by Hedin and Lundqvist, ' whose leading term is just (4).

TABLE I .
Calculated values compared with results of Refs.8 and 10 and empirical values.Parameters u, A, B, C correspond to 8, row 11 reports the results of the calculation byDiVincenzo et a1.(Ref.10),and the last row reports the empirical values, taken from Refs. 13, 12, and 18, respec- tively, for dp, k, and Up. the following values for the interlayer distance, force constant, and exfoliation energy, Uo. .