Assessing the performance of self-consistent hybrid functional for band gap calculation in oxide semiconductors

Oxides are one of the most physically complicated and technologically important systems, exhibiting a variety of structures and interesting properties. Metal oxides have broad applications in many relevant areas including superconductivity [1], ferroelectricity [2], dielectric [3], transparent conducting oxides [4], thermoelectricity [5], photocatalysis [6], and photovoltaic materials [7]. One of the fundamental properties of these materials, which is functional for virtually any type of application, is band gaps. Theoretical understanding of oxides and the accurate estimation of the band gap is, however, a very difficult task. This is due to the nature of metal–oxygen bonding which can vary between nearly ionic to highly covalent or metallic [8], and to the generally high degree of electronic correlation, which poses serious problems to conventional computational methods.

With the rapid development of first principles computational approaches to solve the many-body Schrödinger equation for extended systems, in particular density functional theory (DFT), many materials-specific properties have been calculated without relying on any empirical parameters [9]. However, it is still very difficult to calculate the band gap accurately within the DFT framework using standard semilocal exchange-correlation (XC) functionals such as the local density approximation (LDA) [10] and the generalized gradient approximation (GGA) [11]. This is a particularly pressing problem, considering the fundamental relevance and the technological importance of band gaps. Both LDA and GGA inevitably underestimate the value of band gaps, leading to values which are significantly smaller than the measured ones. The reason for this failure lays on the approximate form of the XC functionals, ultimately related to the lack of derivative discontinuities of XC energy when the electron number increases through an integer point [12, 13]. In fact, in DFT the value of the band gap is obtained from the one-particle eigenvalues and the finite correction arising from this derivative discontinuity is not taken into account. Many different approaches have been proposed for reducing this intrinsic inaccuracy of DFT functionals: the weighted density approximation [14], DFT  +  U approximation [15], self-interaction correction method [16], screened exchange approximation [17], optimized effective potential [18], generalized Kohn–Sham schemes [19] (GKS), various meta-GGA potentials [20], and hybrid functionals [21], to name the most representative ones [22]. Alternative and more advanced ways to go beyond DFT are the GW approximation [23], in which the self energy of a many-body system of electrons is computed explicitly by making use of the single particle Green's function G and the screened Coulomb interaction W.

Similarly to DFT, the Hartree–Fock theory (HF) also does not provide a good prediction of band gap, but in this case this failure arises from the lack of proper correlation (C) effects in the single Slater determinant picture, which leads to many large values of the band gap. In the HF framework, this limitation is cured by incorporating electronic correlation as done in the so-called post-HF methods (Møller–Plesset perturbation theory [24], configuration interaction [25], and coupled cluster [26]).

In this variety of methods, hybrid functionals turned out to be a valid compromise between computational cost and accuracy and have been increasingly used in the computational materials science community [27–37].

Hybrid functionals are typically seen as a suitable mix between HF and DFT, obtained by replacing a portion of the exchange functional with the exact HF exchange (X). The fundamental aspect at the core of the construction of a hybrid functional is the adiabatic connection formula, which continuously transforms the non-interacting particle system to the physical interacting one:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {E_{\rm XC}=\int_0^1 E_{\rm XC,\lambda}{\rm d}\lambda} \label{eq:ad} \nonumber \end{align} \tag{ 1 }$$

where E_XC is the Kohn–Sham (KS) exchange-correlation energy and λ is the coupling parameter that tunes the electron Coulomb potential from the KS ($\lambda=0$) to the fully interacting ($\lambda=1$) limit. The first hybrid functional originally proposed by Becke in 1993 [38] is based on an approximate evaluation of the coupling-constant integral (equation (1)) using a simple linear interpolation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {E_{\rm XC}=\int_0^1 E_{\rm XC,\lambda}{\rm d}\lambda~\rightarrow~E_{\rm XC}^{\rm Hybrid}=\frac{1}{2}E_X^{\rm HF}+\frac{1}{2}E_{\rm XC}^{\rm LDA}.} \label{eq:ad2} \nonumber \end{align} \tag{ 2 }$$

$E_{\rm XC}^{\rm Hybrid}$ is generally referred to as the half-half hybrid. After proposing the half-half hybrid, Becke introduced a parametric hybrid functional including exact exchange and local (LDA) and gradient-corrected (GGA) exchange and the correlation term that has become very popular in the quantum chemistry community with the abbreviation B3LYP [38, 39]. The B3LYP depends on three parameters and incorporates only 20% of the exact HF exchange and have the following form:

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle& {E_{\rm XC}^{\rm B3LYP} = E_{\rm XC}^{\rm LDA} + \alpha_1(E_X^{\rm HF} - E_{X}^{\rm LDA})}\nonumber \\ & {+ \alpha_2(E_{X}^{\rm GGA} - E_{X}^{\rm LDA}) + \alpha_3(E_{C}^{\rm GGA} - E_{C}^{\rm LDA})} \nonumber \end{align*}$$

