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

Ferreira, Luiz G.; Marques, Marcelo; Teles, Lara K.

doi:10.1063/1.3624562

The very old and successful density-functional technique of half-occupation is revisited [J. C. Slater, Adv. Quant. Chem. 6, 1 (1972)]. We use it together with the modern exchange-correlation approximations to calculate atomic ionization energies and band gaps in semiconductors [L. G. Ferreira et al., Phys. Rev. B 78, 125116 (2008)]. Here we enlarge the results of the previous paper, add to its understandability, and show when the technique might fail. Even in this latter circumstance, the calculated band gaps are far better than those of simple LDA or GGA. As before, the difference between the Kohn-Sham ground state one-particle eigenvalues and the half-occupation eigenvalues is simply interpreted as the self-energy (not self-interaction) of the particle excitation. In both cases, that of atomic ionization energies and semiconductor band gaps, the technique is proven to be very worthy, because not only the results can be very precise but the calculations are fast and very simple.

I. INTRODUCTION

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.

II. THE HALF-ION TECHNIQUE

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

\begin{equation}\frac{\partial E}{\partial f_{\alpha }}=e_{\alpha }(f_{\alpha })\end{equation}

\frac{\partial E}{\partial f_{α}} = e_{α} (f_{α})

(1)

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

\begin{equation}E_{\alpha }(0)-E_{\alpha }(-1)=\int _{-1}^0\frac{\partial E}{\partial f_{\alpha }} df_{\alpha }=\int _{-1}^0e_{\alpha }(f_{\alpha })df_{\alpha }=e_{\alpha }(-1/2)\end{equation}

E_{α} (0) - E_{α} (- 1) = \int_{- 1}^{0} \frac{\partial E}{\partial f_{α}} d f_{α} = \int_{- 1}^{0} e_{α} (f_{α}) d f_{α} = e_{α} (- 1 / 2)

(2)

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

\[\frac{e_{\alpha }(0)+e_{\alpha }(-1)}{2}\]

\frac{e_{α} (0) + e_{α} (- 1)}{2}

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

\begin{equation}e_{\alpha }(f_{\alpha })=KE[n_{\alpha }]+\int V_{nucleus}(\vec{r})n_{\alpha } (\vec{r})d^3r+\int \int d^3rd^3r^{\prime }\frac{n_{\alpha }(\vec{r})n (\vec{r}^{\ \prime })}{\big |\vec{r}-\vec{r}^{\ \prime }\big |}+ \int d^3rn_{\alpha }(\vec{r}) \frac{\delta E_{xc}}{\delta n(\vec{r})}\end{equation}

