A Koopmans-compliant screened exchange potential with correct asymptotic behavior for semiconductors

TL;DR

This paper introduces a new screened exchange potential for semiconductors that is physically grounded, Koopmans-compliant, and accurately predicts band gaps, with a tunable parameter suitable for alloys.

Contribution

A novel Koopmans-compliant screened exchange potential with correct asymptotic behavior for semiconductors, adaptable to alloys and based on microscopic screening properties.

Findings

Accurately reproduces a wide range of band gaps.

Maintains a constant tunable parameter across isovalent substitutions.

Demonstrates proper asymptotic behavior and physical grounding.

Abstract

The performance of density functional theory depends largely on the approximation applied for the exchange functional. We propose here a novel screened exchange potential for semiconductors, with parameters based on the physical properties of the underlying microscopic screening and obeying the requirements for proper asymptotic behavior. We demonstrate that this functional is Koopmans-compliant and reproduces a wide range of band gaps. We also show, that the only tunable parameter of the functional can be kept constant upon changing the cation or the anion isovalently, making the approach suitable for treating alloys.

Tables

Table 1. Table 1: Fundamental 0K quasi-particle band gap of various materials. Given are the renormalization factor Z 𝑍 Z , and the resulting band gap with the screened exchange functional presented in this work.

	E${}_{}{}^{\text{target}}{}_G{}^{}$ (GW0)	ssZ${}_{\text{opt}}$ss	EG (present work)
LiF	14.2	0.62	14.2
ZnSe	2.7	0.62	2.8
ZnS	3.7	0.62	3.8
CdS	2.4	0.62	2.4
InAs	0.4	0.72	0.4
GaN	3.6	0.72	3.6
Ga2O3	5.0	0.72	5.1
CuGaS2	2.6	0.65	2.6
CuInS2	1.5	0.65	1.5
Si	1.2	0.53	1.1
3C-SiC	2.7	0.53	2.7
Diamond	5.8	0.53	5.8
SiO2	10.2	0.53	10.1

Table 2. Table 2: Z values from the G ​ W 𝐺 𝑊 GW calculation, averaged over the states filled by the N el subscript 𝑁 el N_{\rm el} electrons, for different materials and in comparison to the optimal values Z opt opt {}_{\text{opt}} used in Eq. ( 3 )

	ssGaNss	ssZnS ss	ss CdSss	ss LiFss	diamond
$\overline{Z_{\text{GW}}}$	0.80	0.77	0.77	0.77	0.81
Z${}_{\text{opt}}$	0.72	0.62	0.62	0.62	0.53

Table 3. Table 3: Fulfillment of the gKT for various materials. Δ Δ \Delta KS (HOMO) and Δ Δ \Delta KS (LUMO) describe the energetic position of the Kohn-Sham levels of the highest occupied molecular orbital for the neutral (HOMO), and of the lowest unoccupied molecular orbital (LUMO) for the +1 charge state of the defect. Δ Δ \Delta SCF is the electron removal energy. For a more detailed explanation, see footnote [46].

	$\Delta$KS (HOMO)-$\Delta$SCF	$\Delta$SCF - $\Delta$KS (LUMO)
GaAs	0.02	0.03
GaN	-0.04	-0.04
Ga2O3	-0.03	-0.04
Diamond	-0.05	0.05

Table 4. Table 4: Computational parameters used in the calculations. Given are the number of bands, the MP set used in the PBE calculation of the dielectric constant, the resulting dielectric constant ε b subscript 𝜀 𝑏 \varepsilon_{b} , and the effective number of electrons per unitcell.

	NBands	xx MP-setxx	xx ${\epsilon }_b$xx	N${}_{\text{el}}$/unitcell