where the three mixing parameters $\alpha_1=0.2$, $\alpha_2=0.72$, and $\alpha_3=0.81$ are determined by fitting experimental atomization energies, electron and proton affinities and ionization potentials of the molecules in Pople's G1 data set. B3LYP has been intensively and successfully adopted for atomic and molecular calculations, but its application to periodic systems is not equally satisfactory, because B3LYP functional does not reproduce the correct exchange correlation energy for free-electron gas. This is particularly problematic for metals and heavier elements beyond the 3d transition metal series [40].

A more appropriate hybrid functional for solid-state applications is the PBEh proposed by Perdew, Burke, and Ernzerhof [41] (also referred to as PBE0) [42], which reproduces the homogeneous electron gas limit and significantly outperforms B3LYP in solids, especially in the case of systems with itinerant character (metals and small gap semiconductors) [40]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pbe0} {E_{\rm XC}^{\rm PBE0} = E_{\rm XC}^{\rm PBE} + \alpha(E_X^{\rm HF} + E_X^{\rm PBE})}. \nonumber \end{align} \tag{ 3 }$$

Here the GGA functional is chosen according to the Perdew, Burke, Ernzerhof (PBE) parameterization [43]. By analyzing the dependence of $E_{X, \lambda}$ on α and by a direct comparison with Møller–Plesset perturbation theory reference energies, Perdew, Ernzerhof and Burke have found that the choice $\alpha=0.25$ yields the best atomization energies of typical molecules [41]. Although this 0.25 choice has become a standard in PBE0 calculations, the same authors have warned that α is a system-specific and property-specific quantity. Indeed, many researchers have treated α as an adjustable parameter to reproduce specific properties, in particular band gaps [27, 31, 44, 45]. It is known that the in hybrid functionals the size of the band gap scales linearly with α [27, 31]; we verify this scaling relation in figure 1(a) for selected semiconductors.

Zoom In Zoom Out Reset image size

The PBE0 is the most representative example of a full-range hybrid, meaning that both the long-range (lr) and short-range (sr) part of the Coulomb potential are treated at HF level. The calculation of the slow-decaying lr part of the exchange integrals and exchange potential is numerically complicated and a large number of k points is needed, leading to a very slow convergence. To solve this issue, Heyd, Scuseria, and Ernzerhof (HSE) have proposed to replace the lr-exchange by the corresponding density functional counterpart [46, 47], and to keep the HF description only in the short-range limit:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rm E_{\rm XC}^{\rm HSE} = {\alpha}E_X^{HF,sr}(\mu) + (1-\alpha)E_X^{\rm PBE,sr}(\mu) + E_X^{\rm PBE,lr}(\mu) + E_C^{\rm PBE}. \label{eq:hse} \nonumber \end{align} \tag{ 4 }$$

Here $\mu^{-1}$ is the critical screening length controlling the range separation, namely the distance at which the long-range Coulomb interactions can be assumed to be negligible. Based on molecular tests the value of μ was set to 0.2 Å$^{-1}$ (corresponding to a screening length $r_s=2/\mu=10$ Å), which is routinely considered as the standard choice for range-separated HSE calculations [47].

Recently there have been successful attempts to build parameter-free PBE0 hybrids, inspired by a comparison between the GKS and the GW quasiparticle equations [30, 33, 48–50]. In Hedin's GW equations, in fact, the XC functional is replaced by the self-energy Σ, evaluated in terms of the single particle Green function G and the screened Coulomb interaction W:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Sigma(\omega) = {\rm i}{\rm G}(\omega){\rm W}(\omega), \nonumber \end{align} \tag{ 5 }$$

where the screened Coulomb interaction

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm W}(\omega) = \epsilon(\omega)^{-1}v \nonumber \end{align} \tag{ 6 }$$

