Slater half-occupation technique revisited: the LDA-1/2 and GGA-1/2 approaches for atomic ionization energies and band gaps in semiconductors

The very old half-occupation technique^1–3 and its formalization through the Janak theorem⁴ have been continuously used since the seventies. Even now we find many papers citing those authors and using their techniques in many computational problems.^5–13 The Janak theorem refers to the derivative of the total energy with respect to the occupation, which equals the Kohn-Sham (KS) eigenvalue.¹⁴ If the derivative is linear with the occupation then the half-occupation technique follows. Part of the recurring publications on the subject refers to verifying the linearity,¹⁵ and this since the early days of the technique.¹⁶ The first applications of the half-occupation technique were for atoms, thus we revisit this problem, now using the modern exchange-correlation approximations.^17,18 In this respect, our work runs parallel to that of Kraisler et al.,¹¹ among others, though they did not use the half-occupation technique, which is a matter of necessity in our case, as we will see in the case of extended systems (crystals).

After atoms, molecules and clusters, the half-occupation technique must face the problem of extended systems such as crystals. In this respect, aside from our work, being reviewed, we only know the work of Liberman.¹⁹ Of course, removing 1/2 electron from an infinite system is no perturbation, so that a special handling of the half-occupation technique is required. Related to this special technique is the definition of “self-energy”, or the “self-interaction”, known to be missing in the Kohn-Sham eigenvalues.²⁰ The self-interaction correction (SIC) also has a long history. It was first defined by Perdew and Zunger¹⁷ and frequently proposed as the difference between the ionization energy (IE) and the KS eigenvalue.^20,21 It happens that whenever the Janak theorem is valid, namely when the energy functional can be differentiated with respect to the occupations, one can prove that the difference between the IE and the KS eigenvalue, is about 1/2 the SIC defined by Ref. 17. This factor of 1/2, apparently first discovered by Trickey,²² has also been numerically verified by Refs. 23. To avoid confusion we are naming the difference between the IE and the KS eigenvalue, following from Janak's theorem, as “self-energy” because its expression is indeed that of a self-energy.

The technique was proposed by Slater's group in the early seventies.^1–3 At that time the density-functional (DF) approximations were unknown and researchers were using the X − α exchange to obtain the best results available at the time. The half-ion technique gained its formalization much latter with Janak⁴ which proved

In this Eq. f_α is the occupation of the Kohn-Sham (KS) eigenstate,¹⁴ not necessarily an integer, e_α is the KS eigenvalue, dependent on the occupation f_α, and E(f_α) is the total energy of the atom (or ion). Since the seminal works of Slater and Janak, the half-ion technique is being intensively used to this day. The proof of Eq. (1) is certainly exact for the local-density approximation (LDA) to DF, and it is likely to be also valid for the generalized-gradient approximation (GGA), though GGA may come in very different formulations.

The half-ion technique begins by integrating Eq. (1) from the ion to the atom

The last equality assumes that the KS eigenvalue is linear with the occupation. The linearity was verified very early in the development of the half-ion technique^15,16 and is now being confirmed by a very large set of atoms. Of course the half-ion technique cannot be as precise as the difference of total energies, yet its precision is typically better than 0.1eV and, sometimes it is an important tool of research.^24,25 In the case of extended systems, as crystals, it is unavoidable. Instead of the approximation of Eq. (2), ref. 13 among other authors, uses

which cannot be interpreted as the energy of a transition state and, as Numerical Analysis says, it is a slightly poorer approximation.

Apparently the authors of the half-ion technique failed to give an explanation of why it works and the physical meaning for the difference between the half-ion eigenvalue e_α( − 1/2) and the eigenvalue e_α(0) at full occupation. An understanding may be reached starting from

where

