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

TL;DR

This study evaluates the self-consistent hybrid functional scPBE0 for calculating band gaps in oxide semiconductors, showing it outperforms some existing methods but has limitations in certain cases.

Contribution

It introduces a self-consistent procedure for determining the mixing parameter in hybrid functionals and assesses its accuracy across diverse oxide semiconductors.

Findings

scPBE0 provides better band gap predictions than PBE0 and HSE06.

Mean absolute percentage error for band gaps is 14%.

Overestimates gaps in materials with weak dielectric screening.

Abstract

In this paper we assess the predictive power of the self-consistent hybrid functional scPBE0 in calculating the band gap of oxide semiconductors. The computational procedure is based on the self-consistent evaluation of the mixing parameter $\alpha$ by means of an iterative calculation of the static dielectric constant using the perturbation expansion after discretization (PEAD) method and making use of the relation $\alpha = 1/\epsilon_{\infty}$. Our materials dataset is formed by 30 compounds covering a wide range of band gaps and dielectric properties, and includes materials with a wide spectrum of application as thermoelectrics, photocatalysis, photovoltaics, transparent conducting oxides, and refractory materials. Our results show that the scPBE0 functional provides better band gaps than the non self-consistent hybrids PBE0 and HSE06, but scPBE0 does not show significant improvement on the description of the static dielectric constants. Overall, the scPBE0 data exhibit a mean absolute percentage error of 14 \% (band gaps) and 10 \% ($\epsilon_\infty$). For materials with weak dielectric screening and large excitonic biding energies scPBE0, unlike PBE0 and HSE06, overestimates the band gaps, but the value of the gap become very close to the experimental value when excitonic effects are included (e.g. for SiO$_2$). However, special caution must be given to the compounds with small band gaps due to the tendency of scPBE0 to overestimate the dielectric constant in proximity of the metallic limit.

Figures

Tables

