Local hybrid density functional for interfaces - ACS Publications

of everyday phenomena in physics, chemistry, geology and biology.1 At a solid interface, the abrupt .... additional point charge and the surrounding p...
0 downloads 11 Views 533KB Size
Subscriber access provided by UNIV OF MISSOURI ST LOUIS

Article

Local hybrid density functional for interfaces Pedro Borlido, Miguel A. L. Marques, and Silvana Botti J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.7b00853 • Publication Date (Web): 11 Dec 2017 Downloaded from http://pubs.acs.org on December 18, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Local hybrid density functional for interfaces Pedro Borlido,† Miguel A. L. Marques,‡ and Silvana Botti∗,† †Institut f¨ ur Festk¨orpertheorie und -optik, Friedrich-Schiller-Universit¨at Jena, Max-Wien-Platz 1, 07743 Jena, Germany ‡Institut f¨ ur Physik, Martin-Luther-Universit¨at Halle-Wittenberg, D-06099 Halle, Germany ¶European Theoretical Spectroscopy Facility E-mail: [email protected]

Abstract Hybrid functionals are by now the state-of-the-art for the calculation of electronic properties of solids within density functional theory. The key to their performance is how a part of Fock exchange is mixed with a semi-local exchange-correlation functional. The choice of the mixing parameter is particularly critical in non-homogeneous systems, such as an interface between two solid phases. In this work we propose a local mixing function that is a functional of the electron density through an estimator of the local dielectric function. Using this mixing function to modify the PBE0 and the HSE06 hybrid functionals, we obtain band gaps and band-edge alignments at interfaces with an accuracy that is comparable to the one of the GW approximation. However, and in contrast to GW and other recent self-consistent schemes for the mixing parameter, our approach does not require the evaluation of the dielectric function and leads to a negligible increase of the computation time with respect to standard PBE0 or HSE06 hybrid calculations.

1

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1

Introduction

Interfaces are ubiquitous in nature and in modern technology. They are the place where exchanges with the surrounding environment occur, and they play a central role in a variety of everyday phenomena in physics, chemistry, geology and biology. 1 At a solid interface, the abrupt structural and chemical inhomogeneity at the nanoscale gives origin to novel functional properties, that differ significantly from the properties of the bulk compounds in contact with one another. 2 This explains why the optimization of many existing technologies (and the development of new ones) often passes by a better understanding and control of interfaces. Transistors, lasers, light-emitting devices, and solar cells all exploit the transport and excitation of electrons in the small number of atomic layers that separate two solids phases. An accurate description of electronic band gaps and band alignments at heterostructure interfaces is hence essential to characterize the functioning and performance of semiconductor devices. These quantities are necessary starting points to further evaluate the quantum efficiency of processes such as the photo-generation of carriers, electron-hole separation or recombination, or the charge and energy transport through the interface. However, such physical processes in multilayer devices are still described to a large extent in terms of the electronic properties of the bulk constituents. Alternatively, empirical or semi-classical approaches, often in the framework of multiscale modeling, are employed to deal with the complex problem of modeling an interface. 3 The attempt to investigate directly interfaces, both theoretically and experimentally, is relatively recent. On the experimental side, access to buried interfaces is more problematic than imaging atoms on exposed surfaces. On the theoretical side, the simulation of internal interfaces requires reliable microscopic models of interface reconstruction, 4 that cannot be obtained directly from experimental data. Finding the most favorable atomic configurations at interfaces, accounting for coordination defects, strain, and impurities is not an easy task, and it leads to large unit cells. These supercells can be characterized using density functional 2

ACS Paragon Plus Environment

Page 2 of 30

Page 3 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

theory 5,6 (DFT) with standard semi-local functionals, or even hybrid functionals. 7 However, when electronic excitations are involved, as in the determination of band alignments and defect energy levels, the quality of the results obtained with standard DFT approximations is often unsatisfactory. On the other hand, a “brute force” application of advanced approaches for excited states, such as GW methods, 8,9 is still prohibitive for realistic interface models. In view of that, there is a clear need for more efficient and more accurate quantum approaches to tackle local electronic excitations at interfaces. Much of the computational effort in GW comes from the calculation of the dynamically screened Coulomb interaction W . However, in cases where non-locality (and not dynamical effects) is the dominating factor, we can approximate W and render the methodology more efficient numerically. A popular strategy consists in replacing the inverse dielectric function by a constant, leading to the so-called (global) hybrid functionals. They can also be seen as non-local exchange-correlation (xc) functionals within the generalized Kohn-Sham picture 10 of DFT. Hybrids were first put forward in theoretical chemistry, 11 but these functionals also became the state-of-the-art in solid state physics, 12,13 thanks notably to the introduction of screened hybrids. 14–16 Global hybrid functionals replace a part of the DFT exchange (x) by the Fock term, while the correlation (c) is accounted for by a semi-local DFT functional. The simplest (one parameter) form of a global hybrid reads hyb Exc = αExFock + (1 − α)ExDFT + EcDFT .

(1)

The mixing parameter α is a constant (in principle the same for all physical systems) and it is usually obtained from theoretical constraints, and/or fitting of theoretical or experimental data. Most hybrid functionals proposed in the literature differ by (i) the percentage of Fock exchange (i.e., the value of α) and (ii) the semi-local DFT functional. The first of these

3

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

turns out to be the most important in several physical situations. For example, the band gap of a semiconductor varies linearly with α, going from the strongly underestimated (semilocal) Kohn-Sham value to the overestimated Hartree-Fock value. A fixed mixing parameter of α ∼0.2–0.3 yields good results for a large class of important semiconductors, even if it becomes less accurate when the gap is very large or very small. In the case of a bulk system the correct experimental gap can always be obtained by using an ad hoc, material dependent, mixing parameter. 17 Unfortunately, the situation becomes much more complicated if the material is inhomogeneous, as in the case of an interface, since the ideal ad hoc parameters for the two materials at the hetero-junction can be substantially different. This is exactly the problem that we address in the following.

2

Theory

Our starting point is the quasi-particle equation of many-body perturbation theory in the GW approximation of Hedin. 8 This equation has the same form as the Kohn-Sham equation if one replaces the xc potential with the non-local and frequency-dependent GW self-energy Σ(r, r′ ; ω). We consider now the Coulomb hole and screened exchange (COHSEX) static approximation to GW : 8

