Predicting structure-dependent Hubbard U parameters via machine learning

1.   Introduction

The complex nature of many-body interactions makes it a long-standing challenge to further improve the exchange-correlation (XC) energy functional for density functional theory (DFT). The semi-local XC functionals, such as the generalized gradient approximation (GGA) of Perdew, Burke, and Ernzerhof (PBE) [1] and many others [2, 3], have well-known self-interaction issues, failing to describe the energy bands of many compounds correctly, especially those ionic compounds with open-shell elements. One viable solution is to add intra-atomic interactions between electrons to mitigate the self-interaction error intrinsic to the local or semi-local XC functionals, namely the Hubbard U correction. The DFT + U method, which was first proposed by Anisimov et al [4] and further developed by Dudarev et al [5], introduces an on-site Coulomb interaction term to penalize partial occupation of the localized orbitals, which can correctly predict behaviors of strongly correlated systems, e.g. Mott insulators [6]. However, finding appropriate U values for a given material system is generally a challenging task. Previously, Wang et al fitted one U value for a given open-shell transition-metal (TM) element based on experimental chemical reaction enthalpies of TM oxides [7], but the one-U-value predicted reaction enthalpies are less satisfactory [8]. On the other hand, first-principles approaches, such as the linear response method [9] and constrained random phase approximation (cRPA) [10-12], were developed to self-consistently determine the U value for a given system, but the computational cost is increased by at least 10 times. Recently, the machine learning (ML) method was employed, specifically the Bayesian optimization (BO), to extract the U value according to higher-level ab initio results [13]. Such a method has been also successfully applied to interface [14, 15] and superlattice [16], with extra computational overhead.

One major challenge to find an appropriate U lies in its strong dependence on local stoichiometry and structure, since by its very nature, the Hubbard U correction represents a short-range electronic interaction which is implicitly associated with local charge density and orbital symmetry. Therefore, the U value has poor transferability. It depends strongly on stoichiometry or d-valence, as reflected, for example, by a value of 6.7 eV in Ce₂O₃ but 5.13 eV in CeO₂ [17]. It also has a strong structural dependence, as demonstrated by its variation in pressure-induced phase transitions of non-magnetic [18] and magnetic structures [19]. This in turn calls for deriving an accurate U value in order to obtain a more reliable potential energy surface [20]. Especially, in general, one has to derive a specific U value, one at a time, for a given system of the specific chemical and structural environment, which is not only a redundant process but also computational costly and time-consuming. Therefore, it is highly desirable and useful to establish an efficient approach, once for all, which allows one to predict an accurate U value for any given system of any chemical and structural environment on the fly.

In this work, we tackle this problem by employing our in-house high-throughput BO-based workflow, and using the Mn-O compounds as our model system. The model is primarily designed to predict the U value for GGA-PBE when utilizing VASP code, aiming to yield electronic structures closely aligned with the HSE06 level. We employ more than 3000 Mn-O configurations whose band gaps and energy bands are fitted to the high-level hybrid functional result (Heyd-Scuseria-Ernzerhof functional [21, 22], HSE) as for the abundant structure availability of the Mn-O system. This step essentially follows what has been done in [13]. In a second, and more important step, we carry on to employ a supervised random forest ML algorithm [23, 24] to train a predictive Hubbard U model, which then predicts the U value with sufficient accuracy and efficiency for any given Mn-O structure, even those not included in the training dataset. The obtained ML model shows a remarkable accuracy and reaches the coefficients of determination R² = 0.97 and mean absolute error (MAE) of 0.128 eV for U values. More importantly, it is unraveled by the regression that the U value is primarily associated with the bond length, which is consistent with the cRPA theory.

2.   Methods

2.1   Vienna Ab initio Simulation Package (VASP)