Table 1. Table 1: Band gaps (E g , in eV) and the averaged ion-clamped static dielectric tensor, ε ∞ ¯ ¯ subscript 𝜀 \overline{\varepsilon_{\infty}} , for our oxides dataset calculated by scPBE0, PBE0, and HSE06. α sc subscript 𝛼 sc \alpha_{\rm{sc}} and α fit subscript 𝛼 fit \alpha_{\rm{fit}} are the self-consistent and fitted value of the mixing parameter α 𝛼 \alpha . ε ∞ ¯ ¯ subscript 𝜀 \overline{\varepsilon_{\infty}} under scPBE0, fitted, HSE06, and PBE0 are high-frequency dielectric constants calculated by using scPBE0, the fitted α 𝛼 \alpha with experiment ( α fit subscript 𝛼 fit \alpha_{\rm{fit}} ), HSE06, and PBE0, respectively. The experiment approaches used for band gap measurement are indicated under Method (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 show 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 2 PdO 2 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			fitted		HSE06		PBE0		Experiment
Compound	Space group	Eg(eV)	$\overline{{\epsilon }_{\infty }}$	${\alpha }_{\mathrm{sc}}$	${\alpha }_{\mathrm{fit}}$	$\overline{{\epsilon }_{\infty }}$	Eg (eV)	$\overline{{\epsilon }_{\infty }}$	Eg (eV)	$\overline{{\epsilon }_{\infty }}$	Eg (eV)	Method	$\overline{{\epsilon }_{\infty }}$
CdO	Fm$\overline{3}$m	0.99	5.35	0.187	0.169	5.45	0.88	5.12	1.52	5.04	0.84Madelung et al.	TR	6.2Van Vechten (1969)
PbO	P4/nmm	2.13	5.86	0.170	0.152	6.09	1.93	5.55	2.57	5.32	2.03Madelung et al.	OA	7.1Van den Broek (1969)
MgO	Fm$\overline{3}$m	8.31	2.78	0.359	0.319	2.82	6.48	2.91	7.20	2.88	7.9Madelung (2004)	OA	2.94Madelung (2004)
ZnO	P63mc	3.38	3.63	0.276	0.283	3.60	2.49	3.74	3.15	3.67	3.44Madelung (2004)	OA	3.72Madelung (2004)
BeO	P63mc	11.31	2.86	0.349	0.282	2.91	9.52	2.93	10.25	2.92	10.59Madelung et al.	OA	2.95Madelung et al.
NiO	Fm$\overline{3}$m	3.51	6.60	0.152	0.198	5.75	4.45	5.13	5.20	5.07	4.30Sawatzky and Allen (1984)	XPS-BIS	5.77Pecharroman and Iglesias (1994)
SnO2	P42/mnm	3.60	3.86	0.269	0.269	3.86	2.71	3.92	3.40	3.23	3.6Madelung et al.	UPS	4.06Madelung et al.
SiO2	P3121	10.49	2.24	0.446	0.306	2.30	7.67	2.32	8.41	2.32	9.0Vella et al. (2011)	OA	2.3Malyi et al. (2016)
GeO2	P42/mnm	4.42	3.77	0.264	0.345	3.63	3.56	3.80	4.26	3.83	5.35Madelung (2004)	OA	4.43Madelung (2004)
HfO2	P42/nmc	6.83	4.46	0.224	0.123	4.70	6.34	4.39	7.07	4.39	5.9Cheynet et al. (2007)	SSD
TiO2	P42/mnm	2.90	7.07	0.142	0.159	6.98	3.12	6.45	3.86	6.41	3.05Robertson (2000); Cronemeyer (1952)	ET	7.37Madelung et al.
Ag2O	Pn$\overline{3}$m	1.27	5.92	0.169	0.160	5.93	1.21	5.70	1.90	5.34	1.20Xu and Schoonen (2000); Fortin and Weichman (1964)	PC
Cu2O	Pn$\overline{3}$m	1.77	7.22	0.139	0.182	6.80	2.07	6.34	2.80	6.28	2.17Madelung (2004)	OA	7.11Madelung (2004)
Al2O3	R$\overline{3}$c	9.63	2.98	0.337	0.258	3.02	7.99	3.03	8.72	3.03	8.8Robertson (2000)	OA	3.4Robertson (2000)
La2O3	P$\overline{3}$m1	6.13	4.03	0.249	0.216	4.10	3.92	4.05	6.14	4.02	5.8Shang et al. (2004)	OA
In2O3	R$\overline{3}$c	3.32	4.10	0.244	0.214	4.01	2.72	3.68	3.38	3.92	3.02King et al. (2009)	OA	3.62Prathap et al. (2006)
CuAlO2	R$\overline{3}$m	4.01	4.45	0.226	0.120	4.82	3.49	4.42	4.24	4.37	2.99Pellicer-Porres et al. (2006)	OA
LiCoO2	R$\overline{3}$m	4.22	4.79	0.209	0.116	5.47	4.12	4.64	4.89	4.63	2.7van Elp et al. (1991)	XPS-BIS
LaAlO3	R$\overline{3}$c	6.27	4.08	0.246	0.254	4.16	5.57	4.05	6.30	4.06	6.33Cicerrella et al. (2005)	OA	4.0Nunley et al. (2016)
LiNbO3	R3c	5.46	4.48	0.223	0.011	5.17	4.97	4.41	5.71	4.32	3.50Xu and Schoonen (2000)	OA	4.87Barker and Loudon (1967)
BiFeO3	R3c	2.87	7.03	0.142	0.125	7.17	3.40	6.12	4.14	6.10	2.67Xu et al. (2009)	OA	5.52Li et al. (2007)
BaTiO3	P4mm	3.11	5.54	0.181	0.197	5.88	3.05	5.54	3.77	5.52	3.26Madelung et al.	OA	5.75Huang et al. (2006)
PbTiO3	P4mm	2.51	7.31	0.138	0.261	6.29	2.62	6.58	3.32	6.43	3.4Robertson (2000); Pandey et al. (2005)	OA	6.25Robertson (2000)
BaSnO3	Pm$\overline{3}$m	3.15	3.98	0.251	0.246	3.34	2.45	4.25	3.14	3.98	3.1Mizoguchi et al. (2004)	OA	3.3Stanislavchuk et al. (2012)
SrTiO3	Pm$\overline{3}$m	3.19	5.75	0.175	0.188	5.59	3.33	5.35	3.85	5.52	3.3Robertson (2000); Van Benthem et al. (2001)	OA	6.1Robertson (2000)
LaMnO3	Pmna	2.12	5.67	0.176	0.141	6.14	2.27	5.16	3.01	5.11	1.7Saitoh et al. (1995)	PE	4.9Arima and Tokura (1995)
BiVO4	C2/c	2.95	6.68	0.149	0.072	7.16	2.97	6.41	3.67	6.34	2.4Dunkle et al. (2009)	DRS
Ag2PdO2	Immm	0.86	8.49	0.115	0.026	11.05	1.19	7.39	1.89	7.34	0.18Schreyer and Jansen (2001)	ET
BiCuSeO	P4/nmm	0.88	12.56	0.080	0.067	12.25	1.31	10.46	1.92	10.18	0.8Hiramatsu et al. (2007)	OA
LaCuSeO	P4/nmm	2.59	6.85	0.147	0.157	6.80	2.68	6.43	3.36	6.36	2.8Hiramatsu et al. (2007)	OA
MAPE (%)		14.3	10.0				18.5	8.8	31.1	9.7

Equations

$${\mathrm{E}}_{\mathrm{XC}}={\int }_0^1{\mathrm{E}}_{\mathrm{XC},\lambda }\mathrm{d}\lambda$$

$${\mathrm{E}}_{\mathrm{XC}}={\int }_0^1{\mathrm{E}}_{\mathrm{XC},\lambda }\mathrm{d}\lambda \to {\mathrm{E}}_{\mathrm{XC}}^{\mathrm{Hybrid}}=\frac{1}{2}{\mathrm{E}}_{\mathrm{X}}^{\mathrm{HF}}+\frac{1}{2}{\mathrm{E}}_{\mathrm{XC}}^{\mathrm{LDA}}$$

$${\mathrm{E}}_{\mathrm{XC}}^{B3LYP}={\mathrm{E}}_{\mathrm{XC}}^{\mathrm{LDA}}+{\alpha }_1({\mathrm{E}}_{\mathrm{X}}^{\mathrm{HF}}-{\mathrm{E}}_{\mathrm{X}}^{\mathrm{LDA}})+{\alpha }_2({\mathrm{E}}_{\mathrm{X}}^{\mathrm{GGA}}-{\mathrm{E}}_{\mathrm{X}}^{\mathrm{LDA}})+{\alpha }_3({\mathrm{E}}_{\mathrm{C}}^{\mathrm{GGA}}-{\mathrm{E}}_{\mathrm{C}}^{\mathrm{LDA}})$$

$${\mathrm{E}}_{\mathrm{XC}}^{PBE0}={\mathrm{E}}_{\mathrm{XC}}^{\mathrm{PBE}}+\alpha ({\mathrm{E}}_{\mathrm{X}}^{\mathrm{HF}}+{\mathrm{E}}_{\mathrm{X}}^{\mathrm{PBE}})$$

$${\mathrm{E}}_{\mathrm{XC}}^{\mathrm{HSE}}=\alpha {\mathrm{E}}_{\mathrm{X}}^{\mathrm{HF},\mathrm{sr}}(\mu )+(1-\alpha ){\mathrm{E}}_{\mathrm{X}}^{\mathrm{PBE},\mathrm{sr}}(\mu )+{\mathrm{E}}_{\mathrm{X}}^{\mathrm{PBE},\mathrm{lr}}(\mu )+{\mathrm{E}}_{\mathrm{C}}^{\mathrm{PBE}},$$

$$\Sigma (\omega )=i\mathrm{G}(\omega )\mathrm{W}(\omega ),$$

$$\mathrm{W}(\omega )=ϵ(\omega {)}^{-1}\mathrm{v}$$

$${\mathrm{W}}^{PBE0}(\omega =0)=\alpha v.$$

$$\frac{1}{r}=\underset{\text{sr}}{\underbrace{\frac{\mathrm{erfc}(\mu r)}{r}}}+\underset{\text{lr}}{\underbrace{\frac{\mathrm{erf}(\mu r)}{r}}}$$

$${\mathrm{W}}^{\mathrm{HSE}(\omega =0)}=\alpha \frac{\mathrm{erfc}(\mu |\mathbf{r}-{\mathbf{r}}^{\prime }|)}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}$$