is expressed in terms of the bare Coulomb interaction $v=1/{\vert {\bf r}-{\bf r}^{\prime}\vert }$ and the frequency dependent dielectric function . As already elaborated by other authors [30, 44, 48], in the GW static approximation (known as the COHSEX approximation) [23] one can replace the screening in W by the static ($\omega=0$) macroscopic dielectric constant $\newcommand{\e}{{\rm e}} \epsilon_{\infty}$, and, by setting $\newcommand{\e}{{\rm e}} \epsilon_{\infty}=1/\alpha$, the COHSEX screened Coulomb potential W^COHSEX becomes analogous to the corresponding PBE0 screened Coulomb potential:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm W}^{\rm PBE0}(\omega=0)={\alpha}v. \nonumber \end{align} \tag{ 7 }$$

Using the relation $\newcommand{\e}{{\rm e}} \alpha = 1/\epsilon_{\infty}$, some authors have recently proposed 'parameter-free' full-range and range-separated hybrids [30, 44, 51–57]. Some of these approaches are based on a fully self-consistent evaluation of $\newcommand{\e}{{\rm e}} \epsilon_{\infty}$, which is computed at the end of each hybrid cycle using different theories and approximations [30, 52, 54, 55, 57]. Among these, the protocol proposed by Skone et al, which adopts the coupled perturbed Kohn–Sham equations including local-field effects (i.e. density response due to external electric field) turned out to provide excellent agreement with experiments for the band gaps of a representative materials dataset composed by mono-atomic and binary compounds with dielectric constants spanning a wide range from 1.23 eV to 15.9 eV [30]. The same authors have also proposed a range-separated version of the same functional [56].

Inspired by the sc-hybrid of Skone and coworkers [30], in this paper we employ a similar sc-hybrid and asses its performance for the prediction of band gaps in a wide class of binary, ternary and quaternary oxides, spanning band gaps between 0 to 12 eV, static dielectric functions within 3 and 15, and diverse structural properties. Our sc-hybrid is based on the vienna ab initio simulation package (VASP) [58, 59] and the dielectric function is computed from the response to finite electric fields using the perturbation expansion after discretization (PEAD) method [60, 61]. During the writing of our manuscript, we became aware of the work of Fritsch et al in which exactly the same type of hybrid is employed for calculations of the structural and electronic properties of three binary oxides: ZnO, MgO and SnO₂ [55]. Consistent with the nomenclature proposed by Fritsch et al we also call this functional scPBE0. To assess the quality of the scPBE0 we have considered a wide material dataset formed by thirty different oxide semiconductors; our oxides dataset includes strongly correlated oxides [62], transparent conducting oxides [63], thermoelectric [64], photocatalysis [6], photovoltaic materials [65], ferroelectrics, refractory materials, and other oxides which are used in the emerging field of oxide-based electronics.

We will also compare the performance of scPBE0 with the standard screened HSE06 and full-range PBE0 results. We would like to underline that the inverse dielectric constant relation for the estimation of the 'optimum' α is not directly applicable to range-separated hybrid functionals such as HSE06, because in screened hybrid functionals screening is already present to some extent as a consequence of the range separation [27]. Considering that in HSE06 the Coulomb kernel is decomposed into a lr and sr part by way of the error function erf and the complementary error function erfc [28, 47]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{1}{r} = \underbrace{\frac{{\rm erfc}(\mu r)}{r} }_{{\rm sr}} + \underbrace{\frac{{\rm erf} (\mu r)}{r}}_{{\rm lr}}, \nonumber \end{align} \tag{ 8 }$$

the HSE06 screened Coulomb potential takes the form:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm W}^{{\rm HSE}(\omega=0)} = \alpha\frac{{\rm erfc}{(\mu{\vert {\bf r}-{\bf r}^{\prime}\vert })}}{\vert {\bf r}-{\bf r}^{\prime}\vert } \nonumber \end{align} \tag{ 9 }$$

which clearly depends on the screening parameter μ. In [27] the effect of the HSE06 screening on α is quantified as a downward shift of about 0.07 with respect to the 'optimum' $\newcommand{\e}{{\rm e}} \alpha = 1/\epsilon_{\infty}$ value.

The scPBE0 functional employed in this study is based on the full-range PBE0 functional as implemented in the VASP [58, 59], and the static dielectric function is evaluated using the perturbation expansion after discretization (PEAD) method [60, 61], by considering the perturbation of a small but finite homogeneous electric field ε on the ground-state of an insulating system. In the VASP this is done by minimizing the total energy functional $E[\{\psi^{(\varepsilon)}, \varepsilon \}]$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {E}[\{\psi^{(\varepsilon)},\varepsilon \}] = {E_0}[\{\psi^{(\varepsilon)}\}] - \Omega\varepsilon\cdot{\bf P}[\{\psi^{(\varepsilon)}\}], \label{eq:electric} \nonumber \end{align} \tag{ 10 }$$