The first-principles electronic structure calculations are done using VASP codes [25] based on DFT with the PBE XC functional [1] and the projector-augmented-wave approach [26, 27]. The energy cutoff of the plane-wave basis in the calculation is set to be 520 eV, which suffices for accurately describing the energetics discussed in this work. The $\mathrm{\Gamma }$-centered K-points grid density of 125 k-points$/{\text{Å}}^{-3}$ is adopted for all the calculations, and materials were modeled as ferromagnetic ordering. All the DFT + U results are performed by VASP using the rotationally invariant DFT + U approach introduced by Dudarev et al [5]. We note that different DFT + U implementations may yield different results, furthermore, as shown in [28], PBE has the closest result between Dudarev approximation and Lichtenstein form [29] than local density approximation (LDA) and PBE revised for solids (PBEsol). We stress that the Hubbard U values are not transferable among different DFT implementations.

2.2   FHI-aims

The hybrid DFT calculation is performed by the full-potential, all-electron, numeric atomic orbital-based FHI-aims code [30-32]. The real-space band structures for the HSE [21, 22] reference have been performed using a standard tier 1 basis set and applying intermediate integration grids. The K-points grid density is set to 5 in units of ${\text{Å}}^{-1}$ for obtaining reasonable results. The high-symmetry points are generated by the HighSymmKpath module in the pymatgen library [33] for all the band structure calculations, including VASP and FHI-aims.

2.3   BO

The hyperparameter $\kappa$ controls the trade-off between exploration and exploitation. Here, we set $\kappa =5$. Based on the HSE06 band structure, we include the top 6 valence bands and the bottom 4 conduction bands in the optimization (the top 4 valence bands and the bottom 4 conduction bands are considered for the calculations when PBE results do not have enough valence bands). More details of the BO method can be seen in [13]. The UCB acquisition function is selected as it converges quicker than the probability of improvement (PI) and the expected Improvement (EI) acquisition functions [35].

2.4   ML model

We employ the Random Forest Regression (RFR) implemented in the scikit-learn library [36] to extract the structure-Hubbard U relationship and rank the relative importance of descriptors. Nearly 600 descriptors are generated for each structure following [37], and the model training process sorts out the best ten descriptors for the final model construction. More details about the descriptors are provided in Supplementary. The model is primarily designed to predict the U value for GGA-PBE when utilizing VASP code, aiming to yield electronic structures closely aligned with the HSE06 level.

3.   Results

3.1   Workflow

The process including data generation, structural distortion, BO, and ML is schematically illustrated in figure 1. The thermodynamic stability of the compounds is evaluated by the physical quantity of energy above hull ($E_{\text{hull}}$), which is the reaction enthalpy required to decompose a material to other stable compounds [38, 39]. We start from selecting the thermodynamically stable ($E_{\text{hull}}<200\text{meV}/\text{atom}$) and Mn-O chemical systems with reasonable size ($N_{\text{atom}}<20$) from 312 possible structures generated from the Atomly [40] materials database and end up with 67 individual compounds. To enlarge the size of our structural space, we apply uniaxial, biaxial, and triaxial strain of $-2\mathrm{\%}<\epsilon <2\mathrm{\%}$ to the selected 67 compounds and obtain 3724 Mn-O structures. The strain deformation step is employed to modify the screened local Coulomb interaction U [41, 42] and ensures that our model surveys among a large enough structural space to yield a good model extrapolation. The BO method, which is developed by Yu et al, is employed to determine the Hubbard U parameter within the PBE + U method by fitting to the HSE06 band gap and band structures for all the structures. In this way, a U value for Mn 3d state can be fitted for each structure. We note that materials were modeled as ferromagnetic ordering. In principle, higher-level methods, such as the GW [43] or coupled cluster singles and doubles [44]-level of the calculation, can be used too if one has enough computational resources, but we use the HSE06 as the ‘ground truth’ for its viable efficiency. Finally, a ML model is constructed by harnessing the RFR to directly predict the optimal Hubbard U parameters for any Mn-O structures. The computational details can be found in supplementary.

Figure  1.  Computational process. The flow chart of the machine learning process to create a structure-dependent Hubbard U model.

DownLoad: Full-Size Img PowerPoint