$$\mathrm{E}[\{{\psi }^{(\text{upvarepsilon})},\text{upvarepsilon}\}]={\mathrm{E}}_0[\{{\psi }^{(\text{upvarepsilon})}\}]-\mathrm{\Omega }\text{upvarepsilon}\cdot \text{P}[\{{\psi }^{(\text{upvarepsilon})}\}],$$

$$\text{P}[\{{\psi }^{(\text{upvarepsilon})}\}]=-\frac{2ie}{(2\pi {)}^3}\sum _n{\int }_{BZ}d\text{k}⟨u_{n\text{k}}^{(\text{upvarepsilon})}|{\nabla }_{\text{k}}|u_{n\text{k}}^{(\text{upvarepsilon})}⟩,$$

$$H|{\psi }_{n\text{k}}^{(\text{upvarepsilon})}⟩=H_0|{\psi }_{n\text{k}}^{(\text{upvarepsilon})}⟩-\Omega \text{upvarepsilon}\cdot \frac{\delta \text{P}[\{{\psi }^{(\text{upvarepsilon})}\}]}{\delta ⟨{\psi }_{n\text{k}}^{(\text{upvarepsilon})}|}$$

$${ϵ}_{ij}={\delta }_{ij}+{\chi }_{ij}={\delta }_{ij}+4\pi \frac{\partial P_i}{\partial {\text{upvarepsilon}}_j}$$

