Thomas — Fermi Calculation of the Interlayer Force in Graphite *

the [100] TA mode. This would imply only a 3% error or so in the corresponding fluorine-ion eigenvectors. As a method for checking these eigenvectors, isotope-induced absorption is very sensitive but allows access to only a few modes. Others are forbidden by symmetry or swamped by two-phonon summation absorption. The frequency of the X, symmetry point has been found to decrease by 2-3% from the room-temperature value measured in Li F, opposite to the shift direction normally expected. Finally, by using the shell-model data in the region above the reststrahlen frequency, it has been shown that the isotope-induced absorption may be responsible for some of the small features in the high values of conductivity measured by other authors.


I. INTRODUCTION
The method of Thomas -Fermi-Dirac (TFD) gives the energy of an electron system by means of the integral' -3~2@2 3 't2j'3 E2 3 2 3 )1/3 U= dv ' - I p"'+ ---I p"' iOm ~8~4 (, where lt is the reduced Planck constant, m is the electron mass and -e its charge, p is the number of electrons per unit volume, and R is the electric field.The first term of the integral gives the kinetic energy of the electrons, the second one gives the potential energy, and the last one gives the exchange energy.The potential energy is given in terms of the electric field for later convenience. The integral (l) must be a minimum subject to the following conditions: the Poisson equation, which relates the electrostatic potential to the charge density, the normalization of the density, and the boundary conditions.The latter are usually stated by giving the positions of the external charges or the electric field at the boundary of the volume occupied by the electrons.The TFD method has been applied with success to the determination of prop- erties of electron syst~~ms when these properties are not very sensitive to the charge distribution, as is the case, for example, in the calculation of the effective potential in atoms.However, the method fails in the calculation o~e interatomic forces because, as 'i'eller has proved, the TFD method without corrections does not explain even the stability of molecules.Nevertheless, there are recent calculations of this type.' As is mell known, the TFD method is a semi- classical approximation of the self -consistent-field Hartree-Fock method, with which the TFD coincides in the case of an homogeneous electron gas.The semiclassical approximation is obtained when the zeroand first-order terms of a series expan- sion in powers of h are retained.Therefore, the obvious corrections to this approximation should be to take into account more terms in the series expansion.In this way, inhomogeneity corrections to the TFD method have been derived mhich are also called quantum corrections.The derivation of the quantum corrections, however, is not as simple a matter because it is known that the series expansion of the density in powers of 5 does not almays converge.Therefore, it is not obvious that the corrections actually refine the TFD meth- .od, and several different inhomogeneity corrections have been proposed to date.
An inhomogeneity correction to the Thomas- Fermi (TF) method was proposed long ago by Weizsacker.He justified, in an intuitive form, a correction to the energy given by the integral U"= (h'/8m) f [(Vp)'/p] dv . (2) The rigorous derivation of the quantum corrections has been made by Kompaneets and Pavlovskii, ' Kirzhnits, and others.' The quantum corrections of lower order have been found to be of the form (3) Actually, only the first term of Eq. ( 3) is usually quoted in the literature because the second one is zero when the external charges are pointlike.In the model of graphite that we will study this is not so, and the second term of Eq. (3) must be taken into account.The existence of the second term is clearly seen in the paper by Kirzhnits, for exam- ple.
Kirzhnits has shown that the correction of Eq.
(3) reduces the error in the energy of atoms from 20 to about 3%%up, while the Weizsacker correction [Eq.(2)] is too large.If the calculation of the en- ergy of an atom is performed correcting also for the lack of validity of the TFD method in the neigh- borhood of the nucleus, the Kirzhnits method [TFD plus Eq. ( 3)] gives a result which is close to the Hartree-Fock energy, differing by about 0. 2%%up.' However good these results are, the value of the inhomogeneity corrections has been questioned ow- ing to the lack of convergence of the series ex- pansion of the electron density in terms of 5.
Actually, Dagens has shown that the series expansion of the energy may be convergent even if the series expansion of the density is not.Another difficulty of the Kirzhnits method is that the in- homogeneity corrections to the exchange energy [not taken into account in Eq. ( 3)] may be as important as the corrections to the kinetic energy when the density is low.' Even the form of the corrections to the kinetic energy is questioned, and some authors believe that Eq. ( 3) is not the best one.' Actually, it has been shown that the corrections of Weizsacker [Eq.( 2)] and Kirzhnits [Eq.( 3)] give upper and lower bounds to the ener- gy, each one being a better approximation in specific cases.This discussion shows that the prob- lem of the inhomogeneity corrections to the TFD method is not solved, so that it is interesting to test the different corrections in some particular electron systems.
A very sensitive test for the correctness of a calculated electron density is the use of this den- sity in the study of interatomic forces.A number of calculations of interatomic forces by the TF method are found in the literature.Most of them refer to noble gas atoms, but the covalent binding of some homonuclear molecules has been studied also by this procedure. 4However, it seems that the calculation of the covalent binding energy has little value because the TF method, without cor- rections taking into account the granularity of electrons, cannot give rise to shell effects, which are essential in the interpretation of the covalent bond.In fact, the method will predict similar bonds for molecules with a similar number of electrons, such as nitrogen, fluorine, and neon, but it is well known that the first two molecules are stable and the last one is not.This shows clearly that the method is not adequate in the study of covalent bond if corrections taking into account the granu- larity of electrons are not included.Some general corrections of this type have been proposed, such as that by Fermi and Amaldi (see Gombds') and by Plaskett et aL ' Without these corrections, which do not seem easy to include in the study of interatomic forces, it is obvious that only forces which depend on the global distribution of electrons can be calculated by means of the TF method.On the other hand, intermolecular long-range forces are also excluded because they depend essentially on the electron correlation, which is not taken into account in the TF method.(These are frequently E. SANTOS AND A. QILLAGHA called van der Waals forces, although this name is sometimes used with a wider meaning.) Therefore, it seems that only some medium-range forces can be calculated with the TF method.This seems to be the case of the interlayer force in graphite, and this assumption is confirmed by the present calculation.
In this paper, we deal with a very simplified problem in which the TFD equation can be analytically solved.Aside from the interest of the calculation of the interlayer force in graphite, the purpose of this paper is to test, in a particular problem, whether the inhomogeneity corrections are an essential refinement to the TFD method and which is the best method to take it into account.We mean that one of the purposes of this paper is to compare the different inhomogeneity corrections (i.e. , Weizsacker, Kirzhnits, etc. ) and to see which is the best one; also, to see how these corrections must be included in the calculation.We have made calculations using the corrections in the form of Weizsacker [Eq.( 2)] and in the form of Kirzhnits [Eq.( 3)].Another aim of the calculation is to make a comparison between the results ob- tained from TF and TFD equations.In the first case, the electron density is found by solving the TF equation, and the exchange (Dirac) and inhomo- geneity corrections are calculated as perturbations.
In the second procedure" the exchange is taken in- to account from the very beginning calculating ihe density from the TFD equation, while the inhomo- geneity correction is considered a perturbation.(Calculations are in progress in which both exchange and inhomogeneity are taken into account in the starting equations.)

II. MODEL OF GRAPHITE
A graphite crystal consists of a set of planar layers, the interlayer distance being 3. 354 A.
The carbon atoms in a layer form a hexagonal lattice, the distance between nearest-neighbor nuclei being 1.421 A. The bonds in every layer are strong, while the forces between layers are much weaker.In order to find the interlayer force, we must calculate the energy per atom of a crystal as a function of the interlayer distance.We will assume that the interatomic distance in a layer does not change when the interlayer distance changes and that the vibrational energy of nuclei is negligi- ble compared with the electronic energy.The en- ergy per atom U will be a function of the interlayer distance R with a minimum at the equilibrium dis- tance R0=3. 354A.
The first step in the calculation is the solution of the TF or TFD equations once the position of every nucleus is given.These equations, depending on three variables and being nonlinear, are extremely difficult to solve accurately.We very much sim- plify the equations by substituting homogeneous, planar distributions of electric charge for the actual charges due to the (pointlike) nuclei.The surface density of charge is so chosen that the total positive charge of a large crystal will coincide with the actual charge of all nuclei in it.
is the interatomic distance, the surface density of atoms in a graphite layer is n= 4/(3v 3 d ) and the charge density g= 6en=0.6405 a. u. [From now on, we will use hartree units (a.u. ) such that 5 = m = e = 1.] Hence, the electric field in the neighborhood of a layer is normal to it and has a magnitude IEOI =E0=2z0=4.024 a. u.In order to justify the use of this model of graphite, it is necessary to show that the change in energy with the interlayer distance is almost independent of the actual distribu- tion of the positive charges of nuclei.It is not needed to prove that the absolute energy is inde- pendent of the model.In fact, the interlayer force depends very much on the electron density in the mean region between layers, but it is probably almost independent of the density in the neighborhood of a layer.It is reasonable to assume that the electron density in the mean interlayer region does not depend very much on the detailed charge distribution at the layers.A test of these hypotheses is the agreement between the results of the calculation and the experimental data.If the agreement is good it is very probable that both the model of graphite and the TF method of calculation are good.Nevertheless, a disagreement may be due to the failure of either.
An estimate of the error involved in our model of graphite can be found by considering the changes in the interlayer force when the layer stacking changes.It is known" that the force constant can change from 2. Ox 10 up to 3. 6x 10 dyn/cm in different specimens of graphite, even when the interlayer distance changes less than 3%, although it is not clear if these changes are in part due to lattice defects.In any case, this does not invali- date our model as a test for the TF method, be- cause the different corrections of the TF theory give rise to changes in the force constant of graph- ite of several orders of magnitude, as will be seen at the end of the paper.
The proposed model of graphite extremely simplifies the starting equations, which become ordinary, instead of partial, differential equations.Indeed, they can be solved analytically.
If x is the distance from a point inside the crystal to the nearest layer, the equations are These equations must be solv a. u. , E -, 'R) = 0 . (6) In order to state cle .a e c early the si ' ns, we also write t 'gn and unit cong, the electri f' ensity p.
ic ield E, and th e electron Itis tobbe noted that the zero chosen so that the Ferm e Fermi energy is zer 4) is e is easily reduced t where C is a c 4 a constant.Takin g (8) The boundar co ion: ry conditions lead to the following rela-2( 5 2)ii5 I'@5 5 -5 f '(E'+C)-"'dz .
This function added to the U»(R) appears also in Fig. 1.
Finally, the inhomogeneity correction was ob- tained from both Eqs. ( 2) and (3).The last one gives the integral U»= (A/120(( ) (60(() i [ f (E + C) i'E dE 0.04 -5Eo(EO+C)' + ~Eo '] 002 where also the energy was redefined so that U»(~) = 0. Integration by parts gives x Ca. u.) which reduces U~to known integrals.From this, the function U»+ UL)+ U~was obtained, which also appears in Fig. 1.The calculation of the inhomo- geneity correction from Eq. ( 2), which is straight- forward, gives a total energy U(R) from which a much stronger interlayer force is obtained, the interlayer equilibrium distance predicted being less than 3 a.u.
The comparison between thefunctions U(R) shows that the exchange and inhomogeneity corrections are essential in the calculation of interlayer forces.It is interesting to note that the stability of graphite is predicted without taking into account the inho- mogeneity correction if the exchange correction is calculated as a perturbation.This is probably an accident because, according to the Teller theorem, the inhomogeneity corrections are essential in explaining the stability.
It has not been necessary to obtain explicitly the electron density in order to find the energy, but it is interesting to perform the calculation.From Eqs. ( 4), (7), and (8) we obtain 4((p= (250/9(( )'i'(E + C) i' .
Eliminating E between this equation and Eq. ( 9) and integrating, we obtain a relation between the den- sity p and the position x in integral form.The in- tegration was performed numerically.
In Fig. 2, the function p(x) is shown for two values of the interlayer distance, R = R0 and R = .The density as a function of x for an infinite interlayer separa- tion (i.e. , C= 0) has the simple analytical form p(x) = ( (( (() [x+ (225(( /2EO)'i ] This is the electron density which the TF method gives for the surface of any metal if the positive charges of nuclei are approximated by an homo- geneous distribution.
FIG. 2. Electron density as a function of the distance to the nearest layer.Dashed lines, from YF equation; continuous lines, from TFD equation.In each case, upper line for equilibrium interlayer separation, lower line for infinite interlayer separation.

EQUATION
The first step is to combine Eqs. ( 5) and (7) in order to obtain the following relation between E and p: E= [ (( (3/(() i p -2(((3/(()' p -C]' (15) where C is an integration constant.On the other hand, from Eq. ( 7) one derives dE EO 1 "' Ep, c 2(( o p 2((p (0) 2 (( p(((i 2)   p' ( 16) where E(p, C) is the right-hand side of (15).Two relations between C, p(-, 'R) and p( 0) can be obtained from Eq. ( 15) taking into account the boundary con- ditions (6), which can be written Hence, the values of p(0) and C were numerically calculated for each value of p(~R).In practice values of p(-, 'R) were taken between 0. 002127 5 and 0. 089 63 a. u.From the values of p(2R), p(0), and C, those of Bwere found by numerical inte gration of Eq. ( 16).The values of R so obtained lie between 8. 60 and 3. 01 a. u.The interlayer distance R = 8. 60 a. u. is the greatest one for which the solution of Eq. ( 5) has no singularities.This value of R corresponds to the lowest possible value of p(2R) (= 125/192m'), which is obtained from Eq. ( 15) taking C = 0 and E=O.If R &8. 6 a. u. the TFD equation gives a discontinuous density, which is zero in all points at a distance greater than 4. 3 a.u. from the nearest layer.This is similar to the result that atoms have finite radii in the TFD theory.
The calculation of the energy is straightforward if Eq. ( 1) is integrated from 0 to 2R in the coordinate x in a similar form to that indicated in Sec.
II of this paper.Going from the variable x to E by means of Eq. ( 7) and then to p by Eq. ( 15), the following expression is obtained [after a redefinition of the energy similar to Eq. ( 12), but such that U(R= 8. 60)= 0]: UT Fn= (A/14411) j[10lT (Sp(0)/'w) 26(3p(0)/7) ]Eo + CR+2(3/iT) J p ~E(p, C)dp} 7. 7123 .P( a/2) The integral was evaluated numerically and the values of the energy so obtained are plotted in Fig. 3.The energy is chosen to be zero for the greatest interlayer separation R = 8. 60 a. u. because the TFD method gives constant energy for R &R The function U»D has no minimum as was expec- ted from Teller theorem.It is seen that the calculated energy is quite different from that obtained starting from the TF equation and adding the ex- change correction as a perturbation.
The inhomogeneity correction is obtained from Eq. (3) which, taking into account Eq. ( 7), gives f'P (0) The derivatives are obtained from Eq. ( 15) and the integral in (17) was calculated by a numerical method.The values of U~, added to that of V»D, give the function U(R) which is shown in Fig. 3.It is not possible to calculate the energy for values of R greater than 8. 6 a. u. because the function p(x) is discontinuous in this case and the integral ( 17) is not defined.
The agreement of the final function U(R) with experimental data is very good, as will be seen in detail in Sec.IV.Nevertheless, a calculation was performed taking into account the inhomogeneity correction also in the form of Weizsacker [Eq.
(2)].The procedure is similar to that described starting with Eq. ( 3).The energy so obtained, U~(R), when added to U»n gives a minimum of the total energy for R &3 a. u. far from the experi- mental equilibrium distance.As a conclusion, it seems that the Kirzhnits correction [Eq.(3)] is much better than the Weizsacher correction [Eq.
The electron density p(x) can be obtained from Eqs. ( 7) and ( 15).The resulting differential equa- tion can be solved analytically if C = 0, which cor- . However, the analytical form is rather complex and it is better to make a numerical integration.The function so obtained is shown in Fig. 2 in comparison with that obtained from the TF equation.A numerical calculation was performed also for the experimental interlayer distance R=RO which is also shown in Fig. 2. It is to be noted that the function p given by the TFD equation is not strictly the electron density, the difference being important when the number of electrons is small.' In a graphite crystal the number of electrons is indefinite and the difference is probably small.

V. DISCUSSION
The experimental data about the interlayer force in graphite are very limited.In fact, only two quantities have been measured with good precision: the interlayer distance and the force constant.' Actually, the force as a function of the interlayer separation is known between R =2. 90 A and Ro = 3. 35 A. ' The interlayer distance corresponds to the position of the minimum of the function U(R) and the force constant to the second derivative of this function at the minimum.
The absolute value of the minimum, i. e. , the binding energy, has not been measured.
In our calculation, the interlayer distance may be very sensitive to the model of graphite, although the agreement with the experimental datum would support both the model and the TF theory.
Until now, there was only one theoretical calculation of the interlayer force in graphite, due to Brennan, who used the molecular-orbital method 20 .07 (G 37)c 0.0114 0.02GO 0.0039 0.00GG 0.001 G9a 0. 0057 0. 01Gg 0. 002;y;j 0. 0019 in the linear combination of atomic orbitals ap- proximation (LCAO).The results of Brennan, how- ever, were very sensitive to the effective charge of the carbon atoms, which became, in fact, an adjustable parameter.
Furthermore, Brennan calculated only the repulsive part of the potential, taking a London-type term -A/R for the attractive part.The parameter A was obtained from the condition that the minimum of U(R) would coincide with the experimental datum.In this way, the Brennan interlayer potential depends on two parameters which were determined from two empirical data (force constant and interlayer separation).A number of interlayer potential functions can be found in the literature, ' obtained empirically or with reference to the paper of Brennan, the most recent being due to Drickamer et al. ' In Table I the experimental data are compared with the results of the present calculations.The comparison clearly shows that good results are ob- tained only when the electron density is calculated from the TFD equation, and inhomogeneity corrections are taken into account by the method of Kirzhnits.As was indicated, the binding energy has not been measured, and we compare our results with that of Brennan for this quantity (see Table I).
The agreement between the result of Brennan and that obtained from the TFD equation plus Kirzhnits correction is also satisfactory, but this is not the case for the other TF calculations.Our calculation of the binding energy as the difference U(~) -U(RO) is straightforward if the electron density is obtained from the TF equation, but it is not possible from the TFD density, which is discontinuous for R = 8. 6 a. u.In this case an extrapolation was performed, which is explained in the Appendix.
The agreement between the function U»D+ U~and the experimental data is better than is usual in TF calculations.
(It is to be noted that force constants are extremely difficult to obtain theoretically with accuracy.) An explanation of this fact can be found From Ref. 19.See, however, comment on this quantity in Sec.II.
"See Appendix for the procedure to determine the parameters.
'From the best empirical data at the time of Brennan calculation.
by considering the conditions of validity of the TF method.The semiclassical approximation which leads to the TF equation is valid when the following inequality holds 'd% dx d I 2&]&(x) I "' dx (18) where the potential Q(x) is measured from the Fermi level and atomic units are used.The values of x for which condition (18) holds depend on the function Q(x) and this function depends on the starting equation (TF or TFD) and the value of the inter- layer distance R. Nevertheless, the difference in the potentials (or in the densities) is not great, as can be seen in Fig. 2. Therefore, we make an estimation of the left-hand side of inequality (18) using the simplest function, which is Eq. ( 14).In this way, the following condition is obtained: x«9a.u.
(19) APPENDIX: ANALYTICAL POTENTIAL FUNCTION As a summary of the calculations, it is useful to obtain a simple analytic function U(R) fitting the It is seen that all calculations have been made un- der this limit because the inequality x & &R always holds and the calculations have been limited to -, 'R &6 a. u. in solving the TF equation.In the solution of the TFD equation, the exchange term imposes the condition x & -'R & 4. 3 a.u. , which is stronger than inequality (19).
The calculation of the electron density in our model of graphite can be compared with that of the free atom.The TF method cannot be applied in the region of the atom which is close to the nucleus because the potential changes rapidly in this region and condition (18) does not hold.In our model of graphite, however, the condition (18) is fulfilled up to the layer of positive charge because the singularity of the nuclear potential has been eliminated.On the other hand, the TF method does not hold in the external region of the atom and this is also true in graphite.However, the planar sym- metry in this caseinstead of the spherical one in the atommakes the charge density decrease more slowly in graphite, so that the region of validity is greater.
As a conclusion of the calculations, it seems that the TF method can be used to determine medium range interatomic forces in quantitative agreement with experimental data.For this, it is necessary to calculate the electron density from the TFD equa- tion (the TF equation being not enough), which must be solved with great precision.Corrections for inhomogeneity must be included in the form of Kirzh nits.numerial values of the best interlayer energy (U»n + Ur).As was indicated, the TFD method does not permit a calculation of the energy for an infinite interlayer separation.
This shows that an unknown constant must be added to every calculated value.On the other hand, the minimum of the function can be sensitive to the model of graphite.For thege reasons we have not identified U»D+ U"with U(R) but rather we have assumed the following relation: U»D(R -&R)+ Ur(R -6R) -C = U(R) . (20) The constant 6R has been determined so that U(R) has a minimum at the experimental interlayer distance, which gives 5R = 0. 264 a. u.At first we chose for U the form used by Drickamer et al. , ' which is U(R) = -a/R +B/R' but the fit obtained was poor.Finally, we found that a good fit is obtained using the function The best values of the parameters are A = 11.868a. u. , n = 1.338 a. u. , and B = 16. 719a. u.The constant C of Eq. ( 20) has the value 0. 0016 a. u. , which was used to estimate the binding energy FIG. 1. Int nterlayer bindin en , YF energy U D, excha e

TABLE I .
Comparison of calculated and experimental quantities.