3.2   Model validation

Random forests are a combination of decision trees that individually make predictions on each input and the overall prediction is determined by a majority voting process. We evaluate the prediction ability of our RFR model by plotting the out-of-bag error, which can be analogous to the conventional cross-validation error but provides a global error estimated for all data points, shown in figure 2(a) along with the detailed data distribution. Our model can predict the Hubbard U values fairly accurately with MAE = 0.128 eV and R² = 0.97, meaning that the predicted U value falls into a small error range of ±0.128 eV statistically. The distribution of Hubbard U values shows two peaks due to the uneven distribution of Mn-O bond length which will be discussed later in the paper. It is noticeable that for some structures, their electronic structures are insensitive to the change of Hubbard U parameters, and hence the heavily-structure-dependent model does not apply so well to these compounds, causing the deficiency of the ML model to some extent, but the model overall has good accuracy. To gain a better insight into the physical connection between the Hubbard U parameters and the materials’ properties, we sort out the 10 most important descriptors for the U value prediction, as shown in figure 2(b). There are CB, $\mathrm{\Delta }\chi$, $P$, $T$, and CN (see supplementary, table S1 for their definitions). Further, these descriptors can be divided into four categories: chemical bond (CB), electronegativity difference ($\mathrm{\Delta }\chi$), line surface angle ($P,T$), and coordinate number (CN). It is obvious that, other than the electronegativity, nearly all the decisive descriptors, the CB, line surface angle, and coordinate number, are primarily a function of the atomistic structure of a compound, indicating that the U is primarily a structure-dependent parameter. The CB length is the most important factor for the predictions (50.7%), which is computed as the (normalized) total reduction of the criterion brought by that feature.

Figure  2.  Machine learning model and its applications. (a) The comparison of Hubbard U values from fitting the HSE06 results (Bayesian optimization) and machine learning. (b) The importance ranking of all descriptors in prediction, highlighting the most significant factor from the bond length. Different kinds of descriptors are colored in different colors, as depicted in the legend. The superscript and subscript of the descriptor denote the inter-site and intra-site mathematical operation, respectively (see supplementary). (c) Pearson correlation coefficient matrix of all descriptors. The radii of circles represent the absolute magnitude of coefficients. (d) The performance of different methods to determine the Hubbard U. The objective function is employed as the evaluation of the performance of different methods. Here five structures are chosen from the Atomly database with their chemical formulas and space groups listed by the horizontal axis.

DownLoad: Full-Size Img PowerPoint

Upon categorizing these 10 factors into four aspects, these parameters have a close connection with each other. Therefore, we look into the Pearson correlation matrix (figure 2(c)) of the 10 descriptors. It can be seen that the choice of the descriptors is fairly orthogonal as those descriptors are weakly coupled with each other, except for the CB descriptor ($C{B}_u^v$) and the electronegativity descriptor (AR $\mathrm{\Delta }{\chi }_u^k$, AR represents Allred-Rockow electronegativity). This is expected considering that the electronegativity descriptors themselves are essentially derived from the crystal structure and the atom species. On the other hand, the strength of a CB should depend on the electronegativity difference ($\mathrm{\Delta }\chi$) between the two bonding atoms. The $\mathrm{\Delta }\chi$ between atoms bonded together will greatly affect the charge density distribution of the local structure, thereby influencing the local electronic screening. Therefore, the correlation between electronegativity and the CB length further suggests it is viable to assign a U value to a compound based on its structure.