$$({ϵ}_{\infty }{)}_{ij}={\delta }_{ij}+4\pi \frac{(\text{P}[\{{\psi }^{(\text{upvarepsilon})}\}]-(\text{P}[\{{\psi }^{(0)}\}]{)}_i}{{\text{upvarepsilon}}_j}$$

$${\mathrm{M}}_{{\mathrm{E}}_{\mathrm{g}}}^{\mathrm{theory}}=\frac{1}{\mathrm{N}}\sum _{\mathrm{i}=1}^{\mathrm{N}}\left|\frac{{\mathrm{E}}_{\mathrm{g}}^{\mathrm{expt}}-{\mathrm{E}}_{\mathrm{g}}^{\mathrm{theory}}}{{\mathrm{E}}_{\mathrm{g}}^{\mathrm{expt}}}\right|\times 100\%$$

$${ϵ}_{\infty }^2=\mathrm{C}\cdot {\mathrm{E}}_{\mathrm{g}}$$

Full text

I Introduction

Oxides are one of the most physical complicated and technological important systems, exhibiting a variety of structures and interesting properties. Metal oxides have broaden applications in many relevant areas including superconductivity Emery (1987), ferroelectricity Cohen (1992), dielectric Shannon (1993), transparent conducting oxides Ginley and Bright (2000), thermoelectricity Fergus (2012), photocatalysis Linsebigler et al. (1995), and photovoltaic materials Fortunato et al. (2007). One of the fundamental properties of materials, which is functional for virtually any type of application is the band gaps. The 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 the metal-oxygen bonding that can very between nearly ionic to highly covalent or metallic Rao (1989), and by 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 the density functional theory (DFT), many materials-specific properties have been calculated without relying on any empirical parameter Lejaeghere et al. (2016). However, it is still very difficult to calculate band gap accurately within the DFT framework using standard semilocal exchange-correlation (XC) functionals such as the local density approximation (LDA) Ceperley and Alder (1980) and the generalized gradient approximation (GGA) Perdew et al. (1992). This is a particularly pressing problem, considering the fundamental relevance and the technological importance of bands 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 the XC energy when the electron number increases through an integer point Perdew and Levy (1983); Sham and Schlüter (1983). 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 to reducing this intrinsic inaccuracy of DFT functionals: the weighted density approximation Alonso and Girifalco (1978), the DFT+U approximation Anisimov et al. (1991), the self-interaction correction method Perdew and Zunger (1981), the screened exchange approximation Bylander and Kleinman (1990), the optimized effective potential Grüning et al. (2006), generalized Kohn-Sham schemes Seidl et al. (1996) (GKS), various meta-GGA potentials Tao et al. (2003), and hybrid functionals Becke (1993a), to name the most representative ones Kümmel and Kronik (2008). Alternative and more advanced ways to go beyond DFT are the GW approximation Hedin (1965), 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, also the Hartree-Fock theory (HF) 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 much 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 Møller and Plesset (1934), Configuration interaction Paldus (1974), and Coupled cluster Cizek (1966)).

In this variety of methods, hybrid functionals turned to be a valid compromise between computational cost and accuracy and have been increasing used in the computational materials science community He and Franchini (2012); Franchini (2014); Marsman et al. (2008); Skone et al. (2014); Alkauskas et al. (2011a); Iori et al. (2012); Henderson et al. (2011); Archer et al. (2011); He et al. (2012); Gerosa et al. (2015); Garza and Scuseria (2016, 2016).

Hybrid functionals are typically seen as a suitable mixing 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:

[TABLE]