with respect to the field-polarized Bloch functions $\{\psi^{(\varepsilon)}\}$; in equation (10) Ω is the cell volume and ${\bf P}[\{\psi^{(\varepsilon)}\}]$ is the macroscopic polarization as defined in the modern theory of polarization [60, 66, 67]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\bf P}[\{\psi^{(\varepsilon)}\}]=-\frac{2{\rm i}e}{(2\pi)^3}\sum_n\int_{\rm BZ}{\rm d}{\bf k} \langle u_{n{\bf k}}^{(\varepsilon)} \vert \nabla_{{\bf k}} \vert u_{n{\bf k}}^{(\varepsilon)} \rangle, \label{eq:mtp} \nonumber \end{align} \tag{ 11 }$$

where $u_{n{\bf k}}^{(\varepsilon)}$ is the cell-periodic part of $\{\psi_{n{\bf k}}^{(\varepsilon)}\}$. To compute the $\{\psi^{(\varepsilon)}\}$ it is necessary to minimize the functional $E[\{\psi^{(\varepsilon)}, \varepsilon \}] = \langle \psi_{n{\bf k}}^{(\varepsilon)} \vert H\vert \psi_{n{\bf k}}^{(\varepsilon)} \rangle$, which results in the following optimization problem:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H\vert \psi_{n{\bf k}}^{(\varepsilon)} \rangle = H_0\vert \psi_{n{\bf k}}^{(\varepsilon)} \rangle - \Omega\varepsilon\cdot \frac{\delta {\bf P}[\{\psi^{(\varepsilon)}\}]} {\delta\langle \psi_{n{\bf k}}^{(\varepsilon)}\vert } \label{eq:min} \nonumber \end{align} \tag{ 12 }$$

which represents the gradient of the functional in equation (10). In the PEAD scheme the dielectric tensor is given by [61]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_{ij} = \delta_{ij} + \chi_{ij} = \delta_{ij} + 4\pi\frac{\partial P_i}{\partial \varepsilon_j} \label{eq:pead} \nonumber \end{align} \tag{ 13 }$$

where $\chi_{ij}=4\pi\frac{\partial P_i}{\partial \varepsilon_j}$ is the susceptibility. If the atoms are kept fixed, as in our case, this formula finally yields the so called ion-clamped static dielectric tensor $\newcommand{\e}{{\rm e}} (\epsilon_\infty)_{ij}$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle (\epsilon_\infty)_{ij} = \delta_{ij} + 4\pi\frac{({\bf P}[\{\psi^{(\varepsilon)}\}]) - ({\bf P}[\{\psi^{(0)}\}])_i}{\varepsilon_j}. \label{eq:epsilon} \nonumber \end{align} \tag{ 14 }$$

In this formalism local-field effects are naturally taken into account through self-consistency.

Technical setup

All the calculations were performed using the projector augmented wave (PAW) pseudopotential method [68, 69]. The crystal structures of the compounds under scrutiny were fully optimized (lattice constants and atom positions) within the PBEsol parametrization for the exchange-correlation functional [70] adopting an energy cutoff for the plane wave expansion of 520 eV and a k-mesh with k-point density per reciprocal atom (KPPRA) of 8000. For all structures considered in our study, we have used k-meshes with KPRRA  >  3000 and an energy cutoff of 400 eV for all hybrid functional calculations. Based on the optimized structures, the dielectric function and the band gap have been computed at PBE0  +  PEAD level using the self-consistent protocol schematized in figure 2.

Zoom In Zoom Out Reset image size

Specifically, we started the self-consistent loop by setting $\alpha_{\rm in}=0.25$, the standard mixing parameter for PBE0 [42]. At the end of this iteration, the new α, $\alpha_{\rm out}$, has been determined as the inverse of the mean value of the diagonal elements of the dielectric function evaluated by PEAD, i.e. $\alpha_{\rm out}=1/\overline{\varepsilon_{\infty}}$. This new value of α was then used to start a new iteration. We have continued this loop until $\alpha_{\rm out} - \alpha_{\rm in} < 0.01$. When the criterion is satisfied, the converged values of α is used to obtain $\newcommand{\e}{{\rm e}} \epsilon_\infty$ as well as the band gap E_g. In figure 3 we show the practical application of this self-consistent loop for the small band gap oxide Ag₂O and for the large band gap oxide MgO: the various quantities (α, $\overline{\varepsilon_{\infty}}$, and E_g) rapidly reach a well-converged value after only four to five iterations. Considering that the computation of the dielectric properties within the PEAD formalism is about four times more expensive than a standard PBE0 cycle, overall the scPBE0 procedure is about 20 times more expensive than standard non-self-consistent PBE0.