ZnSe	600	20$\times$ 20$\times$20	8.29	4
ZnS	600	20$\times$ 20$\times$20	6.35	4
CdS	600	20$\times$ 20$\times$20	6.48	4
LiF	600	20$\times$ 20$\times$20	2.03	5
InAs	380	15$\times$ 15$\times$15	19.8	3
GaN	500	20$\times$ 20$\times$20	5.82	6
Ga2O3	500	20$\times$ 20$\times$20	3.97	24
CuGaS2	300	15$\times$ 15$\times$15	10.5	16
CuInS2	300	15$\times$ 15$\times$15	10.8	16
Si	60	12$\times$ 12$\times$12	13.63	4
3C-SiC	60	12$\times$ 12$\times$12	7.2	2
Diamond	60	12$\times$ 12$\times$12	5.95	4
SiO2	500	15$\times$ 15$\times$15	2.46	24

Equations

$$V^{\text{X}}(q)={\epsilon }^{-1}(q)V_{\text{HF}}^{\text{x}}(q).$$

$${\epsilon }^{-1}(q)=1+\left(\frac{1}{{\epsilon }_b}-1\right)\frac{1}{\cosh (q/\sigma )}.$$

$$\sigma =\frac{Z}{\log (2+\sqrt{3})}\sqrt{k_{\mathrm{TF}}^2\left(\frac{1}{{\epsilon }_b-1}+1\right)}.$$

$$k_{\mathrm{TF}}=4{\left(\frac{3N_{\mathrm{el}}}{\pi V}\right)}^{1/3},$$

Full text

Density functional theory (DFT) is the workhorse of electronic structure calculations in many areas of solid state physics and materials science. The exact form of the exchange-correlation potential is, unfortunately, not known, and the real predictive power of DFT-based methods is mainly limited by the quality of the applied approximate exchange functional. It has been established Perdew:82 that the exact functional provides a total energy, which is a piecewise linear function of the occupation numbers and has a derivative discontinuity at integer values. This is not case for the standard approximations (LDA: local density approximation, or GGA: generalized gradient approximation). As it is well known, these functionals are convex between, and have no derivative discontinuity at integer occupation numbers. Therefore, they underestimate the gap of semiconductors and lead to an artificial delocalization of defect states. In contrast, DFT with unscreened non-local Hartree-Fock (HF) type exchange leads to a strong overestimation of the band gap and to over-localization of defect states. Earlier, self-interaction correction schemes to LDA/GGA have been applied to remedy these problems Perdew:81 ; Baumeier:06 , but in the last decade, screened exchange Clark:10 ; Clark:10a and hybrid functionals Becke:96 ; PBE0 ; HSE ; Krukau:08 ; Zheng:11 ; Alkauskas:11 ; Skone:14 ; Chen:18 (that mix semi-local and non-local exchange) have emerged as useful alternatives.

Hybrid functionals utilize error compensation between semi-local and non-local exchange Deak:19a and their parameters are tuned semi-empirically. While connection between the parameters and the screening properties of the solid were being sought Marques:11 ; Skone:14 ; Chen:18 , the importance of Koopmans-compliance became ever more recognized Stein:10 ; Miceli:18 ; Nguyen:18 . Fulfillment of the generalized Koopmans’ theorem (gKT) is equivalent to the linearity of the total energy as a function of the fractional occupation numbers,Janak:78 ; Perdew:87 ; Perdew:97a and is the precondition for predicting the correct localization of one-electron statesLany:09 . So-called Koopman’s compliant functionals, (see, e.g., Refs. Perdew:81SI ; Dabo:10 ; Nguyen:18 ), fulfill the gKT by construction, which is appealing from a fundamental viewpoint. However, at present their applicability for solids suffers from practical limitations due to, e.g., limited possibility for k-point sampling.

Earlier we have shown that the two parameters of the Heyd-Scuseria-Ernzernhof (HSE) hybrid HSE can be tuned to reproduce the band gap and fulfill the gKT Deak:19a . This has allowed to obtain very accurate results for defects, but the parameters had to be optimized for every material separately Deak:10 ; Han:17 ; Deak:17 ; Deak:19 .