where EXC is the Kohn-Sham (KS) exchange-correlation energy and $\lambda$ 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 Becke (1993b) is based on an approximate evaluation of the coupling-constant integral (Eq.1) using a simple linear interpolation:

[TABLE]

$E_{XC}^{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 correlation terms that has become very popular in the quantum chemistry community with the abbreviation B3LYP Becke (1993c); Stephens et al. (1994). The B3LYP depends on three parameters and incorporates only 20% of the exact HF exchange and have the following form:

[TABLE]

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 the free-electron gas. This is particularly problematic for metals and heavier elements, beyond the 3$d$ transition metal series Paier et al. (2007).

A more appropriate hybrid functional for solid-state applications is the PBEh proposed by Perdew, Burke, and Ernzerhof Perdew et al. (1996a) (also referred to as PBE0) Adamo and Barone (1999a), 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) Paier et al. (2007):

[TABLE]

Here the GGA functional is chosen according to the Perdew, Burke, Ernzerhof (PBE) parameterization Perdew et al. (1996b). By analyzing the dependence of $E_{X,\lambda}$ on $\alpha$ 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 Perdew et al. (1996c). Although this 0.25 choice has become a standard in PBE0 calculations, the same authors have warned that $\alpha$ is a system-specific and property-specific quantity. Indeed many researchers has treated $\alpha$ as an adjustable parameter to reproduce specific properties, in particular band gaps Marques et al. (2011); Alkauskas et al. (2011b); He and Franchini (2012); Varley et al. (2012). It is known that the in hybrid functionals the size of the band gap scales linearly with $\alpha$ He and Franchini (2012); Alkauskas et al. (2011b); we verify this scaling relation in Fig. 1(a) for selected semiconductors.

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 Heyd et al. (2003, 2006), and to keep the HF description only in the short-range limit:

[TABLE]

Here $\mu^{-1}$ is the critical screening length controlling the range separation, namely the distance at which the short-range Coulomb interactions can be assumed to be negligible. Based on molecular tests the value of $\mu$ 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 Heyd et al. (2006).

Recently there has been successful attempts to build parameter-free PBE0 hybrid, inspired by a comparison between the GKS and the GW quasiparticle equations Moussa et al. (2012); Henderson et al. (2011); Del Sole et al. (1994); Gygi and Baldereschi (1989); Skone et al. (2014). In the Hedin’s GW equations, in fact, the XC functional is replaced by the self-energy $\Sigma$, evaluated in terms of the single particle Green function G and the screened Coulomb interaction W:

[TABLE]

where and the screened Coulomb interaction

[TABLE]

is expressed in terms of the bare Coulomb interaction $v=1/{|{\bf r}-{\bf r}^{\prime}|}$ and the frequency dependent dielectric function $\epsilon$. As already elaborated by other authors Moussa et al. (2012); Skone et al. (2014); Marques et al. (2011), in the GW static approximation (known as the COHSEX approximation) Hedin (1965) one can replace the screening in $W$ by the static ($\omega=0$) macroscopic dielectric constant $\epsilon_{\infty}$, and, by setting $\epsilon_{\infty}=1/\alpha$, the COHSEX screened Coulomb potential WCOHSEX become analogous to the corresponding PBE0 screened Coulomb potential:

[TABLE]

Using the relation $\alpha=1/\epsilon_{\infty}$, some authors have recently proposed “parameter-free“ full-range and range-separated hybrids Marques et al. (2011); Shimazaki and Asai (2008, 2009); Refaely-Abramson et al. (2013); Koller et al. (2013); Skone et al. (2014); Fritsch et al. (2017); Skone et al. (2016); Erba (2017). Some of these approaches are based on a fully self-consistent evaluation of $\epsilon_{\infty}$, which is computed at the end of each hybrid cycle using different theories and approximations Shimazaki and Asai (2009); Koller et al. (2013); Skone et al. (2014); Fritsch et al. (2017); Erba (2017). 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.Skone et al. (2014) The same authors have also proposed a range-separated version the same functionalSkone et al. (2016).

