Dielectric dependent hybrid functionals for heterogeneous materials

TL;DR

This paper introduces a parameter-free dielectric-dependent hybrid functional that accurately models electronic properties of heterogeneous interfaces, surfaces, and low-dimensional materials by self-consistently determining a local dielectric function.

Contribution

It generalizes self-consistent hybrid functionals to heterogeneous systems using a spatially varying dielectric, enhancing accuracy without adjustable parameters.

Findings

Accurately predicts band gaps and dielectric constants of 3D and 2D materials.

Provides reliable band offsets for interfaces.

Achieves accuracy comparable to GW calculations.

Abstract

We derive a dielectric-dependent hybrid functional which accurately describes the electronic properties of heterogeneous interfaces and surfaces, as well as those of three- and two-dimensional bulk solids. The functional, which does not contain any adjustable parameter, is a generalization of self-consistent hybrid functionals introduced for homogeneous solids, where the screened Coulomb interaction is defined using a spatially varying, local dielectric function. The latter is determined self-consistently using density functional calculations in finite electric fields. We present results for the band gaps and dielectric constants of 3D and 2D bulk materials, and band offsets for interfaces, showing an accuracy comparable to that of GW calculations.

Tables

Table 1. Table 1: The electronic dielectric constants ( ϵ ∞ subscript italic-ϵ \epsilon_{\infty} ) of three dimensional materials obtained from PBE and spatial-resolved dielectric dependent hybrid functional (DDH) calculations [Eq. ( 10 )], compared with the results of Ref.  Skone  et al. ( 2014 ) and experiment. All calculations (PBE and DDH) were carried out using ONCV pseudopotentials  Schlipf and Gygi ( 2015 ) and by sampling the Brillouin zone with the Γ Γ {\Gamma} point. The number of atoms or units used in the supercell calculations are indicated as subscripts for each solids.

	PBE	DDH	Ref. Skone et al. (2014)	Exp.
Si	12.46	11.80	11.76	11.9 YU and Cardona (2010)
SiC	6.86	6.49	6.50	6.52 YU and Cardona (2010)
AlP	8.08	7.57	7.23	7.54 YU and Cardona (2010)
Diamond	5.77	5.58	5.61	5.70 YU and Cardona (2010)
MgO	3.26	2.99	2.81	2.96 Lide (2009)
LiCl	2.93	2.77	2.77	2.70 Van Vechten (1969)
Ar	1.73	1.66	1.66	1.66 Sinnock and Smith (1969)
Ne	1.29	1.25	1.21	1.23 Schulze and Kolb (1974)

Table 2. Table 2: The fundamental energy gaps (eV) of three dimensional materials obtained from PBE and spatial-resolved dielectric dependent hybrid functional (DDH) calculations [Eq. ( 10 )], compared with the results of Ref.  Skone  et al. ( 2014 ) and experiment. All calculations (PBE and DDH) were carried out using ONCV pseudopotentials  Schlipf and Gygi ( 2015 ) and by sampling the Brillouin zone with the Γ Γ {\Gamma} point. The number of atoms or units used in the supercell calculations are indicated as subscripts for each solids. In columns 4 and 5 we report calculations with a constant mixing fraction (See Eq. ( 8 )), α = 1 / ϵ ¯ 𝛼 1 ¯ italic-ϵ \alpha=1/\bar{\epsilon} , and α = 1 / ϵ ∞ 𝛼 1 subscript italic-ϵ \alpha=1/\epsilon_{\infty} respectively. The zero-phonon renormalization (ZPR) is reported when available from experiment  Chen  et al. ( 2018b ) .

	PBE	DDH	$\alpha =1/\overline{ϵ}$ 111Hybrid functional calculation with $\alpha =1/\overline{ϵ}$ where $\overline{ϵ}$ is the bulk average of $ϵ\left(𝒓\right)$: values reported in Table 1.	$\alpha =1/{ϵ}_{\infty }$222Hybrid functional calculation with $\alpha =1/{ϵ}_{\infty }$ where ${ϵ}_{\infty }$ is from Ref. Skone et al. (2014): values reported in Table 1.	Ref. Skone et al. (2014)	ZPR	Exp.