Current hybrid functionals, used for calculating defects in solids, like PBE0 PBE0 and HSE HSE but also the sX (screened exchange) Clark:10 functional, do not show the physically correct asymptotic behavior of the screening in semiconductors. As an alternative, we present here a novel screened exchange functional with proper asymptotic behavior. We will show that this functional reproduces the relative position of the band edge states and fulfills the gKT, i.e., the total energy shows the correct piecewise linearity as the function of the occupation numbers. These properties make the functional very well suited for defect calculations in semiconductors. The parameters of this functional are derived from the physical properties of the screening, and it contains only one adjustable parameter with a value being constant upon substitution of the cation or the anion in a given class of semiconductors. (As we will show, even that parameter can be obtained approximatively from first-principles.)

We start from the ansatz

[TABLE]

Here $V_{\text{HF}}^{\text{x}}(q)$ is the Hartree-Fock-type non-local exchange potential in wave vector space, and the model screening function is

[TABLE]

Correlations are added on the GGA level in the PBE (Perdew, Burke, Ernznerhof) approximation PBE . The choice of an $1/\cosh$ behavior is guided by a well known result of Green’s function theory, that the $q\rightarrow\infty$ behavior should be exponential Banyai:98 ; Gartner:02 while the $q\rightarrow~{}0$ behavior has to be quadratic Gartner:00 . Our choice for $\varepsilon$ ensures that $V^{\text{X}}(r)$ has the proper 1/r behavior for $\vec{r}\rightarrow\infty$ and at the same time approaches the correct pure Coulomb limit at $\vec{r}\rightarrow 0$.

To determine the screening length $\sigma$, we follow Refs. Cappellini:93 ; Shimazaki:08 in deriving a static approximation to the dielectric function in random phase approximation (RPA). This results in a Lorentzian dependence of $\varepsilon^{-1}(q)$ with a screening length $\kappa$. We than choose the same FWHM in our $\cosh$-ansatz as the RPA would give for its Lorentzian screening. This procedure results in

[TABLE]

The TF wave vector $k_{\rm TF}$ can be expressed as Cappellini:93 ; Shimazaki:08 ; Shimazaki:10

[TABLE]

with the cell volume $V$. We will elaborate on the choice of the effective number of electrons $N_{\rm el}$ per unit cell and on the renormalization factor $Z$ below. Our exchange potential can be seen as a static approximation to a $GW$ calculation with a model screening function. In it’s idea, the method is similar to the sX functional Clark:10 , however, the inclusion of the correct limits and q-dependence marks a substantial improvement. In Fig. 1, we compare our model screening function to the ones used in the HSE and sX functionals. The sX clearly neglects the background screening $\varepsilon_{b}$ and uses a Lorentzian q-dependence, while the HSE approach differs substantially and compensates by mixing with PBE exchange.

To apply this theoretical model to semiconductor materials we have to determine the dielectric constant $\varepsilon_{b}$, the effective number of electrons $N_{\rm el}$ and the renormalization factor $Z$. The value of $\varepsilon_{b}$ is taken from an independent particle calculation of the optical properties at the GGA (PBE) level. This is in the spirit of the GW0 approximation Shishkin:07 , which in most cases reproduces the band gaps better than fully self-consistent GW schemes (without vertex corrections).

It is well known, that energetically deep lying electronic states contribute insignificantly to the screening, as they are strongly localized, their spatial overlap is weak and their energetic distance to the top of the VB is high. Therefore, as an approximation, we disregard them in the screening completely. The ”effective number of electrons”, $N_{\rm el}$, is hence determined as follows. We evaluate which states primarily form the highest valence band and choose the corresponding number of electrons from the constituting atomic states 111technically from Vasp’s POTCAR file. As an example, in GaN, the top of the VB is made up of nitrogen 2p orbitals. In a nitrogen atom, these hold 3 electrons, hence for the primitive cell, containing 2 nitrogen atoms, $N_{\rm el}=6$. $N_{\rm el}$ should, of course, be increased for a supercell, according to the number of atoms involved. It should be noted, that neither $Z$ nor $N_{\rm el}$ enter $\sigma$ independently, but only in combination $ZN_{\rm el}^{1/3}$. Therefore the choices for $N_{\rm el}$ and $Z$ are not unique, as they lead to the same functional if the same $\sigma$ is generated.