ΣCOHSEX (r, r′ ; ω = 0) = ΣsX + δ (r − r′ ) ΣCOH (r, r′ ) =

occ X i

1 φ∗i (r)φ(r′ )W (r, r′ , ω = 0) + δ (r − r′ ) [W (r, r′ , ω = 0) − v(r, r′ )] , (2) 2

where the first term is the statically-screened exchange ΣsX , which has the same structure as the Hartree-Fock self-energy. Similarly to the Hartree-Fock self-energy, ΣsX involves only occupied orbitals, however in ΣsX the bare Coulomb potential v is replaced by the screened Coulomb interaction W = ǫ−1 (r, r′ ; ω = 0)/|r − r′ |. The second term in Eq. (2) is the local Coulomb hole ΣCOH (r, r′ ). The COH term represents the classical interaction between an

4

ACS Paragon Plus Environment

Page 4 of 30

Page 5 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

additional point charge and the surrounding polarization cloud, induced by the same point charge. A global hybrid functional is obtained from (2) when the screening that enters in the definition of W is approximated by a constant, 18 and the COH term is replaced by the static and local part of an hybrid functional. Physically, this constant should be α = 1/ǫ∞ , where ǫ∞ is the low-frequency dielectric constant of the material. This is, in fact, the approach taken in the so-called screened-exchange methods. 10,19,20 These ideas have been widely explored in the past to calculate electronic states of solids: some authors propose to calculate the dielectric constant self-consistently, 21–25 eventually using model dielectric functions to lower the computational burden, while others combine a dielectric-dependent mixing with range-separation. 26,27 There is also a recent attempt to extend the formalism from solids to molecules. 28 A position-dependent dielectric constant, obtained by partitioning the dielectric constant into atomic contributions was proposed by Shimazaki and Nakajima, 29 however the method was only applied to bulk materials. The procedure to obtain the mixing parameter from the dielectric function is not very appealing, as the calculation of this physical quantity is significantly more involved than a standard DFT calculation. That is why some of us have recently proposed 30 to bypass this calculation by using an estimator of the inverse dielectric function. This estimator g¯ is readily available from ground-state calculations, and leads to a mixing function written in terms of the electron density n. It is obtained by averaging a density functional over the whole space (i.e., over the unit cell volume Vcell for a periodic solid): 30

g¯ =

1 Vcell

Z

d3 r

cell

s

|∇n(r)| . n(r)

(3)

The mixing parameter built using Eq. (3) is therefore meaningful for bulk crystals, as the dielectric-dependent mixing parameters based on the macroscopic dielectric constant. We explore here an interesting path to access electronic properties of nanostructured

5

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 30

systems, which involve interfaces between two different materials. More specifically, at the interface between two solids with very different screening properties, it is too restrictive to select a single mixing parameter for the composed material (for example, the average of the two bulk mixing parameters). To extend the idea of a density-dependent mixing, we start by considering a general coordinate-dependent mixing function α(r, r′ ), leading to the (non-)local hybrid functional: occ

Exc

1X =− 2 i,j

Z Z

d3 rd3 r′ φ∗i (r)φ∗j (r′ )

α(r, r′ ) φj (r)φi (r′ ) |r − r′ | Z  + d3 r n(r) [1 − α(r, r)] eDFT (r) + eDFT (r) , (4) x c

where eDFT and eDFT are the exchange energy and correlation energy per unit particle in x c a semi-local approximation, and φi are Kohn-Sham orbitals . Note that a full off-diagonal mixing may be used also in the second term, provided that a non-local function is employed for the DFT exchange energy. 31 We emphasize that Eq. (4) is more general than the usual local hybrids that have the form proposed by Jaramillo et al. : 32

Exc =

Z

n o DFT DFT d3 r n(r) α(r)eFock (r) + [1 − α(r)] e (r) + e (r) . x x c

(5)

Although the application of local hybrids to extended systems has not been explored yet, this class of functional has already gained considerable attention for molecular calculations. After the local hybrid by Jaramillo et al., 32 many local mixing functions have been proposed and tested on molecular systems. 33–38 These local mixing functions differ by their form and by the choice of the exchange and correlation functionals to be used in their construction. Different prescriptions were also considered to satisfy known constraints and fix free parameters. To construct the non-local mixing function of the form α(r, r′ ) we decided to use the separable expression: α(r, r′ ) =

p 6

a(r)a(r′ ) .

ACS Paragon Plus Environment

(6)

Page 7 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

This choice leads to a symmetric α(r, r′ ) and therefore to a Hermitian Hamiltonian. Furthermore, this form is separable, which greatly reduces the computational effort. The local mixing function a(r; σ) is obtained by convoluting the local estimator of the screening 30 of Eq. (3) with a Gaussian of variance σ: 1 g¯(r; σ) = (2πσ 2 )3/2

Z

3 ′

dr

s

  |∇n(r′ )| |r − r′ |2 . exp − n(r′ ) 2σ 2

(7)

This quantity can be easily computed by using fast Fourier transforms (FFT). Since we want to generalize the hybrid functional of Ref., 30 we impose that the new and the old functional give the same results when they are used in bulk crystals. We build therefore a(r; σ) using directly the two sets of parameters m, a1 and a2 fitted in Ref., 30 so that the global mixing parameter g¯ of Eq. (2) of Ref. 30 is automatically recovered when g¯(r; σ) does not vary significantly in the unit cell: a(r; σ) = a1 + a2 [¯ g (r; σ)]m ,

(8)

with m = 1, a1 = −1.00778, and a2 = 1.10507 for a PBE0 form of the hybrid functional 12 and m = 4, a1 = 0.121983 and a2 = 0.130711 for a HSE form. 14,15 The equations above define the hybrid functionals that we propose in this work, and that we call PBE0-local and HSElocal. These names emphasize that the only change with respect to the original functionals HSE06 and PBE0 is the replacement of the constant mixing with the coordinate-dependent mixing function defined by Eq. (6), Eq. (7) and Eq. (8). A well known problem of local hybrids is their gauge freedom, 35,39,40 originating from the non-unique definition of energy densities. Gauge inconsistencies between the exact exchange and PBE exchange energy densities lead to calibration errors for inhomogeneous systems, that have been shown to hamper further refinement of the mixing function for molecular systems. As accurate thermochemical calculations are not our priority, we ignore at this stage this issue. 7

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