Zoom In Zoom Out Reset image size

To systematically evaluate the performance of the PEAD-based scPBE0 functional, we have computed the bandgap and ion-clamped dielectric constant of thirty different oxide semiconductors using scPBE0, PBE0 and HSE06, and compared the results with available experimental data. The data are collected in table 1, where we also provide the list of all compounds scrutinized in our study; a graphical summary of the calculated data is provided in figure 4.

Zoom In Zoom Out Reset image size

Table 1. Band gaps (E_g, in eV) and the averaged ion-clamped static dielectric tensor, $\overline{\varepsilon_{\infty}}$, for our oxides dataset calculated by scPBE0, PBE0 (using both $\alpha_{\rm{fit}}$ and the standard choice $\alpha=0.25$), and HSE06. $\alpha_{\rm{sc}}$ and $\alpha_{\rm{fit}}$ are the self-consistent and fitted value of the mixing parameter α. The experimental methods used for band gap measurement are also indicated (TR: thermoreflectance; OA: optical absorption; XPS-BIS: x-ray photoelectron spectroscopy and bremsstrahlung isochromat spectroscopy; UPS: ultraviolet photoelectron spectroscopy; SSD: single scattering distribute spectra; PC: photoconductivity; PE: photoemission; DRS: diffuse reflectance spectra; ET: electrical transport). The last row shows the mean absolute percentage errors (MAPE) for the calculated band gaps and dielectric constants. For the evaluation of the bandgap-MAPE the data relative to Ag₂PdO₂ have been excluded due to the exceedingly large relative errors (statistically not relevant). The materials are ordered by their chemical formula, from binary to quaternary compounds.

		scPBE0			PBE0-$\alpha_{\rm{fit}}$		HSE06		PBE0		Experiment
