Bayesian Optimization of Hubbard U’s for Investigating InGaN Superlattices

: In this study, we undertake a Bayesian optimization of the Hubbard U parameters of wurtzite GaN and InN. The optimized Us are then tested within the Hubbard-corrected local density approximation (LDA+U) approach against standard density functional theory, as well as a hybrid functional (HSE06). We present the electronic band structures of wurtzite GaN, InN, and (1:1) InGaN superlattice. In addition, we demonstrate the outstanding performance of the new parametrization, when computing the internal electric-ﬁelds in a series of [InN] 1 –[GaN] n superlattices (n = 2–5) stacked up along the c -axis.


Introduction
InGaN multi-layers form multiple-quantum wells (MQWs) that are widely used for producing light-emitting diodes with high efficiency, low power consumption, and the ability to tune the wavelength of the emitted light by varying the structure of the MQWs [1][2][3][4]. The system has received a lot of attention from experiments [5][6][7][8], and theory [9][10][11][12][13][14][15]. The latter focuses on the electronic structure, such as band gap and band alignments, and the internal (built-in) electric fields. Modeling the electronic structure of InGaN multi-layers or superlattices (SL) is complicated by the shortcomings of the standard density functional theory (DFT), when it comes to describing the electronic band gaps. The root cause of this problem is the self-interaction error [16], which leads to charge delocalization and fractional orbital occupations [17]. The problem is especially severe in the case of InN, which comes out as a metal under this approach. Hybrid functionals, like HSE06 or PBE0, are known to improve the description of the electronic properties of semiconductors significantly, as compared to the standard DFT [18,19]. Many body methods, like GW and BSE, offer even better fidelity in terms of reproducing the electronic structure of the real materials. However, the computational cost of the latter is extremely high [19]. Therefore, hybrids are widely used in ab initio calculations of semiconductors. Unfortunately, even hybrids are quite expensive, hence the attempt to develop a much cheaper surrogate, able to reproduce properties of interest with reasonable precision. More affordable methods introduce corrections on top of the standard DFT [20]. One such method is LDA+C [21] applied to correct the electronic structure by adding extra potential terms that are sharply peaked on the atomic sites after obtaining the equilibrium geometry with the standard DFT. Gorczyca and co-workers have used this method for InGaN superlattices [11][12][13]. Another alternative is LDA+U [22][23][24]. It amends the local density approximation (LDA) energy by an additive term U eff 2 ∑ σ (n m,σ − n 2 m,σ ) [23], where U eff is the effective Hubbard U parameter, n m,σ is the occupation number of the given orbital, while m and σ are the magnetic quantum number and the spin index, respectively. Essentially, it attempts to correct the self-interaction error by adding an energy penalty to the orbitals with non-integer occupations. In what follows, we will refer to U eff simply as U. Effectively, in this approach one enforces integer occupations of the orbitals to improve the description of the electronic structure. LDA+U has been used by Terentjevs et al. [25,26] for studying InN nanowires. They proposed to apply the Hubbard U of 6 eV to the 4d orbital of In and 1.5 eV to the 2p orbitals of N. Additionally, Carvalho and colleagues [14] studied the optical properties of In x Ga 1−x N alloys using LDA+U combined with a scissor correction, derived from their advanced HSE+G 0 W 0 study [27]. The Hubbard U parameters they employed were: U(Ga 3d) = 5.7 eV, U(In 4d) = 3.7 eV. This scatter in the Us is unfortunate, since LDA+U offers an attractive route to mitigate the deficiencies of the standard DFT at virtually no computational overhead. Reliable LDA+U parametrization can enable computing systems that are inherently large, e.g., diluted alloys, long-period superlattices, nanowires, or quantum dots. These considerations stimulated us to search for the optimized LDA+U parametrizations for GaN and InN and apply a method for Bayesian optimization of Hubbard U parameters, recently proposed by Yu et al [28]. Bayesian optimization [29,30] is a machine learning method that is used for constrained global optimization of black box functions, that are expensive to evaluate. The black box function can involve DFT calculations, which are quite computationally demanding. The method relies upon creating an approximation to the investigated (objective) function, that is iteratively improved by intelligently choosing where to evaluate the objective function in the next step so as to maximize the obtained insight (utility) [29]. This way, Bayesian optimization requires lesser function evaluations as compared to, e.g., grid-based methods, and thus speeds up the search for an optimal U value. Having obtained the optimized Us, we evaluate their performance in predicting the electronic band structures of bulk GaN, InN, and a 1:1 InGaN superlattice and compare them to the results of the standard DFT and a hybrid. Furthermore, we also compute the values of the internal electric fields in a series of short-period InGaN superlattices.