Inspired by the sc-hybrid of Skone and coworkers Skone et al. (2014), 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) Kresse and Hafner (1993); Kresse and Furthmüler (1996) and the dielectric function is computed from the response to finite electric fields using the perturbation expansion after discretization (PEAD) method Nunes and Gonze (2001); Souza et al. (2002). 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 SnO2 Fritsch et al. (2017). 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 Dagotto and Tokura (2008), transparent conducting oxides Zhang et al. (2016), thermoelectric Li et al. (2012), photocatalysis Linsebigler et al. (1995), photovoltaic materials Sullivan et al. (2016), ferroelectrics, refractory materials, and other oxides which are sued 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’ $\alpha$ 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 He and Franchini (2012). 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 Heyd et al. (2006); Franchini (2014):

[TABLE]

the HSE06 screened Coulomb potential takes the form:

[TABLE]

which clearly depends on the screening parameter $\mu$. In Ref. He and Franchini (2012) the effect of the HSE06 screening on $\alpha$ is quantified as a downward shift of about 0.07 with respect to the ’optimum’ $\alpha=1/\epsilon_{\infty}$ value.

II Methodology and Computational Details

The scPBE0 functional employed in this study is based on the full-range PBE0 functional as implemented in the VASP Kresse and Hafner (1993); Kresse and Furthmüler (1996), and the static dielectric function is evaluated using the perturbation expansion after discretization (PEAD) method Nunes and Gonze (2001); Souza et al. (2002), by considering the perturbation of a small but finite homogeneous electric field $\upvarepsilon$ on the ground-state of an insulating system. In the VASP this is done by minimizing the total energy functional $E[\{\psi^{(\upvarepsilon)},\upvarepsilon\}]$ :

[TABLE]

with respect to the field-polarized Bloch functions $\{\psi^{(\upvarepsilon)}\}$; in Eq. 10 $\Omega$ is the cell volume and $\textbf{P}[\{\psi^{(\upvarepsilon)}\}]$ is the macroscopic polarization as defined in the modern theory of polarization Nunes and Gonze (2001); King-Smith and Vanderbilt (1993); Resta (1994):

[TABLE]

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

[TABLE]

which represents the gradient of the functional in Eq. 10. In the PEAD scheme the dielectric tensor is given by Souza et al. (2002):

[TABLE]

where $\chi_{ij}=4\pi\frac{\partial P_{i}}{\partial\upvarepsilon_{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 $(\epsilon_{\infty})_{ij}$:

[TABLE]

In this formalism local-field effects are naturally taken into accounted through self-consistently.

II.1 Technical Setup

All the calculations were performed using the projector augmented wave (PAW) pseudopotential method Blöchl (1994); Kresse and Joubert (1999). 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 Perdew et al. (2008) adopting an energy cutoff for the plane wave expansion of 520 eV and a k-mesh with KPPRA (k-point density per reciprocal atom) of 8000. For all structures considered in our study we have used k-meshes with KPRRA >3000 and 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 Fig. 2.

Specifically, we started the self-consistent loop by setting $\alpha_{\rm in}$=0.25, the standard mixing parameter for PBE0 Adamo and Barone (1999b). At the end of this iteration the new $\alpha$, $\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 $\alpha$ 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 $\alpha$, $\epsilon_{\infty}$ as well as the band gap Eg are obtained. In Fig. 3 we show the practical application of this self-consistent loop for the small band gap oxide Ag2O and for the large band gap oxide MgO: the various quantities ($\alpha$, $\overline{\varepsilon_{\infty}}$, and Eg) rapidly reach a well-converged value after only 4-5 iterations. Considering that the computation of the dielectric properties within the PEAD formalism is about 4 times more expensive than a standard PBE0 cycle, overall the scPBE0 procedure is about 20 times more expensive than standard non-self-consistent PBE0.

III Results and discussion

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 Tab. 1, where we also provide the list of all compounds scrutinized in our study; a graphical summary of the calculated data is provided in Fig. 4.

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 Fig. 5, as well as the mean absolute percentage error defined by the formula:

[TABLE]

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}}$ ($\rm{M^{theory}_{\overline{\varepsilon_{\infty}}}}$). The value of the MAPE, collected in the last row of Tab. 1 indicate that scPBE0 is the most accurate method for the evaluation of the band gaps, with an overall MAPE M${}^{\rm{scPBE0}}_{\rm E_{g}}$=14.3 %, which is smaller than the corresponding HSE06 (M${}^{\rm{HSE06}}_{\rm E_{g}}$= 18.5%), and PBE0 (M${}^{\rm PBE0}_{\rm E_{g}}$= 31.1%) values. On the other side, for what concern the ion-clamped dielectric function, the three methods are almost equivalent with 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%.