Si	0.603	1.00	1.01	1.01	0.99	0.06	1.17 Kittel (2004)
SiC	1.38	2.35	2.35	2.35	2.29	0.11	2.39 Choyke et al. (1964)
AlP	1.56	2.27	2.28	2.32	2.37	0.02	2.51Monemar (1973)
Diamond	4.17	5.48	5.54	5.53	5.42	0.37	5.48 Clark C. D. et al. (1964)
MgO	4.78	7.70	8.08	8.30	8.33	0.53	7.83 Whited et al. (1973)
LiCl	6.47	9.38	9.56	9.56	9.62	0.17	9.40 Baldini and Bosacchi (1970)
Ar	8.70	13.93	14.34	14.34	14.67		14.2 Schwentner et al. (1975)
Ne	11.62	20.60	22.38	22.72	23.67		21.7 Schwentner et al. (1975)

Table 3. Table 3: Band offsets (eV) computed at different levels of theory (using the PBE functional and the functional of Eq. ( 10 ), with the procedure of Fig.  5 ) for the silicon-silicon nitride interface, compared with experiment (from Anh Pham  et al. ( 2013 ) and the references therein)

Si/Si3N4	PBE	DDH	Exp.
Conduction band offset	1.2	1.9	1.83 $-$ 2.83
Valence band offset	0.7	1.3	1.5 $-$ 1.78

Equations

$$\overline{ϵ}$$

$${\bar{\chi }}^0(\mathbf{r},{\mathbf{r}}^{\prime })=\sum _n{\lambda }_n{\varphi }_n(\mathbf{r}){\varphi }_n^{\ast }({\mathbf{r}}^{\prime }),$$

$${\mathcal{M}}_{ij}:={\int }_{\mathbf{r}\in {\Omega }_S}d\mathbf{r}{\varphi }_i^{\ast }(\mathbf{r}){\varphi }_j(\mathbf{r}),$$

$$\mathcal{M}\cdot V_m=w_m{V}_m,w_m\in [0,1],$$

$$w_m=\frac{{\int }_{\mathbf{r}\in {\Omega }_S}d\mathbf{r}|{\varphi }_m^{\text{loc}}(\mathbf{r}){|}^2}{\int d\mathbf{r}|{\varphi }_m^{\text{loc}}(\mathbf{r}){|}^2}.$$

$${\mathcal{F}}_A$$

$${\mathcal{F}}_B$$

$${\mathcal{F}}_I$$

$$\text{eig}({\bar{\chi }}_A^0)\cup \text{eig}({\bar{\chi }}_B^0)\cup \text{eig}({\bar{\chi }}_I^0)\simeq \text{eig}({\bar{\chi }}^0).$$

$$W(\mathbf{r},{\mathbf{r}}^{\prime })=\alpha (\mathbf{r},{\mathbf{r}}^{\prime })\frac{1}{|\mathbf{r}-{\mathbf{r}}^{\prime }|},$$

$$E_x$$

$$\alpha (\mathbf{r},{\mathbf{r}}^{\prime })\simeq \frac{1}{\sqrt{ϵ(\mathbf{r})ϵ({\mathbf{r}}^{\prime })}}.$$

$$E_{xc}$$

$$F(\mathbf{E},[\rho ])$$

$$U(\mathbf{D},[\rho ])=E_{\text{KS}}[\rho ]+\frac{1}{8\pi }\int d\mathbf{r}(\mathbf{D}-4\pi \mathbf{P}{)}^2,$$

$${ϵ}_{\alpha \beta }={\delta }_{\alpha \beta }+4\pi \frac{\partial P_{\alpha }}{\partial E_{\beta }},({ϵ}^{-1}{)}_{\alpha \beta }={\delta }_{\alpha \beta }-4\pi \frac{\partial P_{\alpha }}{\partial D_{\beta }},$$

$$\Delta \mathbf{P}=\frac{-e}{\Omega }\sum _{i=1}^{N_{\text{s}}}\Delta {\mathbf{r}}_c^i.$$

$$\Delta \mathbf{P}(\mathbf{r})=-e\sum _{i=1}^{N_{\text{c}}}{N}_i\Delta {\mathbf{R}}_c^i\delta (\mathbf{r}-{\mathbf{R}}_c^i),$$

$$\Delta \mathbf{P}=\frac{1}{\Omega }{\int }_{\Omega }\Delta \mathbf{P}(\mathbf{r})d\mathbf{r}.$$