Density Functional Calculations
All ab initio calculations presented in this paper were performed using DFT as implemented in the VASP [31][32][33][34] code, within the framework of the projector augmentedwave method (PAW) [35]. We employed an energy cut-off of 520 eV. The Brillouin zone integration was performed over a Γ-centered 8 × 8 × 5 k-mesh in the case of bulk nitrides and adjusted for superlattices to take into account the different lengths along the c-axis. Electronic band structures were computed along the 'Γ-M-K-Γ-A-L-H-A|L-M|K-H' path, where the special points are defined according to Setyawan and Curtarolo [36]. The exchange-correlation (XC) interactions were treated at LDA [37], PBEsol [38], rotationally invariant LDA+U [23,39,40], and the Heyd-Scuseria-Ernzerhof hybrid functional (HSE06) [41] levels of theory. As discussed in the introduction, the Hubbard U corrections have been applied to different orbitals in the literature. This is likely due to the fact that the frontier orbitals (those close to the Fermi level) are formed by hybridized p-and d-orbitals of Ga/In, and p-orbitals of N. The d-orbitals are expected to suffer the most from the delocalization problem, and, hence, they are often subject to the LDA+U treatment. In this work, the Hubbard U-corrections was applied to the 3d orbitals of Ga and 4d orbitals of In. We also considered simultaneously applying the U-terms to the d-orbitals of Ga/In and the p-orbitals of nitrogen, as well as applying Us only to the p-orbitals of both Ga/In and N, however, less satisfactory results were obtained (see Supplementary Information for more details). In the case of bulk wurtzite GaN [42] and InN [42], we fully relaxed the ionic positions, as well as the cell degrees of freedom for every XC functional. The considered superlattices were constructed by stacking 1 monolayer (ML) of InN and 1 to 5 MLs of GaN along the c-axis, while GaN acting as a template (i.e., InN layer assuming the in-plane lattice parameters of GaN). The superlattices were always relaxed keeping the a and b lattice parameters fixed. All structural optimizations were performed so as to reach the residual forces not exceeding 10 −3 eV/Å.