The first term of S_α is typically the self-energy of the charge distribution n_α. In many instances this first term will be the most important and that is why we are naming S_α as self-energy. A numerical comparison between the first term and the whole self-energy is in the footnote.²⁶ The first two terms of Eq. (5) coincide with Trickey's post-hoc correction,²² except for the unit systems being used. The self-energy is always positive and it is larger if the number-density n_α is denser (core states) and smaller if it is diffuse (valence states). In deriving Eq. (3) to obtain Eq. (5) we apparently miss the derivative of n_α with respect to f_α, except when it is contained in the total density n. The reason is that the eigenvalue is stationary with respect to variations in the normalized wavefunctions ψ_α, thus the derivative is zero (Hellman-Feynman theorem).²⁷ It must be stressed that the self-energy will never be calculated with Eq. (5) but in a much simpler way, even simpler than by using any known self-interaction expression.^17,22,28 The simplicity makes the method very useful.

Integrating Eq. (4) obtain

the last equality again depending on linearity, which again may be assumed because it is usually verified. Truly, the self-energy tends to be independent of the occupations.

Now inserting Eq. (6) into Eq. (2) obtain

Meaning that the ionization of electron α requires not only the negative of the eigenvalue energy but the self-energy as well.

We calculated many ionization energies (IE) with the half-occupation technique, using both the LDA approximation of Perdew and Zunger¹⁷ to the results of Ceperley and Alder²⁹ and the GGA of Perdew, Burke, and Ernzerhof.¹⁸ Our calculations were all with spin polarization,³⁰ thus for the heavier atoms we had to make a perturbation correction to include relativistic effects.

Our results are compared with experiment and with calculated results by other methods. The NIST site³¹ presents very large tables of IE data. Of course we cannot compare our results with coupled-clusters. Coupled-cluster is known to reach a precision of sub-meV but only for very light atoms,³² for heavier atoms their precision is not better than ours.³¹ We are mainly interested in extended systems, specially in semiconductors, where heavier atoms are important. Thus we are comparing our results calculated with spin-polarization, radial grid integration, LDA and GGA half-ion technique with basis functions results, Hartree-Fock (HF) and Density-Functional (DF) results presented in the NIST site. In the NIST site we chose the results for the largest basis (aug-cc-pVQZ, when available), and for the honorable B3LYP functional. The comparison is in Fig. 1. As expected, only HF has important deviations from experimental values. The results presented in Ref. 11 also seem to coincide with ours and NIST’s, barring HF, though a detailed comparison is not possible because those authors do not tabulate their results. Ref. 11 also obtains the ionization energies as the difference between the total energies of atom and ion, not from half-occupation. Looking more closely, at the rms of the difference between calculated and experimental values, our GGA-1/2¹⁸ calculation is outstanding. Its superiority with respect to LDA-1/2¹⁷ only proves that, in searching for better exchange- correlation approximations, scientists have not been idle. The apparent superiority of GGA-1/2 with respect to B3LYP might be a problem of method (numerical against bases), not necessarily a question of quality of the approximation.

One important advantage of the half-ion technique is its ability to calculate heavier atoms and higher order IE. To show this ability we plot in Fig. 2 the first and second IE's for the atoms tabulated in reference 31 and some heavier atoms. In the calculation of the second IE one must delete 3/2 electrons from the atom. The inclusion of the heavier atoms increased the rms error to some extent. The second IE has higher rms error than the first. In the case of Cs and Tl we had to apply a perturbation scheme to include the relativistic corrections into the spin-polarized calculation.

Following Eq. (6), the IE, which equals minus the eigenvalue at half occupation, is obtained if we subtract the self-energy S_α from the eigenvalue at full occupation. The self-energy is the quantum-mechanical average

of the “self-energy potential” V_S defined as follows

The self-energy potential keeps a simple relationship with the KS potential for the half ion and the atom. Indeed consider

Assuming linearity and comparing with Eq, (9) obtain