Compound	Space group	$E_{\rm{g}}$(eV)	$\overline{\varepsilon_{\infty}}$	$\alpha_{\rm{sc}}$	$\alpha_{\rm{fit}}$	$\overline{\varepsilon_{\infty}}$	$E_{\rm{g}}$ (eV)	$\overline{\varepsilon_{\infty}}$	$E_{\rm{g}}$ (eV)	$\overline{\varepsilon_{\infty}}$	$E_{\rm{g}}$ (eV)	Method	$\overline{\varepsilon_{\infty}}$
CdO	Fm$\bar{3}$m	0.99	5.35	0.187	0.169	5.45	0.88	5.12	1.52	5.04	0.84 [71]	TR	6.2 [72]
PbO	P4/nmm	2.13	5.86	0.170	0.152	6.09	1.93	5.55	2.57	5.32	2.03 [71]	OA	7.1 [73]
MgO	Fm$\bar{3}$m	8.31	2.78	0.359	0.319	2.82	6.48	2.91	7.20	2.88	7.9 [74]	OA	2.94 [74]
ZnO	P63mc	3.38	3.63	0.276	0.283	3.60	2.49	3.74	3.15	3.67	3.44 [74]	OA	3.72 [74]
BeO	P63mc	11.31	2.86	0.349	0.282	2.91	9.52	2.93	10.25	2.92	10.59 [71]	OA	2.95 [71]
NiO	Fm$\bar{3}$m	3.51	6.60	0.152	0.198	5.75	4.45	5.13	5.20	5.07	4.30 [75]	XPS-BIS	5.77 [76]
SnO2	P42/mnm	3.60	3.86	0.269	0.269	3.86	2.71	3.92	3.40	3.23	3.6 [71]	UPS	4.06 [71]
SiO2	P3121	10.49	2.24	0.446	0.306	2.30	7.67	2.32	8.41	2.32	9.0 [77]	OA	2.3 [78]
GeO2	P42/mnm	4.42	3.77	0.264	0.345	3.63	3.56	3.80	4.26	3.83	5.35 [74]	OA	4.43 [74]
HfO2	P42/nmc	6.83	4.46	0.224	0.123	4.70	6.34	4.39	7.07	4.39	5.9 [79]	SSD
TiO2	P42/mnm	2.90	7.07	0.142	0.159	6.98	3.12	6.45	3.86	6.41	3.05 [80, 81]	ET	7.37 [71]
Ag2O	Pn$\bar{3}$m	1.27	5.92	0.169	0.160	5.93	1.21	5.70	1.90	5.34	1.20 [82, 83]	PC
Cu2O	Pn$\bar{3}$m	1.77	7.22	0.139	0.182	6.80	2.07	6.34	2.80	6.28	2.17 [74]	OA	7.11 [74]
Al2O3	R$\bar{3}$c	9.63	2.98	0.337	0.258	3.02	7.99	3.03	8.72	3.03	8.8 [80]	OA	3.4 [80]
La2O3	P$\bar{3}$m1	6.13	4.03	0.249	0.216	4.10	3.92	4.05	6.14	4.02	5.8 [84]	OA
In2O3	R$\bar{3}$c	3.32	4.10	0.244	0.214	4.01	2.72	3.68	3.38	3.92	3.02 [85]	OA	3.62 [86]
CuAlO2	R$\bar{3}$m	4.01	4.45	0.226	0.120	4.82	3.49	4.42	4.24	4.37	2.99 [87]	OA
LiCoO2	R$\bar{3}$m	4.22	4.79	0.209	0.116	5.47	4.12	4.64	4.89	4.63	2.7 [88]	XPS-BIS
LaAlO3	R$\bar{3}$c	6.27	4.08	0.246	0.254	4.16	5.57	4.05	6.30	4.06	6.33 [89]	OA	4.0 [90]
LiNbO3	R3c	5.46	4.48	0.223	0.011	5.17	4.97	4.41	5.71	4.32	3.50 [82]	OA	4.87 [91]
BiFeO3	R3c	2.87	7.03	0.142	0.125	7.17	3.40	6.12	4.14	6.10	2.67 [92]	OA	5.52 [93]
BaTiO3	P4mm	3.11	5.54	0.181	0.197	5.88	3.05	5.54	3.77	5.52	3.26 [71]	OA	5.75 [94]
PbTiO3	P4mm	2.51	7.31	0.138	0.261	6.29	2.62	6.58	3.32	6.43	3.4 [80, 95]	OA	6.25 [80]
BaSnO3	Pm$\bar{3}$m	3.15	3.98	0.251	0.246	3.34	2.45	4.25	3.14	3.98	3.1 [96]	OA	3.3 [97]
SrTiO3	Pm$\bar{3}$m	3.19	5.75	0.175	0.188	5.59	3.33	5.35	3.85	5.52	3.3 [80, 98]	OA	6.1 [80]
LaMnO3	Pmna	2.12	5.67	0.176	0.141		2.27	5.16	3.01	5.11	1.7 [99]	PE	4.9 [100]
BiVO4	C2/c	2.95	6.68	0.149	0.072	7.16	2.97	6.41	3.67	6.34	2.4 [101]	DRS
Ag2PdO2	Immm	0.86	8.49	0.115	0.026	11.05	1.19	7.39	1.89	7.34	0.18 [102]	ET
BiCuSeO	P4/nmm	0.88	12.56	0.080	0.067	12.25	1.31	10.46	1.92	10.18	0.8 [103]	OA
LaCuSeO	P4/nmm	2.59	6.85	0.147	0.157	6.80	2.68	6.43	3.36	6.36	2.8 [103]	OA
MAPE (%)		14.3	10.0				18.5	8.8	31.1	9.7

To quantify the relative predictive power of the employed computational schemes we have also computed the relative errors for all considered semiconductors, which are shown as histograms in figure 5, as well as the mean absolute percentage error defined by the formula:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {M^{\rm theory}_{E_{\rm g}}} = \frac{1}{N} \sum_{i=1}^N \left\vert \frac{E_{\rm g}^{\rm expt}-E_{\rm g}^{\rm theory}}{E_{\rm g}^{\rm expt}} \right\vert \times 100\% \label{eq:mare} \nonumber \end{align} \tag{ 15 }$$

where $N=30$ is the total number of compounds considered, and theory refers to the specific type of functional used in the calculation. A similar formula has been used for evaluating the MAPE relative to $\overline{\varepsilon_{\infty}}$ (${M^{\rm theory}_{\overline{\varepsilon_{\infty}}}}$). The value of the MAPE, collected in the last row of table 1 indicates that scPBE0 is the most accurate method for the evaluation of the band gaps, with an overall MAPE $M^{\rm {scPBE0}}_{E_{\rm g}}=14.3\%$, which is smaller than the corresponding HSE06 ($M^{\rm HSE06}_{E_{\rm g}}= 18.5\%$), and PBE0 ($M^{\rm PBE0}_{E_{\rm g}}= 31.1\%$) values. On the other side, for what concerns the ion-clamped dielectric function, the three methods are almost equivalent to PBE0, delivering slightly (1%) better estimations: $M^{\rm scPBE0}_{\overline{\varepsilon_{\infty}}}=10\%$, $M^{\rm PBE0}_{\overline{\varepsilon_{\infty}}}= 8.8\%$), and $M^{\rm HSE06}_{\overline{\varepsilon_{\infty}}}= 9.7\%$.