We investigate now the impact of the variance σ in Eq. (7) on the electronic properties of bulk crystals and interfaces. From physical arguments we assume that σ should be large enough such that averages of the electronic states are meaningful, but small enough such that microscopic inhomogeneities (such as interfaces) are still accounted for. We expect that this parameter should be of the order of magnitude of the inverse of the screening parameter (1/ω) which defines the separation range 41 in screened hybrids. For the HSE06 functional, for example, ω = 0.11 Bohr−1 , leading to σ = 9.09 Bohr = 4.81 ˚ A. As a last comment, we consider how our functional behaves in some relevant limit cases. In the limit of a uniform electron gas, when the gradient of the density is zero, we can deduce that g¯(r; σ) of Eq. (7) vanishes, implying that a(r; σ) and α(r, r′ ; σ) become constant and equal to the parameter a1 . In Eq. (4) we can see that Exc reduces then correctly to its value for the homogeneous electron gas (HEG), as eFock = eHEG = ePBE when the gradient x x x of the density is negligible. At the positions of the nuclei the gradient of the density diverges. We never face this situation as we use pseudopotentials, but we can anyway consider what would happen if our q functional was implemented in an all-electron DFT code. The function |∇n(r)| would still n(r) be integrable over the unit cell, 42 yielding finite values for the mixing function. Note that, to

reproduce our calculations in an all-electron code, it is likely that the local mixing function a(r; σ) should include only contributions from valence electrons. For a surface (i.e., an interface with the vacuum), in the asymptotic region where the electron density vanishes,

|∇n(r)| n(r)

does not diverge, but we do not recover the correct limit

α → 1. Unfortunately, it is not possible to impose the correct asymptotic limit without modifying the analytic form of our hybrid functional. Further work along this line is in progress. Despite its limitations for surfaces, the functional can already be reliably employed for a variety of systems that contain solid/solid interfaces, such as core-shell nanostrucutres, nanostructures in a matrix, grain boundaries and internal hetero-interfaces.

8

ACS Paragon Plus Environment

Page 8 of 30

Page 9 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

3

Interface models and calculation details

We chose to investigate the interface of Si with SiO2 , a typical building block of silicon (nano)electronics and photovoltaics. This interface has been widely investigated both experimentally and theoretically. One can indeed find in the literature DFT calculations with semi-local functionals, 43 hybrid functionals 44 and different GW approaches, 45 as well as experimental measurements. 46 We use the interface model generated and validated by Giustino and Pasquarello in Ref., 43 which was also used in Ref. 45 In fact, the choice of the supercell requires particular attention, as previous density-functional studies revealed that the band offsets are very sensitivity to the model adopted for the interfacial bonding pattern. 44,47,48 √ √ The model of the Si/SiO2 supercell was built by matching a 8× 8 Si (100) surface with a β-cristobalite SiO2 (100) surface, obtaining a perfect bonding network that includes 11 Si monolayers and 10 SiO2 molecular layers. We used for this supercell a k-point grid restricted to the Γ point. Note that the selected supercell model differs from the one constructed by Alkauskas et al., 44 which includes 9 monolayers of Si and a 17 ˚ A layer of SiO2 . We can therefore expect a variation of few hundreds of meV in the calculated band offsets when comparing to this work. Another group of extensively studied heterostructures are semiconductor superlattices. 49 We will use these systems for further validation of our functional. We computed band alignments for almost lattice matched GaAs/AlAs (100), AlP/GaP (100) and Si/GaP (110) hetero-interfaces, using the same supercells used for the PBE, PBE0 and GW calculations of Ref. 49 For the interfaces with (100) orientation, the supercells were built using four monolayers of each constituent in a tetragonal cell. We adopted in that case a 8 × 8 × 1 k-point grid. For the Si/GaP (110) interface, seven monolayers of each constituent were joined in an orthorombic supercell and a 6 × 4 × 1 k-point grid was selected. In both cases, q-point grids with half the number of points were used. We implemented the calculation of the mixing function α(r, r′ ; σ), required to build the PBE0-local and HSE-local hybrid functionals, in the quantum espresso 50 electronic struc9

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ture code. In order to simplify the functional and preserve computational efficiency, we neglected the derivatives of the mixing function with respect to the density and its gradient, i.e., we applied the mixing directly to the non-local Fock potential. The full functional derivatives are given in Appendix A and the details of the implementation of the functional are discussed in Appendix B. Note that the choice to neglect all derivatives of the mixing function leads to the consequence that the xc potential is not the functional derivative of the energy. This implies that other properties, like forces and dipole moments, stop also being exact derivatives of the energy and a number of exact conditions is violated. 51,52 Many other approaches (see, e.g., the LB94 functional 53 or the modified Becke-Johnson meta-GGA functional 42 ) face similar consistency problems. In our case, this inconsistency can be solved by including the missing functional derivatives. The price to pay in terms of complexity of the implementation and computational burden is however large, as we can deduce from the form of the neglected terms. More details can be found in Appendix A. All calculations were performed with norm-conserving pseudopotentials within a spin unpolarized formalism. We remark that the pseudopotential of Ga did not include Ga 3d electrons in the valence. The use of pseudopotentials implies that our hybrid xc energy and potentials are functionals of the density of valence electrons only. A cutoff energy of 30 Ha was used for the Si/SiO2 interface, while a cutoff energy of 55 Ha was set for the semiconducting superlattices. Beside the new local hybrid functionals, we used for comparison other approximations to the xc functional, namely PBE, 54 HSE06, 14,15 and PBE0. 12 Given an A/B interface, the calculation of band offsets can be done in different ways. Following Van de Walle and Martin, 55 three calculations are required: one for the bulk crystal A, one for the bulk crystal B, and one for the supercell model including an interface interface. From the calculation of the interface, the potential offset ∆V between the two components is obtained using a macroscopic average of the electrostatic potential. This value is then used to align the energy levels extracted from the two separate calculations for the

10

ACS Paragon Plus Environment

Page 10 of 30

Page 11 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