Therefore, instead obtaining IE from the eigenvalue of the half ion, we may as well subtract the self-energy potential from the potential of the atom and run a single non self-consistent iteration to determine the eigenvalue at half occupation. After the subtraction the potential is V( − 1/2, r), no longer V(0, r), and the eigenvalue of a single iteration is e_α( − 1/2). This seems an undue complication in the case of atoms, but it is important in the case of many-atom systems. In those cases we will want to use the self-energy potential defined in the atom and subtract it from the potentials of the atoms (of the same species) composing the system.

First we recall that the KS eigenstates in a crystal are Bloch states whose wavefunctions extend through the infinite crystal. As such, the Bloch state has an exactly null self-energy. Now, when removing or exciting an electron from an occupied band, say the topmost valence band, the excitation or hole does behave as being in a localized wavefunction, therefore with non null self-energy. In fact, it would be impossible to remove an electron infinitely delocalized in any finite time. Further, theories of excitons, potential barriers, etc, would be in trouble if holes were Bloch states. Though these arguments seem to be very simple, there is no calculation of a localized hole in a solid that we know of, specially in the case when such localized hole could scatter into the continuum of band states, or, in other words, when the localized hole is degenerate with the band continuum.³³ 

Thus, to proceed with the half-occupation for crystals, we have to assume that holes in a valence band, and electrons in a conduction band are localized with non-null self-energy. Generalizing Eqs. (2) and (6), to describe the excitation of an electron from the valence to the conduction band, we write

This last equation is frequently writen as

where Δ_xc is known as the “derivative discontinuity” that is caused by an exchange-correlation beyond LDA or GGA.^34–36 The half-occupation technique allows an entirely different, and perhaps simpler, interpretation of this discontinuity: It is a self-energy caused by the localization of the particle excitation and it is not the effect of exotic exchange-correlation functionals.

Since we do not know how to calculate localized holes and electrons in solids, which are not even stationary states, we assume that their self-energies are again quantum mechanical averages of the atomic self-energy potential V_S. In a crystal such as NaCl we would add V_S(Cl) to the chlorine atoms and V_S(Na) to the sodiums. We add to all atoms of the crystal because the calculation with just one or a few perturbed atoms is much too complicate.

There is a problem with adding V_S to all the atoms of an infinite crystal: the potential will diverge. A V_S is a potential of an excess charge of 1/2 proton (see Eq. (9)) and has a tail of 0.5/r that cannot be summed in an infinite lattice. Therefore the tail has to be trimmed by a step function, and we have been using

Thus we are making

In this Eq., the exponent m is as large as possible and it is limited only by the band-calculation code. We have been using always m = 8. Thus we are introducing a single parameter in the calculation, CUT, which is determined variationally, as explained in ref. 37.

Thus our recipe, named LDA-1/2, to determine band gaps, consists of adding the atomic self-energy potential to the atoms of the crystal. Now, in all our calculations with crystals, and we made many, we never found an instance when the self-energy of the electron in the conduction band had any importance in Eq. (12). The valence band is usually made of the valence state of the anion, while the conduction band is usually a mixture of many atomic states, there included the valence state of the cation. Then, usually the self-energy of the cation could be neglected, or, which is the same, the extremizing solution usually was CUT_cation = 0. In this respect we found the calculation of the work functions of Al, Cu, Na very symptomatic. A work function is a band gap between the Fermi level at the conduction band and the free electron band at infinity. In those cases we always found CUT = 0, that is the LDA-1/2 method reducing to standard LDA.

Fig. 3 shows the situation of the exceptional case ZnO. This material is exceptional in two ways: first the pure LDA or GGA calculations give a band gap smaller than 1.0 eV, while our LDA-1/2 is able to bring the gap to a value close to the experimental 3.4 eV. A major shift and a major success of our technique. Secondly, it is exceptional because the self-energy of the cation Zn cannot be neglected. Fig. 3 explains how the calculation was made. We began by including the self-energy potential for the anion

and by maximizing the gap with respect to CUT_O. Then with the maximizing CUT_O fixed, we included the self-energy potential for Zn

A new maximization was made, resulting in a gap near the experimental value.