Full text

I Introduction

Density Functional Theory (DFT) was first applied to compute the structural and electronic properties of condensed systems more than 35 years ago Yin (1982); Parr and Yang (1995); Martin (2004), using the local density approximation Ceperley and Alder (1980); Perdew and Zunger (1981) of the exchange and correlation (xc) energy functional. Approximately ten years later, when gradient corrected approximations (GGA) Langreth and Perdew (1980); Langreth and Mehl (1983); Hu and Langreth (1985); Perdew et al. (1992, 1996a); Zhang and Yang (1998) for the xc energy were derived, DFT was adopted for some molecular investigations by the quantum chemistry community. Shortly after the first GGA molecular calculations, hybrid functionals were proposed Becke (1993a, b, 1996); Perdew et al. (1996b); Jaramillo et al. (2003); Chai and Head-Gordon (2008); Maier et al. (2018) and most DFT applications for finite systems, which use localized basis sets, have been carried out with hybrid functionals Bauschlicher (1995); Ren et al. (2012), most notably B3LYP Lee et al. (1988); Becke (1993b); Stephens et al. (1994). These are functionals where the exchange energy is defined as a linear combination of exact (Hartree-Fock) and local exchange Ghosh et al. (2018). The condensed matter physics community adopted hybrid functionals later than the quantum chemistry community, due to computational difficulties in evaluating the Hartree-Fock (HF) exchange energy using plane wave (PW) basis sets; these are the basis set of choice in most of the codes used for materials Gonze et al. (2002, 2009); Gygi (2008); Giannozzi et al. (2009); VAS (2018); CP2 (2018); CPM (2018); Lejaeghere et al. (2016), although periodic DFT codes using localized basis sets are also in use Soler et al. (2002); Dovesi et al. (2018); Frisch et al. (2016); Blum et al. (2009). The difficulties in evaluating HF exchange in PW basis sets have now been largely overcome, with the advent of fast algorithms based on bisection techniques Gygi (2009); Gygi and Duchemin (2013); Dawson and Gygi (2015); Lin (2016) or maximally localized Wannier functions Wu et al. (2009); DiStasio et al. (2014). Nevertheless periodic DFT calculations with hybrid functionals and PW basis sets remain substantially heavier, from a computational standpoint, than local or semi-local DFT calculations. The functionals PBE0 Adamo and Barone (1999) and HSE Heyd et al. (2003); Heyd and Scuseria (2004); Krukau et al. (2006) are among the most popular hybrid functionals used for condensed systems, and lately dielectric dependent hybrid functionals Skone et al. (2014, 2016); Brawand et al. (2017); Gerosa et al. (2018a) have been increasingly used to predict structural and electronic properties of solids Skone et al. (2014); Ferrari et al. (2015); Skone et al. (2016); Seo et al. (2016); Chen et al. (2018a, b); Gerosa et al. (2015, 2018b, 2018a); Kronik and Kümmel (2018); Manna et al. (2018) and liquid Gaiduk et al. (2016); Pham et al. (2017); Gaiduk et al. (2018) and of several molecules Skone et al. (2016); Brawand et al. (2016, 2017). Another category of orbital dependent functionals recently proposed is that of Koopmans-compliant functionals, used for both molecules and solids Stein et al. (2010); Borghi et al. (2014); Nguyen et al. (2018).

A drawback of most of the functionals mentioned above is that while they work well for certain classes of homogeneous systems, e.g. solids, they are usually not as accurate for heterogeneous systems, e.g. surfaces and interfaces, where the dielectric screening of different portions of the system differ substantially. For heterogeneous semiconductors, Shimazaki et al. Shimazaki and Nakajima (2015) introduced an estimator of the electrostatic environment surrounding the atoms in a semiconductor leading to the definition of position-dependent atomic dielectric constants. For solid/solid interfaces, Borlido et al. Borlido et al. (2018) introduced a non-local mixing fraction, based on an estimator of a local dielectric function that contains parameters to be evaluated with system-dependent fitting procedures.