To understand the role of $Z$, let us recall that in $GW$ theory the screened potential $W$ is usually calculated in RPA. The polarizability $P$ that controls $W$ via the Dyson equation $W=V+VPW$ (V is the unscreened Coulomb potential), is given by the product of two Green’s functions, $P=GG$ Baym:62 . We base our approach on the same framework, with the role of $W$ taken by $V_{X}$. One can define a renormalization factor $Z$ that describes the fraction of spectral weight attributed to the main quasi-particle peak within the full Green’s function (GF) $G$. In our approximation only this fraction $Z$ of the full GF $G$ is considered in determining the polarizability $P$. Via the RPA ($P=GG$ ), this leads (with the usual steps used to derive static screening Schaefer:02 ) to the occurence of $Z$ in Eq. (3). Physically, this means that we assume significant contributions to the screening only by the main quasi-particle peak of the GF G.

Obviously, this value of $Z$ varies for the states around the valence band edge and between different bands $\lambda$ in the $GW$ approximation, i.e. $Z_{\text{GW}}$=$Z_{\lambda}(\vec{k})$. As an approximation we average Z over the bands that hold the $N_{\rm el}$ electrons to get an average value $\overline{Z_{\text{GW}}}$. However, as effects beyond a static approximation and beyond the RPA clearly influence the screening, we employ a pragmatic approach, using $\overline{Z_{\text{GW}}}$ as a starting point, and modifying it to reproduce the band gap, yielding Z${}_{\text{opt}}$. As we will show below, both the average value $\overline{Z_{\text{GW}}}$ and hence also the optimal value Z${}_{\text{opt}}$ are constant for isovalent substitution of atoms.This allows to tackle semiconductor alloys with a common $Z$ value.

We have implemented the exchange functional as discussed above into the Vienna Ab initio Simulation Package, Vasp 5.4.4 VASP:3 ; VASP:4 , using the projector augmented wave method and treating the semi-core d-states as part of the valence shell. The modified Vasp source code can be made available to certified owners of a Vasp user license. Calculations on the unit cell of the bulk materials where performed using a 6$\times$6$\times$6 $\Gamma$-centered Monkhorst-Pack Monkhorst:76 grid. For defect calculations, 64 (GaAs), 96 (GaN), 160 (Ga2O3), and 512 (diamond) atom supercells were used, applying the $\Gamma$-point approximation. The defect geometries were fully relaxed. A 450eV (900eV) cutoff was applied for the expansion of the wave functions (charge density). The computational cost for practical calculations of defects and adsobates is on par with that of DFT calculations with the PBE0 exchange potential. Charge corrections for the total energy were performed by the method of Freysoldt, Neugebauer, and van der Walle Freysoldt:09 , while localized defect levels were corrected using the formula derived by Chen and Pasquarello Chen:13 . GW0 calculations were performed on top of PBE calculations with 1000 bands and a 6$\times$6$\times$6 $\Gamma$-centered Monkhorst-Pack grid to provide a reference for the band gap and a starting value for $Z$. In the present calculations, the experimental lattice constants are used for consistency with the literature.

For the purpose of this letter, we study 13 semiconductor compounds with band gaps ranging between 0.5 and 14 eV. The band gap of most of these materials is systematically underestimated by standard HSE06Deak:19a . These compounds represent different lattice structures and vary in composition and in the degree of ionicity, therefore, they can serve as a representative set for testing the method.

Table 1 shows the calculated band gaps, in comparison with the GW0 results. The optimal value of the scaling factor Z${}_{\text{opt}}$ is also given. We note that our value for the quasi-particle band gap, as given in Table 1, is determined to include neither excitonic nor polaronic renormalizations. Both the electron-hole Coulomb interaction and the electron-phonon interaction lower the band gap (resulting in the measurable optical gap) and are not included here Deak:17 ; Deak:19 . An overview of the comparison between the band gaps of the proposed method and the GW0 ones is given in Fig. 2. We find excellent agreement, with a mean average error of 0.05eV and a mean average percentage error of only 3%.