We have been naming our method³⁷ of calculating band gaps as the LDA-1/2, though sometimes we use the Generalized Gradient Approximation, for instance of Ref. 18, instead of the Local-Density Approximation of Refs. 17 or.38 The NSCF version is that explained in the last paragraph of section III. In other words,

1. The (model) number-density is chosen that of the ground state unperturbed by V_S.

2. For a single interaction we choose the potential as the sum of the ground state potential plus the CUT-trimmed self-energy potential.

3. Solve the Schroedinger equation.

4. Then CUT is chosen so as to make the band gap extreme.

The procedure works very well for highly ionic crystals, as shown in Fig 4. There we compare calculated band gaps with measured numbers. The pink balls are calculations according to the present procedure. They fall very near the straight line of perfection. Simple LDA, without any correction, gives band gaps smaller, as shown by the blue crosses. Then, using the procedure of the next subsection, that is the SCF LDA-1/2 we obtain the red squares that are half way between the pure LDA results and the NSCF LDA-1/2 results.

This is the original LDA-1/2 method.³⁷ The atomic self-energy potential, trimmed by the function Θ, Eq. (14), is added to the atoms of the crystal, as an external potential. Then self-consistent iterations are run. The parameter CUT is chosen according to make the gap extreme. Since the 2008 paper, many new semiconductors were calculated. The new plot of Calculated Band Gaps versus Experimental Band Gaps (Fig. 5) is very impressive, probably the best set of calculated band gaps available in the literature.

It is very difficult to calculate the self-energy in the crystal itself and we have been recurring to atomic physics to calculate an approximation. We found that for most semiconductors the self-energy potential of the atom was transferable to crystals after trimming its 0.5/r tail. For the highly ionic crystal we found that no self-consistent loops could be made. For the crystals in the Fig. 6 we are finding that neither procedure, self-consistent or non self-consistent is very good. In the case of these oxides we found that a working procedure was to calculate the self-energy potential of Oxygen as

instead of Eq. (11) for the neutral atoms.

We are paying a price for using an atomic self-energy potential instead of a self-energy potential defined in the crystal: we do not know a priori whether to run self-consistent iterations or NSCF. On the other hand, it is easy to verify that all covalent semiconductors and insulators require full SCF calculations, while the highly ionic insulators require treating the self-energy potential in a single NSCF iteration. Further, the localized valence band hole energy, e_v − S_v, falls within the band continuum for those cases when we make the standard SCF LDA-1/2, and falls outside for the cases when we use NSCF. Thus there is a simple explanation for the difference between the two procedures: when the hole energy falls within the band continuum the atomic self-energy potential is screened by the other band electrons and the self-consistent iterations must be completed to account for that screening. Contrarily, when the hole energy falls outside the band, no screening by the electrons of the other atoms must be considered and the situation is all similar to that of atoms.

In this paper we show the power of Slater's half-ion technique (transition state technique). Aside from its known usefulness to finite systems, such as atoms, molecules, and cluster, we now extend its use to crystals and semiconductors. The extension is attained by taking a second derivative of the energy functional, one step beyond Janak's first derivative, which allowed us to define a self-energy, means to calculate ionization energies and band gaps. The set of band gaps we are presenting have a GW precision,⁴² though it is obtained with minimal computational effort.

This paper resulted from our experience with the half-occupation technique after its proposal in Ref. [37]. In the case of atoms we compared our ionization energies with those of other methods and verified that we are performing better, except for the coupled cluster results of the first line of the Periodic Table. In the case of crystals we now included highly ionic insulators, for which our technique had to be modified. We are also presenting band gap results for many new covalent semiconductors. And to be complete, we found some solids (oxides) for which the technique did not perform very well.

The authors are thankful to the Brazilian Funding Agencies FAPESP and CNPq for continuous support. The authors are indebted to Profs. Marcio H. F. Bettega, Maria C. Santos and Mauro F. S. Ribeiro Jr. for their suggestions.