In this work, we propose a hybrid functional that describes equally accurately three- and two-dimensional solids, as well as surfaces and interfaces, and which is derived entirely from first principles, with no need to define any adjustable parameter. The functional is based on an approximation of the screened Coulomb interaction using a local dielectric function, which is derived from first principles by minimizing a dielectric enthalpy functional. We first discuss (Section II) the foundation of dielectric disentanglement by showing that the dielectric screening of a system composed of two subsystems interfaced with each other, may be decomposed into the screening of the two subsystems plus an interfacial contribution. The disentanglement is carried out using a localized representation of the eigenvectors of the dielectric matrix, obtained using bisection techniques originally proposed for the eigenfunctions of Kohn-Sham Hamiltonians Gygi (2009). Our results on dielectric decomposition are used to justify the definition of a local, spatial dependent dielectric function (Section II), which in the bulk portion of the subsystems coincides with their respective dielectric constants. We then use this local dielectric function to define a dielectric hybrid functional for heterogeneous systems (Section III); the functional is derived from first principles, without any adjustable parameter, by carrying out calculations in finite electric field. Finally we present results for 3D and 2D solids in Section V.1 and for surfaces and interfaces in Section V.2, with focus on the calculations of band gaps, dielectric constants and band-offsets.

II Spatial disentanglement of dielectric spectra

In this section, we address the following question: can the dielectric matrix of an heterogeneous system (composed, e.g. of two solids or a liquid and a solid) be expressed in terms of the dielectric matrices of the subsystems? For simplicity we restrict our attention to a system of volume $\Omega$ composed of two subsystems, A and B interfaced with each other and we consider a single interface between A and B. We address the question by writing a spectral decomposition of the dielectric matrices of the heterogeneous system and of A and B, and then we use bisection techniques Gygi (2009) to localize the eigenvectors of the dielectric matrices in desired regions of space.

According to linear response theory, the density-density response function $\chi$ and the irreducible polarizability $\chi^{0}$ are related to the dielectric matrix ($\epsilon$) of the system by the following equation:

[TABLE]

where the bar in Eq. (1) indicates that the functions have been symmetrized with respect to the Coulomb potential (see, e.g. Ref. Govoni and Galli (2015)). We represent $\bar{\chi}^{0}$ using its spectral decomposition,

[TABLE]

where $\phi_{n}$ and $\lambda_{n}$ are eigenvectors and eigenvalues, respectively. In the following we focus on static dielectric responses.

Fig. 1 shows the eigenvalues of $\bar{\chi}^{0}$ (left panels) for two representative interfaces, H-Si/H2O and Si/Si3N4, one where the two subsystems are non-covalently bonded and one where there are covalent bonds at the interface (the geometry of the model slabs and how they were obtained are described in the SM SM ). The square moduli of selected eigenvectors projected in the direction perpendicular to the interface ($z$) and their corresponding eigenvalues (dots) are shown on the right and left panels of Fig. 1, respectively. For both surfaces, we see that some eigenmodes are predominantly localized on one side of the slabs while other modes, especially those corresponding to $|\lambda_{i}|\to 0$ (green and red curves) are localized over the entire slab.

In order to express response functions of the entire system in terms of those of the subsystems, we represent the dielectric matrix in terms of localized functions, instead of eigenfunctions.

We first define two subsystems using the projection of the electronic charge density on the $z$-axis perpendicular to the interface, as illustrated in Fig. 2. The use of the charge density to define regions A and B introduces a certain degree of arbitrariness, as a criterion is required to determine charge density minima, in correspondence of which interface planes are defined. Such a criterion is system dependent. While the charge density is used in this section to define interfacial planes for the purpose of illustrating the concept of disentanglement of the dielectric response, it will not be used in practical calculations. As we will see in Section IV, a general, system independent procedure can be defined to compute local dielectric functions.

After partitioning the full system into subsystems using the charge density, we obtain a set of localized functions from the set of eigenvectors $\phi_{i}$ by constructing and diagonalizing the filtered overlap matrix $\mathcal{M}$,

[TABLE]

where $w_{m}$ and $V_{m}$ are eigenvalues and eigenvectors of $\mathcal{M}$, and $\Omega_{S}$ is the volume of either subsystem A or B as defined using the electronic charge density (see Fig. 2). The set of eigenvectors of $\mathcal{M}$ provides the transformation matrix from the set of $\phi_{i}(\bm{r})$’s to a set of localized orbitals. The eigenvalues $w_{m}$ represent the weights of the localized orbital $\phi^{\text{loc}}_{m}(\bm{r})$ within the subspace $\Omega_{S}$:

[TABLE]

If $w_{i}\simeq 1$, $\phi_{i}^{\text{loc}}$ is localized on $\Omega_{S}$; if $w_{i}\simeq 0$, $\phi_{i}^{\text{loc}}$ is localized on $\Omega-\Omega_{S}$. We classify the $\phi_{i}^{\text{loc}}(\bm{r})$’s into three subsets:

[TABLE]

where $w_{thr}$ is a chosen localization threshold that can be systematically varied to verify the robustness and convergence of the localization procedure (it was chosen to be 0.01 in the examples shown in the figures).

Fig. 3 displays the weights $w$ and the square moduli of localized basis functions for the H-Si/H2O and Si/Si3N4 interfaces: we found that most of the basis functions are localized in one of the two subsystems, with the rest of them localized near the interface.

After obtaining the localized basis set $\mathcal{F}$ ($\mathcal{F}_{A}\cup\mathcal{F}_{B}\cup\mathcal{F}_{I}$), we expressed the matrix elements of $\bar{\chi}^{0}$ as $\bar{\chi}^{0}=\bar{\chi}^{0}_{A}+\bar{\chi}^{0}_{B}+\bar{\chi}^{0}_{I}+\bar{\chi}^{0}_{\text{off-diag.}}$, where $\chi^{0}_{\text{off-diag.}}$ includes all the off-diagonal blocks representing the coupling between the two subsystems. By diagonalizing $\bar{\chi}^{0}_{A}$, $\bar{\chi}^{0}_{B}$ and $\bar{\chi}^{0}_{I}$ in the respective subspaces $\mathcal{F}_{A}$, $\mathcal{F}_{B}$, and $\mathcal{F}_{I}$ defined in Eq. (II), we found that the response of the whole system can be disentangled into contributions from the subsystems, i.e. we found that for all systems studied here:

[TABLE]

Fig. 4 shows decomposed spectra [Eq. (6)] compared with the spectrum of the whole system. It is seen that [$\text{eig}(\bar{\chi}^{0}_{A})\cup\text{eig}(\bar{\chi}^{0}_{B})\cup\text{eig}(\bar{\chi}^{0}_{I})$] and $\text{eig}(\bar{\chi}^{0})$ give very similar results, with small differences in the low eigenvalue regions, due to the neglect of the elements of $\chi^{0}_{\text{off-diag.}}$. As expected neglecting these elements is a better approximation for aqueous interfaces than for the Si-Si3N5 interface, where covalent bonds are formed. Therefore, we conclude that the dielectric screening of the whole slab may be approximated as the sum of contributions from the subsystems plus an interfacial dielectric screening contribution.

The results of this section indicate that it is reasonable to approximate the screening of the entire slab by a local dielectric function $\epsilon(\bf r)$, a smooth function expected to describe accurately the screening of the two separate subsystems in their respective bulk regions. We will see in the next section that these assumptions lead to a definition of a generalized dielectric hybrid functional which yields accurate band gaps and dielectric constants for 2D and 3D systems and band offsets for complex interfaces.

We now turn to describing a procedure to obtain $\epsilon(\bf r)$ which does not rely on the definition of an interface plane based on the electronic charge density, nor on any parameters defining subsystems A and B.

III Dielectric-dependent hybrid functionals

The results on dielectric disentanglement described in the previous section led to the idea of defining a local dielectric function whose limiting values in the two subsystems is expected to coincide with the dielectric constants of the respective bulk subsystems. Such a local dielectric function can then be used to generalize the hybrid functionals introduced in Ref. Skone et al. (2014). In the following, we define the local dielectric function from first principles.

In Ref. Skone et al. (2014), the screened Coulomb interaction in a homogeneous system is approximated as $W(\bm{r},\bm{r}^{\prime})=\frac{1}{\epsilon_{\infty}|\bm{r}-\bm{r}^{\prime}|}$, where $\epsilon_{\infty}$ is the macroscopic static dielectric constant. This approximation is used in the definition of a hybrid functional similar to PBE0 but with mixing fraction $\alpha=1/\epsilon_{\infty}$ instead of 0.25. Several authors have suggested using $\alpha$ as an adjustable parameter to reproduce the experimental band gap of solids Alkauskas et al. (2008, 2011); Shimazaki and Asai (2008); Marques et al. (2011); Refaely-Abramson et al. (2013); Conesa (2013); Gerosa et al. (2018a).