A striking feature of these results is that, in a particular sub-class of materials, the optimal value Z${}_{\text{opt}}$ remains the same when the cation or the anion is replaced. The reason for that is a similar behavior of the average $GW_{0}$ value, $\overline{Z_{\text{GW}}}$, as shown in Table 2. We can infer that Z${}_{\text{opt}}$ remains constant for all II/VI, III/V, and IV/IV compounds. Deviations between the optimal and average $GW$ values most probably stems from the use of a static approximation to the GW self energy in Eq. (2) and hence they reflect the influence of frequency dependent screening. The biggest deviation is found for silicon (1.23eV vs. 1.05eV), where it is known Thygesen:13 that static approximations tend to underestimate the band gap. This is also the case for our approximation, where the agreement for silicon is worse than for SiC or diamond.

The results above mean that our method allows to treat alloys of isovalent elements without parameter retuning. It is important to note that, as we will discuss below, not only the band gap is reproduced, independent of the alloy composition, but also the gKT is satisfied.

We want to point out, that, in an alloy, $\varepsilon_{b}$ has to be recalculated for the given composition, while $N_{\rm el}$ is determined by the prevailing crystal structure. In principle, the common Z value would allow to treat heterostructures between compounds of such sub classes but presently the spatial variation of the dielectric constant is not taken into account. Work in that direction is in progress.

As discussed in previous works Deak:17 ; Deak:19 , not only the minimum band gap, but the band edges positions over the entire Brillouin-zone have to fit the results of a GW0 calculation, in order to provide a proper description of defect levels, as they can be generated by a superposition of all band edge states. This condition is also satisfied with our approach, as shown on the example of GaN in Fig. 3.

To correctly describe the energetic position and the localization of defect states, the gKT has to be fulfilled,Deak:19a ; Lany:09 i.e., the total energy has to be linear with respect to the fractional occupation numberNote2 . We have tested this criterion on the (0/+1) transition levels of the antisite pair in GaAs, the nitrogen vacancy in GaN, the oxygen vacancy in Ga2O3, and the B${}_{\text{C}}$ substitutional in diamond. As shown in Table 3, we find an excellent fulfillment of the gKT.

It has been shown Marini:03 that issues in the description of optical properties within linear response TDDFT (with the Cassida equation) are closely connected to the description of screening in the exchange kernel.

Hence, we investigate the performance of our proposed method in comparison to state-of-the-art many-body perturbation theory calculations with the GW0+BSE scheme in Fig. 4 for the cases of GaN and LiF with moderate and strong excitonic effects in the optical response, respectively. Apart from a slight underestimation of the oscillator strength, very good agreement is found for a fraction of the computational effort.

In conclusion, we have presented a functional based on exact exchange with improved screening for semiconductors. As it is based on the correct asymptotic limits of the exchange potential, its parameters can be derived from physical principles, to a large extent eliminating the need for ”tuning” which is a common procedure with current hybrid functionals. Because of the correct piecewise linearity of the total energy, this functional reproduces the relative position of the band edge states and fulfills the generalized Koopmans’ theorem. This is of high significance in various application areas. For example, in studying photo-assisted reactions on semiconducting catalysts, it is critical to reproduce the position of the band edges, which measure the chemical potential of photo-generated holes and electron Zhao:14 ; Deak:16 . In general, the correct gap is the starting point to determine the optical and transport properties as well. Besides the reproduction of the gap, the fulfillment of the gKT is a condition to accurately predict the localization and energy of defects states in semiconductors, which is a prerequisite in the successful identification of defects, which influence device behavior in micro/opto-electronics and photovoltaics Deak:19a ; Deak:10 ; Deak:17 ; Deak:19 . A great advantage of our approach is that it can be used to describe alloys without any re-tuning, as the single parameter of this approach is transferable between semiconductors of a similar type, when the cation or anion is replaced. In addition, it is also suitable for the description of excitonic effects within linear response TDDFT, for a fraction of the cost of GW0+BSE.

Appendix A Parameters

In this appendix, we present computational details about the calculation of the optical properties (MP sets and number of unoccupied states), and provide the effective number of electrons $N_{el}.$. These results are summarized for all materials in Table 4. It should be noted that the results for the dielectric constant is strongly dependent on the number of bands, as also known from GW-type calculations.