Figure 3 demonstrates the performance of PBE with predicted U for MnO₂ and MnO compounds. It can be found that the PBE functional is unable to fully capture the electronic structure of the Mn-O compounds (and in fact also other transition metal oxides) due to the incomplete self-interaction error cancellation [46-49], and yields an overestimation of Coulomb repulsion. For MnO₂ (figure 3(a)), the PBE band gap of 0.91 eV (figure 3(b)) is considerably underestimated compared to the HSE06 result of 2.96 eV (figure 3(c)). Moreover, the locations of the conduction band minimum (CBM) and valence band maximum (VBM), as well as the spin-up and spin-down channels, from PBE are drastically different from those from HSE06. Using the Hubbard U predicted in this work, we obtain the spin-polarized band structure that matches the HSE06 result well (figure 3(d)). For MnO, the usage of our predicted Hubbard U successfully corrects the band position closer to the HSE06 values, and more importantly results in a band gap opening for this compound. This material is a conductor according to PBE (figure 3(f)) and a semiconductor with a band gap of 0.39 eV in PBE + U (figure 3(g)), amending the PBE result significantly. We note that adding a proper U correction still underestimates the band gap to some extent, e.g. in MnO₂, the ML U correction increases the gap to 1.63 eV, whereas the HSE06 band gap is 2.96 eV; in MnO, the ML U correction increases the gap to 0.39 eV, whereas the HSE06 band gap is 1.09 eV (figure 3(h)). It reflects that the Hubbard U correction cannot completely capture the features of exchange interactions in HSE06. Overall, our ML model predicts reliable Hubbard U values that apply well to the PBE + U calculations and reproduces the qualitative features of HSE06 band structures.

Figure  3.  Crystal structures and band structures of MnO and MnO₂ obtained using different methods. (a) The crystal structure of MnO₂. Band structures of MnO₂ obtained in different methods: (b) PBE; (c) HSE06; (d) PBE with ML predicted U. (e) The crystal structure of MnO. Band structures of MnO were obtained in different methods: (f) PBE; (g) HSE06; (h) PBE with ML predicted U. The red line and the cyan line represent the spin-up and spin-down bandstructure, respectively and materials were modeled as ferromagnetic ordering.

DownLoad: Full-Size Img PowerPoint

3.3   Structural dependence

As discussed in the previous session, our ML model indicates that the Hubbard U parameter is greatly structure-dependent. In order to investigate the correlation between the Hubbard U parameter and Mn-O bond length, their distributions among different valence states are plotted in figure 4. It can be seen that the distribution of Mn-O band length is shape-wise similar to that of the U parameters. For example, Mn³⁺ exhibits a wide distribution of bond lengths, while its Hubbard U spreads over a wider range (from 0 to 10 eV) compared to the Mn²⁺ and Mn⁴⁺ cases. Furthermore, the overall larger Mn-O bond lengths in Mn²⁺ compounds correspond to their overall larger U values, whereas the smaller Mn-O bond lengths in Mn⁴⁺ compounds lead to smaller U values. The dashed line in figure 4(b) represents the single U value used in the Materials Project (MP) [50] as obtained by fitting experimental data, showing an apparent discrepancy with the distribution of Hubbard U values from this work. Given that, it implies that the Mn-O bond length is a sensible parameter to describe the variation in U, suggesting the necessity of further investigating the relationship between the Mn-O bond length and the Hubbard U.

Figure  4.  Distribution of bond length and Hubbard U. The (a) bond length and (b) Hubbard U distributions with the contribution of different valence states (Mn²⁺, Mn^3+, and Mn⁴⁺) are colored in red, blue, and green, respectively. The dashed line denotes the U value used in the Materials Project.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

This work showcases a ML model for predicting Hubbard U to skip the expensive first-principles U value calculation process without sacrificing accuracy. Although the Mn-O system is selected as the model system, the out-of-box models can be created for the community for all the open-shell elements as the U value is essentially local-structure dependent, which is consistent with the findings reported in previous work [17, 20]. This method has the advantage that the model allows one to assign an appropriate Hubbard U parameter to a system prior to the DFT calculation and yields improved results that are close to the higher-level methods such as HSE06 or GW. The GW-level of accuracy is also tangible once such a GW dataset is available. When applying the DFT + U method to structural relaxations or molecular dynamics simulations, it would be ideal to adjust the U parameter to appropriate values on the fly as the structure evolves [18, 51, 52]. However, this will become prohibitively expensive if the U value is determined using the conventional first-principles approaches, such as the linear response or the cRPA schemes. The pre-trained ML model, as demonstrated in the present work, will make all this readily happen.

Another advantage is that this approach can be extended to several other properties of systems other than the energy band difference. For example, the model can also calibrate the adhesive energy of the system by including the energies in the objective function. Also, the intersite interaction parameter, V [53], can be also incorporated into the model to further improve its predictive power, which we hope to spark a future investigation.

Moreover, the accuracy and robustness of our ML model can be further enhanced with the reinforced dataset. The purpose of this paper is to showcase this approach, while we are aware that with the hybrid-functional-level treatment adopted, only a small dataset is produced (3724 data points) due to computational cost. If the size of the dataset is presumably enlarged by one or two orders of magnitude, the ML model could evolve into a deep neural network, meanwhile, the model accuracy, extrapolation, and generalization can be greatly enhanced.

Finally, our work demonstrates that the Hubbard U parameter is local-structure dependent to some extent. However, the U value we used in this work is a kind of global measure of the electronic screening effect, which may be not sensitive to the change in local structure. One solution can be to assign the U values for every inequivalent site, which requires a tremendous amount of calculation resources and a fairly large dataset for model training. To this end, we hope efforts can be made by the entire community to collaboratively carry forward this method to generally reliable and efficient models to predict Hubbard U values for all open shell elements.

In summary, we developed a data-driven method for predicting the value of Hubbard U for DFT calculations. Specifically, a ML model is constructed to predict the Hubbard U for Mn-O systems, which can accurately assess the U value of a system without running costly first-principles calculations. It is also demonstrated that the predicted Hubbard U can reproduce hybrid functional-level band gap and band structures without actual hybrid functional-level runs. In addition, our ML model reveals the bond length which shares similar distribution with the Hubbard U is the most decisive factor in determining the U value, which can be justified by cRPA theory. Developing a ML model that can accurately yield appropriate U values for a given structure, without actually running expensive and sophisticated electronic-structure calculations, is a long-sought goal. We demonstrate in this work that this is indeed possible, at least in a given type of system. More work is needed to extend the present model from Mn-O systems to general atomic species and structures, but we do not expect essential difficulties that prevent us from eventually achieving this goal. Our ML model not only opens up a new avenue to calculate Hubbard U values for all open-shell elements, but also provides insights into the physical correlation between the U parameters and local structure in condensed matter, which is relevant to many other important physical questions, such as metal-insulator transitions, superconductivity, magnetic phase transition, etc.

As mentioned by the editor, we noticed that a recent paper [54] introduced a methodology for fitting the U + V parameters within a specific system as an improvement to Yu et al [13]. We would point out that this paper is distinctively different from these two papers [13, 54] as they developed methods to fit U or U+ V parameters for the studying system, but our method does not need the extra fitting once the ML model is trained for Mn-O as demonstrated. Thus the two papers [13, 54] can be employed as a tool to generate data for our ML model.

4.1   Future perspectives

Accurately predicting Hubbard U parameters has been a longstanding pursuit due to its relevance to a wide range of physical inquiries concerning the fractional, magnetic, lattice, and charge excitations fundamental to quantum materials and devices. Conventional methods for fitting the Hubbard U parameters, such as the linear response method and cRPA method, requires circumventing the cumbersome steps of first-principles calculations. This work has provided significant physical insights, revealing the close association between the U value and the local environment of the open-shell cations, including factors such as bond lengths, coordination numbers, and more. Consequently, it is proposed that a ML model can be developed to reliably predict structure-specific U values for any given structure, if one have enough data to establish the relationship between U and local structure. The study, then, has successfully demonstrated the development of an ML model for the Mn-O system, capable of predicting Hubbard U parameters without relying on expensive linear response method or cRPA method based on a structure-specific U dataset with 3000 datapoints. Importantly, this work suggests the potential to develop a universal U-predicting model directly from the atomistic structures of compounds, which could significantly reshape the ways of DFT + U calculations.