Following the definition of exchange in Ref. Skone et al. (2014); Borlido et al. (2018), if we write the screened Coulomb potential as

[TABLE]

the exchange energy of the entire system takes the following form Borlido et al. (2018)

[TABLE]

We assume that the function $\alpha(\bm{r},\bm{r}^{\prime})$ is a simple separable function of $\epsilon(\bm{r})$ and $\epsilon(\bm{r}^{\prime})$, with $\alpha(\bm{r},\bm{r})=\epsilon(\bm{r})$ and we write:

[TABLE]

We then arrive at the following ansatz for the exchange and correlation energy:

[TABLE]

where we have chosen the PBE approximation to represent the local part of the exchange and correlation energy.

The exchange-correlation functional defined in Eq. (10) is similar, in spirit, to the local functional proposed in Ref. Borlido et al. (2018). However we emphasize two important conceptual and practical differences: we have provided a theoretical justification of Eq. (10) based on the decomposition of the screened Coulomb interaction into that of subsystems and an interfacial region. Next we show that $\epsilon(\bm{r})$ may be obtained from first principles by carrying out calculations in finite field, eliminating the need to tune any arbitrary parameter, or adopt any fitting, system-dependent procedure, which are necessary instead in the formalism of Ref. Borlido et al. (2018).

IV Self-consistent determination of local dielectric functions using a finite field approach

Here we describe a finite field approach to compute $\epsilon(\bm{r})$. In general, the macroscopic dielectric tensor of any condensed system can be obtained by carrying out calculations in a finite electric field and by minimizing the functional Umari and Pasquarello (2002); Souza et al. (2002); Stengel et al. (2009):

[TABLE]

where $\int\bm{E}\cdot\bm{r}\rho(\bm{r})d\bm{r}$ is called the electric enthalpy, and $E_{\text{KS}}$ is the Kohn-Sham energy of the system. Alternatively one could minimize the functional:

[TABLE]

where $\bm{D}=\bm{E}+4\pi\bm{P}=\bm{\epsilon}\cdot\bm{E}$, and $\bm{P}$ is the polarization of the system; the components of the dielectric tensor $\epsilon$ are:

[TABLE]

where $\alpha$ and $\beta$ are Cartesian coordinates. In periodic systems, the induced polarization can be computed from the shift of the centers of the Wannier functions ($\Delta\bm{r}^{i}_{c}$) of the unperturbed system when an electric field is applied King-Smith and Vanderbilt (1993); Stengel and Spaldin (2006). For a homogeneous system of $N_{\text{s}}$ occupied states, the average change in macroscopic polarization is given by

[TABLE]

This allows us to define a spatial dependent polarization for heterogeneous systems (e.g., 2D materials, surfaces and interfaces):

[TABLE]

where $N_{\text{s}}$ Wannier centers have been grouped in $N_{\text{c}}$ clusters: $N_{i}$ is the number of Wannier centers in the $i$-th cluster, $\bm{\Delta}R_{c}^{i}=\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}\bm{\Delta}r_{c}^{j}$ is the shift of the center of the $i$-th cluster induced by the applied electric field. In practical calculations the $\delta$-function is replaced by a Gaussian function of finite width equal to the average of the spreads of the corresponding Wannier orbitals belonging to the same cluster. We note that $\Delta\bm{P}$ entering Eq. (14) can be obtained from $\Delta\bm{P}(\bm{r})$ using the following relation:

[TABLE]

The spatial dependence of $\epsilon$ is then defined by the spatial dependence of the polarization, as given in Eq. (15).

We computed the local dielectric function $\epsilon(\bf r)$ by minimizing the electric enthalpy [Eq. (11)] with the Kohn-Sham energy defined using the exchange correlation functional of Eq. (10). The minimization is carried out using a finite field approach, as implemented in the Qbox code Gygi (2008, 2018). The function $\epsilon(\bf r)$ is computed self-consistently. The whole procedure is schematically shown in Fig. 5. At the first iteration we perform a DFT calculation at the PBE level [$\alpha(\bm{r},\bm{r}^{\prime})=0$]. At the second iteration we set $\epsilon(\bm{r})=\epsilon^{\text{PBE}}(\bm{r})$ in Eq. (9) and repeat the process until $\epsilon(\bm{r})$ and the total energy are converged.