A and B bulk compounds. Results for the Si/SiO2 interface obtained using this procedure are available for comparison in the Supporting Information. Alternatively, one can extract directly the band profiles from a single interface calculation using, for example, the method of Ref. 56 In practice, this method divides the supercell in slabs parallel to the interface of thickness ∆z (we set ∆z = 0.6 ˚ A). A local density of states is defined as the contribution to the total density of states coming from the electron density in the i-th slab of volume Ωi : Z 1 X i d3 k ωkν δ (ǫ − ǫkν ) , Di (ǫ) = VBZ ν

(9)

where we have introduced the weight function i ωkν

=

Z

d3 r|φkν (r)|2 .

(10)

Ωi

The integral in Eq. (9) is over the k-points in the Brillouin zone of volume VBZ , φkν is the Kohn-Sham wavefunction with momentum k, band index ν, and eigenvalue ǫkν . From these local density of states the local conduction band minimum and the local valence band maximum can then be extracted as a function of the value of z associated to each slab.

4

Results

In the upper panel of Fig. 1 we plot the planar average of the mixing parameter α ¯ xy (z; σ) = R A−1 dx dy α(r, r; σ) as a function of the coordinate z perpendicular to the interface, for

different values of the variance σ. The results shown refer to the local hybrid based on the HSE06 screened hybrid. The same qualitative behavior is obtained for the local hybrid based on PBE0.

We observe that for large values of σ the quantity α ¯ xy (z; σ) approaches a constant equal to the average of the mixing parameters of the two bulk materials composing the interface. 11

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