Bayesian Optimization
In the method of Yu et al [28], one minimizes the deviation of the electronic band structure from a reference one. We use HSE06 for computing the reference, since HSE06 is known to reproduce well the electronic properties of GaN and InN [9,[43][44][45]. The objective function formulated in [28] is defined as: Here, E HSE06 g and E DFT+U g are the Kohn-Sham (KS) band gaps, i.e., E CBM -E VBM , computed at the HSE06 and DFT+U levels of theory, α 1 and α 2 are control parameters, which allow to tune the relative importance of the band gap and the overall resemblance of the band structure. ∆Band is the mean-square deviation of the DFT+U band structure from the HSE06 one [28]: where N E is the total number of eigenvalues ( ), N k is the number of k-points, and N b is the number of bands taken into account. Following [28], we include only 10 bands around the gap-6 below and 4 above-into optimization. In addition, we put focus on reproducing the band gap as closely as possible by setting α 1 = 0.5 and α 2 = 0.5. To fit the objective function, the Gaussian process [29,46] is used as implemented in the Bayesian optimization package [30]. We used the upper confidence bound (UCB) [29] as the utility/acquisition function. At each iteration of the optimization process, the next sampling point proposed by the algorithm corresponds to the maximum of the utility function. The parameter κ that decides the balance between exploitation (search near the already known maxima) and exploration (search in less explored areas) was kept equal to the default value (2.576).
Though it was demonstrated in [28] that Bayesian optimization often yields reliable results already after 10-20 iterations, we decided to be conservative and performed 50 iterations for each considered material for the sake of performance analysis. The search for the optimized Us was performed within the range of 1 to 10 eV.

Bayesian Optimization of U in GaN
We chose to first test the approach of Yu et al. [28] on GaN which is correctly predicted to be a semiconductor already at the simplest level of theory (LDA). Bayesian optimization of the Hubbard U value applied to the 3d orbitals of Ga results in a U value equal to 5.61 eV, in good agreement to 5.7 eV previously proposed in [14]. The details of the performance of the Bayesian optimization are presented in Figure 1.
It is evident that the maximum of the objective function is located within a ±0.001 margin proposed in [28] already after 15 steps. To converge the U value to ±0.01 eV, the procedure required 43 steps. It should be noted that we set up the Bayesian optimization rather conservatively, for instance, we used the first 10 steps to bootstrap the method. Using the optimized U, we compute the electronic band structure and compare it to other XC functionals considered in this work (see Figure 2).  The band structures computed along the high-symmetry paths proposed by Setyawan and Curtarolo [36] enable a comprehensive overview of the electronic structure and its change as a function of the employed XC-functional. All considered levels of theory yield qualitatively the same result, i.e., a semiconductor. Although the shape of the semi-core bands (below −10 eV) is reproduced satisfactory, none of the approximations employed is able to match their position along the energy axis, as compared to HSE06. The best agreement corresponds to LDA+U. The valence bands (−10 to 0) are well described at all levels of theory. Most of the features exhibited by HSE06 are captured by LDA, PBEsol, and LDA+U. Again, the overall best match to the hybrid is scored by LDA+U. The predicted Kohn-Sham band gap varies substantially: 1.9 eV for PBEsol, 2.11 eV in the case of LDA, 3.38, and 3.35 eV for LDA+U with the U-parameter of Carvalho et al. [14] and the one obtained in this work, respectively. The HSE06 KS band gap amounts to 3.49 eV, in good agreement to the HSE06 gap reported in [44] (3.51 eV) and the experimental band gap of 3.51 eV [47]. The Kohn-Sham band gaps computed with LDA and GGAs, to which PBEsol belongs, are known to often underestimate the experimental band gaps [18,48]. The nearly identical performance of the considered LDA+U approximations is also not surprising, given similar values of the Hubbard U parameters. Overall, LDA+U exhibits the best match to HSE06 and experiment in reproducing the electronic band gap.

Bayesian Optimization of U in InN
Bayesian optimization of the Hubbard U value applied to the 4d orbitals of In results in U = 7.14 eV. The performance of Bayesian optimization is detailed in Figure 3. In the case of InN, the maximum of the objective function is located to within ±0.001 even faster than for GaN, i.e., after 13 steps. The U-term maximizing the objective function is converged to ±0.01 eV already after 20 steps, significantly faster than for GaN. The electronic band structure computed with the optimized U and other XC approximations are presented in Figure 4. Similar to GaN, the shape of the semi-core states (below −10 eV) is reproduced satisfactory in all considered approximations. Although LDA and PBEsol better describe the shape of these bands, especially, of the lowest energy ones, the position of these bands along the energy axis is best matched by LDA, as compared to HSE06. The valence states (−10 to 0 eV) are described relatively well by all XC functionals, except for the Γ-point close to 0 eV, where LDA and PBEsol fail badly. Here, standard DFT (LDA and PBEsol) erroneously predicts InN to be a metal. HSE06 exhibits a band gap of 0.84 eV, in agreement to 0.7-0.77 eV reported in literature [9,45]. Though the optimized U (7.14 eV) is much bigger than the one proposed in [14] (3.7 eV), it allows achieving a better description of the KS band gap without resorting to the scissor correction. Thus, LDA+U with U taken from [14] yields the KS band gap of 0.23 eV, whereas the U obtained in this study gives 0.74 eV and compares more favorably to the HSE06 and experimental band gaps of 0.69-0.7 eV [49,50].

Transferability Check: 1:1 InN/GaN Superlattice
The optimized Hubbard Us for bulk GaN and InN were used to investigate InGaN superlattices. Figure 5 shows the electronic band structure of the 1:1 superlattice, with 1ML of InN and 1ML of GaN stacked up along the c-axis, assuming GaN as a template. All considered XC approximations predict a non-zero band gap. However, the values of the KS band gaps strongly depend on the XC functional used: 0.25 eV for LDA, 0.11 eV for PBEsol, 0.64 eV and 1.21 eV for LDA+U with U taken from [14] and this work, respectively. For a comparison, Gorczyca et al. [11][12][13] reported a band gap of 1.5-1.6 eV for the (1:1) SL, obtained within LDA + correction approach (LDA+C). Sun and co-workers [15] estimated a band gap of 1.64 eV based on the bowing parameters from literature. Tsai and Bayram [51] reported a KS gap of 1.57 eV for the ternary wurtzite-based nitride In 0.5 Ga 0.5 N, which has the same chemical composition as the (1:1) SL. Again, the best match to HSE06 (1.46 eV) among the considered XC functionals is achieved while using the optimized U-values obtained in this work. Though the agreement is worse than for bulk GaN and InN, it is still quite remarkable, given that InN experiences substantial compressive in-plane and tensile out-of-plane strains in this SL due to the significant lattice mismatch between InN and GaN. Similar performance is revealed when computing the band gaps of longer period superlattices, as shown in Figure 6.
To sum up, LDA+U with the optimized Us clearly outperforms LDA, PBEsol, and LDA+U with the Us from [14] in predicting the band gap of the considered superlattice.

Internal Electric Fields in 1:2-5 InN/GaN Superlattices
In this section, the performance of the optimized Us in describing the electric fields present in InGaN superlattices is investigated. To this end, we consider SL's made of 1ML InN + 2 to 5 ML of GaN stacked up along a polar direction (c-axis), again assuming GaN acting as a template.  [14] and obtained in this work (U * ). For the comparison, LDA+C values from [13] are also shown.
To evaluate the E-fields, we employ the total potential averaged over the plane perpendicular to the c-axis ( Figure 7a). A closer look into the GaN region of the superlattice (see Figure 7b) reveals a linear change of the total potential as a function of the distance-a clear signature of the E-field. The numerical value of this field is obtained by locating the positions of the potential minima (red dots), fitting a straight line to these data-points, and evaluating the slope. The results for all considered SL's and XC approximations are depicted in Figure 7c,d. As compared to HSE06, LDA+U with U from [14] exhibits the worst performance. Surprisingly, LDA and PBEsol perform quite well, delivering values in a good agreement to HSE06. LDA+U with optimized Us demonstrate a significant improvement over [14], and reveals a deviations from HSE06 on par with LDA and PBEsol. Compared to values reported in the literature [13]: for the (1:3) SL an E-field of ∼7 MV/cm is predicted, which agrees well with our results for HSE, LDA, and PBEsol, but not for LDA+U, which deviates somewhat more, especially when the Us from [14] are employed; for the (1:5) SL an E-field slightly below 5 MV/cm is reported, which is a bit less than our HSE06 result of ∼6 MV/cm. Overall, most of the considered XC approximations provide good results as compared to HSE06, except for LDA+U with Us from [14].

Conclusions
In this paper, we undertook Bayesian optimization of the Hubbard U-parameters of bulk GaN and InN to improve the description of the electronic (Kohn-Sham) band gap as compared to the standard DFT (LDA and PBEsol GGA), as well as to the already existing U-parameters presented by Carvalho et al. [14]. We used a hybrid functional (HSE06) as a reference for the optimization. It resulted in U (Ga 3d) = 5.61 eV and U (In 4d) = 7.14 eV. The new Us improve the prediction of the KS band gaps in the bulk GaN and InN without using the scissor correction, and also provide a significantly improved description of the band gaps of the considered InGaN superlattices. In addition, LDA+U with the optimized Us exhibited a good performance in evaluating the E-fields present in these superlattices. Interestingly, standard LDA and GGA, in spite of their limitations to describe the conducting character of the systems considered, provide E-fields within the same accuracy as LDA+U. It must, however, be admitted that this agreement is fortuitous and results from a cancellation of errors. Hence, we are confident that LDA+U with the optimized Us is an attractive approach offering performance similar to that of HSE06 at a tiny fraction of its cost when computing band gaps and E-fields of InGaN superlattices, e.g., in digital alloys [15]. This performance gain might prove helpful when studying, engineering, and optimizing the electronic properties of InGaN nanostructures for light-emitting diodes [52][53][54], gas sensing [55], electrochemical devices [56], solar energy harvesting, and conversion [57][58][59][60], such as InGaN nanowires, core-shell structures, and quantum dots. Data Availability Statement: The data presented in this work are available from corresponding author upon reasonable request.