V Validation of self-consistent hybrid functionals for 3D & 2D materials, surface & interfaces

V.1 Three-dimensional and two-dimensional materials

Fig. 6 shows the band gap at each iteration for bulk Si and a 3C-(SiC) computed using supercells with 512 atoms and the $\Gamma$ point to sample the Brillouin zone (the corresponding $\epsilon(z)$ [average of $\epsilon(\bm{r})$ in the xy plane] are shown in the supplementary information [SI]). In both cases, calculations rapidly converge and the computed band gap agrees with the experimental one within $\sim 0.1$ eV (see Table. 1). The results for dielectric constants and band gaps of several solids, including covalently, ionic and van der Waals bonded systems, are shown in Tables. 1 and  2, respectively. Our results for the dielectric constants are all close to those of self-consistent hybrid calculations reported in Ref. Skone et al. (2014) [using the functional of Eq. (10) with $\epsilon_{\infty}$ replacing $\epsilon(\bm{r})$]. The use of the microscopically averaged $\epsilon$ over the whole cell appear to yield results in slightly better agreement with experiments. Part of the small differences between column 3 and 4 in Table I is due to the use of pseudopotentials (this work) versus all electron calculations (Ref. Skone et al. (2014)).

Table II shows band gaps obtained with the functional of Eq. (10) and the procedure shown in Fig. 5 (column 3) and those obtained with the global hybrid functional defined in Ref. Skone et al. (2014) with two different values of $\epsilon_{\infty}$: the bulk average of $\epsilon(\bm{r})$ computed in this work (column 4), and the $\epsilon_{\infty}$ from Ref. Skone et al. (2014) (column 5). Considering that the all-electron results of Ref. Skone et al. (2014) (reported in column 6) are obtained with all electrons and a localized basis set, the comparison between columns 5 and 6 shows differences arising from the use of pseudopotentials and the plane-wave basis set. The comparison between column 4 and 5 shows the sensitivity of the band gaps to slightly different values of $\alpha$. The most interesting comparison is between column 3 and 4 which shows that the spatial variations of $\epsilon(\bm{r})$ hardly affect the band gap of covalently bonded systems; however they do influence the computed gap for ionic and especially van der Waals bonded solids.

Fig. 6 shows the band gap for monolayer MoS2 and h-BN. The dielectric hybrid hybrid functional (DDH) of Eq. (10) predicts a fundamental gap of 3.1 eV for MoS2. The effect of spin-orbit coupling, known to lead to a splitting of the degenerate valence bands of about 0.1 eV Cheiwchanchamnangij and Lambrecht (2012), was neglected in our calculations. Therefore, we conclude that our quasiparticle gap is in reasonable agreement with the experimental value of 2.78(2) eV Yao et al. (2017).

The self-consistent hybrid functional of Eq. (10) predicts a gap of 8.2 eV for h-BN. This is consistent with that obtained with self-consistent GW calculations ($\sim 8.4$ eV) in Ref. Berseneva et al. (2013). The Kohn-Sham gap obtained in PBE calculations is about 4 eV smaller, and G0W0 and GW0 results using PBE wavefunctions also underestimate the quasiparticle gap by $\sim 2$ and $\sim 1$ eV, respectively Berseneva et al. (2013); Smart et al. (2018).

The dielectric function $\epsilon(z)$ of the 2D systems studied here turns out to be localized at the monolayers (see Fig. S3 in Ref. SM ). This provides a physical measure of the “dielectric thickness” of the 2D layers, which we define as $w_{\epsilon}=\frac{\int dz(z-z_{0})^{2}\chi(z)}{\int dz\chi(z)}$ where $\chi(z):=\epsilon(z)-1\,$. We obtain a thickness of 3.4 and 1.6 Bohr for MoS2 and h-BN respectively. The spreads of the charge density (see Fig. S3 in Ref. SM ) are 2.3 and 1.3 Bohr respectively, slightly smaller than those of the respective dielectric functions, but comparable.

V.2 Surfaces and interfaces