In this sense, applying large values of σ is equivalent to the procedure followed by Alkauskas et al., 44 i.e. to use a constant fraction of Fock exchange equal to the average of the optimal opt + αBopt )/2. mixing values of the two bulk constituents A and B: α = (αA

Smaller values of σ make the variation of α ¯ xy (z; σ) more pronounced, until we observe a saturation of α ¯ xy (z; σ) in the interior of the bulk-like regions. This is visible already for values of the variance parameter smaller than about σ = 4 ˚ A. The height of the plateau corresponds to the value of α in the bulk crystals. We conclude therefore that 2 ˚ A. σ . 4˚ A is a perfectly justified choice. In the following we decided to use a value in the lower end of this range, namely σ = 2.12 ˚ A, as this allows to use smaller supercells to model the interface. For bulk crystals (see also the results reported in the Supporting Information), it turns out that g¯(r; σ) already saturates to a constant for values of σ larger than 1 ˚ A. Consequently, α(r, r′ ; σ) ≈ α assumes the optimal value of the constant mixing of Ref. 30 if σ . 1 ˚ A. This means that our local estimator of the dielectric function is already reasonably defined in a volume that includes only a few atoms. Considering that the interface region extends typically over a few atomic layers, we expect that a small value of about 2–3 ˚ A is a reasonable compromise. This choice allows us to obtain the bulk limit in the middle of the bulk layers, while maintaining enough spatial resolution to sample the interface region. We will test this assumption in the following. Reversing the perspective, the curves of α ¯ xy (z; σ) in Fig. 1 can also be used as a measure of the quality of the model of the interface. After setting the value of σ, can say that the Si/SiO2 supercell is adequate if we recover a plateau in the middle of the slabs. The height of this plateau should, of course, be the value of the mixing parameter α of the bulk crystals. To further assess the effect of the parameter σ on the local electronic properties at the Si/SiO2 interface, we plot in the lower panel of Fig. 1 the band diagram across the interface, i.e. the dependence on the coordinate z of the local valence band maximum (VBM) and the local conduction band minimum (CBM). We can deduce from the profile of the local

12

ACS Paragon Plus Environment

Page 12 of 30

Page 13 of 30

0.55

σ σ σ σ σ σ

α ¯ xy (z)

0.50 0.45 0.40

= 2.12 Å = 2.65 Å = 3.18 Å = 4.81 Å = 7.41 Å = 13.76 Å

0.35 0.30 0.25

CBM/VBM (eV)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

4 3 2 1 0 -1 -2 -3 -4 -5

σ = 2.12 Å σ = 4.81 Å σ = 13.76 Å 0

5

10

15

20

25

30

35

z (Å)

Figure 1: Upper panel: planar averaged mixing parameter α ¯ xy (z; σ). Middle panel: Profile view of the Si/SiO2 structure. Si atoms are blue and O atoms are red. Lower panel: band profiles of the CBM and VBM. The calculations are performed with the HSE-local functional using different values of σ, as indicated in the legend.

13

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation

5 4 3 2 1 0 −1 −2 −3 −4 −5 −6

E (eV)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

PBE0-local HSE-local PBE0 HSE PBE 0

5

10

15

20

25

30

z (Å)

Figure 2: Local CBM and VBM obtained with different xc functionals: PBE, PBE0, HSE, PBE0-local, HSE-local. band edges that the overall width of the interface is around 6 ˚ A. Inside the slabs (for HSElocal and σ =2.12 ˚ A) we obtain local gaps of 1.50 eV for Si and 8.82 eV for SiO2 . This compares well with the values of the band gaps calculated for bulk Si and bulk SiO2 with the same functional (1.17 eV for Si and 8.67 eV for SiO2 ; see the Supporting Information for more details). In Fig. 1 we can also observe that the variations of the frontier states extend significantly more in the SiO2 layer rather than in the Si layer, in agreement with previous PBE calculations. 43 We conclude that the thickness of the SiO2 layer cannot not be further reduced without compromising the reliability of the model. Thanks to our local description of the screening, we gain in fact a direct insight on the variation of local electronic states with the coordinate z. In Fig. 2 we show the local band edges calculated using different xc functionals: PBE, ˙ PBE0, HSE, and our functionals PBE0-local and HSE-local. Also in this case σ = 2.12 ˚ AIn Table 1 we present the band offsets obtained by direct measure from the band edges (i.e., requiring only one calculation for the interface system), together with experimental data 46 and other calculations for the same supercell that we can find in the literature. 45 We see that HSE-local and PBE0-local perform better than standard PBE0 and HSE hybrids, leading to results which are of the same quality as GW . Furthermore, HSE-local performs better than

14

ACS Paragon Plus Environment

Page 14 of 30

Page 15 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Table 1: Band offsets at the Si/SiO2 interface. These values were directly computed from the local band edges. Experimental offsets were obtained from Ref. 46 while the GW values are from Ref. 45 Functional PBE PBE0 HSE GW QPGW PBE0-local HSE-local Expt.

∆Ev 2.30 3.03 2.96 4.1 4.0 4.73 4.27 4.44

∆Ec 2.06 2.27 2.26 2.9 2.7 3.59 3.05 3.38

its PBE0-local counterpart for the case of Si/SiO2 . This was expected since HSE-local gives a better estimate of the bulk SiO2 band gap. We extended then our study to the other semiconductor heterointerfaces specified above. In this case, the supercells used in literature 49 are relatively small and the atoms at the interface are strained to obtain a perfect matching of the two layers. In fact, in the plots of α ¯ xy (z; σ) (not shown here) we can see that it is impossible to reach a constant value of the mixing inside the semiconducting layers within physical values of σ. This indicates that our approach does not describe bulk-like states in the center of the layers and we are therefore rather studying a thin superlattice. Because of this, for these interfaces we resorted to the method of Van de Walle and Martin 55 to calculate the band offsets. While a direct comparison with experiments should be done based on a more realistic interface model, the comparison with published GW calculations is fully consistent, as we are using the same supercells. The resulting values of the band offsets can be found in Table 2. The qualitative description of the band offsets is always similar for both PBE0-local and HSE-local functionals: we determine, in agreement with GW calculations, 49 that AlAs/GaAs and Si/GaP present a type I band alignment while AlP/GaP has a type II staggered configuration. From a quantitative point of view, HSE-local is also in this case closer to the GW results 49 and to experiment.

15

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 2: Band offsets of selected hetero-interfaces. Experimental, and GW values are from Ref. 49 Interface

Functional PBE PBE0 HSE GaP/Si PBE0-local HSE-local G 0 W0 Expt. PBE PBE0 HSE AlP/GaP PBE0-local HSE-local G 0 W0 Expt. PBE PBE0 HSE AlAs/GaAs PBE0-local HSE-local G 0 W0 Expt.

5

∆Ev 0.36 0.41 0.40 0.46 0.40 0.53 0.80 0.42 0.53 0.51 0.44 0.50 0.67 0.55 0.29 0.53 0.52 0.45 0.51 0.60 0.53

∆Ec 0.68 0.74 0.74 0.79 0.77 0.83 0.38 -0.48 -0.52 -0.52 -0.50 -0.52 -0.76 -0.39 0.46 0.21 0.17 0.26 0.17 0.17 0.18

Conclusions

In conclusion, we presented two local hybrid functionals, based on the HSE06 and the PBE0 hybrid functionals, which are defined through a coordinate-dependent extension of the density-dependent mixing of Ref. 30 These local hybrids are built as approximations to the GW self-energy, and they include an accurate description of the local screening. Our functionals allow to access directly, and with accuracy comparable to more involved methods, band offsets and local electronic states at an interface. Their implementation is relatively straightforward and their computational cost is comparable with the one of a standard HSE06 or PBE0 hybrid functional calculation. Our new method allows therefore to calculate efficiently and accurately quasiparticle states at interfaces, using larger and more realistic supercell models, opening the way to a better understanding of the physical processes at the 16

ACS Paragon Plus Environment

Page 16 of 30

Page 17 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

heart of electronic devices. As a future perspective, it will be interesting to test the performance of this functional in the context of time-dependent DFT calculations for electronic excitations in heterostructures.

Appendix A

Functional derivative of the local hybrid functional

We derive here the xc potential of our PBE0-local functional with the non-local mixing function α(r, r′ ) of Eq. (6). We consider here also the terms depending on the functional derivatives of the mixing function α [n, ∇n], that are neglected in the present implementation. This will allow to analyze what would be the additional computational burden caused by their inclusion. A strictly analogous derivation can be done in the case of the HSE06-local functional, remembering that the parameter ω, which controls the separation of the interaction in a long-range and a short-range contribution, is unaffected by our modifications. α,hyb , which we separate in three We consider the energy functional of our local hybrid Exc

terms: the exact-exchange energy (Exα,Fock ), the semi-local PBE exchange energy (Exα,PBE ) and the PBE correlation energy (EcPBE ), α,hyb Exc = Exα,Fock + Exα,PBE + EcPBE ,

(11)

where each component is given by Exα,Fock Exα,PBE EcPBE

occ Z Z ′ 1X 3 3 ′ ∗ ∗ ′ α(r, r ) =− d rd r φi (r)φj (r ) φj (r)φi (r′ ) , ′ 2 i,j |r − r | Z   = d3 r n(r) (1 − α(r, r)) eDFT (r) , x Z = d3 r n(r)eDFT (r) . c

17

ACS Paragon Plus Environment

(12a) (12b) (12c)

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 30

We want to compute within the generalized Kohn-Sham scheme

ui (r) =

1 φ∗i (r)

δExc δφi (r)

(13)

for each of the three terms that constitute the xc energy functional. p The mixing function α(r′ , r) = a(r)a(r′ ) of Eq. (6) is a separable function of a(r) = a1 + a2 g¯m (r), where the parameters a1 , a2 and m are those reported in Ref. 30 The local

estimator of the dielectric function g¯(r) of Eq. (7) is written in terms of the density n(r) and its gradient: 1 g¯(r) = (2πσ 2 )3/2

Z

3 ′

dr

s

  |r − r′ |2 |∇n(r′ )| exp − . n(r′ ) 2σ 2

(14)

We calculate now the functional derivative of g¯(r; σ), as this is the essential ingredient to determine the xc potential.

(  δ¯ g (r′ ) 1 |r − r′ |2 ∂ = exp − 2 3/2 2 δn(r) (2πσ ) 2σ ∂n(r)

s

|∇n(r)| n(r)

!

∂ − ∇· ∂∇n(r)

s

|∇n(r)| n(r)

!)

. (15)

The first term of Eq.(15) is simply ∂ ∂n(r)

s

|∇n(r)| n(r)

For the second derivative we need ∂

p

!

1 =− 2

s

|∇n(r)| . n3 (r)

(16)

|∇n(r)|/∂∇n(r), which we can evalute for each

18

ACS Paragon Plus Environment

Page 19 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

vector component: p p ∂ |∇n(r)| ∂ 4 (∂k n(r)) · (∂k n(r)) = ∂(∂j n(r)) ∂(∂j n(r)) ∂[(∂i n(r)) · (∂i n(r))] 1 = [(∂k n(r)) · (∂k n(r))]−3/4 · 4 ∂(∂j n(r)) 1 2(∂i n(r)) ∂(∂i n(r)) = 4 |∇n(r)|3/2 ∂(∂j n(r)) 1 ∂j n(r) = 2 |∇n(r)|3/2 p ∂ |∇n(r)| 1 ∇n(r) ⇒ = . ∂|∇n(r)| 2 |∇n(r)|3/2

(17)

We calculate now the functional derivative of the semi-local exchange term Exα,PBE , which we rewrite as Exα,PBE

=

Z

d

3

r n(r)ePBE (r) x



Z

d3 r n(r)α(r, r)eDFT (r) . x

(18)

The first term is trivially identified as the traditional PBE exchange energy, with the corresponding functional derivative, δExPBE δE PBE = φ∗i (r) x δφi (r) δn(r) = φ∗i (r)vxPBE (r) .

(19)

To treat the second term, we remember that we need only the diagonal terms of α(r, r) = a(r). Taking the functional derivative we obtain two terms: δ δφi (r)

Z

3 ′



d r n(r )a(r



)ePBE (r′ ) x

=

Z

δePBE (r′ ) ′ d r n(r ) x a(r ) + δφi (r) 3 ′



19

ACS Paragon Plus Environment

Z

d3 r′ n(r′ )ePBE (r′ ) x

δa(r′ ) , δφi (r) (20)

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 30

with Z

  Z δ n(r′ )exPBE (r′ ) δ n(r′ )ePBE (r′ ) x ′ 3 ′ ∗ dr a(r ) = d r φi (r) a(r′ ) δφi (r) δn(r) 3 ′

= φ∗i (r)a(r)vxPBE (r)

(21)

and δa(r′ ) δa(r′ ) = φ∗i (r) δφi (r) δn(r) δ¯ g m (r′ ) = bφ∗i (r) δn(r) = mb φ∗i (r) g¯m−1 (r′ )

δ¯ g (r′ ) . δn(r)

(22)

Gathering the different terms we conclude that 1 δExα,PBE φ∗i (r) δφi (r)

=

vxPBE (r) [1

− a(r)] − mb

Z

d3 r′ n(r′ )ePBE (r′ ) g¯m−1 (r′ ) x

δ¯ g (r′ ) . (23) δn(r)

α For the derivative of the exact exchange energy Exx we obtain again two terms:

Z X α(r′ , r) δExα,Fock ∗ =− φm (r) d3 r′ φ∗i (r′ )φm (r′ ) ′ δφi (r) |r − r| m Z Z φ∗ (y)φ∗m (y)φk (x)φm (x) δα(x, y) 1X d3 xd3 y k . − 2 m,k |x − y| δφi (r)

(24)

The first term corresponds to the usual exchange correlation term of a hybrid functional, while the second term accounts for all derivatives of the mixing function α(r, r′ ). We need therefore to calculate also the functional derivative of the non-diagonal part of α: δα(x, y) a2 m ∗ = φ (r) δφi (r) 2 i

"

a(x) a(y)

1/2

δ¯ g (y) + g¯m−1 (y) δn(r)

20



a(y) a(x)

ACS Paragon Plus Environment

1/2

δ¯ g (x) g¯m−1 (x) δn(r)

#

,

(25)

Page 21 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

which reduces to Eq.(22) for x = y. Since this expression is symmetric, we can replace it in Eq. (24) and exchange the integration indices, obtaining: X δExα,Fock =− φ∗m (r) δφi (r) m

Z

d3 r′ φ∗i (r′ )φm (r′ )

a2 m ∗ X φ (r) − 2 i j,k

=−

X m

a2 m ∗ X φ (r) − 2 i j,k

φ∗k (y)φ∗j (y)φk (x)φj (x) d xd y |x − y|

ZZ

φ∗m (r)a1/2 (r)

3

Z

ZZ

α(r′ , r) |r′ − r|

3

d3 r′ φ∗i (r′ )φm (r′ )

a(x) a(y)

1/2

g¯m−1 (y)

δ¯ g (y) δn(r)



a(x) a(y)

1/2

g¯m−1 (y)

δ¯ g (y) δn(r)

a1/2 (r′ ) |r′ − r|

φ∗k (y)φ∗j (y)φk (x)φj (x) d xd y |x − y| 3



3

(26) It is easy to deduce from Eq. (26) that the extra terms involving the derivatives of the mixing function are not trivial to implement, and they would bring a significant increase of the computational burden, related to the calculation of additional Coulomb matrix elements. On the other hand, if the derivatives of the mixing function are neglected, the computational cost remains essentially unchanged with respect to a standard calculation with hybrid functionals.

Appendix B

Implementation details

The implementation of our local hybrid functional was carried out in quantum espresso, taking advantage of the existing routines to calculate hybrid functionals, described in Appendix A.5 of Ref. 50 The exact exchange integrals contained in Eq. (4),

Exc

occ Z Z 1 X α(r, r′ ) =− φk′ ,v′ (r)φk,v (r′ ) , d3 rd3 r′ φ∗k,v (r)φ∗k′ ,v′ (r′ ) 2 k,v;k′ ,v′ |r − r′ |

(27)

are evaluated exploiting the dual-space formalism. Auxiliary codensities, defined as ρk′ ,v′ ;k,v (r) = 21

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ψk∗ ′ ,v′ (r)ψk,v (r), are computed in real space and transformed by fast Fourier transform (FFT) to the reciprocal space, where the Fock exchange energy is calculated. The application of the Fock-exchange operator to a wavefunction involves additional FFTs and real-space array multiplications, that need to be repeated for every point k of the Brillouin zone grid and for every occupied band. An auxiliary coarser grid of points q in the BZ, centered at the Γ point, is introduced to reduce the computational effort by restricting the sums over k′ to the points k′ = k = q. The integrable divergence that appears in the q → 0 limit is treated using the procedure proposed by Gygi and Baldereschi. 57 Since our mixing is separable in r and r′ (see the form of α in Eq. (6)) and we are neglecting the functional derivatives of α [n, ∇n], the modifications that we needed to include in the existing code were very limited: they consisted in multiplying wavefunctions and codensities with the functions a(r; σ) given in Eq. (8). As such, the only additional effort with respect to standard PBE0 and HSE06 calculations comes from some extra multiplications of arrays and the calculation of g¯(r) defined in Eq. (7), involving a simple convolution handled by a FFT. This explains why there is essentially no increase in computational time in comparison with standard HSE06 or PBE0 calculations. This is also the main reason why we favored a fully non-local and separable form of α(r, r′ ).

Associated Content The Supporting Information is available free of charge on the ACS Publications website at DOI: It contains calculations for bulk systems and further data on the Si/SiO2 interface model.

22

ACS Paragon Plus Environment

Page 22 of 30

Page 23 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Acknowledgments We thank F. Giustino for providing the supercell geometry of the Si/SiO2 interface. M.A.L.M acknowledges partial support from the DFG though projects SFB-762 and MA-6786/1. Computational resources were provide by the Leibniz Supercomputing Centre through the SuperMUC projects p1841a and pr48je.

References (1) Allara, D. L. A perspective on surfaces and interfaces. Nature 2005, 437, 638–639. (2) Brillson, L. J. Surfaces and interfaces of electronic materials; Wiley, 2010. (3) Steinhauser, M. O. Computational multiscale modeling of fluids and solids: theory and simulation; Springer, 2008. (4) Nakagawa, N.; Hwang, H. Y.; Muller, D. A. Why some interfaces cannot be sharp. Nature 2006, 5, 204–209. (5) Hohenberg, P.; Kohn, W. Inhomogeneous electron gas. Phys. Rev. 1964, 136, B864– B871. (6) Kohn, W.; Sham, L. J. Self-consistent equations including exchange and correlation effects. Phys. Rev. 1965, 140, A1133–A1138. (7) Scuseria, G. E.; Staroverov, V. N. Progress in the development of exchange-correlation functionals. In Theory and applications of computational chemistry: The first forty years; Dykstra, C. E., Frenking, G., Kim, K. S., Scuseria, G. E., Eds.; Elsevier: Amsterdam, 2005; pp 669 – 724. (8) Hedin, L. New method for calculating the one-particle Green’s function with application to the electron-gas problem. Phys. Rev. 1965, 139, A796. 23

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(9) Onida, G.; Reining, L.; Rubio, A. Electronic excitations: density-functional versus many-body Green’s-function approaches. Rev. Mod. Phys. 2002, 74, 601. (10) Seidl, A.; G¨orling, A.; Vogl, P.; Majewski, J. A.; Levy, M. Generalized Kohn-Sham schemes and the band-gap problem. Phys. Rev. B 1996, 53, 3764–3774. (11) Becke, A. D. Perspective: Fifty years of density-functional theory in chemical physics. J. Chem. Phys. 2014, 140, 18A301. (12) Ernzerhof, M.; Scuseria, G. E. Assessment of the Perdew-Burke-Ernzerhof exchangecorrelation functional. J. Chem. Phys. 1999, 110, 5029–5036. (13) Brothers, E. N.; Izmaylov, A. F.; Normand, J. O.; Barone, V.; Scuseria, G. E. Accurate solid-state band gaps via screened hybrid electronic structure calculations. J. Chem. Phys. 2008, 129 . (14) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207–8215. (15) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Erratum: “Hybrid functionals based on a screened Coulomb potential” [J. Chem. Phys. 118, 8207 (2003)]. J. Chem. Phys. 2006, 124 . (16) Paier, J.; Marsman, M.; Kresse, G. Dielectric properties and excitons for extended systems from hybrid functionals. Phys. Rev. B 2008, 78, 121201. (17) Alkauskas, A.; Broqvist, P.; Pasquarello, A. Defect levels through hybrid density functionals: Insights and applications. Phys. Status Solidi B 2011, 248, 775–789. (18) Gygi, F.; Baldereschi, A. Quasiparticle energies in semiconductors: Self-energy correction to the local-density approximation. Phys. Rev. Lett. 1989, 62, 2160–2163.

24

ACS Paragon Plus Environment

Page 24 of 30

Page 25 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(19) Maksimov, E. G.; Maxin, I. I.; Savrasov, S. Y.; Uspenski, Y. A. Excitation spectra of semiconductors and insulators: a density-functional approach to many-body theory. J. Phys.: Condens. Matter 1989, 1, 2493. (20) Bylander, D. M.; Kleinman, L. Good semiconductor band gaps with a modified localdensity approximation. Phys. Rev. B 1990, 41, 7868–7871. (21) Shimazaki, T.; Asai, Y. First principles band structure calculations based on selfconsistent screened Hartree–Fock exchange potential. J. Chem. Phys. 2009, 130 . (22) Shimazaki, T.; Nakajima, T. Dielectric-dependent screened Hartree–Fock exchange potential and Slater-formula with Coulomb-hole interaction for energy band structure calculations. J. Chem. Phys. 2014, 141 . (23) Skone, J. H.; Govoni, M.; Galli, G. Self-consistent hybrid functional for condensed systems. Phys. Rev. B 2014, 89, 195112. (24) Gerosa, M.; Bottani, C. E.; Caramella, L.; Onida, G.; Di Valentin, C.; Pacchioni, G. Electronic structure and phase stability of oxide semiconductors: Performance of dielectric-dependent hybrid functional DFT, benchmarked against GW band structure calculations and experiments. Phys. Rev. B 2015, 91, 155201. (25) Koller, D.; Blaha, P.; Tran, F. Hybrid functionals for solids with an optimized Hartree–Fock mixing parameter. J. Phys.: Condens. Matter 2013, 25, 435503. (26) Kronik, L.; Stein, T.; Refaely-Abramson, S.; Baer, R. Excitation gaps of finite-sized systems from optimally tuned range-separated hybrid functionals. J. Chem. Theory Comput. 2012, 8, 1515–1531. (27) Skone, J. H.; Govoni, M.; Galli, G. Nonempirical range-separated hybrid functionals for solids and molecules. Phys. Rev. B 2016, 93, 235106.

25

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(28) Brawand, N. P.; V¨or¨os, M.; Govoni, M.; Galli, G. Generalization of dielectric-dependent hybrid functionals to finite systems. Phys. Rev. X 2016, 6, 041002. (29) Shimazaki, T.; Nakajima, T. Theoretical study of a screened Hartree–Fock exchange potential using position-dependent atomic dielectric constants. J. Chem. Phys. 2015, 142 . (30) Marques, M. A. L.; Vidal, J.; Oliveira, M. J. T.; Reining, L.; Botti, S. Density-based mixing parameter for hybrid functionals. Phys. Rev. B 2011, 83, 035119. (31) Krukau, A. V.; Scuseria, G. E.; Perdew, J. P.; Savin, A. Hybrid functionals with local range separation. The Journal of Chemical Physics 2008, 129, 124103. (32) Jaramillo, J.; Scuseria, G. E.; Ernzerhof, M. Local hybrid functionals. J. Chem. Phys. 2003, 118, 1068–1073. (33) Arbuznikov, A. V.; Kaupp, M.; Bahmann, H. From local hybrid functionals to “localized local hybrid” potentials: Formalism and thermochemical tests. J. Chem. Phys. 2006, 124, 204102. (34) Arbuznikov, A. V.; Kaupp, M. Local hybrid exchange-correlation functionals based on the dimensionless density gradient. Chem. Phys. Lett. 2007, 440, 160 – 168. (35) Perdew, J. P.; Staroverov, V. N.; Tao, J.; Scuseria, G. E. Density functional with full exact exchange, balanced nonlocality of correlation, and constraint satisfaction. Phys. Rev. A 2008, 78, 052513. (36) Schmidt, T.; Kraisler, E.; Makmal, A.; Kronik, L.; K¨ ummel, S. A self-interactionfree local hybrid functional: Accurate binding energies vis-`a-vis accurate ionization potentials from Kohn-Sham eigenvalues. J. Chem. Phys. 2014, 140, 18A510. (37) Johnson, E. R. Local-hybrid functional based on the correlation length. J. Chem. Phys. 2014, 141, 124120. 26

ACS Paragon Plus Environment

Page 26 of 30

Page 27 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(38) de Silva, P.; Corminboeuf, C. Local hybrid functionals with orbital-free mixing functions and balanced elimination of self-interaction error. J. Chem. Phys. 2015, 142, 074112. (39) Arbuznikov, A. V.; Bahmann, H.; Kaupp, M. Local Hybrid Functionals with an Explicit Dependence on Spin Polarization. The Journal of Physical Chemistry A 2009, 113, 11898–11906, PMID: 19848429. (40) Maier, T. M.; Haasler, M.; Arbuznikov, A. V.; Kaupp, M. New approaches for the calibration of exchange-energy densities in local hybrid functionals. Phys. Chem. Chem. Phys. 2016, 18, 21133–21144. (41) Krukau, A. V.; Vydrov, O. A.; Izmaylov, A. F.; Scuseria, G. E. Influence of the exchange screening parameter on the performance of screened hybrid functionals. J. Chem. Phys. 2006, 125 . (42) Tran, F.; Blaha, P. Accurate Band Gaps of Semiconductors and Insulators with a Semilocal Exchange-Correlation Potential. Phys. Rev. Lett. 2009, 102, 226401. (43) Giustino, F.; Pasquarello, A. Electronic and dielectric properties of a suboxide interlayer at the silicon–oxide interface in MOS devices. Surf. Sci. 2005, 586, 183 – 191. (44) Alkauskas, A.;

Broqvist, P.;

Devynck, F.;

Pasquarello, A. Band offsets at

semiconductor-oxide interfaces from hybrid density-functional calculations. Phys. Rev. Lett. 2008, 101, 106802. (45) Shaltaf, R.; Rignanese, G.-M.; Gonze, X.; Giustino, F.; Pasquarello, A. Band offsets at the Si/SiO2 interface from many-body perturbation theory. Phys. Rev. Lett. 2008, 100, 186401. (46) Keister, J. W.; Rowe, J. E.; Kolodziej, J. J.; Niimi, H.; Madey, T. E.; Lucovsky, G. Band offsets for ultrathin SiO2 and Si3 N4 films on Si(111) and Si(100) from photoemission

27

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

spectroscopy. J. Vac. Sci. Technol., B: Nanotechnol. Microelectron.: Mater., Process., Meas., Phenom. 1999, 17, 1831–1835. (47) Peacock, P. W.; Robertson, J. Bonding, energies, and band offsets of Si−ZrO2 and HfO2 gate oxide interfaces. Phys. Rev. Lett. 2004, 92, 057601. (48) Sharia, O.; Demkov, A. A.; Bersuker, G.; Lee, B. H. Theoretical study of the insulator/insulator interface: Band alignment at the SiO2 /HfO2 junction. Phys. Rev. B 2007, 75, 035306. (49) Steiner, K.; Chen, W.; Pasquarello, A. Band offsets of lattice-matched semiconductor heterojunctions through hybrid functionals and G0 W0 . Phys. Rev. B 2014, 89, 205309. (50) Giannozzi, P. QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. J. Phys.: Condens. Matter 2009, 21, 395502. (51) Gaiduk, A. P.; Staroverov, V. N. How to tell when a model Kohn–Sham potential is not a functional derivative. J. Chem. Phys. 2009, 131, 044107. (52) Gaiduk, A. P.; Staroverov, V. N. Construction of integrable model Kohn-Sham potentials by analysis of the structure of functional derivatives. Phys. Rev. A 2011, 83, 012509. (53) van Leeuwen, R.; Baerends, E. J. Exchange-correlation potential with correct asymptotic behavior. Phys. Rev. A 1994, 49, 2421–2431. (54) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 3865–3868. (55) Van de Walle, C. G.; Martin, R. M. Theoretical study of band offsets at semiconductor interfaces. Phys. Rev. B 1987, 35, 8154–8165. (56) Yamasaki, T.; Kaneta, C.; Uchiyama, T.; Uda, T.; Terakura, K. Geometric and electronic structures of SiO2 /Si(001) interfaces. Phys. Rev. B 2001, 63, 115314. 28

ACS Paragon Plus Environment

Page 28 of 30

Page 29 of 30 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(57) Gygi, F.; Baldereschi, A. Self-consistent Hartree-Fock and screened-exchange calculations in solids: Application to silicon. Phys. Rev. B 1986, 34, 4405–4408.

29

ACS Paragon Plus Environment

Journal of Chemical Theory and Computation 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

ACS Paragon Plus Environment

Page 30 of 30