Zoom In Zoom Out Reset image size

The calculated data indicate that scPBE0 performs better than HSE06 and PBE0 for band gaps smaller than  ≈8 eV, whereas for larger gaps PBE0 appears to be a better choice. The best performance of scPBE0 with respect to the fixed-α PBE0 functional is exemplified in figure 3: for small band gap materials such as Ag₂O (1.20 eV), α decreases monotonically with increasing numbers of iterations and finally reaches a converged value of 0.169, significantly smaller than the standard PBE0 value, 0.25; in the larger band gap limit, represented by MgO ($E_{\rm g} = 7.9$ eV), the self-consistent value of α, 0.359, is, on the other hand, larger than 0.25.

Apparently, there are exceptions to this general trend: for BeO ($E_{\rm g} = 10.59$ eV), SiO₂ ($E_{\rm g} = 9.0$ eV), and Al₂O₃ ($E_{\rm g} = 8.8$ eV) scPBE0 predicts larger band gaps than PBE0 due to the larger value of the self-consistent α as compared to the standard 0.25 choice. However, for these low dielectric constant materials (≈3) excitonic effects, not taken into account in the present study, can be very strong and thus unavoidably affect the comparison between calculated and measured data. For example, in SiO₂, the excitonic binding energy is as large as 1.2 eV [104] and the measured band gap is 9.0 eV; by taking into account this excitonic shift in our calculations, scPBE0 turns out to provide the best prediction of the band gap, 9.29 eV, whereas PBE0 and HSE are off by more than 2 eV ($E_{\rm g}^{\rm PBE0} = 7.21$ eV and $E_{\rm g}^{\rm HSE06} = 6.47$ eV). Therefore, the fact that scPBE0 gives band gaps larger than experiments for materials with a weak dielectric screening (i.e. typically large excitonic effects) should be considered a positive aspect. The inclusion of excitonic effects would lead to a decrease of the band gap thus improving the comparison with the experiment. Unfortunately, measuring or calculating excitonic binding energies is an extremely complicated task [105] and only very few values are available in literature for the compounds included in our materials dataset. The materials for which scPBE0 overestimates the band gaps with respect to the measured values are: BeO, SiO₂, HfO₂, Al₂O₃, CuAlO₂, LiCoO₂, NaSbO₃, and LiNbO₃. Among these materials, to our knowledge, the excitonic binding energy is accurately reported only for SiO₂ [105].

In the small band gap limit (usually large $\newcommand{\e}{{\rm e}} \epsilon_{\infty}$) excitonic effects are expected to be weak (from few tens to few hundreds of meV), scPBE0 performs much better than PEB0. It is worth noting that for small band gaps in the range 0–3 eV, the screened HSE06 functionals performs relatively well due to the intrinsic screening incorporated in this type of functional via range-separation which allows for a better description of the metallic limit. Considering that no extra screening is adopted in scPBE0, a comparable performance with HSE06 for these compounds indicates a success of the scPBE0 functional for small band gap oxides.

Lastly, we would like to mention that we found that scPBE0 breaks down for two oxides with small band gaps, i.e. PdO ($E_{\rm{g}}=1.0$ eV [82]) and AgBiO₃ ($E_{\rm{g}}=0.8$ eV [106]): the self-consistent procedure in these two cases ends up with a metallic solution due to the a considerable overestimation of $\newcommand{\e}{{\rm e}} \epsilon_{\infty}$.