e_{α} (f_{α}) = K E [n_{α}] + \int V_{n u c l e u s} (\vec{r}) n_{α} (\vec{r}) d^{3} r + \int \int d^{3} r d^{3} r^{'} \frac{n_{α} (\vec{r}) n ({\vec{r}}^{'})}{| \vec{r} - {\vec{r}}^{'} |} + \int d^{3} r n_{α} (\vec{r}) \frac{δ E_{x c}}{δ n (\vec{r})}

(3)

where

$n_{\alpha }=\psi ^*_{\alpha }\psi _{\alpha }$

n_{α} = ψ_{α}^{*} ψ_{α}

is the number density of the KS state α, n = ∑_if_in_i is the number density of the whole electron cloud, and KE[] is the one-electron kinetic energy functional. Taking the derivative of Eq. with respect to f_α obtain

\begin{equation}\frac{\partial e_{\alpha }}{\partial f_{\alpha }}=2S_{\alpha },\end{equation}

\frac{\partial e_{α}}{\partial f_{α}} = 2 S_{α},

(4)

\begin{eqnarray}S_{\alpha }&=&\frac{1}{2} \int \int d^3rd^3r^{\prime }\frac{n_{\alpha }(\vec{r})n_{\alpha } (\vec{r}^{\ \prime })}{\big |\vec{r}-\vec{r}^{\ \prime }\big |}+ \frac{1}{2}\int \int d^3rd^3r^{\prime }n_{\alpha }(\vec{r}) \frac{\delta ^2E_{xc}}{\delta n(\vec{r})\delta n(\vec{r}^{\ \prime })} n_{\alpha }(\vec{r}^{\ \prime })\\&+&\frac{1}{2}\int \int d^3rd^3r^{\prime }\frac{n_{\alpha }(\vec{r})}{\big |\vec{r}-\vec{r}^{\ \prime }\big |}\sum _{\beta }f_{\beta } \frac{\partial n_{\beta } (\vec{r}^{\ \prime })}{\partial f_{\alpha }} +\frac{1}{2}\int \int d^3rd^3r^{\prime }n_{\alpha }(\vec{r}) \frac{\delta ^2E_{xc}}{\delta n(\vec{r})\delta n(\vec{r}^{\ \prime })} \sum _{\beta }f_{\beta }\frac{\partial n_{\beta }(\vec{r}^{\ \prime })}{\partial f_{\alpha }} \nonumber\end{eqnarray}

\begin{matrix} S_{α} & = & \frac{1}{2} \int \int d^{3} r d^{3} r^{'} \frac{n_{α} (\vec{r}) n_{α} ({\vec{r}}^{'})}{| \vec{r} - {\vec{r}}^{'} |} + \frac{1}{2} \int \int d^{3} r d^{3} r^{'} n_{α} (\vec{r}) \frac{δ^{2} E_{x c}}{δ n (\vec{r}) δ n ({\vec{r}}^{'})} n_{α} ({\vec{r}}^{'}) \\ + & \frac{1}{2} \int \int d^{3} r d^{3} r^{'} \frac{n_{α} (\vec{r})}{| \vec{r} - {\vec{r}}^{'} |} \sum_{β} f_{β} \frac{\partial n_{β} ({\vec{r}}^{'})}{\partial f_{α}} + \frac{1}{2} \int \int d^{3} r d^{3} r^{'} n_{α} (\vec{r}) \frac{δ^{2} E_{x c}}{δ n (\vec{r}) δ n ({\vec{r}}^{'})} \sum_{β} f_{β} \frac{\partial n_{β} ({\vec{r}}^{'})}{\partial f_{α}} \end{matrix}

(5)

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

\begin{equation}e_{\alpha }(-1/2)-e_{\alpha }(0)=\int _0^{-1/2}\frac{\partial e_{\alpha }}{\partial f_{\alpha }} =-S_{\alpha }\end{equation}

e_{α} (- 1 / 2) - e_{α} (0) = \int_{0}^{- 1 / 2} \frac{\partial e_{α}}{\partial f_{α}} = - S_{α}

(6)

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

\begin{equation}E(-1)=E(0)-e_{\alpha }(0)+S_{\alpha }\end{equation}

E (- 1) = E (0) - e_{α} (0) + S_{α}

(7)

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.

FIG. 1.

View large Download slide

(color online) Ionization energies: Comparison of calculated results with experiment for the atoms tabulated in the figure. The experimental data, HF, and B3LYP are taken from the NIST site.³¹ The LDA-1/2 uses the exchange-correlation approximation of ref. 17. The GGA-1/2 uses ref. 18.

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.

FIG. 2.

View large Download slide

(color online) Ionization energies: Comparison of calculated results for the first and second IE with experiment for the atoms tabulated in the Fig. The exchange-correlation used was GGA of 18.

III. NON SELF-CONSISTENT IE CALCULATIONS

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

\begin{equation}S_{\alpha }=\int d^3r n_{\alpha }(\vec{r})V_S(\vec{r})\end{equation}

S_{α} = \int d^{3} r n_{α} (\vec{r}) V_{S} (\vec{r})

(8)

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

\begin{eqnarray}V_S(\vec{r})&=&\frac{1}{2}\int d^3r^{\ \prime } \frac{n_{\alpha }(\vec{r}^{\ \prime })}{|\vec{r}- \vec{r}^{\ \prime }|}+\frac{1}{2}\int d^3r^{\ \prime }\frac{\delta ^2 E_{xc}}{\delta n(\vec{r})\delta n(\vec{r}^{\ \prime })}n_{\alpha }(\vec{r}^{\ \prime }) \\&+&\frac{1}{2}\int d^3r^{\ \prime }\frac{\sum _{\beta } f_{\beta }\frac{\partial n_{\beta } (\vec{r}^{\ \prime })}{\partial f_{\alpha }}}{|\vec{r}- \vec{r}^{\ \prime }|}+\frac{1}{2}\int d^3r^{\ \prime }\frac{\delta ^2 E_{xc}}{\delta n(\vec{r})\delta n(\vec{r}^{\ \prime })}\sum _{\beta }f_{\beta } \frac{\partial n_{\beta } (\vec{r}^{\ \prime })}{\partial f_{\alpha }}\nonumber\end{eqnarray}

\begin{matrix} V_{S} (\vec{r}) & = & \frac{1}{2} \int d^{3} r^{'} \frac{n_{α} ({\vec{r}}^{'})}{| \vec{r} - {\vec{r}}^{'} |} + \frac{1}{2} \int d^{3} r^{'} \frac{δ^{2} E_{x c}}{δ n (\vec{r}) δ n ({\vec{r}}^{'})} n_{α} ({\vec{r}}^{'}) \\ + & \frac{1}{2} \int d^{3} r^{'} \frac{\sum_{β} f_{β} \frac{\partial n_{β} ({\vec{r}}^{'})}{\partial f_{α}}}{| \vec{r} - {\vec{r}}^{'} |} + \frac{1}{2} \int d^{3} r^{'} \frac{δ^{2} E_{x c}}{δ n (\vec{r}) δ n ({\vec{r}}^{'})} \sum_{β} f_{β} \frac{\partial n_{β} ({\vec{r}}^{'})}{\partial f_{α}} \end{matrix}

(9)

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

\begin{eqnarray}V(0,r)&-&V(-1/2,r)=\int _{-1/2}^0df_{\alpha }\frac{\partial }{\partial f_{\alpha }} \left\lbrace -\frac{Z}{r}+\int d^3r^{\ \prime }\frac{n(\vec{r}^{\ \prime })}{|\vec{r}- \vec{r}^{\ \prime }|}+\frac{\delta E_{xc}}{\delta n(\vec{r})}\right\rbrace \nonumber \\&=&\int _{-1/2}^0df_{\alpha }\left\lbrace \int d^3r^{\ \prime }\frac{n_{\alpha } (\vec{r}^{\ \prime })}{|\vec{r}- \vec{r}^{\ \prime }|}+\int d^3r^{\ \prime }\frac{\delta ^2 E_{xc}}{\delta n(\vec{r})\delta n(\vec{r}^{\ \prime })}n_{\alpha }(\vec{r}^{\ \prime }) \right\rbrace \\&+&\int _{-1/2}^0df_{\alpha }\left\lbrace \int d^3r^{\ \prime }\frac{\sum _{\beta } f_{\beta }\frac{\partial n_{\beta } (\vec{r}^{\ \prime })}{\partial f_{\alpha }}}{|\vec{r}- \vec{r}^{\ \prime }|}+\int d^3r^{\ \prime }\frac{\delta ^2 E_{xc}}{\delta n(\vec{r})\delta n(\vec{r}^{\ \prime })}\sum _{\beta }f_{\beta } \frac{\partial n_{\beta } (\vec{r}^{\ \prime })}{\partial f_{\alpha }}\right\rbrace \nonumber\end{eqnarray}

\begin{matrix} V (0, r) & - & V (- 1 / 2, r) = \int_{- 1 / 2}^{0} d f_{α} \frac{\partial}{\partial f_{α}} \{- \frac{Z}{r} + \int d^{3} r^{'} \frac{n ({\vec{r}}^{'})}{| \vec{r} - {\vec{r}}^{'} |} + \frac{δ E_{x c}}{δ n (\vec{r})}\} \\ = & \int_{- 1 / 2}^{0} d f_{α} \{\int d^{3} r^{'} \frac{n_{α} ({\vec{r}}^{'})}{| \vec{r} - {\vec{r}}^{'} |} + \int d^{3} r^{'} \frac{δ^{2} E_{x c}}{δ n (\vec{r}) δ n ({\vec{r}}^{'})} n_{α} ({\vec{r}}^{'})\} \\ + & \int_{- 1 / 2}^{0} d f_{α} \{\int d^{3} r^{'} \frac{\sum_{β} f_{β} \frac{\partial n_{β} ({\vec{r}}^{'})}{\partial f_{α}}}{| \vec{r} - {\vec{r}}^{'} |} + \int d^{3} r^{'} \frac{δ^{2} E_{x c}}{δ n (\vec{r}) δ n ({\vec{r}}^{'})} \sum_{β} f_{β} \frac{\partial n_{β} ({\vec{r}}^{'})}{\partial f_{α}}\} \end{matrix}

(10)

Assuming linearity and comparing with Eq, (9) obtain

\begin{equation}V(0,r)-V(-1/2,r)=V_S(r)\end{equation}

V (0, r) - V (- 1 / 2, r) = V_{S} (r)

(11)

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.

IV. EXTENDED SYSTEMS: BAND GAPS IN CRYSTALS

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

\begin{eqnarray}\text{True Band Gap}&=&E_{excited}(val-1,cond+1)-E_{ground}(0,0)= e_c(1/2)-e_v(-1/2)\nonumber \\&=&e_c(0)-e_v(0)-S_c+S_v\nonumber \\&=&\text{Kohn-Sham Band Gap}+S_v-S_c\end{eqnarray}

\begin{matrix} True Band Gap & = & E_{e x c i t e d} (v a l - 1, c o n d + 1) - E_{g r o u n d} (0, 0) = e_{c} (1 / 2) - e_{v} (- 1 / 2) \\ = & e_{c} (0) - e_{v} (0) - S_{c} + S_{v} \\ = & Kohn-Sham Band Gap + S_{v} - S_{c} \end{matrix}

(12)

This last equation is frequently writen as

\begin{equation*}\text{True Band Gap}=\text{Kohn-Sham Band Gap}+\Delta _{xc} \nonumber\end{equation*}

True Band Gap = Kohn-Sham Band Gap + Δ_{x c}

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

\begin{equation}\Theta (r)=\left\lbrace \begin{array}{lr}\left[1-\big (\frac{r}{CUT}\big )^m\right] ^3&r\le CUT\\0&r>CUT\end{array}\right.\end{equation}

Θ (r) = \{\begin{matrix} {[1 - {(\frac{r}{C U T})}^{m}]}^{3} & r \leq C U T \\ 0 & r > C U T \end{matrix})

(13)

Thus we are making

\begin{equation}V_S(r)\longrightarrow V_S^{\prime }(r)=\Theta (r)V_S(r)\end{equation}

V_{S} (r) \to V_{S}^{'} (r) = Θ (r) V_{S} (r)

(14)

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

\begin{equation}V_S^O(r)=V(2s^22p^4,r)-V(2s^22p^{3.5},r)\end{equation}

V_{S}^{O} (r) = V (2 s^{2} 2 p^{4}, r) - V (2 s^{2} 2 p^{3.5}, r)

(15)

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

\begin{equation}V_S^{Zn}(r)=V(3d^{10}4s^2,r)-V(3d^{9.5}4s^2,r)\end{equation}

V_{S}^{Z n} (r) = V (3 d^{10} 4 s^{2}, r) - V (3 d^{9.5} 4 s^{2}, r)

(16)

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

FIG. 3.

View large Download slide

(color online) The figure depicts the choice of CUT for O and Zn in the case of ZnO. First we maximize the gap varying CUT^O, then we vary CUT^Zn to reach a band gap value near the correct value.

A. The non self-consistent (NSCF) LDA-1/2 method

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,

The (model) number-density is chosen that of the ground state unperturbed by V_S.
For a single interaction we choose the potential as the sum of the ground state potential plus the CUT-trimmed self-energy potential.
Solve the Schroedinger equation.
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.

FIG. 4.

View large Download slide

(color online) Comparison of calculated band gaps with experiment. The best results were obtained with the NSCF LDA-1/2 technique (pink balls). The worst were obtained with pure LDA (crosses). The standard LDA-1/2 technique³⁷ (red squares), here being renamed SCF LDA-1/2, gives results between the two. The band structure calculations were made with the codes VASP^39,40 and WIEN2k.⁴¹

B. The self-consistent (SCF) LDA-1/2 method

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.

FIG. 5.

View large Download slide

(color online) Comparison of calculated band gaps with experiment. The red square are the SCF LDA-1/2 (standard LDA-1/2). The crosses are standard LDA. The small gap semiconductors are metals (negative gaps), when calculated with LDA. LDA-1/2 corrects the situation. The band structure calculations were made with the codes VASP^39,40 and WIEN2k.⁴¹

C. Intermediate cases

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

\begin{equation}V(O^-,r)-V(O^{-0.5-},r)=V_S(r)\end{equation}

V (O^{-}, r) - V (O^{- 0.5 -}, r) = V_{S} (r)

(17)

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

FIG. 6.

View large Download slide

(color online) Comparison of calculated band gaps with experiment. The red squares are the SCF LDA-1/2 (standard LDA-1/2). The crosses are standard LDA, with no self-energy potential. The pink balls are NSCF LDA-1/2 calculations. The band structure calculations were made with the codes VASP^39,40 and WIEN2k.⁴¹

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.

V. CONCLUSIONS

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.

ACKNOWLEDGMENTS

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.

REFERENCES

1.

J. C.

Slater

and

K. H.

Johnson

,

Phys. Rev. B

5

,

844

(

1972

).

https://doi.org/10.1103/PhysRevB.5.844

Google Scholar

Crossref

2.

J. C.

Slater

,

Adv. Quantum Chem.

6

,

1

(

1972

).

https://doi.org/10.1016/S0065-3276(08)60541-9

Google Scholar

Crossref

3.

J. C.

Slater

and

J. H.

Wood

,

Int. J. Quant. Chem. Suppl.

4

,

3

(

1971

).

Google Scholar

4.

J. F.

Janak

,

Phys. Rev. B

18

,

7165

(

1978

).

https://doi.org/10.1103/PhysRevB.18.7165

Google Scholar

Crossref

5.

Y.

Imamura

,

R.

Kobayashi

, and

H.

Nakai

,

J. Chem. Phys.

134

,

124113

(

2011

).

https://doi.org/10.1063/1.3569030

Google Scholar

Crossref

PubMed

6.

D. H. E.

abd

E. R.

Johnson

,

X.

Hu

, and

W.

Yang

,

J. Phys. Chem. A

115

,

76

(

2011

).

https://doi.org/10.1021/jp109280y

Google Scholar

Crossref

PubMed

7.

S.

Sanna

,

W. G.

Schmidt

,

T.

Frauenheim

, and

U.

Gerstmann

,

Phys. Rev. B

80

,

104120

(

2009

).

https://doi.org/10.1103/PhysRevB.80.104120

Google Scholar

Crossref

8.

J.-W.

Song

,

M. A.

Watson

, and

K.

Hirao

,

J. Chem. Phys.

131

,

144108

(

2009

).

https://doi.org/10.1063/1.3243819

Google Scholar

Crossref

PubMed

9.

S.

Sanna

,

T.

Frauenheim

, and

U.

Gerstmann

,

Phys. Rev. B

78

,

085201

(

2008

).

https://doi.org/10.1103/PhysRevB.78.085201

Google Scholar

Crossref

10.

T.

Körzdörfer

and

S.

Kümmel

,

Phys. Rev. B

82

,

155206

(

2010

).

https://doi.org/10.1103/PhysRevB.82.155206

Google Scholar

Crossref

11.

E.

Kraisler

,

G.

Makov

, and

I.

Kelson

,

Phys. Rev. A

82

,

042516

(

2010

).

https://doi.org/10.1103/PhysRevA.82.042516

Google Scholar

Crossref

12.

E. R.

Johnson

,

W.

Yang

, and

E. R.

Davidson

,

J. Chem. Phys.

133

,

164107

(

2010

).

https://doi.org/10.1063/1.3497190

Google Scholar

Crossref

PubMed

13.

F.

Gallino

,

G.

Pacchioni

, and

C. D.

Valentin

,

J. Chem. Phys.

133

,

144512

(

2010

).

https://doi.org/10.1063/1.3491271

Google Scholar

Crossref

PubMed

14.

W.

Kohn

and

L. J.

Sham

,

Phys. Rev.

140

,

A1133

(

1965

).

https://doi.org/10.1103/PhysRev.140.A1133

Google Scholar

Crossref

15.

C.

Göransson

,

W.

Olovsson

, and

I. A.

Abrikosov

,

Phys. Rev. B

72

,

134203

(

2005

).

https://doi.org/10.1103/PhysRevB.72.134203

Google Scholar

Crossref

16.

J. R.

Leite

and

L. G.

Ferreira

,

Phys. Rev. A

3

,

1224

(

1971

).

https://doi.org/10.1103/PhysRevA.3.1224

Google Scholar

Crossref

17.

J. P.

Perdew

and

A.

Zunger

,

Phys. Rev. B

23

,

5048

(

1981

).

https://doi.org/10.1103/PhysRevB.23.5048

Google Scholar

Crossref

18.

J. P.

Perdew

,

K.

Burke

, and

M.

Ernzerhof

,

Phys. Rev. Lett.

77

,

3865

(

1996

).

https://doi.org/10.1103/PhysRevLett.77.3865

Google Scholar

Crossref

PubMed

19.

D. A.

Liberman

,

Phys. Rev. B

62

,

6851

(

2000

).

https://doi.org/10.1103/PhysRevB.62.6851

Google Scholar

Crossref

20.

T.

Körzdörfer

,

J. Chem. Phys.

134

,

094111

(

2011

).

https://doi.org/10.1063/1.3556979

Google Scholar

Crossref

PubMed

21.

I.

Dabo

,

A.

Ferretti

,

N.

Poilvert

,

Y.

Li

,

N.

Marzari

, and

M.

Cococcioni

,

Phys. Rev. B

82

,

115121

(

2010

).

https://doi.org/10.1103/PhysRevB.82.115121

Google Scholar

Crossref

22.

S. B.

Trickey

,

Phys. Rev. Letters

56

,

881

(

1986

).

https://doi.org/10.1103/PhysRevLett.56.881

Google Scholar

Crossref

23.

C. D.

Pemmaraju

,

T.

Archer

,

D.

Sánchez-Portal

, and

S.

Sanvito

,

Phys. Rev. B

75

,

045101

(

2007

).

https://doi.org/10.1103/PhysRevB.75.045101

Google Scholar

Crossref

24.

J. M.

García-Lastra

,

P. L.

Cook

,

F. J.

Himpsel

, and

A.

Rubio

,

J. Chem. Phys.

133

,

151103

(

2010

).

https://doi.org/10.1063/1.3497188

Google Scholar

Crossref

PubMed

25.

R. R.

Zope

,

T.

Baruah

,

S. L.

Richardson

,

M. R.

Pederson

, and

B. I.

Dunlap

,

J. Chem. Phys.

133

,

034301

(

2010

).

https://doi.org/10.1063/1.3459056

Google Scholar

Crossref

PubMed

26.

In the table below the Electrostatic self-energy is the first term of Eq. , Self-energy is calculated according to Eq. . Entries are in Rydbergs.

\[\begin{array}{l@{\quad}c@{\quad}c@{\quad}c@{\quad}c@{\quad}c}\displaystyle &{\rm Be} &{\rm Ne} &{\rm Mg} &{\rm Ar} &{\rm Ca}\\[3pt]\displaystyle {\rm Electrostatic\ Self\hbox{-}Energy} &0.349 &0.960 &0.287 &0.547 &0.224.\\[3pt]\displaystyle {\rm Self\hbox{-}energy} &0.252 &0.624 &0.214 &0.402 &0.169\\\end{array}\]

\begin{matrix} Be & Ne & Mg & Ar & Ca \\ Electrostatic Self - Energy & 0.349 & 0.960 & 0.287 & 0.547 & 0.224 . \\ Self - energy & 0.252 & 0.624 & 0.214 & 0.402 & 0.169 \end{matrix}

27.

We rewrite Eq. as

\[e_{\alpha }=\frac{<\psi _{\alpha } H \psi _{\alpha }>}{<\psi _{\alpha } \psi _{\alpha }>}\]

e_{α} = \frac{< ψ_{α} H ψ_{α} >}{< ψ_{α} ψ_{α} >}

where H is the part of the Eq. operating on n_α. Then, as ψ_α is a normalized eigenfuntion of H, the derivative with respect to any of its parameters is null.

28.

T.

Körzdörfer

,

S.

Kümmel

,

N.

Marom

, and

L.

Kronik

,

Phys. Rev. B

79

,

201205

–

R

(

2009

).

https://doi.org/10.1103/PhysRevB.79.201205

Google Scholar

Crossref

29.

D. M.

Ceperley

and

B. J.

Alder

,

Phys. Rev. Lett.

45

,

566

(

1980

).

https://doi.org/10.1103/PhysRevLett.45.566

Google Scholar

Crossref

30.

S.

Froyen

,

N.

Troullier

,

J. L.

Martins

, and

A.

Garcia

, Code ATOM inside the SIESTA package (

2006

).

31.

cccbdb.nist.gov, Calculated data D.6, Ionization Energies.

32.

W.

Klopper

,

R. A.

Bachorz

,

D. P.

Tew

, and

C.

Hättig

,

Phys. Rev. A

81

,

022503

(

2010

).

https://doi.org/10.1103/PhysRevA.81.022503

Google Scholar

Crossref

33.

The reader should not confuse localized quasi-stationary states in crystals, such as the holes we are talking about, with localized basis functions as the Wannier functions.

34.

J. P.

Perdew

and

M.

Levy

,

Phys. Rev. Lett.

51

,

1884

(

1983

).

https://doi.org/10.1103/PhysRevLett.51.1884

Google Scholar

Crossref

35.

L. J.

Sham

and

M.

Schlüter

,

Phys. Rev. Lett.

51

,

1888

(

1983

).

https://doi.org/10.1103/PhysRevLett.51.1888

Google Scholar

Crossref

36.

N.

Helbig

,

N. N.

Lathiotakis

, and

E. K. U.

Gross

,

Phys. Rev. A

79

,

022504

(

2009

).

https://doi.org/10.1103/PhysRevA.79.022504

Google Scholar

Crossref

37.

L. G.

Ferreira

,

M.

Marques

, and

L. K.

Teles

,

Phys. Rev. B

78

,

125116

(

2008

).

https://doi.org/10.1103/PhysRevB.78.125116

Google Scholar

Crossref

38.

J. P.

Perdew

and

Y.

Wang

,

Phys. Rev. B

45

,

13244

(

1992

).

https://doi.org/10.1103/PhysRevB.45.13244

Google Scholar

Crossref

39.

G.

Kresse

and

J.

Hafner

,

Phys. Rev. B

47

,

R558

(

1993

).

https://doi.org/10.1103/PhysRevB.47.558

Google Scholar

Crossref

40.

G.

Kresse

and

J.

Furthmüller

,

Comput. Mater. Sci.

6

,

15

(

1996

).

https://doi.org/10.1016/0927-0256(96)00008-0

Google Scholar

Crossref

41.

P.

Blaha

,

K.

Schwarz

,

P.

Sorantin

, and

S. B.

Trickey

,

Comput. Phys. Commun.

59

,

399

(

1990

), see www.wien2k.at.

https://doi.org/10.1016/0010-4655(90)90187-6

Google Scholar

Crossref

42.

G.

Onida

,

L.

Reining

, and

A.

Rubio

,

Rev. Mod. Phys.

74

,

601

(

2002

).

https://doi.org/10.1103/RevModPhys.74.601

Google Scholar

Crossref

2011

Author(s)

This article is distributed under a Creative Commons Attribution 3.0 Unported License.

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

I. INTRODUCTION

II. THE HALF-ION TECHNIQUE

III. NON SELF-CONSISTENT IE CALCULATIONS

IV. EXTENDED SYSTEMS: BAND GAPS IN CRYSTALS

A. The non self-consistent (NSCF) LDA-1/2 method

B. The self-consistent (SCF) LDA-1/2 method

C. Intermediate cases

V. CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

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

I. INTRODUCTION

II. THE HALF-ION TECHNIQUE

III. NON SELF-CONSISTENT IE CALCULATIONS

IV. EXTENDED SYSTEMS: BAND GAPS IN CRYSTALS

A. The non self-consistent (NSCF) LDA-1/2 method

B. The self-consistent (SCF) LDA-1/2 method

C. Intermediate cases

V. CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

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