The calculated data indicate that scPBE0 performs better than HSE06 and PBE0 for band gaps smaller than $\approx$ 8 eV, whereas for larger gaps PBE0 appears to be a better choice. The best performance of scPBE0 with respect to the fixed-$\alpha$ PBE0 functional is exemplified in Fig. 3: for small band gap materials such as Ag2O (1.20 eV) $\alpha$ decreases monotonically with increasing number of iteration 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 (Eg = 7.9 eV), the self-consistent value of $\alpha$, 0.359, is, on the other hand, larger than 0.25.

Apparently, there are exceptions to this general trend: for BeO (Eg = 10.59 eV), SiO2 (Eg = 9.0 eV), and Al2O3 ($E_{g}$ = 8.8 eV) scPBE0 predicts larger band gaps than PBE0 due to the larger value of the self-consistent $\alpha$ as compared to the standard 0.25 choice. However, for these low dielectric constant materials ($\approx$ 3) excitonic effects, not taken into account in the present study, can be very strong and this unavoidably affect the comparison between calculated and measured data. For example in SiO2 the excitonic binding energy is as larger as 1.2 eV Kresse et al. (2012) 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 as a positive aspect. The inclusion of excitonic effects would lead to a decrease of the band gap thus improving the comparison with experiment. Unfortunately, measuring or calculating excitonic binding energies is an extremely complicated task Shishkin et al. (2007) 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, SiO2, HfO2, Al2O3, CuAlO2, LiCoO2, NaSbO3, and LiNbO3. Among these materials, to our knowledge, the excitonic binding energy is accurately reported only for SiO2 Shishkin et al. (2007).

In the small band gap limit (usually large $\epsilon_{\infty}$) excitonic effects are expected to be week (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 $\sim$ 3 eV, the screened HSE06 functionals performs relatively well, due to the intrinsic screening incorporated in this type of functionals via the range-separation which allow 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 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 (Eg=1.0 eV Xu and Schoonen (2000)) and AgBiO3 (Eg=0.8 eV Mizoguchi and Woodward (2004)): the self-consistent procedure in these two cases end up with a metallic solution due to the a considerable overestimation of $\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 $\alpha$ as a fitting parameter, for instance by setting $\alpha$ to the value that reproduce the measured band gap, $\alpha_{\rm fit}$. This can be done by using the linear relation between Eg and $\alpha$ depicted in Fig. 1. The resulting values of $\alpha_{\rm fit}$, also reported in Tab. 1, ranges from 0.011 (LiNbO3) to 0.345 (GeO2); 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_{fit}$ and the corresponding self-consistent value adopted in scPBE0, $\alpha_{\rm sc}$, shown in Fig. 6(a), show strong deviations, with an overall tendency of scPBE0 to yield larger values of $\alpha$. This difference between scPBE0 and fitted-PBE0 is due to the fact that $\alpha$ and ${\varepsilon_{\infty}}$ does not scale linearly Alkauskas et al. (2011a); Moss (1950, 1985). This is graphically shown in Fig. 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$, $\alpha_{sc}$. As a result, by fitting $\alpha$ with respect to the band gap the MAPE associated with the band gaps is minimized but the description of the dielectric properties get 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 $\alpha$ obtained through the relation $\alpha=1/\epsilon_{\infty}$ allows for a generally good prediction of the band gap.

The non-linearity between Eg and $\epsilon_{\infty}$ has been pointed out first by Moss Moss (1950, 1985), who proposed the formula, known as Moss relation:

[TABLE]

which has been verified for many elements and semiconducting compounds, including oxidesReddy and Ahammed (1995), though deviations with respect to the measured values have been reported for band gap smaller than 1.4 eV Hervé and Vandamme (1995).

In Fig. 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$\approx$0.012.

IV Summary and Conclusions

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 $\alpha$ is determined by a self-consistent evaluation of ion-clamped dielectric constants using the PEAD method and by making use of the relation $\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 wit ha MAPE of about 10 %.

Importantly, for materials characterized by a 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 experiment only for scPBE0; this has been verified for SiO2, for which the excitonic binding energy is as large as 1.2 eV. On the other side, in the very small band gap limit, we found that scPBE0 could lead to metallic solutions. This was the case for PdO and AgBiO3, which are better described at PBE0 and HSE06 level.

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

In conclusion, scPBE0 represents an 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 error 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.

V ACKNOWLEDGMENTS

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