As mentioned in the introduction, a practical remedy to improve the agreement with experiment within a non-self-consistent hybrid framework is to treat α as a fitting parameter, for instance by setting α to the value that reproduces the measured band gap, $\alpha_{\rm fit}$. This can be done by using the linear relation between $E_{\rm{g}}$ and α depicted in figure 1. The resulting values of $\alpha_{\rm fit}$, also reported in table 1, range from 0.011 (LiNbO₃) to 0.345 (GeO₂); this wide range of variation of $\alpha_{\rm fit}$ is a further demonstration that this quantity is indeed material-dependent. A direct comparison between $\alpha_{\rm fit}$ and the corresponding self-consistent value adopted in scPBE0, $\alpha_{\rm sc}$, shown in figure 6(a), shows strong deviations, with an overall tendency of scPBE0 to yield larger values of α. This difference between scPBE0 and fitted-PBE0 is due to the fact that α and ${\varepsilon_{\infty}}$ do not scale linearly [31, 107, 108]. This is graphically shown in figure 6(b), where we report the values of $\overline{\varepsilon_{\infty}}$ according to the relation $\alpha=1/\overline{\varepsilon_{\infty}}$ used for the self-consistent procedure, as well as the values of $\overline{\varepsilon_{\infty}}$ obtained at PBE0-PEAD level using the self-consistent value of α, $\alpha_{\rm sc}$. As a result, by fitting α with respect to the band gap the MAPE associated with the band gaps is minimized but the description of the dielectric properties gets worse. On the other hand we emphasize once more that scPBE0 is based on a self-consistent evaluation of the dielectric constant, which guarantees a good description of the dielectric screening and, at the same time, the self-consistent α obtained through the relation $\newcommand{\e}{{\rm e}} \alpha = 1/\epsilon_{\infty}$ allows for a generally good prediction of the band gap.

Zoom In Zoom Out Reset image size

The non-linearity between E_g and $\newcommand{\e}{{\rm e}} \epsilon_{\infty}$ has been pointed out first by Moss [107, 108], who proposed the formula, known as Moss relation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon_{\infty}^2 = C \cdot {E_{\rm g}} \nonumber \end{align} \tag{ 16 }$$

which has been verified for many elements and semiconducting compounds, including oxides [109], though deviations with respect to the measured values have been reported for band gap smaller than 1.4 eV [110].

In figure 7 we show that the Moss relation is sufficiently well fulfilled for our material dataset, not only experimentally but also based on the scPBE0 results: a linear fit of the measured and scPBE0 data leads to two almost overlapping lines with $C\approx0.012$.

Zoom In Zoom Out Reset image size

In summary, in this paper we have assessed the performance of the unscreened self-consistent hybrid functional scPBE0 for a wide set of oxide materials with different dielectric, optical and structural characteristics. In the scPBE0 scheme the HF mixing parameter α is determined by a self-consistent evaluation of ion-clamped dielectric constants using the PEAD method and by making use of the relation $\newcommand{\e}{{\rm e}} \alpha = 1/\epsilon_{\infty}$. Therefore this method is de facto an ab initio approach, without adjustable or empirical parameters. We have compared the scPBE0 results with those obtained using the standard PBE0 functional (with a fix value of $\alpha=0.25$) and the screened (range-separated) HSE06 functional. Our results shows that scPBE0 outperforms PBE0 and, to a lesser extent, HSE06 for the prediction of band gaps, with an overall MAPE of 14.3%. As for the prediction of the ion-clamped dielectric constant all methods deliver similar results, associated with a MAPE of about 10%.

Importantly, for materials characterized by weak dielectric screening and therefore prone to exhibit large excitonic effects, scPBE0 (unlike PBE0 and HSE06) furnishes band gaps larger than the measured data. The inclusion of excitonic effects (neglected in this study) would therefore improve the agreement with experiments only for scPBE0; this has been verified for SiO₂, for which the excitonic binding energy is as large as 1.2 eV. On the other hand, in the very small band gap limit, we found that scPBE0 could lead to metallic solutions. This was the case for PdO and AgBiO₃, which are better described at PBE0 and HSE06 level.

Moreover, we have shown that by taking advantage of the linear relation between α and E_g it is possible to obtain a set of optimally fitted α that guarantees an excellent agreement with experiments at PBE0 level. However, due to the non-linearity between E_g and $\newcommand{\e}{{\rm e}} \epsilon_{\infty}$, with these fitted values of α the dielectric properties are not well described. Finally, we have verified that the Moss relation ($\newcommand{\e}{{\rm e}} \epsilon_{\infty}^2 = C \cdot E_{\rm g}$) is fulfilled for both the experimental and the scPBE0 sets of data.

In conclusion, scPBE0 represents a valuable scheme to obtain a satisfactory description of band gaps and dielectric properties in materials fully ab initio, thus representing a step forward with respect to parameter-dependent hybrid functionals. Clearly, to achieve a more accurate account of the optical and screening properties of materials (with errors smaller than 2–3%) it is necessary to go beyond the Hartree–Fock picture. A natural option is the GW approach with the inclusions of excitonic effects, which is however computationally much more demanding than hybrid functional schemes.

This work was sponsored by the FWF project INDOX (Grant No. I1490-N19). All calculations were performed on the Vienna Scientific Cluster (VSC).