In the case of surfaces and interfaces, we carried out calculations with the scheme outlined in Fig. 5, applying the $\bm{E}$ field parallel to the surface/interface, insuring that the tangential part of the $\bm{E}$ field is continuous across the interface. (If a constant $\bm{D}$ field were applied, when minimizing the functional Eq. (12), the $\bm{D}$ field would be instead perpendicular to the interface).

Fig. 7 shows the dielectric function and band offsets for an unreconstructed, hydrogen terminated silicon (111) surface (H-Si). We find that the dielectric constant in the silicon bulk regions is $\sim$ 9, which is smaller than that reported in Table 1, due to finite size effects. Indeed, the silicon slab has only 72 Si atoms, a size insufficient to converge the dielectric constant to the bulk value.

The band gap of the silicon portion of the slab and the band offsets between the surface and vacuum obtained from DDH calculations are in good agreement with those of G0W0@PBE calculations; we note that there is a slight difference in the spatial variation of the conduction band at the interface, which is sharper in the case of the hybrid functional calculations, possibly indicating differences between the PBE wavefunctions and charge density (not updated in the GW calculations) and the respective quantities computed self-consistently at the hybrid level.

Calculations for representative interfaces (H-Si/H2O, CH3-Si/H2O, COOH-Si/H2O, and Si/Si3N4) are shown in Fig. 8. We again observe that the calculation of $\epsilon(\bm{r})$ converges rapidly, after 3-4 iterations (see Fig. 8, upper panels). We can clearly see that there are two distinct average values of $\epsilon$ in the two bulk regions where $\epsilon$ oscillates around a constant value. The transition regions in the four interfaces, defined as the region where $\epsilon(\bm{r})$ changes sharply, have a thickness of approximately 5 Bohr for aqueous interfaces and 10 Bohr for the silicon-silicon nitride interface.

As already found for the hydrogenated Si-surface, the DDH functional of Eq. (10) predicts the band gap in the silicon bulk regions (Fig. 8) in agreement with G0W0@PBE. In the water region of the aqueous interfaces however, the VBM and CBM are substantially different from those predicted by G0W0@PBE calculations; this is understandable since the PBE wavefunctions are not a good approximation of the band edges of water, as shown in Ref. Gaiduk et al. (2018). The DDH calculations are instead in good agreement with the values reported in Ref. Gaiduk et al. (2018) and obtained at the G0W0@sc-hybrid level, where the mixing fraction was taken equal to the electronic dielectric constant of water. The band gap ($10.5$ eV) is also in good agreement with that found in Ref. Gaiduk et al. (2018). In the case of Si/Si3N4 (Table 3) we compare our DDH results with experiment, and we find good agreement (the band gap of silicon is again larger than in experiment, due to finite size effects, i.e. to the small slab chosen in our calculations).

VI Conclusions

We introduced a general dielectric-dependent functional, which is applicable to any semiconductor and insulator and does not contain any adjustable parameter. The functional is a generalization of the self-consistent hybrid functional for homogeneous solids introduced in Ref. Skone et al. (2014), and it is defined using a local, spatially dependent dielectric function. We justified the definition of the functional and the spatial variation of the dielectric function using the disentanglement of the dielectric spectra of heterogeneous systems in terms of the spectra of subsystems; such a disentanglement was achieved using linear combinations of dielectric eigenvectors localized in real space. The local dielectric function was then computed self-consistently by carrying our density functional calculations in finite electric fields.

We showed that the dielectric hybrid functional introduced here predicts the band gaps and dielectric constants of three- and two-dimensional solids, as well as band offsets of surfaces and interfaces, with an accuracy comparable to that of GW calculations, thus paving the way to efficient and accurate calculations of the electronic properties of complex heterogeneous systems.

Finally we note that the formulation introduced in our work provides a definition of the dielectric thickness of interfaces and 2D systems, and a physical interpretation of the spatial variations of single particle energy levels upon the formation of interfaces.

Acknowledgements.

We thank Francois Gygi, Christopher Knight, and Jonathan Skone for numerous discussions. This work was supported by Argonne Leadership Computing Facility (ALCF) Theta Early Science program and the Midwest Integrated Center for Computational Materials (MICCoM). This work used computing resources of Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility under Contract DE-AC02-06CH11357. We also gratefully acknowledge the computing resources provided on Bebop, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory.