---
author:
- 'C. E. Ekuma'
- 'M. Jarrell'
- 'J. Moreno'
- 'D. Bagayoko'
title: |
    \[Supplementary Information\]\
    Re-examining the electronic structure of germanium: A first-principle study
---

DFT and Progress in band gap in materials
=========================================

Despite the great progress made possible by density functional theory (DFT), from 1964 to present, problems associated with obtaining theoretically the measured energy or band gaps, for finite and crystalline semiconductors, respectively, have persisted. Specifically, most DFT calculations, with emphasis on those utilizing local density approximation (LDA) and semi-local potentials, have led to semiconductor band gaps that are 30 – 50% smaller than their corresponding, measured values. Much effort has been deployed to find explanations of and remedies to this recalcitrant band gap problem. Perdew and Zunger [@PhysRevB.23.5048] introduced the self interaction correction (SIC) to local spin density (LSD) approximation calculations. While the exact functional for the ground state is self interaction free, these authors discussed corrections that appear to be needed for the description of finite systems, beginning with atoms, and of localized states in solids. This self interaction is argued to contribute to the underestimation of the band gaps of insulators by DFT calculations [@PhysRevB.23.5048]. Consequently, self interaction corrections (SIC) are expected to improve the agreement between calculated band gaps and measured ones, in addition to improving binding energies and bringing orbital energies closer to removal energies [@PhysRevB.23.5048; @PhysRevB.26.5445]. While self interaction corrections have led to some improvements in band gap calculations, they have not totally resolved the problem. Applications of SIC have mostly overestimated the band gap of semiconductors [@PhysRevB.52.R14316; @PhysRevB.54.5495]. According to Cohen [@Cohen2008], self interaction is well-defined only for one-electron systems.

According to the literature, a major source of the theoretical underestimation of band gaps consists of the derivative discontinuity of the exchange correlation energy, Exc [@PhysRevLett.51.1888; @PhysRevLett.49.1691; @PhysRevLett.51.1884; @PhysRevB.32.3883]. Perdew [@PhysRevLett.49.1691], following a thought experiment on a diatomic molecule, established the existence of a derivative discontinuity of the exchange correlation energy, i.e., a discontinuity in the exchange correlation potential, V$_{xc}$. Perdew and Levy [@PhysRevLett.51.1884] generalized this discontinuity to the case of semiconductors. They showed that the exchange correlation potential may jump by the discontinuity, $\Delta_{xc}$, when the number of electrons in the system under study increases by one. Band gaps calculated with a local density approximation (LDA) potentials, according to their findings, are to be augmented by this discontinuity in order to reproduce the corresponding, measured values. The authors suggested, without claiming to have a proof of it, that this discontinuity is a non zero (and positive) in real semiconductors and insulators. Sham and Schlüter [@PhysRevLett.51.1888] also found a derivative discontinuity of E$_{xc}$ in insulators. These authors, however, asserted that their work does not show whether or not this discontinuity is non zero in real insulators. Subsequent work by Sham and Schlüter [@PhysRevB.32.3883] derived the discontinuity of the functional derivative of E$_{xc}$ in insulators by considering an increase of the number of electrons by one. Cautiously, these authors concluded that the discrepancy between calculated and measured band gap is a measure of the discontinuity $\Delta_{xc}$ - given the results from several calculations – if the employed LDA potentials are assumed to be good approximations. The description of our method below indicates the strong possibility that some current LDA and GGA potentials may be very good approximations.

Despite its popular use to explain the disagreement between calculated and measured band gaps, the above discontinuity has not yet been established to be non zero in real semiconductors or insulators. Further, Sham and Schlüter [@PhysRevB.32.3883] underscored the fact that, in principle, DFT and Kohn Sham LDA hold only if the number of particle is kept constant. The question could arise whether or not the discontinuity, derived by considering a change of the number of particle, is strictly applicable to DFT or LDA calculations. From the preceding, it has not yet been established that DFT or LDA calculations cannot obtain the correct band gaps, despite the fact that presently known LDA potentials do not have a discontinuity and that most of the numerous, previous ab-initio DFT and LDA calculations did not.

Another presumed contributor to the band gap underestimation by theory stems from the use of local (LDA) and semi-local (GGA) potentials. The question naturally arises as to what extent the local and semi-local potentials fail to capture key feature of the exact one. We are aware of no definitive answer, given that the exact one is not known. We would have had to delve into this matter further if we were dealing with molecules or their dissociation. The solid state systems of interest to us, to judge by previous results obtained with our method [@Bagayoko2004; @Bagayoko2008], possible errors due to the use of local and semi-local potentials appear to be very small.

There exist several approaches that have been introduced to address the band gap problem. Review articles and books are the best sources for discussing these approaches and for examples of the many DFT calculations that led to band gaps much smaller than their corresponding, experimental counterparts. In contrast, a summary of results from BZW LDA calculations for over 10 materials show agreement between theory and experiment. Illustrative examples of discrepancies between theory and experiment follow. The case of Ge is summarized in this article. Some previous LDA, GGA, and GW calculations did not yield the measured band gap, from first principle. A table provided by Ekuma and Bagayoko [@Ekumab2011] shows a multitude of DFT calculations with vastly different band gaps for titanium dioxide. With the computational method described here, Ekuma and Bagayoko obtained the measured, direct gap and predicted an indirect one. For elemental silicon, Grüning [@Myrta2006] reported an LDA band gap of 0.7 eV, much smaller than the 1.25 eV they reported as the measured value. These authors also performed calculations with the exact exchange (EXX), EXX plus LDA, EXX plus the random phase approximation (RPA). The last approach or scheme yielded 0.6 eV, a gap smaller than the above LDA gap, while the first two led to 1.5 and 1.6 eV respectively, values much larger than the experimental one. With the original version of our method, Zhao et al. utilized an LDA potential to obtain a gap of 1.02 for Si, much closer to the experimental one. Generalized gradient approximation (GGA) calculations have led to improvements of calculated properties of materials, including lattice parameters. Specifically, Hao [@PhysRevB.85.014111] reported revised Tao-Perdew-Staroverov-Scuseria (revTPSS) meta-GGA calculated lattice parameters that are in agreement with experiment following a zero-point phonon correction, for over 50 materials. Despite this very significant success, most GGA and meta-GGA calculations, including the previous ones discussed here for Ge, have not produced band gaps in agreement with experiment.

From the above summary, historical overview of the band gap problem, it appears that the scientific community believes that the derivative discontinuity of the exchange correlation energy is the main source of the disagreement between DFT calculated energy and band gaps and their corresponding, measured ones. This belief led to the development of several schemes aimed at resolving the band gap problem. Except for the few, most of these schemes are ad hoc as they include adjustable parameters that vary with the material under study. The continuing growth in the number of these schemes seems be a problem in itself, the ad hoc nature of most of them does not lend itself to predictive capabilities from first principle, the aim of theory to inform and to guide experiment. The only exception to the above trend consists of the work of our group. This work has not yet gotten the attention of the community at large, presumably due to the strength of the above belief, on the one hand, and the preponderance of results that are explained with the discontinuity, on the other hand. As we previously noted [@Bagayoko2005], the situation resembles that of the Ptolemaic model of the solar system where epicycles were continually introduced to explain its disagreement with observations. The quintessential point in support of the our method, described below, is the following: For all DFT calculations of energy bands, the *minima of the occupied energies, which add up to yield the ground state energy of the electron system, are obtained from the theory if the “correct” ground state charge density is utilized, subject to the constraint that the number of particle is kept constant* [@PhysRev.136.B864; @PhysRev.140.A1133].

Most of the previous DFT calculations, including those with GGA and LDA potentials, have consisted of judicious selecting large basis set and of performing iterations to obtain self consistent eigenvalues of the Kohn-Sham type equation. It is assumed that the single basis set in question leads to the correct representation of the electronic cloud in the system under study, a system that can be drastically different from an atomic or ionic one. *In particular, as we recently pointed out, polarization ($p$, $d$, and $f$ orbitals) has primacy over spherical symmetry ($s$ orbital) for systems varying from diatomic molecules to solids.* Hence, utilizing basis sets derived from calculations of properties of atoms for the study of solids is potentially problematic. Indeed, the angular symmetries in these systems are vastly different from those for atoms and ions. The need for the method described below becomes apparent with the realization that, irrespective of the degree of convergence of the iterations, a single trial basis set that has a major symmetry inadequacy for the description of the system is not going to lead to physically valid DFT eigenvalues as the implacable condition of using the “correct” ground state density will not be met.

Description of the Bagayoko-Zhao-Williams-Ekuma-Franklin Method
===============================================================

The original version of our method, named after Bagayoko, Zhao, and Williams (BZW) was introduced in 1998 [@Bagayoko1998] and further explained in 1999 [@Zhao1999]. The method consists of implementing the linear combination of atomic orbitals (LCAO) formalism by methodically searching for the smallest basis set, called the optimal basis set, that leads to the minima of the occupied energies. This search begins with a deliberately small basis set that is not smaller than the minimum basis set, i.e., the smallest one needed to account for all the electrons in the system. The self consistent calculation with this basis set is followed by another whose basis set uses the previous one plus one additional orbital. The occupied energies from the two calculations are compared numerically and graphically. In the more than 20 systems we have studied, these occupied energies from these first two calculations have been different, with those of calculation II generally lower than the corresponding one from Calculation I. Calculation III is then carried out, using the basis set in Calculation II plus an additional orbital. The occupied, self consistent energies from Calculations II and III are also compared. This process of augmenting a basis set, performing new self-consistent calculations, and comparing its results with those of the one immediately preceding it, continues until a calculation is found, say N, to have exactly the same occupied energies as the one immediately following it. This perfect superposition establishes that the minima of the occupied energies have been reached and that the corresponding basis set give the best representation of the ground state charge density of the system.

Before elaborating further on the physical content of the method, we note that adding an orbital means increasing the size of the basis set (and hence the dimension of the Hamiltonian) by 2, 6, 10, or 14 depending on the $s$, $p$, $d$, or $f$ character, respectively, of the orbital in question. In the original BZW, we added orbitals in the order of their energies resulting from the atomic calculations, i.e., the orbitals corresponding to the lowest laying excited, atomic state were successively added. As apparent from our earlier results, the BZW method practically led to the measured band gap of semiconductors we studied. Further, our predictions of the band gaps and other properties for cubic Si$_\mathrm{3}$N$_\mathrm{4}$ [@Bagayoko2001] and cubic InN [@Bagayoko2004] were totally confirmed by experiment [@PhysRevB.65.161202; @Egdell2003; @Sch2006]. Following the works of Ekuma and Franklin [@Ekuma2011; @ekuma:012189; @Ekumab2011; @ibid2011], we realized that valence electrons in multi-atomic systems simply do not follow the symmetry landscape that prevail for isolated atoms or ions. The aim is to obtain a better representation of the electronic cloud (ground state charge density) of the system under study. Hence, in most of our subsequent calculations, for a given principal quantum number, we add $p$, $d$, and $f$ (if applicable) orbitals before the $s$ orbital for that principal quantum number. This counter-intuitive ordering, for isolated atoms, is simply needed for ‘some’ multi-atomic systems. The initials of Ekuma and Franklin (EF) are added to the name of the enhanced method \[BZW-EF\] in recognition of their extensive calculations whose results led Bagayoko to see the necessity for this new order. While the BZW method led to band gaps that were sometimes smaller by 0.1 – 0.3 eV, insignificant differences between BZW-EF calculated gaps and corresponding experiment ones are in the second decimal place for the systems studied to date. For ZnO [@Franklin2013729], the BZW-EF method led to an upper valence band width more than 1.0 eV larger than was obtained with the BZW, with a significant improvement in agreement with experiment.

The origin of the changes in the band structure and the band gap when the basis set increases toward the optimal one consists of the progressively better representation of the ground state charge density. As per the derivation of DFT, the minima of the occupied energies are obtained if the “correct” charge density for the ground state is employed. These changes are due to physical interactions, given that the charge density, the potential and hence the Hamiltonian change from one calculation to the next. We should underscore here that while our focus is on occupied energies (i.e., DFT is a fundamentally ground state theory), when these energies reach their minima, so do the low laying unoccupied energies, up to 9 – 10 eV for the materials studied to date with the BZW-EF method. For the many systems with the BZW, most low laying unoccupied energies also converged up to 5 – 6 eV, as was the case for wurtzite indium nitride [@Jin2007]. For a few materials, this convergence of the lowest unoccupied energies was not achieved with the original BZW method. For metals, as shown by Zhao [@Zhao1999], the low-laying unoccupied energies converge when the occupied ones do, due to the fact that at least one band crosses the Fermi level. This fact may partly explain the early successes of DFT in describing metals as compared to semiconductors.

The description of what occurs when basis sets much larger than the optimal ones are employed will complete the description of the BZW-EF method. We first recall that the basis set immediately following the optimal one leads to the occupied energies obtained with the optimal one and to the same unoccupied energies up to 9 – 10 eV. So, by much larger basis sets, we mean the ones that are larger than the basis sets immediately following the optimal ones. Earlier works by Bagayoko and Co-workers [@Zhao1999; @Bagayoko1998] verified that basis sets larger than the optimal one do not change the charge density, the potential, and the Hamiltonian, nor do they change the occupied energies. In the absence of changes in the Hamiltonian, i.e., the physics of the study, the additional lowering of unoccupied energies with these much larger basis set cannot be ascribed to DFT. However, the Rayleigh theorem provides an explanation of the unphysical lowering of these energies. The theorem states that when the same eigenvalue equation is solved with two basis sets of different sizes, such that the larger one includes the smaller one, then the eigenvalue obtained with the larger basis set are lower than or equal to their corresponding one obtained with the smaller basis set.

Clearly, after the optimal basis set is reached, and that the occupied energies are no longer changed from their values obtained with the optimal one, the lowering of unoccupied energies can be ascribed to a mathematical artifact that is the manifestation of the above theorem. We therefore contend that this extra lowering, a variational basis set effect [@Zhao1999; @Bagayoko1998], is a major source of discrepancies between many previous DFT calculations and between these calculations and experiment, as far as band gaps are concerned. *The lowest laying conduction bands, with full physical meaning when they result from the optimal basis set, partly lose their physical content due to the above effect.* It is important to note that this is the case for any basis set that is not the optimal one, whether it is smaller or larger than the optimal, or simply lacks orbital or angular features of the optimal one. For basis sets that do not totally include their corresponding optimal ones, even the occupied energies are not totally DFT results. The preceding lines in this paragraph point to the great difficulty in obtaining physically meaning DFT occupied and low laying unoccupied energies with a single trial basis set, irrespective of the degree of convergence of the applicable iterations for self consistency.

[28]{} natexlab\#1[\#1]{}bibnamefont \#1[\#1]{}bibfnamefont \#1[\#1]{}citenamefont \#1[\#1]{}url \#1[`#1`]{}urlprefix\[2\][\#2]{} \[2\]\[\][[\#2](#2)]{}

, ****, ().

, ****, ().

, , , ****, ().

, , , ****, ().

, , , ****, ().

, , , , ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , , ****, ().

, , , ****, ().

, ****, ().

, , , ****, ().

, , , , , , ****, ().

, ****, ().

, ****, ().

, ****, ().

, , , , ****, ().

, , , ****, ().

, ****, ().

, , , , , , , , , , , ****, ().

, , , , ****, ().

, , , , , , , , ****, ().

, , , , , ****, ().

, , , , ****, ().

, , , , , ****, ().

, , , , ****, ().

, , , ****, ().
