Machine learning generalised DFT+U projectors in a numerical atom-centred orbital framework

Abstract

Accurate electronic structure simulations of strongly correlated metal oxides are crucial for the atomic level understanding of heterogeneous catalysts, batteries and photovoltaics, but remain challenging to perform in a computationally tractable manner. Hubbard corrected density functional theory (DFT+U) in a numerical atom-centred orbital framework has been shown to address this challenge but is susceptible to numerical instability when simulating common transition metal oxides (TMOs), e.g., TiO₂ and rare-earth metal oxides (REOs), e.g., CeO₂, necessitating the development of advanced DFT+U parameterisation strategies. In this work, the numerical instabilities of DFT+U are traced to the default atomic Hubbard projector, which we refine for Ti 3d orbitals in TiO₂ using Bayesian optimisation, with a cost function and constraints defined using symbolic regression (SR) and support vector machines, respectively. The optimised Ti 3d Hubbard projector enables the numerically stable simulation of electron polarons at intrinsic and extrinsic defects in both anatase and rutile TiO₂, with comparable accuracy to hybrid-DFT at several orders of magnitude lower computational cost. We extend the method by defining a general first-principles approach for optimising Hubbard projectors, based on reproducing orbital occupancies calculated using hybrid-DFT. Using a hierarchical SR-defined cost function that depends on DFT-predicted orbital occupancies, basis set parameters and atomic material descriptors, a generalised workflow for the one-shot computation of Hubbard U values and projectors is presented. The method transferability is shown for 10 prototypical TMOs and REOs, with demonstrable accuracy for unseen materials that extends to complex battery cathode materials like LiCo_1−xMg_xO_2−x. The work highlights the integration of advanced machine learning algorithms to develop cost-effective and transferable workflows for DFT+U parameterisation, enabling more accurate and efficient simulations of strongly correlated metal oxides.

---

1 Introduction

A knowledge of the electronic structure of transition metal oxides (TMOs) and rare-earth metal oxides (REOs) is essential for the development of energy materials for net zero. For example, TMOs and REOs are widely used as support materials in heterogeneous catalysis to enhance the reactivity of metal nanoparticles via the formation of catalytically active electron traps (i.e., electron polarons) at interfacial and defect sites. By controlling the formation of electron polarons, catalytic activity and selectivity can be tuned to achieve controllable CO₂ conversion using CeO₂-based catalysts,¹ H₂ production using MoO₃-based catalysts² and the production of sustainable aviation fuel via the Fischer Tropsch process using TiO₂-based catalysts.³ TMOs and REOs are also widely applicable for technologies in the generation and storage of electricity, due to the ability to tune the mobility of electrons and ions. Tuneable charge carrier mobility can be exploited for the development of photovoltaics and lithium ion battery cathodes using dopants such as Nb and Mg to optimise the electrical performance of TiO₂ and LiCoO₂, respectively.^4,5

Despite extensive experimental research across a broad range of TMO- and REO-based energy materials, realising a step change in their design and optimisation requires first-principles simulation methods, such as density functional theory (DFT),^6,7 to bridge theory and experiment; however, accurate DFT simulations of the ground state electronic structures of TMOs and REOs can be very challenging. The self-interaction error (SIE), which arises when using local and semi-local exchange–correlation density functionals to describe the correlation between electrons occupying localised d or f orbitals, results in electronic repulsion that causes systematic errors in predicted material properties, e.g., underestimated insulator band gaps, inaccurate lattice parameters and inaccurate formation energies of point defects and electron polarons.^8–10 To overcome the limitations of DFT for strongly correlated materials, “beyond-DFT” methods can be applied that add corrective schemes to combat the SIE. Hubbard corrected density functional theory (DFT+U) is a suitable and popular beyond-DFT method, involving a tuneable on-site Coulomb repulsion, i.e., an energy penalty against electron delocalisation that is applied selectively to correlated orbitals in the system.¹¹ DFT+U is popular because it adds minimal computational cost compared to standalone DFT,¹² whilst achieving the accuracy of higher levels of theory such as hybrid-DFT if parameterised carefully.¹³

To accurately parameterise DFT+U, one must account for both the magnitude and basis of the Hubbard correction, which are defined using the Hubbard U value and the Hubbard projector function (or Hubbard projector) respectively. These parameters are used to correct the DFT-predicted total energy (E_DFT) with a corrective Hubbard term that treats localised states only (E⁰_U) and a double counting correction (E^dc_U) that prevents the double counting of localised states in both E_DFT and E⁰_U:

E_DFT+U[ρ(r), n_I,m^σ] = E_DFT[ρ(r)] + E⁰_U[n_I,m^σ] − E^dc_U[n_I,m^σ] (1)

where ρ(r) is the electron density and n_I,m^σ is the occupation matrix, whose diagonal elements correspond to orbital occupation numbers for all atoms (I), orbital magnetic quantum numbers (m) and spin channels (σ). According to the rotationally invariant, spherically averaged implementation proposed by Dudarev et al.,¹⁴ the corrective Hubbard term is calculated using the trace (Tr) of the occupation matrix and its square:

The occupation matrix is calculated by projecting all DFT-predicted Kohn–Sham states onto reference orbitals defined by the Hubbard projector, i.e., calculating the overlap between Kohn–Sham states and spatially localised orbitals, before an overlap-dependent assignment of the occupancy of each Kohn–Sham state to the localised orbitals.¹⁵ After evaluating E⁰_U, the corresponding Hubbard potential, i.e., the added correction to the standard Kohn–Sham effective potential, is then obtained by taking the functional derivative of the corrective Hubbard energy with respect to the occupation matrix, yielding an orbital-dependent potential that acts only on the subspace defined by the chosen projector.¹⁴ The Hubbard correction therefore acts as an occupancy-based bias potential that corrects the total energy using the Hubbard U value and the occupation matrix, which necessitates careful choice of both the Hubbard U value and the projector. Choosing an appropriate Hubbard projector is particularly important for accurate simulations of materials with strong covalent character,^15–18 and their specific representation prevents the transferability of Hubbard parameters across electronic structure codes that employ different types of Hubbard projector, e.g., atomic orbitals,¹⁵ Wannier functions,^19–21 projector augmented wave (PAW) projectors²² and muffin-tin orbitals (MTOs).¹⁶

Typically, the Hubbard U value is parameterised in DFT+U calculations with no consideration of the Hubbard projector, using semi-empirical benchmarking of DFT+U-predicted electronic, geometric and energetic properties against reference data from experiments or higher levels of theory.^23,24 These semi-empirical approaches are reliant on accurate and available reference data, which prevents the high-throughput optimisation of Hubbard U values for vast numbers of materials. Several first-principles approaches for computing Hubbard U values have also been developed, including the linear response approach based on constrained DFT (LR-cDFT),²⁵ density functional perturbation theory (DFPT),²⁶ constrained random phase approximation (cRPA)²⁷ and Hartree Fock based approaches (e.g., UHF and ACBN0).^28,29 Of these methods, LR-cDFT is most popular and has been applied in a planewave basis for the high-throughput optimisation of Hubbard U values for over 1000 transition metal oxides.³⁰ However, LR-cDFT can yield unphysical U values, suffers from numerical instability for closed-shell systems^31,32 and is computationally expensive due to the requirement of calculations using large supercells.¹⁸ In addition, Hubbard U values computed using LR-cDFT are site-specific and therefore not transferable from stoichiometric systems with symmetry-equivalent sites to defective systems with broken symmetry.³³ The lack of transferability can prevent the simulation of experimentally observed electrical conductivities of defective TMOs, including the Li-ion battery cathode material LiCo_1−xMg_xO₂, where deep defect states are predicted using Hubbard U values from LR-cDFT in a planewave basis, thus incorrectly suggesting material resistivity.⁵

To minimise the cost and instability of first-principles methods for computing Hubbard U values, active learning methods have been combined with global optimisation algorithms, such as Bayesian optimisation (BO) and Monte Carlo sampling, to minimise a cost function based on reproducing geometric and electronic properties predicted using higher levels of theory.^34–36 In active learning, the cost function is refined by comparison with the output of successive DFT+U calculations, which means there is no a priori knowledge of the DFT+U potential energy surface and the entire approach must be repeated for different systems. Supervised machine and deep learning approaches have been used to improve the transferability of active learning methods for computing Hubbard U values, by attempting to learn the DFT+U potential energy surface for different materials. For example, BO and Random Forest Regression have been used to determine the structure-dependence of the Mn 3d Hubbard U value required to reproduce the electronic band structures of various MnO_x polymorphs computed using hybrid-DFT.³⁷ Similarly, equivariant neural networks have been used to estimate DFPT-predicted Hubbard U values using ground state atomic occupation matrices and interatomic distances for a range of materials, including Li_xFePO₄ and MnO₂.³⁸ However, none of these methods address the parameterisation of Hubbard projectors, as well as differences in the numerical stability of DFT+U calculations with different Hubbard parameters.

The aforementioned challenges in DFT+U parameterisation (accuracy, numerical instability, computational cost and transferability across materials) are considered in this work in the context of simulations of stoichiometric and defective TMOs and REOs in a numerical atom-centred orbital (NAO) framework. In a NAO framework, there is a strong dependence of the accuracy and numerical stability of DFT+U calculations on both the Hubbard U value and the projector. Relying on semi-empirically derived Hubbard U values and the default atomic Hubbard projector can result in inaccurate, unphysical (calculations terminate due to excessive polaron localisation) and unstable (SCF cycle does not converge) simulations of common TMOs, e.g., TiO₂ and REOs, e.g., CeO₂. Thus, we demonstrate simultaneous optimisation of the Hubbard U value and the projector for Ti 3d orbitals in anatase TiO₂ using BO with a cost function defined using symbolic regression (SR), to minimise the errors of target properties relative to experimental references. BO is also subject to constraints on the DFT+U-predicted covalency, to ensure the numerical stability of point defect calculations, as determined using support vector machines (SVMs). Combining SR, SVMs and BO in this manner avoids the need for multiple successive DFT+U calculations, which significantly reduces the overall computational cost for Hubbard parameter optimisation compared to the existing first-principles and active learning approaches. We then extend the method across different materials by expanding the primary feature space for SR to include DFT-predicted orbital occupancies, basis set parameters and atomic material descriptors for a diverse training set of TMOs and REOs; application of hierarchical SR³⁹ allows us to optimise Hubbard U values and projectors from first-principles, by targeting orbital occupancies calculated using hybrid-DFT. The outcome is a transferable approach for the one-shot computation of Hubbard U values and projectors, with good accuracy that is also achieved for unseen materials. Overall, the work demonstrates the development of cost-effective and transferable machine learning-based workflows for more complete DFT+U parameterisation, enabling more accurate and efficient simulations of complex energy materials.

2 Methodology

2.1 Electronic structure calculations

2.1.1 DFT. All electronic structure calculations were performed using the Fritz-Haber Institute ab initio materials simulation (FHI-aims) software,⁴⁰ which uses an all electron numerical atom-centred orbital (NAO) basis set, interfaced with the Python based Atomic Simulation Environment (ASE).⁴¹ The standard light basis set (2020) was used, with equivalent accuracy to the TZVP Gaussian-type orbital basis set,⁴² as decided after benchmarking of the TiO₂ formation energy (see SI, Section 1.1). Relativistic effects were accounted for using the zeroth order regular approximation (ZORA)⁴⁰ as a scalar correction, whilst the system charge and spin was set to zero. Periodic boundary conditions were applied using converged k-point sampling (see SI Section S1.1), whilst the mBEEF meta-GGA exchange–correlation density functional was used,^43,44 as defined in Libxc,⁴⁵ providing the best balance of accuracy and cost compared to other local and semi-local functionals (see SI Section S1.2). Hybrid-DFT reference calculations were performed using the PBE0 exchange–correlation density functional.^46,47 Self-consistent field (SCF) optimisation of the electronic structure was achieved using a convergence criteria of 1 × 10⁻⁶ eV for the change in total energy, 1 × 10⁻⁴ eV for the change in the sum of eigenvalues and 1 × 10⁻⁶ e a₀⁻³ for the change in charge density. Geometry optimisation was performed using the quasi-Newton BFGS algorithm^48–51 with a force convergence criteria of 0.01 eV Å⁻¹. Unit cell equilibrium volumes (V₀) were calculated by fitting to the Birch–Murnaghan equation of state using ASE.⁵² Where presented, formation energies (ΔE_form) for TiO₂, CeO₂ and LiCoO₂ were calculated using the energies of bulk Ti (in the hexagonal close packed, HCP, crystal structure), bulk Ce, Li and Co (all in the cubic crystal structure) and an isolated O₂ molecule using:where i denotes the metal species index in each compound and n_i (n_O) is the number of metal (oxygen) atoms in the formula unit.

2.1.2 DFT+U parameterisation. All DFT+U calculations were performed using the on-site definition of the occupation matrix and the Fully Localised Limit (FLL) double counting correction.¹⁵ In Section 3.1, a Hubbard correction is applied to correct for the Coulomb self-interaction of Ti 3d and Ce 4f orbital electrons in tetragonal TiO₂ and cubic CeO₂, respectively, using the corresponding default Hubbard projector, which is the atomic NAO basis function in the minimal basis set (i.e., the solution of the non spin-polarised single atom Schrödinger equation). The use of localised atomic basis functions provides a reasonable initial guess for constructing the correlated subspace for the DFT+U correction, but is known to overestimate orbital occupation numbers leading to inaccurate predictions of oxidation states and energetic properties of complex oxides.^53–55 Alternative definitions of the Hubbard projector can provide more accurate predictions of orbital occupancies, such as maximally-localised Wannier functions;⁵⁶ however, evaluating these projectors can introduce significant computational overhead. To maximise computational efficiency, FHI-aims allows the Hubbard projector (Φ_Im) to be defined as a linear combination of NAO basis functions (Φ_Imⁱ), specifically the atomic basis function in the mimimal basis set and auxiliary hydrogenic basis functions for the same atomic site I and orbital magnetic quantum number m:where c_i denotes the linear expansion coefficients for a maximum of four basis functions per projector.⁵⁷

Where presented, modified atomic-like Hubbard projectors were therefore constructed as a linear combination of atomic and hydrogenic auxiliary basis functions, given the available basis functions in the light basis set (SI Section S1.3, Table S2). With this definition of the Hubbard projector, the use of basis sets with multiple radial functions for the same localised orbital can reportedly lead to erroneous ground state predictions due to the occupation of electronic states outside the correlated subspace;^15,57 therefore, we restrict our study to the unmodified light basis set, which contains at most one auxiliary basis function per localised orbital as listed in the SI Section S1.3, Table S2. The Hubbard projector is therefore defined using the linear expansion coefficients c₁ and c₂, which correspond to the atomic and auxiliary basis functions, respectively. Positive values of c₂ were avoided to prevent the mixing of the metal d or f auxiliary basis functions with O 2p states in an unphysical manner, as determined by Kick et al. from comparisons with hybrid-DFT.¹⁵ Before constructing modified projectors, the auxiliary basis functions are subject to a Gram–Schmidt orthogonalisation with respect to the corresponding atomic functions, which avoids double counting of the Hubbard correction from the overlap of basis functions with long tails of radial decay.^15,58

2.1.3 Simulating defects and polarons. Constrained DFT+U calculations were performed using the “occupation matrix control” (OMC) method to fix polarons at specific atomic orbitals, by modifying the corresponding occupation matrices (which remain fixed for the entire calculation), e.g., fixing the 3d_z² orbital occupation number as 1 for a nearest neighbour Ti atom relative to an substitutional W dopant to simulate the formation of Ti³⁺ in W-doped TiO₂ (Fig. 1).^15,59,60

Fig. 1 Simulating an electron polaron in W-doped TiO₂ at a nearest neighbour Ti atom, denoted Ti³⁺ by setting the 3d_z² orbital occupation number to 1. The electron polaron can be fixed using the occupation matrix control (OMC) method or initialised using the occupation matrix release (OMR) method. The diagonal elements of the occupation matrix correspond to orbital occupancies for a given magnetic quantum number (3d_m) in the order (from top left to bottom right) 3d₋₂, 3d₋₁, 3d₀, 3d₁ and 3d₂ corresponding to the 3d_xy, 3d_yz, 3d_z², 3d_xz and 3d_x²−y² orbitals, respectively.⁶⁰ Off-diagonals elements in the occupation matrix reflect orbital hybridisation. In this work, the diagonal elements of the occupation matrix are used as a quantitative measure of local chemical bonding environments and to construct workflows for Hubbard parameter optimisation by assessing how the occupation matrix varies with the chosen simulation method (i.e., DFT, DFT+U or hybrid-DFT, detailed in the SI Section S1.3), Hubbard parameters and atomic properties of different materials.

In Section 3.2, self-consistent DFT+U calculations are performed for TiO₂ using a modified atomic-like Ti 3d Hubbard projector, determined using the semi-empirical machine learning approach detailed in Section 2.2. This is expanded in Section 3.3, where self-consistent DFT+U calculations are performed for all materials in Table S2 in the SI Section S1.3, using modified atomic-like Hubbard projectors that were determined using the first-principles machine learning approach detailed in Section 2.3. Self-consistent defect calculations were performed using the “occupation matrix release” (OMR) method to initialise polaron(s); the DFT+U-predicted total energy (E) is then pre-converged using OMC until ΔE ≤ 0.001 eV, at which point the OMC constraint is relaxed and the orbital occupancies are calculated self-consistently.¹⁵ All point defect calculations in this work were performed in a 3 × 3 × 3 supercell, which avoids spurious long-range defect–defect interactions between periodic images. The oxygen vacancy formation energy (ΔE_OV) and the defect energies (ΔE_defect) following substitution of a host metal atom (Ti in TiO₂ and Co in LiCoO₂) with a Nb, W, Co, Mn, Pt, Au, Pd or Mg atom, are calculated as:

ΔE_OV = E_{oxygen deficient bulk} + µ_O − E_{stoichiometric bulk} (5)

ΔE_defect = E_{doped bulk} + µ_host − E_{stoichiometric bulk} − µ_dopant (6)

where the chemical potentials were calculated using the energy of half an isolated O₂ molecule (µ_O), bulk Ti in the HCP crystal structure (µ_Ti) and the bulk dopant species (µ_dopant) in the cubic crystal structure except Mn (tetragonal) and Mg (HCP). No Hubbard correction was applied to the dopant atoms.

2.2 Semi-empirical machine learning approach

A semi-empirical approach was adopted to optimise the Ti 3d Hubbard U value and projector to enable accurate and numerically stable simulations of anatase TiO₂ using DFT+U-calculated properties of the TiO₂ unit cell (targets for regression) and the classified outcome of point defect calculations in a TiO₂ supercell using the OMR method (targets for classification), as illustrated in Fig. 2.

Fig. 2 Semi-empirical approach for simultaneously optimising the Ti 3d Hubbard U value and projector for anatase TiO₂, using the DFT+U-predicted band gap (E_bg), unit cell equilibrium volume (V₀), occupation matrix trace for Ti 3d (Tr[n(Ti 3d)]) and O 2p (Tr[n(O 2p)]) orbitals, total energy (E) and the classified results of bulk oxygen vacancy calculations using OMR.

2.2.1 Symbolic regression. SR was performed using the Sure Independence Screening and Sparsifying Operator (SISSO) algorithm,⁶¹ as implemented in the SISSO++ package,^62,63 to fit empirical correlations for target properties in terms of the primary features U, c₁ and c₂. Empirical correlations were constructed by searching a non-linear secondary feature space by recursively combining the primary features using the algebraic operators , before using sparse regression techniques to select a minimal set of secondary features and linear regression to optimise the coefficients of the final expression. Empirical correlations were fitted with up to three terms, using a recursive depth of three, yielding a linear combination of non-linear terms, e.g., F₁, F₂ and F₃ using the constants a₀, a₁, a₂ and a₃ for a three term correlation:

Target = a₀ + a₁ × F₁ + a₁ × F₂ + a₃ × F₃ (7)

Empirical correlations were constructed to estimate the DFT+U-predicted band gap (E_bg), V₀ and traces of the Ti 3d (Tr[n(Ti 3d)]) and O 2p (Tr[n(O 2p)]) occupation matrices, with the accuracy of each correlation evaluated using the Pearson's coefficient of determination (R²) and root mean squared error (RMSE). The SISSO correlations were then used to evaluate a unitless, regularised cost function (J^SE), which is defined as the Euclidean norm of two terms J^SE₁ and J^SE₂. J^SE₁ is itself the Euclidean norm of the percentage errors of the DFT+U-predicted E_bg and V₀ of bulk anatase TiO₂versus experimental references from the literature (E^Exp_bg and V^Exp₀, respectively).^64,65J^SE₂ is an additional regularisation term to bias J^SE towards larger values of U and c₁, which favours stronger polaron localisation at point defects that is consistent with the experimentally observed formation of mid-gap states within the TiO₂ band gap.^15,66,67J^SE therefore takes the form:

where J^SE₂ involves a constant α = 1 eV⁻¹ to ensure dimensional consistency, whilst being scaled by a factor of 1000 to ensure normalisation with respect to J^SE₁.

2.2.2 Support vector machines. To investigate numerical instability in self-consistent point defect calculations in TiO₂, classification was performed with linear support vector machines (SVMs) using the Scikit-learn Python library.⁶⁸ The SVMs were used to determine the equations of the boundaries S₁ and S₂ separating regions in the feature space U, Tr[n(Ti 3d)] and Tr[n(O 2p)], where calculations rapidly terminated due to unphysical predictions or did not converge due to “charge sloshing” when simulating a bulk oxygen vacancy using the OMR method. SVM classification was also performed to investigate the relationship between the Hubbard parameters and erroneous oxygen vacancy formation energies, which ranged from −7.11 eV to 14.35 eV depending on the choice of U, c₁ and c₂. Here, classification was performed to determine the equation of the linear boundary S₃ separating regions of “physical” (4 eV ≤ ΔE_OV ≤ 6 eV) and “unphysical” (ΔE_OV < 4 eV or ΔE_OV > 6 eV) oxygen vacancy formation energies, in terms of the partial derivatives of the SISSO-predicted total energy (E^SISSO) with respect to U, c₁ and c₂. Numerical partial derivatives of E^SISSO with respect to each Hubbard parameter , and were calculated using the forward finite difference method with a step size of 0.01. Both SR and SVM classification were performed using a training set of ≤60 geometry optimised unit cells and bulk oxygen vacancy calculations, where all computational settings are kept constant except the Ti 3d Hubbard parameters. Feature scaling was performed for all primary features and target properties (denoted using x̄) for the SISSO correlation for V₀ and SVM boundaries S₁ and S₂. The normalisation constants are listed with the non-linear terms, constants and accuracy metrics for all SISSO correlations and SVM boundaries in the SI Section S2.1.

2.2.3 Bayesian optimisation. BO was performed using the GPyOpt Python library⁶⁹ to probabilistically search the Hubbard parameter space by minimising J^SE whilst satisfying three constraints, including those derived from S₁ and S₂. The Hubbard parameter space was searched using the standard Expected Improvement acquisition function,⁷⁰ with a sampling weight of 0.01 to encourage parameter exploration within bounds for U between 0.5 eV and 5 eV, c₁ between 0 and 1.3 and c₂ between −0.6 and 0; beyond these values, numerical instability including calculation termination and non-convergence was observed when simulating the TiO₂ unit cell. An initial population of 1000 combinations of Hubbard parameters were defined using Latin Hypercube Sampling,⁷¹ using the pyDOE Python library,⁷² to build a surrogate model of J^SE that was minimised further using BO for 350 iterations. Sampled Hubbard parameters that violated constraints were penalised by assigning them a high value of J^SE = 1000, therefore discouraging the exploration of the Hubbard parameter space that leads to numerically unstable calculations. The effectiveness of BO in identifying the global minimum J^SE was assessed by comparing the BO-sampled landscape of J^SE with that evaluated using random sampling for 50 000 iterations.

2.3 First-principles machine learning approach

To extend the semi-empirical approach across different materials, a first-principles approach was adopted to parameterise Hubbard U values and projectors by targeting the O 2p orbital occupancies calculated using hybrid-DFT, which was identified as crucial to enable self-consistency of DFT+U simulations of intrinsic and extrinsic defects in TiO₂ using the OMR method. The first-principles approach was similar to the semi-empirical approach outlined in Section 2.2, by aiming to ensure the accuracy and numerical stability of DFT+U simulations of stoichiometric and defective TMOs and REOs using generalised symbolic regression and classification, as illustrated in Fig. 3.

Fig. 3 First-principles approach for optimising the Hubbard U value and the projector. Generalised symbolic regression is used to target the hybrid-DFT-predicted O 2p occupation matrix. Generalised symbolic classification is used to determine constraints on the Hubbard parameter space to ensure numerically stable point defect calculations.

2.3.1 Hierarchical symbolic regression. To generalise SR across different materials, the primary feature space was expanded beyond the Hubbard parameters U, c₁ and c₂, to include data that can either be determined from a single reference DFT calculation or is widely available in the literature. The expanded primary feature space included (1) basis set parameters for the correlated subshell subject to the Hubbard correction, such as the type (Z_type, i.e., hydrogenic or ionic) and effective core charge (Z_val) of the auxiliary basis function, (2) DFT-predicted metal d or f and O 2p orbital occupancies averaged over all atoms in the unit cell and (3) atomic material descriptors, including the metal atom electronegativity (χ) and atomic radius (r), as well as the outer subshell type (S, encoding d or f subshells), principal quantum number (Q) and number of electrons in the ion (e^ion). Both Z_type and S were label-encoded as integers to enable their use in SR.

The expanded set of 20 primary features makes the SISSO method computationally intractable due to the attempted exhaustive search of the secondary feature space, which scales exponentially with the number of primary features, therefore, a two-step hierarchical SISSO (HI-SISSO)³⁹ approach was adopted where the output from a first step using SISSO, with 13/20 of all primary features (Hubbard parameters and DFT-predicted orbital occupancies), was used as an input for a second step using HI-SISSO where the remaining primary features are included as inputs. HI-SISSO was used to fit ten empirical correlations to estimate DFT+U-predicted orbital occupancies across the full range of magnetic quantum numbers (m), with seven correlations for m = −3 to +3 (grouping together d and f orbitals) and three correlations for m = −1 to 1 for O 2p orbitals. The training set contained 197 sets of DFT+U-calculated orbital occupancies from optimised unit cells of anatase and rutile TiO₂ (46%), CeO₂ (13%), Cu₂O (11%), Y₂O₃ (11%), ZrO₂ (11%), WO₃ (7%) and MoO₃ (1%), where the percentages correspond to the proportion of each material in the dataset. LiCoO₂ and LiFePO₄ were used as blind test cases with no training data. The constants, accuracy metrics and non-linear terms F₁, F₂ and F₃ for all HI-SISSO correlations are listed in the SI Section S2.2.1. As illustrated by the parity plots in Fig. 4 for O 2p orbitals, the HI-SISSO approach (including basis set parameters and atomic descriptors) improved the predictive accuracy of all three correlations compared to the single step SISSO approach, which is equivalent to Δ-machine learning with respect to DFT-predicted orbital occupancies.⁷³

Fig. 4 Parity plots for the DFT+U- and SR-predicted O 2p orbital occupancies for (a) n_m = −1, (b) n_m = 0 and (c) n_m = 1. Blue markers show the predictions using a single step SISSO fitting using the primary features U, c₁, c₂ and all DFT predicted metal d or f and O 2p orbital occupancies. Orange markers show the predictions after a second HI-SISSO fitting, where the outputs from the first step are included within a new set of primary features including S, Z_type, Z_val, Q, χ, e^ion and r.

The HI-SISSO-predicted O 2p occupancies were then used to construct a cost function for optimising Hubbard U values and projectors from first-principles, J^FP, by targeting O 2p orbital occupancies calculated using hybrid-DFT, as outlined in Section 3.3.1.

2.3.2 Symbolic classification. The observed numerical instability of self-consistent DFT+U simulations of point defects in TiO₂, including excessive polaron localisation that causes calculations to terminate and charge sloshing preventing SCF convergence, was observed when modelling point defects across the broader range of TMOs and REOs in Table S2 in the SI Section S1.3. Therefore, the use of SVMs for classifying the stability of point defect calculations in Section 2.2.2 was generalised across materials. Generalised classification was achieved by defining new primary features to account for the different degrees of covalent character across the materials in Table S2 in the SI Section S1.3, which is key for the generalisation of predictive models across complex oxides.⁷⁴ The primary features for generalised classification were U, J^FP, the average error in the DFT+U-predicted metal d or f orbital occupancies relative to hybrid-DFT (E^metal) and the ratio of the traces of the metal d or f and O 2p occupation matrices predicted using hybrid-DFT (R^hybrid), as outlined in Section 3.3.2.

Generalised constraints on the Hubbard parameter space were determined using two linear SVMs to determine the equations of the boundaries S₄ and S₅ that separated regions in the primary feature space leading to the termination and non-convergence of bulk oxygen vacancy calculations using the OMR method. The training set consisted of 86 DFT+U simulations across anatase and rutile TiO₂ (76%), CeO₂ (6%), ZrO₂ (6%), MoO₃ (6%), WO₃ (5%) and Cu₂O (2%), where the percentages correspond to the proportion of each material in the dataset. To reduce the number of misclassified data points associated with S₄ and S₅, classification was performed symbolically using the SISSO algorithm to recursively combine the primary features using the same algebraic operators as in SR, but with the objective of minimising the number of data points in the overlap region of a two-dimensional convex hull.⁶³ The secondary features generated from SISSO were then used as inputs for two linear SVMs, which identified simple boundaries to perform binary classifications in a highly non-linear feature space. The associated constants and accuracy metrics for S₄ and S₅ are listed in the SI Section S2.2.1.

2.3.3 One-shot optimisation of Hubbard parameters. Optimisation of U, c₁ and c₂ from first-principles was performed by a linear search of the landscape of the HI-SISSO-predicted cost function (J^FP_predicted) for each material, over the range of U between 0 eV and 5 eV, c₁ between 0.5 and 1 and c₂ between −0.6 and 0, with each feature split into 50 intervals. The output of each linear search was a family of candidate solutions for a material, with several combinations of U, c₁ and c₂ optimised to minimise J^FP_predicted (see SI Section S2.2.2). The accuracy of the HI-SISSO correlations to predict O 2p orbital occupancies was evaluated by comparing J^FP_predicted to the corresponding DFT+U-computed cost function (J^FP_validated) for up to ten different combinations of Hubbard U values and projectors per material (94 in total), before calculating the mean absolute error (MAE), which is averaged across all tested Hubbard parameters (N) for each material:

The relationship between the accuracy of the one-shot approach for minimising J^FP_predicted and the training set size for each material was investigated using the Mahalanobis distance (D_M)⁷⁵ to quantify the distance of the primary feature vector for each material (x), in a reduced two dimensional feature space determined using principal component analysis (PCA) with the Scikit-learn Python library,⁶⁸ from the mean vector of the training data (κ), using the inverse covariance matrix of the training data (C⁻¹):

An integrated one-shot approach for optimising Hubbard U values and projectors from first-principles was tested for computing the Co 3d Hubbard U value and projector for the simulation of stoichiometric, Mg-doped and oxygen defective LiCoO₂ (i.e., LiCo_1−xMg_xO_2−x), which is unseen by any of the trained regression or classification models. Here, the three HI-SISSO correlations to estimate the DFT+U-predicted O 2p orbital occupancies (for m = −1, 0 and 1) were used to screen the landscape of J^FP_predicted for stoichiometric LiCoO₂, before all Hubbard parameters that violate the generalised constraints S₄ and S₅, which are evaluated using all ten HI-SISSO correlations, are removed from the landscape. The remaining Hubbard parameters were reduced to a smaller set of candidates with K-means clustering, using the Scikit-learn Python library,⁶⁸ which uses unsupervised learning to partition the landscape of J^FP_predicted into smaller clusters by minimising intra-cluster variance. The centroids of these clusters were then used as screened Hubbard parameters for the simulation of stoichiometric and defective LiCoO₂, using the OMR method.

3 Results and discussion

3.1 Projector sensitivities in DFT+U simulations

3.1.1 Stoichiometric oxides. DFT+U is known to be non-trivial when simulating the ground state character (i.e., metallic, semiconducting or insulating) of TMOs and REOs. For example, in a planewave basis, material-dependent transitions in the DFT+U-predicted ground state from metallic to insulating can occur upon increasing the Hubbard U value, which can restore the experimentally observed electronic structures of Mott insulators such as NiO and Ce₂O₃ (ref. 76–79) and point defects in CeO₂ surfaces.⁷⁸ Furthermore, erroneous changes in the DFT+U-predicted hybridisation between metal d and O 2p orbitals can drive the predictions of ground state crystal structures and magnetic properties away from experimental observations, as is reported for BaTiO₃ and layered AMoO₂ (A = Li, Na, K).^80,81 The existence of metastable states on the DFT+U potential energy surface, near integer orbital occupations, can also result in erroneous trapping in local energy minima using local optimisation algorithms, resulting in incorrect ground state predictions for point defects in actinide oxides such as UO₂.⁵⁹

In a NAO framework, similar observations are made when modelling stoichiometric REOs using the default atomic Hubbard projector. For example, applying a Hubbard correction to Ce 4f electrons in stoichiometric CeO₂, using the atomic Ce 4f Hubbard projector, results in an insulator-metal transition (IMT) in the predicted ground state. Upon increasing the Ce 4f Hubbard U value, there is a monotonic change in the DFT+U-predicted E_bg in Fig. 5(a) and ΔE_form in Fig. 5(b) between U values of 0 eV to 9 eV. Beyond U = 9 eV, which would be required if adopting the standard approach of benchmarking against the experimental E_bg of 3.2 eV,⁸² there is a sudden deviation in these trends and DFT+U predicts strong electron localisation in the Ce 4f m = −2 orbital, corresponding to metallic behaviour with no band gap. To trace the root cause of the observed IMT at U = 9.5 eV, constrained DFT+U calculations were performed using the ground state Ce 4f occupation matrix calculated using U = 9.5 eV, but with the m = −2 and m = −3 orbital occupancies systematically varied and all other occupancies fixed. As illustrated in Fig. 5(c), there is a metallic global energy minimum and an insulating low-lying metastable state in the potential energy surface, therefore the observed IMT is caused by the small energy differences between these energy minima, which decreases as the Hubbard U value increases. For U = 9.5 eV, self-consistent determination of the Ce 4f occupation matrix leads to a metallic solution, irrespective of the initial occupation matrix; therefore, we conclude that there are no other insulating solutions with lower energy than the metallic solution and that this represents an IMT in the potential energy surface.

Fig. 5 Overview of errors when modelling stoichiometric CeO₂ using the default Ce 4f atomic Hubbard projector, including the variation in the DFT+U-predicted (a) band gap and (b) formation energy with respect to the Ce 4f Hubbard U value, relative to experimental references denoted by the red dashed lines.^82,83 Blue markers correspond to insulating ground states, yellow markers correspond to metallic ground states and green markers correspond to insulating metastable states. (c) Contour plot of the constrained DFT+U-predicted total energy relative to the ground state energy at U = 9.5 eV, after constraining the m = −2 and m = −3 orbital occupancies. The two regions in dark red correspond to global and low-lying local minima in the potential energy surface with respect to Ce 4f orbital occupancies. The metallic global minimum is 0.273 eV more stable than the insulating local minimum. (d) The radial functions corresponding to the atomic Ce 4f (blue) and hydrogenic auxiliary (orange) basis functions available for constructing a modified atomic-like Hubbard projector. The green and red radial functions correspond to modified projectors that do not include any contribution from the hydrogenic auxiliary function (i.e., c₂ = 0) and are noted with the corresponding shift of the observed IMT.

Beyond U = 11 eV, the insulating ground state character is seemingly restored in Fig. 5(a) and (b); however, these were found to be metastable states as confirmed using OMR with a modified initial Ce 4f occupation matrix, which enabled convergence to the true metallic ground state that exists at all U values beyond the IMT at U = 9.5 eV. The IMT was also found to vary with the definition of the Hubbard projector, as shown in Fig. 5(d), where the more localised projector in green brings the IMT forward to U > 5.5 eV and the more diffuse projector defers the IMT beyond U = 12 eV. The results exemplify the importance of developing new strategies for parameterising the Hubbard U value and projector to enable the accurate simulation of insulating ground states of REOs like CeO₂, whilst avoiding erroneous IMTs and metastable states; this goal cannot be achieved using the standard approach of semi-empirical benchmarking of the Hubbard U value to reproduce the experimental E_bg.

3.1.2 Defective oxides. Whilst DFT+U simulations using semi-empirically-derived Hubbard U values and the default atomic Hubbard projector can accurately model point defects in TMOs such as Li₄Ti₅O₁₂,^84,85 the numerical stability of point defect calculations is generally highly sensitive to the choice of the Hubbard U value, as illustrated in Fig. 6 for anatase TiO₂. Upon increasing the Ti 3d Hubbard U value, there is a monotonic change in the DFT+U-predicted E_bg in Fig. 6(a) and ΔE_form in Fig. 6(b) of anatase TiO₂. An appropriate Hubbard U value could be chosen by considering the compromise in accuracy of these properties versus the experimental references denoted by the red dashed lines in each plot; however, the numerical stability of DFT+U simulations of a bulk oxygen vacancy in anatase TiO₂ were found to vary strongly with both the Hubbard U value and projector. For example, Fig. 6(c) shows the radial functions of two modified atomic-like Ti 3d Hubbard projectors (green and red functions), which are examples of a series of systematically tested Hubbard projectors for the simulation of a bulk oxygen vacancy using the OMR method. In each simulation, two Ti³⁺ polarons were initialised at Ti atoms A and B in Fig. 6(d), which are nearest neighbours relative to the oxygen vacancy. With some combinations of Hubbard U values and projectors, full geometry optimisation successfully completed, whilst in other cases, the simulations were terminated within 2 to 5 iterations of self-consistent optimisation of the system occupation matrices, due to excessive polaron localistion at Ti atom A resulting in the predicted 3d_z² orbital occupancy increasing beyond 4 (which is the condition for termination) in Fig. 6(e). With other combinations of Hubbard U values and projectors, the OMR calculations did not converge due to strong oscillations in the charge density (Fig. 6(f)), which is known as charge sloshing between partially filled, degenerate orbitals and is often associated with metallic systems.⁸⁶

Fig. 6 Overview of errors when modelling stoichiometric and defective TiO₂, including the variation of the DFT+U-predicted (a) band gap and (b) formation energy with respect to the Ti 3d Hubbard U value, using the default atomic Ti 3d Hubbard projector, relative to experimental references denoted by the red dashed lines.^64,87 (c) The radial functions corresponding to the atomic Ti 3d (blue) and hydrogenic auxiliary (orange) basis functions available for constructing a modified atomic-like Hubbard projector. The green and red radial functions correspond to modified projectors that incorporate a contribution from hydrogenic auxiliary function given by the linear expansion coefficient c₂. (d) The nearest neighbour Ti atoms surrounding a bulk oxygen vacancy in anatase TiO₂. (e) The evolution of Tr[n(Ti 3d)] for Ti atoms A, B and C in (d) during an oxygen vacancy calculation using U = 3 eV, c₁ = 1, c₂ = −0.1, which leads to calculation termination due to excessive polaron localisation at atom A, after 3 SCF iterations of OMR (which begins after 23 SCF iterations of OMC). (f) The evolution of the change in charge density during SCF optimisation, during the 1st geometry optimisation step for an oxygen vacancy calculation using U = 3 eV, c₁ = 1, c₂ = −0.5. The calculation does not achieve the convergence criteria of changes in charge density below 1 × 10⁻⁶ e a₀⁻³, denoted by the black dashed line, due to charge sloshing.

3.1.3 Tracing numerical instability to the Hubbard projector. The termination of defect calculations in Fig. 6(e) was observed for Hubbard U values >1 eV using the default atomic Ti 3d Hubbard projector and occurred irrespective of initialising Ti³⁺ polarons further away from the oxygen vacancy as well as tuning available parameters such as SCF mixing parameters and initial occupation matrices. We also tested initialising Ti³⁺ polarons in different Ti 3d orbitals (other than the 3d_z² orbital) and performing local symmetry breaking via targeted bond distortions; however, these did not alleviate the termination of defect calculations in a NAO framework, despite their reported success for aiding SCF covergence in planewave codes.^60,88,89 The SCF non-convergence due to charge sloshing in Fig. 6(f) was also very difficult to alleviate with common strategies. After extensive testing, we could find no robust strategy to mitigate this sloshing, despite tuning available parameters including the basis set size, initial geometries and occupation numbers, SCF mixing parameters, Gaussian broadening parameters, Kerker preconditioning,⁹⁰ the OMR pre-convergence criteria, k-point spacing and using wavefunction restarts. Whilst the charge sloshing in Fig. 6(f) appears insensitive to the aforementioned strategies to improve SCF convergence, other cases of SCF non-convergence such as plateauing of the electron density are more straightforward to address, as observed for Mn-doped rutile TiO₂ in Section 3.2.4, which requires an increased Pulay mixing history beyond the default value. The issue of charge sloshing therefore appears to be specific to DFT+U in the NAO framework, although planewave implementations are not without their own projector-related issues, as highlighted by Warda et al., who discuss the challenge of spurious Hubbard forces and the resulting inaccuracies in the predicted phase stabilities of AUO₄ (A = Ni, Mg, Co, Mn) compounds when using atomic Hubbard projectors.⁹¹

We note three approaches that are observed to overcome the numerical instabilities when simulating point defects using DFT+U in a NAO framework. Kick et al. overcame SCF non-convergence when simulating Ti³⁺ polarons at oxygen vacancies in a rutile (110) TiO₂ surface by increasing the effective core charge (Z_val) of the Tier 1 hydrogenic auxiliary basis function from 2.7e to 4.4e, before using an atomic Ti 3d Hubbard projector to define the basis for the Hubbard correction.¹⁵ Using very small Hubbard U values, or employing constrained DFT+U using the OMC method, can also overcome the observed numerical instability; however, we have previously shown that this comes at the expense of accuracy when modelling polarons in Nb- and W-doped bulk anatase and rutile TiO₂.⁶⁷ For these systems, self-consistent resolution of the system occupation matrices was achieved using a modified Ti 3d Hubbard projector, defined as a linear combination of the atomic and Tier 1 hydrogenic basis functions with the expansion coefficients c₁ = 0.828 and c₂ = −0.561 determined by Jakob and Oberhofer using a first-principles generalised LR-cDFT method in FHI-aims.⁵⁷ With a modified Ti 3d Hubbard projector, DFT+U successfully simulated the polymorph-dependent formation of Nb⁴⁺ and W⁵⁺ polarons as observed in electron paramagnetic resonance (EPR) spectroscopy, which cannot be rationalised by the current planewave DFT+U studies in the literature.⁶⁷ It is therefore essential to understand how to define an appropriate Hubbard projector to enable DFT+U simulations with appreciable Hubbard U values without sacrificing self-consistency and remaining robust with respect to the default basis sets, and we investigate this challenge herein for TiO₂ before establishing a generalised understanding for a broader range of TMOs and REOs.

3.2 Bayesian optimisation of the Ti 3d Hubbard projector

3.2.1 Method configuration. The semi-empirical cost function (J^SE) defined in Section 2.2.1 favours regions in the Hubbard parameter space that achieve a compromise between modelling localised polarons at point defects and the accurate geometric and electronic structure of stoichiometric anatase TiO₂. To rapidly sample the Hubbard parameter space without requiring multiple DFT+U calculations, the terms in J^SE were evaluted using SISSO-computed empirical correlations for the DFT+U-predicted E_bg, V̄₀, Tr[n(Ti 3d)] and Tr[n(O 2p)] in terms of the primary features U, c₁ and c₂, as listed in the SI Section S2.1. Three constraints for sampling the Hubbard parameter space were also defined to ensure physicality of the model and avoid numerical instability during point defect calculations using OMR, which result in either termination or non-convergence in Fig. 6. Hubbard parameters that are predicted by the SISSO models to give a negative occupation matrix trace for Ti 3d or O 2p orbitals, or occupation matrix traces that deviate from the respective hybrid-DFT-predicted occupation matrix trace by over 50%, were excluded. The third constraint was applied using two linear SVMs that classify the validity of a bulk oxygen vacancy calculation in terms of the normalised Ū, and , as illustrated in Fig. 7, which shows the boundaries S₁ and S₂ that separate regions in the Hubbard parameter space that lead to unphysical (termination), charge sloshing (preventing SCF convergence) and ground state (stable convergence) calculation outcomes.

Fig. 7 Illustration of the linear boundaries used to classify simulations of a bulk oxygen vacancy in anatase TiO₂. The boundaries separate successful convergence (green markers), termination due to an unphysical ground state (red markers) and charge sloshing preventing SCF convergence (orange markers). The convex hull associated with each binary classification S₁ and S₂ is shown to illustrate the basis for constructing the constraint in eqn (14).

Fig. 7 illustrates how the numerical stability of defect calculations in TiO₂ is sensitive to the DFT+U-predicted covalency, with calculation termination occurring at high covalency, i.e., and charge sloshing occurring at high metallicity, i.e., . To ensure the successful convergence of point defect calculations whilst mitigating the effects of the changing DFT+U-predicted covalency, the equations of the decision boundaries S₁ and S₂ (trained on DFT+U data, with all constants listed in the SI Section S2.1) were used as constraints on the sampled Hubbard parameter space, based on the SISSO-computed cost function J^SE, with the critical values for S₁ and S₂ chosen based on the convex hull plots in Fig. 7:

S₁ > −9.84 ∧ S₂ > −9.16 (14)

3.2.2 Bayesian optimisation and defect energy corrections. BO was used to search the Hubbard parameter space by minimising J^SE whilst satisfying the constraints on the SISSO-predicted traces of the Ti 3d and O 2p occupation matrices, Tr[n(Ti 3d)]^SISSO and Tr[n(O 2p)]^SISSO, respectively and the SVM-derived constraints S₁ and S₂, as illustrated in Fig. 8.

Fig. 8 The sampled Hubbard parameter space for anatase TiO₂ using (a) BO and (b) random sampling, with markers coloured according to their value of the cost function J^SE. Hubbard parameters that violate the constraints on Tr[n(Ti 3d)]^SISSO, Tr[n(O 2p)]^SISSO, S₁ and S₂ are excluded. (c) The distribution of values of J^SE corresponding to the 1350 sampled Hubbard parameters using BO (red and purple markers) and the results of the first 1350 iterations using random sampling (pink markers). In BO, the prior is trained using evaluations of 1000 randomly sampled Hubbard parameters selected using Latin Hypercube Sampling. During BO, any sampled Hubbard parameters that result in constraint violation are assigned a value of J^SE = 1000. After 1000 iterations, BO is performed for 350 iterations to efficiently optimise U, c₁ and c₂.

Fig. 8(a) and (b) show the results of BO and random sampling, respectively, where BO is able to efficiently optimise U, c₁ and c₂, yielding an almost equivalent set of optimised Hubbard parameters, to those obtained using random sampling, with minimal values of J^SE. This is further supported by Fig. 8(c), which shows the values of J^SE corresponding to the sampled Hubbard parameters. The markers in pink show the first 1350 randomly sampled Hubbard parameters, where constraint violations do not inform subsequent sampling, making the approach highly inefficient. The markers in purple correspond to BO, where the 1000 randomly sampled points selected using Latin Hypercube Sampling are used to train the BO prior and enable efficient optimisation of U, c₁ and c₂. The Hubbard parameters U = 2.749 eV, c₁ = 0.758 and c₂ = −0.354, from the region of lowest J^SE in Fig. 8(a) and (b), were tested for the simulation of both anatase and rutile TiO₂ polymorphs. With the refined atomic-like Ti 3d Hubbard projector, oxygen vacancies and the substitutional dopants Nb, W, Co, Mn, Pt, Au and Pd all successfully converged to the ground state, whilst 11 of these 16 calculations failed when using the same Hubbard U value with the default atomic Ti 3d Hubbard projector (listed in Section 3.2.3, Table 2).

Refining the Ti 3d Hubbard projector therefore enables numerically stable self-consistent defect calculations; however, the predictions for anatase TiO₂ have unphysical defect energies, ranging from −11.59 eV and −1.28 eV; the same is not observed for defective rutile TiO₂, with defect energies ranging from −0.97 eV to 8.22 eV. These unphysical defect energies for anatase TiO₂ necessitate further study into how tuning the Ti 3d Hubbard projector affects the DFT+U-predicted total energy (E). As discussed in the SI Section S2.1.2, we construct an additional SVM constraint S₃ to classify whether a particular combination of U, c₁, and c₂ leads to “physical” (4 eV ≤ ΔE_OV ≤ 6 eV) or “unphysical” (ΔE_OV < 4 eV or ΔE_OV > 6 eV) oxygen vacancy formation energies. Filtering out the data points in Fig. 8(b) that do not satisfy S₃ yields the optimal Hubbard parameters U = 2.575 eV, c₁ = 0.752 and c₂ = −0.486. The approximate similarity is noted of our optimised Ti 3d Hubbard projector with that determined by Jakob and Oberhofer for bulk rutile TiO₂, defined by c₁ = 0.828 and c₂ = −0.561, using a first-principles generalised LR-cDFT method in FHI-aims.⁵⁷

3.2.3 Defect energies and computational cost vs. hybrid-DFT. Simulations of stoichiometric bulk anatase and rutile TiO₂, as listed in Table 1, show DFT+U simultaneously predicts both target properties, E_bg and V₀, with superior accuracy than DFT, relative to experimental references.

Table 1 Geometric, electronic and energetic properties of bulk anatase and rutile TiO₂, predicted using DFT (mBEEF density functional), DFT+U (mBEEF density functional, U = 2.575 eV, c₁ = 0.752 and c₂ = −0.486) and hybrid-DFT (PBE0 density functional), presented alongside experimental references: band gap (E_bg, eV), unit cell equilibrium volume (V₀, ų), formation energy (ΔE_form, eV per atom), Ti 3d occupation matrix trace (Tr[n(Ti 3d)]), O 2p occupation matrix trace (Tr[n(O 2p)])

TiO2 polymorph	Method	E bg (eV)	V 0 (Å^3)	ΔEform (eV per atom)	Tr[n(Ti 3d)]	Tr[n(O 2p)]
Anatase	DFT	2.54	137.26	−3.28	1.18	4.91
Anatase	DFT+U	2.75	137.09	−3.02	0.69	4.64
Anatase	Hybrid-DFT	4.29	135.49	−3.21	1.25	4.66
Anatase	Experiment	3.20 (ref. 64)	136.28 (ref. 65)	−3.24 (ref. 87)	N/A	N/A
Rutile	DFT	2.21	63.50	−3.26	1.17	4.95
Rutile	DFT+U	2.42	62.51	−3.01	0.68	4.63
Rutile	Hybrid-DFT	4.05	62.39	−3.20	1.23	4.69
Rutile	Experiment	3.00 (ref. 64)	62.44 (ref. 65)	−3.26 (ref. 87)	N/A	N/A

---

The hybrid-DFT-predicted V₀ is closest to the experimental value, although hybrid-DFT significantly overestimates E_bg. DFT+U predicts a smaller ΔE_form than DFT, hybrid-DFT and experiment, although all computational methods maintain a similar relative difference between ΔE_form for anatase and rutile TiO₂, with ΔE^anatase_form − ΔE^rutile_form values of −0.02 eV, −0.01 eV and −0.01 eV for DFT, DFT+U and hybrid-DFT, respectively. The DFT+U-predicted covalency of TiO₂ is also suppressed compared to hybrid-DFT, as evidenced by the reduced total Ti 3d orbital occupancy, i.e., Tr[n(Ti 3d)], whilst the DFT+U-predicted total O 2p orbital occupancy, given by Tr[n(O 2p)], remains close to that of hybrid-DFT.

Next, self-consistent DFT+U simulations were performed for anatase and rutile TiO₂ containing an oxygen vacancy and the substitutional dopants Nb, W, Au, Pd, Pt, Co and Mn. Table 2 summarises the success or failure of these simulations, both with or without a refined Ti 3d Hubbard projector and with or without the Hubbard parameters satisfying the error correction constraint derived from S₃ in the SI Section S2.1.2.

Table 2 The effect of Hubbard U value and projector modification on the numerical stability of point defect calculations in TiO₂. Ticks (crosses) correspond to successful convergence (calculation termination) of self-consistent calculations using the OMR method. The satisfaction of constraints derived from the SVM boundaries S₁, S₂ and/or S₃ affects the predicted defect energies corresponding to each set of Hubbard parameters, as shown in Fig. 9(a)

Polymorph	U (eV)	c 1	c 2	Satisfied constraints	Point defect
Anatase	2.749	1	0	None	✗	✗	✗	✗	✗	✗	✗	✗
Anatase	2.575	1	0	None	✗	✗	✗	✗	✓	✓	✓	✓
Anatase	2.749	0.758	−0.354	S 1, S2	✓	✓	✓	✓	✓	✓	✓	✓
Anatase	2.575	0.752	−0.486	S 1, S2, S3	✓	✓	✓	✓	✓	✓	✓	✓
Rutile	2.749	1	0	None	✗	✗	✗	✓	✓	✓	✓	✓
Rutile	2.575	1	0	None	✗	✗	✗	✓	✓	✓	✓	✓
Rutile	2.749	0.758	−0.354	S 1, S2	✓	✓	✓	✓	✓	✓	✓	✓
Rutile	2.575	0.752	−0.486	S 1, S2, S3	✓	✓	✓	✓	✓	✓	✓	✓

---

Comparing the rows in Table 2 that use the same Ti 3d Hubbard U value but a different Ti 3d Hubbard projector shows that refining the Hubbard projector has a direct effect on enabling self-consistent defect simulations. All calculation failures were due to unphysical ground states causing the termination of calculations due to excessive polaron localisation as seen in Fig. 6(e). Refining the Ti 3d Hubbard parameters further from U = 2.749 eV, c₁ = 0.758 and c₂ = −0.354 to U = 2.575 eV, c₁ = 0.752 and c₂ = −0.486, using the constraint derived from the SVM boundary S₃, gave physical defect energies for all studied defects in anatase and rutile TiO₂, as illustrated in Fig. 9(a). Comparing the DFT-, DFT+U- (U = 2.575 eV, c₁ = 0.752 and c₂ = −0.486) and hybrid-DFT predicted defect energies, the DFT+U-predicted ΔE_OV in anatase (5.42 eV) and rutile (5.48 eV) are both much closer to the corresponding values computed using hybrid-DFT (5.29 and 5.76 eV, respectively) than those computed using DFT (4.44 and 4.43 eV, respectively). Without the application of a Hubbard correction to the dopant atom d orbitals, there was no obvious trend in the raw DFT+U-predicted substitutional defect energies when compared to hybrid-DFT. However, when normalised with respect to ΔE_OV in the same TiO₂ polymorph, the DFT+U-predicted values of ΔE_defect were closer to the corresponding hybrid-DFT-predicted values then DFT, as illustrated in Fig. 9(b), which shows the relative defect energies for Pt- and Pd-doped anatase and rutile TiO₂ computed using DFT, DFT+U and hybrid-DFT (single point calculations using the DFT+U geometry). The blue columns in Fig. 9(c) plot the cost of the hybrid-DFT single-point calculation in core-hours per atom, relative to the cost of the DFT+U geometry optimisation. For the four systems presented in Fig. 9(b), DFT+U with geometry optimisation is between 79–189× faster than the hybrid-DFT single-point calculations.

Fig. 9 (a) DFT+U-predicted values of ΔE_OV and ΔE_defect using a refined set of Hubbard parameters that satisfy the constraints derived from the SVM boundaries S₁ and S₂, calculated using the mBEEF density functional, U = 2.749 eV, c₁ = 0.758 and c₂ = −0.354 vs. Hubbard parameters that satisfy the constraints derived from the SVM constraints S₁, S₂ and S₃, calculated using the mBEEF density functional, U = 2.575 eV, c₁ = 0.752 and c₂ = −0.486. (b) ΔE_defect relative to ΔE_OV for Pt- and Pd-doped anatase and rutile TiO₂, calculated using geometry optimisation with DFT and DFT+U (mBEEF density functional, U = 2.575 eV, c₁ = 0.752 and c₂ = −0.486) and single-point calculations using hybrid-DFT (PBE0 density functional, using the DFT+U-optimised geometry). (c) The cost of the hybrid-DFT single point calculations relative to the DFT+U geometry optimisations in core-hours per atom.

3.2.4 Electron polarons in anatase vs. rutile TiO₂. DFT+U using U = 2.575 eV, c₁ = 0.752 and c₂ = −0.486 can be used to simulate experimentally observed differences in electron polaron formation in anatase vs. rutile TiO₂, demonstrating the impact of our parameterisation work. For example, DFT+U predicts the polymorph-dependent formation of small Holstein electron polarons in W-doped TiO₂, characterised by strong electron–lattice interactions that results in a localised state within the TiO₂ band gap.^10,67 This is shown in Fig. 10, where the formation of a localised W 5d_z² state within the rutile TiO₂ band gap is increasingly formed by incrementing the level of theory from DFT in Fig. 10(a), to DFT+U in Fig. 10(b) and hybrid-DFT in Fig. 10(c).

Fig. 10 The elemental species projected density of states for W-doped rutile TiO₂, calculated using (a) DFT (mBEEF density functional), (b) DFT+U (mBEEF density functional, U = 2.575 eV, c₁ = 0.752 and c₂ = −0.486) and (c) hybrid-DFT (PBE0 density functional single point calculation using the DFT+U optimised geometry). The Fermi level is denoted by the red dashed line. The corresponding charge density isosurfaces for the highest occupied molecular orbital (HOMO), at the 0.025 e Å⁻³ level, are shown for (d) DFT and (e) DFT+U.

The degree of defect state localisation in W-doped rutile TiO₂, given by the energy gap between the localised state and the TiO₂ conduction band minimum, is predicted as 0.46 eV using DFT+U and 1.76 eV using hybrid-DFT. As illustrated by the charge density isosurfaces of the highest occupied molecular orbital (HOMO) predicted by DFT in Fig. 10(d) and DFT+U in Fig. 10(e), DFT predicts both W 5d and Ti 3d states to be delocalised, i.e., predicting W⁶⁺ formation, whilst DFT+U predicts a hybridised defect state of W 5d and Ti 3d character, indicating the predicted formation of W⁵⁺. Similarly, DFT+U simulations of Nb-doped rutile TiO₂ predict a Nb 4d signature at the Fermi level that is 5.5× larger than that in Nb-doped anatase TiO₂. The results support previous work using EPR spectroscopy to characterise the formation of paramagnetic Nb⁴⁺ and W⁵⁺ in doped rutile TiO₂, which is recoverable using self-consistent DFT+U as opposed to standalone DFT or constrained DFT+U using OMC.⁶⁷

More generally, the DFT+U simulations reveal clear differences in the electronic structures of substitutionally doped anatase and rutile TiO₂, particularly in the occupancies of the dopant atom d orbitals, as illustrated in Fig. 11(a). In rutile TiO₂, there are greater occupancies of the dopant atom e_g orbitals (m = 0 and 2)⁶⁰ than in anatase, indicating a structural and electronic environment that favours the filling of orbitals aligned along the metal–oxygen bonds. Conversely, substitutionally doped anatase TiO₂ shows increased occupancies of the dopant atom t_2g orbitals (m = −2, −1 and 1)⁶⁰ reflecting a local electronic environment that favours the filling of orbitals aligned between the metal–oxygen bonds. The formation of localised vs. delocalised defects states is directly affected by the filling of the e_g and t_2g orbitals, respectively, leading to polymorph-sensitivity in the electrical conductivity and chemical reactivity of TiO₂-based materials.⁶⁷ Localised defect states within the TiO₂ band gap (or present at the valence or conduction band edges), as observed for Nb- and W-doped rutile in Fig. 11(c) and (d), respectively, and both TiO₂ polymorphs containing an oxygen vacancy and Au, Pd, Pt, Co and Mn dopants in Fig. 11(b), (e)–(i), respectively, give rise to electronic conduction via thermally activated polaron hopping (conductivity increases with temperature) and provide readily available sites for activating reactant molecules.^92,93 Conversely, delocalised defect states located deeper in the valence or conduction bands, as observed for Nb- and W-doped anatase in Fig. 11(c) and (d), respectively, give rise to improved electronic conductivity according to a band-like model and can both reduce the degrading recombination of electrons and holes in solar cells and lower the rate determining step of the oxygen evolution reaction on rutile TiO₂ surfaces.^94,95

Fig. 11 The DFT+U-predicted (a) occupancies of the dopant atom t_2g and e_g orbitals, in the 3d, 4d or 5d subshell, for all extrinsic defects in anatase (circles) and rutile (triangles), calculated using the mBEEF density functional, U = 2.575 eV, c₁ = 0.752 and c₂ = −0.486. The corresponding elemental species projected density of states for (b) , (c) , (d) , (e) , (f) , (g) , (h) , (i) in bulk anatase (solid lines) and rutile (dashed lines) TiO₂, are normalised with respect to the different defect concentrations in the anatase and rutile simulation supercells. The Fermi level is denoted by the red dashed line.

3.3 Optimising Hubbard U values and projectors from first-principles

3.3.1 Motivation and method reconfiguration. Following the success of the semi-empirical machine learning approach in mitigating numerical instability during simulations of defects in TiO₂ (Section 3.2), we now turn to the development of a fully first-principles strategy for optimising Hubbard U values and projectors. Simultaneous optimisation of Hubbard U values and projectors is particularly attractive as a step towards the long-term goal of formulating DFT+U as a true density functional, in which both the correlated subspace and the on-site Coulomb repulsion are optimally defined by the electron density and orbitals of a material. However, the reformulation of the semi-empirical machine learning approach is contingent on two major modifications. Firstly, it is necessary to redefine the semi-empirical cost function (J^SE) by removing all experimental reference data, but in a manner that will ultimately yield similar results leading to the numerically stable self-consistent defect calculations in Section 3.2.3. Secondly, the termination and SCF non-convergence of defect calculations observed for TiO₂ were also observed for a broad range of TMOs and REOs; therefore, it is necessary to generalise the SVM constraints in Section 3.2 in order to filter out Hubbard parameters that lead to unstable defect calculations for systems beyond TiO₂.

To investigate how best to construct a first-principles cost function (J^FP), the DFT+U-predicted Ti 3d and O 2p orbital occupancies in anatase and rutile TiO₂, using a refined and atomic Ti 3d Hubbard projector, were compared with the occupancies calculated using hybrid-DFT, as illustrated in Fig. 12. Fig. 12(a) and (c) show that DFT+U using an atomic Ti 3d Hubbard projector predicts Ti 3d orbital occupancies closer to those calculated using hybrid-DFT, although this is expected as the hybrid-DFT-calculated occupancies are derived from the atomic Ti 3d Hubbard projector and so cannot be compared in a like-for-like manner to those derived from a modified Hubbard projector (detailed in the SI Section S1.3). Fig. 12(b) and (d) show that DFT+U using the refined Ti 3d Hubbard projector predicts O 2p orbital occupancies closer to those predict using hybrid-DFT, with errors <5% for all values of m, which is a like-for-like comparison as the O 2p orbital occupancies are derived from the atomic O 2p Hubbard projector across all methods (detailed in the SI Section S1.3). The first-principles cost function was therefore defined as the average error in the DFT+U-predicted O 2p orbital occupancies vs. hybrid-DFT (defined in Section 2.3.1), which is obtained from a single reference calculation of a geometry optimised unit cell. To investigate the requirement for SVM constraints in a first-principles approach for optimising Hubbard U values and projectors, the dataset for classifying the validity of bulk oxygen vacancy calculations in anatase TiO₂ (Fig. 12(e)) was considered again in terms of percentage errors of Tr[n(Ti 3d)] and Tr[n(O 2p)] vs. hybrid-DFT (Fig. 12(f)). Fig. 12(f) shows that minimising J^FP towards zero is necessary but not a sufficient condition to ensure the numerical stability of point defect calculations in anatase TiO₂, highlighting the requirement of satisfying SVM-derived constraints as well as minimising J^FP for accurate Hubbard parameter optimisation.

Fig. 12 Comparing the percentage errors of the DFT+U-predicted Ti 3d and O 2p orbital occupancies in anatase and rutile TiO₂vs. hybrid-DFT (PBE0 density functional); (a) Ti 3d in anatase TiO₂, (b) O 2p in anatase TiO₂, (c) Ti 3d in rutile TiO₂ and (d) O 2p in rutile TiO₂. Blue bars correspond to DFT+U using the mBEEF density functional, U = 2.575, c₁ = 0.752 and c₂ = −0.486. Red bars correspond to DFT+U using the mBEEF density functional, U = 2.575, c₁ = 1 and c₂ = 0. The same outcomes of bulk oxygen vacancy calculations plotted in Fig. 7 are plot in (e) and (f), where (e) is plotted in terms of the raw U, Tr[n(Ti 3d)] and Tr[n(O 2p)] values, whilst (f) is plotted in terms of U and the percentage errors of Tr[n(Ti 3d)] and Tr[n(O 2p)] vs. hybrid-DFT.

In the following sections, we outline the generalisation of both the regression and classification methods adopted in Section 3.2, extending them towards the development of a single unified model that can automatically parameterise Hubbard U values and projectors from material-dependent inputs, i.e., in the spirit of a foundation model for DFT+U parameterisation. This serves as a test of how well the insights gained from TiO₂ transfer across chemical space, with the ultimate goal of enabling accurate, self-consistent simulations of defects and polarons across a broad range of TMOs and REOs. We note that such a framework could equally be formulated as an iterative active learning scheme for systems that are difficult to incorporate within a single unified model, as outlined in Section 4.2.

3.3.2 Generalised symbolic regression. Estimating the DFT+U-predicted O 2p occupancies was achieved using three HI-SISSO correlations, which were constructed using the expanded primary feature set of Hubbard parameters, basis set parameters, DFT predicted orbital occupancies and atomic material descriptors (listed in the SI Section S2.2.1), before searching the landscape of J^FP_predicted for the ten materials listed in Table S2 in the SI Section S1.3, using the corresponding material-dependent descriptors and hybrid-DFT reference orbital occupancies. The outputs of the linear search were families of solutions for each material (discussed in the SI Section S2.2.2), of which 94 combinations of U, c₁ and c₂ were validated using DFT+U to compare the accuracy of the J^FP_predicted values (using the HI-SISSO correlations) versus J^FP_validated values (from DFT+U calculations), achieving an average MAE across all materials of 5.02% in Fig. 13(a).

Fig. 13 (a) Parity plot of the HI-SISSO-predicted (J^FP_predicted) and DFT+U-validated (J^FP_validated) cost function across ten TMOs and REOs using the generalised approach. (b) Comparison of the MAE of J^FP_predicted for each material in (a) vs. the corresponding Mahalanobis distance (D_M), averaged over all tested combinations of Hubbard parameters, to visualise the dependence of the accuracy of the generalised approach on the training set size for each material.

To investigate the relationship between the residuals in Fig. 13(a) and the size of the training set per material, the MAE was compared with D_M (both averaged per material) to quantify whether Hubbard projector optimisation is interpolative or extrapolative. As illustrated in Fig. 13(b), the generalised approach is more interpolative for anatase TiO₂, rutile TiO₂ and ZrO₂, which have the three smallest values of D_M. Large values of D_M, as seen for CeO₂, Cu₂O and LiCoO₂, indicate the generalised approach is comparatively extrapolative for these materials. As there is no correlation between the MAE and D_M, the larger residuals in Fig. 13(a) for ZrO₂ and MoO₃ do not appear to be related to the training set size and are likely due to the choice of primary features used in HI-SISSO, such as the lack of structure-dependent features, and requires further study to enhance the method transferability.

3.3.3 Generalised symbolic classification. Given that the minimisation of J^FP is necessary but not sufficient for numerically stable point defect calculations in TiO₂, as illustrated in Fig. 12(f), the constraints on the DFT+U-predicted covalency established in Fig. 7 needed to be generalised across materials. Generalised constraints were determined using the primary features U and J^FP, as well as two additional material-dependent descriptors of covalency, including the average error in the DFT+U-predicted metal d or f orbital occupancies compared to hybrid-DFT (E^metal) and the ratio of the traces of the metal d or f and O 2p occupation matrices predicted using hybrid-DFT (R^hybrid), as defined in Section 2.3.2. With the primary features U, J^FP, E^metal and R^hybrid, generalised classification was performed to predict the same outcomes shown in Fig. 7 for bulk oxygen vacancy calculations, but across TiO₂, CeO₂, ZrO₂, MoO₃, WO₃ and Cu₂O, as illustrated in Fig. 14(a) and (b), where the unitless SVM boundaries S₄ and S₅ are defined in the SI Section S2.2.1 and the constant α = 1 eV⁻¹ is introduced to ensure dimensional consistency. S₄ and S₅ are then combined to form a constraint on the generalised Hubbard parameter space, to ensure the numerical stability of point defect calculations using the OMR method:

S₄ ≥ 0 ∧ S₅ ≥ 0 (15)

Fig. 14 The linear boundaries (a) S₄ and (b) S₅ that classify the numerical stability of DFT+U simulations of a bulk oxygen vacancy in TiO₂, CeO₂, ZrO₂, MoO₃, WO₃ and Cu₂O, separating regions in the DFT+U-computed feature space that lead to successful convergence, termination due to an unphysical ground state, and charge sloshing that prevents SCF convergence.

3.3.4 Electron and hole polarons in LiCo_1−xMg_xO_2−x. The generalised regression and classification models were then integrated to perform a one-shot screening of the Hubbard parameter space for a previously unseen material, to test the accuracy and numerical stability of DFT+U simulations of the stoichiometric and defective electronic structure. We chose the common Li-ion battery cathode material LiCoO₂, which is the focus of wide-ranging studies to enhance its electrochemical properties, such as electrical conductivity and long-term cycling stability, via the introduction of low-valence ions resulting in charge compensation from the formation of point defects.^96–98 For example, substitution of Co³⁺ in LiCoO₂ with Mg²⁺ is proposed to result in the formation of Co⁴⁺ holes, that exist delocalised^99–101 and are attributed to the 100× enhanced electrical conductivity of LiCo_1−xMg_xO₂ compared to pristine LiCoO₂.¹⁰² Other reports suggest that charge compensation and enhancements in electrical conductivity occur via the formation of oxygen vacancies, whilst Mg doping does not cause a change in the oxidation state of Co atoms.^99,103

Literature reported DFT+U studies, in a planewave basis, have yet to provide an unambiguous rationale for the increased electrical conductivity of Mg-doped LiCoO₂, with significant sensitivity of the DFT+U-predicted localisation of Co⁴⁺ holes with respect to the choice of Hubbard parameters. For example, DFT+U using a Co 3d Hubbard U value determined from first-principles using LR-cDFT (between 4.91 eV to 5.62 eV) predicts deep Co 3d defect states within the LiCoO₂ band gap, which is inconsistent with the experimentally observed high electrical conductivity of LiCo_1−xMg_xO₂, whilst DFT+U using smaller Hubbard U values predicts shallow defect states, supporting the experimental observations.^5,104 Whilst no explanation for the sensitivity of DFT+U simulations of LiCo_1−xMg_xO₂ was given, Hoang and Johannes used DFT+U to simulate the self-trapping of Ni⁴⁺ holes in LiNiO₂, where the DFT+U-predicted instability of Ni⁴⁺ was attributed to the predicted O 2p character of the LiNiO₂ valence band edge, highlighting the importance of considering band edge effects when choosing appropriate Hubbard parameters.¹⁰⁵

To test the DFT+U-predicted localisation of hole polarons in LiCo_1−xMg_xO₂ in a NAO framework, the one-shot approach was first applied to determine an appropriate Co 3d Hubbard U value and projector for stoichiometric LiCoO₂. The parameters were determined using a linear search of the HI-SISSO-predicted cost function (J^FP_predicted) for U between 0 eV and 5 eV, c₁ between 0.5 and 1 and c₂ between 0 and -0.6, before all Hubbard parameters that violate eqn (15) were discarded. The output of this initial screening is a region of the Hubbard parameter space (U, c₁ and c₂) that minimises J^FP_predicted and satisfies eqn (15) but contains multiple possible solutions. We therefore narrow this region down to three candidate Hubbard parameters for validation using K-means clustering of a reduced subset of the screened Hubbard parameter space that corresponds to the lowest 10% of J^FP_predicted (illustrated in the SI Section S2.2.2).

From the three candidate Hubbard parameters determined using K-means clustering, DFT+U using U = 3.342 eV, c₁ = 0.792 and c₂ = −0.506, which is henceforth referred to as “refined-DFT+U”, provided the best accuracy in the predicted E_bg (4.18 eV), V₀ (97.38 ų) and ΔE_form (−2.40 eV per atom) after validating against experimental references (listed in the SI Section S2.2.2). Comparing the DFT+U-predicted electronic structure of stoichiometric LiCoO₂ using refined-DFT+U and U = 3.342 eV, c₁ = 1 and c₂ = 0 (i.e., the same Co 3d Hubbard U value with an atomic projector, which is henceforth referred to as “atomic-DFT+U”), projector refinement was found to be essential to predict the mixed Co 3d and O 2p valence band edge character of LiCoO₂ in Fig. 15(a), which is also predicted using hybrid-DFT and reported experimentally.¹⁰⁶ Conversely, DFT+U using atomic-DFT+U, or a much larger Hubbard U value than in Fig. 15(a), both predict the LiCoO₂ valence band edge to be dominated by O 2p states, as illustrated in Fig. 15(b).

Fig. 15 Charge density isosurface at the 0.05 e Å⁻³ level for the eigenstate corresponding to the HOMO and the corresponding Mulliken-projected band structure for stoichiometric LiCoO₂ along the high-symmetry k-point path Γ–M–K–Γ–A–L–H–A–L–M–K–H, calculated using (a) DFT+U with the mBEEF density functional, U = 3.342 eV, c₁ = 0.792 and c₂ = −0.506 and (b) DFT+U with the mBEEF density functional, U = 3.342 eV, c₁ = 1 and c₂ = 0. Marker sizes and colours in the band structure plots correspond to the relative contribution for that species to the band. The valence band edge character is either (a) a mixture of Co 3d and O 2p states or (b) dominated by O 2p states. The band structure plots are centred with respect to the Fermi level.

The defective LiCoO₂ bulk containing was simulated using both refined-DFT+U and atomic-DFT+U, using the OMR method with the 3d_z² orbital occupancy of the Co atom nearest the Mg dopant reduced by 1 to initialise Co⁴⁺. In both simulations, SCF convergence and geometry optimisation completed successfully without termination or charge sloshing. Both simulations also predict the formation of shallow defect states at the top of the valence band in Fig. 16(b), corresponding to a delocalised hole polaron that exists hybridised between Co 3d and O 2p states using refined-DFT+U in Fig. 16(d), compared to non-hybridised across only O atoms using atomic-DFT+U in Fig. 16(e). These differences in the DFT+U-predicted character of the hole polaron also directly affect their stability, with ΔE_defect = 1.48 eV using refined-DFT+U and ΔE_defect = 10.89 eV using atomic-DFT+U. Given that compositions with a 1 : 1 ratio of Co and Mg, i.e., LiCo_0.5Mg_0.5O₂, have been synthesised experimentally,¹⁰⁷ DFT+U with the default atomic Co 3d Hubbard projector therefore incorrectly predicts the total insolubility of Mg in LiCoO₂.

Fig. 16 Elemental species projected density of states for (a) stoichiometric LiCoO₂, (b) defective LiCoO₂ containing and (c) defective LiCoO₂ containing both and , calculated using DFT+U with the mBEEF density functional, U = 3.342 eV, c₁ = 0.792 and c₂ = −0.506 (refined-DFT+U, solid lines) and DFT+U with the mBEEF density functional, U = 3.342 eV, c₁ = 1 and c₂ = 0 (atomic-DFT+U, dotted lines). All plots are relative to the Fermi level (red dashed line). The corresponding charge density isosurfaces at the 0.025 e Å⁻³ level for the eigenstate corresponding to the HOMO are shown for defective LiCoO₂ containing , calculated using (d) refined-DFT+U and (e) atomic-DFT+U, as well as (f) defective LiCoO₂ containing both and , calculated using refined-DFT+U.

Both refined- and atomic-DFT+U were further used to simulate a bulk oxygen vacancy in Mg-doped LiCoO₂ using the OMR method. Here, the nearest O atom to the Mg dopant was removed and the two nearest Co atoms to the vacancy were initialised as Co²⁺ by increasing their 3d_yz orbital occupancies by 1. With the atomic Co 3d Hubbard projector, the simulation did not converge due to charge sloshing, i.e., oscillations in the charge density illustrated in the SI Section S2.2.2, similar to TiO₂ in Section 3.1.2. Conversely, DFT+U using the refined Co 3d Hubbard projector successfully converged without termination or charge sloshing, predicting a deep defect state of 3d_x²−y² character and an associated oxygen vacancy formation energy of ΔE_OV = 3.61 eV. These results suggest that both Mg-doping and oxygen vacancy formation could enhance the electrical conductivity of LiCoO₂via the transport of delocalised hole and localised electron polarons, respectively, whilst the generalised approach for optimising the Co 3d Hubbard U value and projector is robust with respect to numerical instability that is often observed when using the default atomic Hubbard projector. These results give promise to the development of the first-principles approach by extending the size and diversity of the training set used for generalised regression and classification.

4 Future work

Whilst we have demonstrated how supervised machine learning can be applied to simultaneously optimise Hubbard U values and projectors to improve the accuracy and numerical stability of DFT+U in a NAO framework, we note clear avenues for further research that are not accounted for in the current work.

4.1 Extension to challenging open-shell systems

The methods and results presented generally hold well for closed-shell systems or those with ions in a low-spin state. In contrast, we encountered significant challenges applying DFT+U to open-shell systems containing ions in a high spin state, e.g., Co²⁺ in CoO, Mn⁴⁺ in MnO₂, Fe³⁺ in Fe₂O₃ and Cr³⁺ in Cr₂O₃. For these materials, DFT+U was found to consistently suffer from SCF non-convergence as a result of charge sloshing, i.e., excessive oscillations in the charge density, which even occurred during single point calculations of small unit cells, irrespective of the chosen Hubbard U values and projectors. Upon analysing the DFT-predicted ground state electronic structures for these materials, all cases are incorrectly predicted to be metallic with zero band gap, therefore we hypothesise that the starting electronic structure from DFT makes SCF convergence for DFT+U very challenging in a NAO framework. After further testing, we find that switching to spin polarisation in combination with meta-GGA exchange–correlation density functionals and carefully selected SCF mixing parameters can improve the DFT-predicted electronic structure to restore a bandgap (albeit which is underestimated compared to experiment). These changes make the convergence of DFT+U with refined Hubbard projectors feasible and will be the focus of a future study translating this work to defects in magnetic TMOs and REOs. Whilst this remains an unresolved challenge at the moment, we note that the restoration of semiconducting or insulating ground states from an initial DFT-predicted metallic state is routinely achievable in planewave implementations of DFT+U,¹⁰⁸ therefore a more detailed investigation into exact differences between implementations may help.

4.2 Universal models vs. active learning

In Section 3.3, we outline the development of a single universal model for parameterising Hubbard U values and projectors across materials. This approach served as a test for how well the insights gained from TiO₂ transfer across chemical space, with the ultimate goal of enabling accurate, self-consistent simulations of defects and polarons across a broad range of TMOs and REOs. A valid alternative approach could be an iterative active learning scheme, where U, c₁ and c₂ are progressively refined so that the DFT+U-predicted electronic structure resembles that from hybrid-DFT, e.g., via comparison of orbital occupation numbers or band structures. This active learning approach may work well for systems that are not accurately incorporated within a single universal model, but would require repeating for different materials and incur a potentially significant cost for the total number of DFT+U calculations required. We have shown that Bayesian optimisation can efficiently sample a SISSO-derived semi-empirical cost function in Section 3.2.2, therefore this could be repurposed to iteratively perform DFT+U geometry optimisation calculations for a small unit cell, sampling Hubbard parameters and targeting orbital occupancies calculated using hybrid-DFT. If the cost function landscape is very difficult to optimise with an iterative approach, both the universal and active learning schemes could be combined, i.e., the outputs from the one-shot approach used as initial positions for active learning, which may improve convergence to good solutions whilst requiring fewer DFT+U iterations.

4.3 Compatibility with high-throughput simulations

The current framework outputs several combinations of Hubbard U values and projectors, which serve as optimised candidate parameters for validation with the aim of enabling self-consistent defect simulations. Further development into a high-throughput optimisation requires an improvement in accuracy of the generalised symbolic regression model on out-of-training data, which should be feasible since this work considers a very small training set of 197 sets of DFT+U-calculated orbital occupancies from geometry optimised unit cells, of which nearly half correspond to TiO₂. With a better extrapolative accuracy, the generalised approach could be used to output a single Hubbard parameter, rather than relying on K-means clustering to output several candidate parameters for further validation. The reliance on small unit cells further reduces the computational overhead associated with generating corresponding hybrid-DFT reference data.

4.4 Choice of high-level reference and extension to other codes

In the current work, we chose hybrid-DFT as the high-level reference method to obtain target orbital occupancies for metal d or f and O 2p states. This choice was made in order to establish a consistent and computationally tractable benchmark against which machine learning-optimised U values and projectors could be compared. While this introduces a dependency on the choice of reference, it is entirely possible within our framework to select another level of theory depending on the target material and available computational resources. The framework could equally be extended to optimise against alternative physics-informed quantities, e.g., minimising errors in band structure plots, enforcing piecewise linearity of the total energy with respect to fractional occupancy of the correlated subspace, or maximising occupation matrix idempotency. In the case of magnetic systems, for which identifying the true ground state is complicated by the presence of metastable states in the potential energy surface, one could instead target experimentally observed magnetic orderings as the optimisation criteria.

The methods presented could also be extended to other electronic structure codes, for which there already exist several schemes that require simultaneous parameter optimisation, e.g., DFT+U+V and orbital resolved DFT+U.^32,109 In such frameworks, rather than tuning linear expansion coefficients as in our NAO-based approach, one would instead adjust the definition of the Hubbard projector via the projector augmented wave (PAW) augmentation radius in a planewave basis,¹⁷ or the muffin-tin radius in a linear augmented planewave (LAPW) basis, both of which have been reported to significantly influence the predicted geometric, energetic and electronic properties of complex oxides.^16,110 Direct adaptation of the framework for other electronic structure codes will be essential, as Hubbard parameters are generally not transferable across codes with different basis sets and definitions of the Hubbard projector. As discussed by Kick et al., DFT+U in a NAO framework typically requires U values ∼1–2 eV smaller than other formalisms to achieve comparable charge localisation.¹⁵ The observation agrees with the results of this work concerning Ti 3d Hubbard parameters for simulating TiO₂ (U = ∼ 2.5–3 eV) vs. literature reported values computed in a real-space formalism (U = ∼ 5–6 eV).¹¹¹

5 Conclusions

To navigate the numerical instability of self-consistent DFT+U simulations of TMOs and REOs, including the prediction of erroneous metallic ground states and the termination or non-convergence of point defect calculations, it is essential to carefully define an appropriate Hubbard projector. Simultaneous optimisation of the Hubbard U value and projector has been demonstrated semi-empirically using SR and SVMs, before using BO to minimise the errors of target properties relative to experimental references. The Ti 3d Hubbard U value and projector have been semi-empirically refined to enable self-consistent DFT+U simulations of intrinsic and extrinsic defects in both anatase and rutile TiO₂. The outcome is DFT+U simulations with comparable accuracy to hybrid-DFT in terms of the relative stabilities of point defects and the formation of localised Holstein polarons, but at orders of magnitude lower cost. The DFT+U-predicted occupation matrices reproduce the hybrid-DFT O 2p occupation matrix, which can be an effective cost function for a first-principles strategy for Hubbard U value and projector optimisation.

The semi-empirical approach was therefore defined as a first-principles approach and generalised across materials using hierarchical SR to screen the Hubbard parameter space using empirical correlations to learn the DFT+U potential energy surface in terms of orbital occupancies. Predictions of metal d or f orbital and O 2p orbital occupancies were made in terms of Hubbard parameters, basis set parameters, DFT-predicted orbital occupancies and atomic material descriptors. The first-principles approach enables the development of a generalised workflow for the one-shot computation of Hubbard U values and projectors for different materials that requires no successive DFT+U calculations, as in active learning schemes. The method transferability is demonstrated for 10 prototypical TMOs and REOs (anatase and rutile TiO₂, Cu₂O, MoO₃, WO₃, Y₂O₃, ZrO₂, CeO₂, LiCoO₂ and LiFePO₄), which each require one reference DFT and hybrid-DFT calculation as inputs, whilst generating families of solutions for each material, i.e., optimised Hubbard projectors for a given Hubbard U value. Upon validating a subset of these solutions, a MAE of 5.02% for the DFT+U-predicted O 2p orbital occupancies was achieved, with demonstrated accuracy for materials unseen from model training (LiCoO₂ and LiFePO₄).

Predicting the numerical stability of point defect calculations can also be generalised across materials using symbolic classification, using Hubbard U values and material-dependent descriptors of covalency, enabling the determination of Hubbard U values and projectors that are robust against numerical instability. The validity of Hubbard U values and projectors determined from first-principles has been investigated for the self-consistent simulation of Mg-doped and oxygen deficient LiCoO₂, where refining the Co 3d Hubbard projector enables the numerically stable simulation of experimentally reported charge compensation mechanisms driving the material's high electrical conductivity. The same results were not possible using an atomic Co 3d Hubbard projector and did not require any prior testing of suitable Co 3d Hubbard U values or projectors, which gives promise for the development of a foundational tool for simultaneously determining multiple Hubbard parameters in a NAO framework and beyond. The work demonstrates how advanced machine learning algorithms can aid the development of inexpensive and transferable workflows for DFT+U parameterisation, achieving extrapolative accuracy beyond the limits of small training sets, for more accurate and efficient simulations of complex energy materials.

Author contributions

All authors contributed to the conceptualisation of the project and to software choices and method development in this work. AC performed the electronic structure calculations presented in this work and performed the analysis of results. All authors contributed to the preparation of the manuscript.

Conflicts of interest

There are no conflicts to declare.

Data availability

All Python scripts for global optimisation, datasets for regression/classification and input/output files for electronic structure calculations have been uploaded to the GitHub repository https://github.com/amitmc1/Hubbardprojectors, with an archived version accessible at the DOI: https://doi.org/10.5281/zenodo.17454878.

Supplementary information: information on the determination of DFT and DFT+U parameters, settings and outputs for machine learning methods, and combined application in a first-principles screening approach for the DFT+U parameters. See DOI: https://doi.org/10.1039/d5dd00292c.

Acknowledgements

We thank Harald Oberhofer, Matthias Kick, and Maximilian Brand for valuable scientific discussions regarding the implementation of DFT+U in FHI-aims. We also acknowledge funding by the Prosperity Partnership project Sustainable Catalysis for Clean Growth, funded by the UK Engineering and Physical Sciences Research Council (EPSRC), bp through the bp International Centre for Advanced Materials (bp-ICAM) and Johnson Matthey plc in collaboration with Cardiff University and The University of Manchester (EPSRC grant number EP/V056565/1). A. J. L. acknowledges funding by the UKRI Future Leaders Fellowship program (MR/T018372/1 and MR/Y034279/1). We acknowledge computational resources and support from the Supercomputing Wales project, which is part-funded by the European Regional Development Fund (ERDF) via the Welsh Government; and the UK National Supercomputing Services ARCHER and ARCHER2, accessed via membership of the Materials Chemistry Consortium, which is funded by Engineering and Physical Sciences Research Council (EP/L000202/1, EP/R029431/1 and EP/T022213/1).

Notes and references

1. K. Chang, H. Zhang, M.-j. Cheng and Q. Lu, Application of Ceria in CO₂ Conversion Catalysis, ACS Catal., 2020, 10, 613–631.
2. A. Avani and E. Anila, Recent advances of MoO₃ based materials in energy catalysis: Applications in hydrogen evolution and oxygen evolution reactions, Int. J. Hydrogen Energy, 2022, 47, 20475–20493.
3. M. Lindley, P. Stishenko, J. W. M. Crawley, F. Tinkamanyire, M. Smith, J. Paterson, M. Peacock, Z. Xu, C. Hardacre, A. S. Walton, A. J. Logsdail and S. J. Haigh, Tuning the Size of TiO₂-Supported Co Nanoparticle Fischer–Tropsch Catalysts Using Mn Additions, ACS Catal., 2024, 14, 10648–10657.
4. A. Folli, J. Z. Bloh, A. Lecaplain, R. Walker and D. E. Macphee, Properties and photochemistry of valence-induced-Ti³⁺ enriched (Nb,N)-codoped anatase TiO₂ semiconductors, Phys. Chem. Chem. Phys., 2015, 17, 4849–4853.
5. J. A. Santana, J. Kim, P. R. C. Kent and F. A. Reboredo, Successes and failures of Hubbard-corrected density functional theory: The case of Mg doped LiCoO₂, J. Chem. Phys., 2014, 141, 164706.
6. W. Kohn and L. J. Sham, Self-Consistent Equations Including Exchange and Correlation Effects, Phys. Rev., 1965, 140, A1133–A1138.
7. W. Kohn, A. D. Becke and R. G. Parr, Density Functional Theory of Electronic Structure, J. Phys. Chem., 1996, 100, 12974–12980.
8. N. L. Nguyen, N. Colonna, A. Ferretti and N. Marzari, Koopmans-Compliant Spectral Functionals for Extended Systems, Phys. Rev. X, 2018, 8, 021051.
9. J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 5048–5079.
10. M. Reticcioli, U. Diebold and C. Franchini, Modeling polarons in density functional theory: lessons learned from TiO₂, J. Condens. Matter Phys., 2022, 34, 204006.
11. H. J. Kulik, Perspective: Treating electron over-delocalization with the DFT+U method, J. Chem. Phys., 2015, 142, 240901.
12. M. Capdevila-Cortada, Z. Łodziana and N. López, Performance of DFT+U Approaches in the Study of Catalytic Materials, ACS Catal., 2016, 6, 8370–8379.
13. D. S. Lambert and D. D. O'Regan, Use of DFT + U + J with linear response parameters to predict non-magnetic oxide band gaps with hybrid-functional accuracy, Phys. Rev. Res., 2023, 5, 013160.
14. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Electron-energy-loss spectra and the structural stability of nickel oxide: An LSDA+U study, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
15. M. Kick, K. Reuter and H. Oberhofer, Intricacies of DFT+U, Not Only in a Numeric Atom Centered Orbital Framework, J. Chem. Theory Comput., 2019, 15, 1705–1718.
16. K. Park, M. Raman, A.-J. Olatunbosun and J. Pohlmann, Revisiting DFT+U calculations of TiO₂ and the effect of the local-projection size, AIP Adv., 2024, 14, 065114.
17. Z. Wang, C. Brock, A. Matt and K. H. Bevan, Implications of the DFT + U method on polaron properties in energy materials, Phys. Rev. B, 2017, 96, 125150.
18. I. Timrov, N. Marzari and M. Cococcioni, Self-consistent Hubbard parameters from density-functional perturbation theory in the ultrasoft and projector-augmented wave formulations, Phys. Rev. B, 2021, 103, 045141.
19. D. D. O'Regan, N. D. M. Hine, M. C. Payne and A. A. Mostofi, Projector self-consistent DFT + U using nonorthogonal generalized Wannier functions, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 081102.
20. Y.-Y. Ting and P. M. Kowalski, Refined DFT+U method for computation of layered oxide cathode materials, Electrochim. Acta, 2023, 443, 141912.
21. G. L. Murphy, Z. Zhang, R. Tesch, P. M. Kowalski, M. Avdeev, E. Y. Kuo, D. J. Gregg, P. Kegler, E. V. Alekseev and B. J. Kennedy, Tilting and Distortion in Rutile-Related Mixed Metal Ternary Uranium Oxides: A Structural, Spectroscopic, and Theoretical Investigation, Inorg. Chem., 2021, 60, 2246–2260.
22. I. Timrov, F. Aquilante, L. Binci, M. Cococcioni and N. Marzari, Pulay forces in density-functional theory with extended Hubbard functionals: From nonorthogonalized to orthogonalized manifolds, Phys. Rev. B, 2020, 102, 235159.
23. O. Y. Long, G. Sai Gautam and E. A. Carter, Evaluating optimal U for 3d transition-metal oxides within the SCAN+U framework, Phys. Rev. Mater., 2020, 4, 045401.
24. L. Wang, T. Maxisch and G. Ceder, Oxidation energies of transition metal oxides within the GGA + U framework, Phys. Rev. B: Condens. Matter Mater. Phys., 2006, 73, 195107.
25. M. Cococcioni and S. de Gironcoli, Linear response approach to the calculation of the effective interaction parameters in the LDA + U method, Phys. Rev. B: Condens. Matter Mater. Phys., 2005, 71, 035105.
26. I. Timrov, N. Marzari and M. Cococcioni, Hubbard parameters from density-functional perturbation theory, Phys. Rev. B, 2018, 98, 085127.
27. M. Springer and F. Aryasetiawan, Frequency-dependent screened interaction in Ni within the random-phase approximation, Phys. Rev. B: Condens. Matter Mater. Phys., 1998, 57, 4364–4368.
28. N. J. Mosey, P. Liao and E. A. Carter, Rotationally invariant ab initio evaluation of Coulomb and exchange parameters for DFT+U calculations, J. Chem. Phys., 2008, 129, 014103.
29. L. A. Agapito, S. Curtarolo and M. Buongiorno Nardelli, Reformulation of DFT + U as a Pseudohybrid Hubbard Density Functional for Accelerated Materials Discovery, Phys. Rev. X, 2015, 5, 011006.
30. G. C. Moore, M. K. Horton, E. Linscott, A. M. Ganose, M. Siron, D. D. O'Regan and K. A. Persson, High-throughput determination of Hubbard U and Hund J values for transition metal oxides via the linear response formalism, Phys. Rev. Mater., 2024, 8, 014409.
31. K. Yu and E. A. Carter, Communication: Comparing ab initio methods of obtaining effective U parameters for closed-shell materials, J. Chem. Phys., 2014, 140, 121105.
32. E. Macke, I. Timrov, N. Marzari and L. C. Ciacchi, Orbital-Resolved DFT+U for Molecules and Solids, J. Chem. Theory Comput., 2024, 20, 4824–4843.
33. C. Ricca, I. Timrov, M. Cococcioni, N. Marzari and U. Aschauer, Self-consistent site-dependent DFT+U study of stoichiometric and defective SrMnO₃, Phys. Rev. B, 2019, 99, 094102.
34. M. Yu, S. Yang, C. Wu and N. Marom, Machine learning the Hubbard U parameter in DFT+ U using Bayesian optimization, npj Comput. Mater., 2020, 6, 180.
35. W. Yu, Z. Zhang, X. Wan, H. Guo, Q. Gui, Y. Peng, Y. Li, W. Fu, D. Lu, Y. Ye and Y. Guo, Active Learning the High-Dimensional Transferable Hubbard U and V Parameters in the DFT+U+V Scheme, J. Chem. Theory Comput., 2023, 19, 6425–6433.
36. P. Tavadze, R. Boucher, G. Avendaño-Franco, K. X. Kocan, S. Singh, V. Dovale-Farelo, W. Ibarra-Hernández, M. B. Johnson, D. S. Mebane and A. H. Romero, Exploring DFT+U parameter space with a Bayesian calibration assisted by Markov chain Monte Carlo sampling, npj Comput. Mater., 2021, 7, 182.
37. G. Cai, Z. Cao, F. Xie, H. Jia, W. Liu, Y. Wang, F. Liu, X. Ren, S. Meng and M. Liu, Predicting structure-dependent Hubbard U parameters via machine learning, Mater. Futures, 2024, 3, 025601.
38. M. Uhrin, A. Zadoks, L. Binci, N. Marzari and I. Timrov, Machine learning Hubbard parameters with equivariant neural networks, npj Comput. Mater., 2025, 11, 19.
39. L. Foppa, T. A. R. Purcell, S. V. Levchenko, M. Scheffler and L. M. Ghiringhelli, Hierarchical Symbolic Regression for Identifying Key Physical Parameters Correlated with Bulk Properties of Perovskites, Phys. Rev. Lett., 2022, 129, 055301.
40. V. Blum, R. Gehrke, F. Hanke, P. Havu, V. Havu, X. Ren, K. Reuter and M. Scheffler, Ab initio molecular simulations with numeric atom-centered orbitals, Comput. Phys. Commun., 2009, 180, 2175–2196.
41. A. H. Larsen, J. J. Mortensen, J. Blomqvist, I. E. Castelli, R. Christensen, M. Dułak, J. Friis, M. N. Groves, B. Hammer, C. Hargus, E. D. Hermes, P. C. Jennings, P. B. Jensen, J. Kermode, J. R. Kitchin, E. L. Kolsbjerg, J. Kubal, K. Kaasbjerg, S. Lysgaard, J. B. Maronsson, T. Maxson, T. Olsen, L. Pastewka, A. Peterson, C. Rostgaard, J. Schiøtz, O. Schütt, M. Strange, K. S. Thygesen, T. Vegge, L. Vilhelmsen, M. Walter, Z. Zeng and K. W. Jacobsen, The atomic simulation environment—a Python library for working with atoms, J. Phys.: Condens. Matter, 2017, 29, 273002.
42. O. Lamiel-Garcia, K. C. Ko, J. Y. Lee, S. T. Bromley and F. Illas, When Anatase Nanoparticles Become Bulklike: Properties of Realistic TiO₂ Nanoparticles in the 1–6 nm Size Range from All Electron Relativistic Density Functional Theory Based Calculations, J. Chem. Theory Comput., 2017, 13, 1785–1793.
43. J. Wellendorff, K. T. Lundgaard, K. W. Jacobsen and T. Bligaard, mBEEF: An accurate semi-local Bayesian error estimation density functional, J. Chem. Phys., 2014, 140, 144107.
44. J. P. Perdew, A. Ruzsinszky, G. I. Csonka, O. A. Vydrov, G. E. Scuseria, L. A. Constantin, X. Zhou and K. Burke, Restoring the Density-Gradient Expansion for Exchange in Solids and Surfaces, Phys. Rev. Lett., 2008, 100, 136406.
45. S. Lehtola, C. Steigemann, M. J. Oliveira and M. A. Marques, Recent developments in libxc — A comprehensive library of functionals for density functional theory, SoftwareX, 2018, 7, 1–5.
46. C. Adamo and V. Barone, Toward reliable density functional methods without adjustable parameters: The PBE0 model, J. Chem. Phys., 1999, 110, 6158–6170.
47. X. Ren, P. Rinke, V. Blum, J. Wieferink, A. Tkatchenko, A. Sanfilippo, K. Reuter and M. Scheffler, Resolution-of-identity approach to Hartree–Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions, New J. Phys., 2012, 14, 053020.
48. C. G. Broyden, The Convergence of a Class of Double-rank Minimization Algorithms 1. General Considerations, IMA J. Appl., 1970, 6, 76–90.
49. R. Fletcher, A new approach to variable metric algorithms, Comput. J., 1970, 13, 317–322.
50. D. F. Shanno, Conditioning of Quasi-Newton Methods for Function Minimization, Math. Comput., 1970, 24, 647–656.
51. D. Goldfarb, A Family of Variable-Metric Methods Derived by Variational Means, Math. Comput., 1970, 24, 23–26.
52. F. Birch, Finite Elastic Strain of Cubic Crystals, Phys. Rev., 1947, 71, 809–824.
53. G. L. Murphy, Z. Zhang, R. Tesch, P. M. Kowalski, M. Avdeev, E. Y. Kuo, D. J. Gregg, P. Kegler, E. V. Alekseev and B. J. Kennedy, Tilting and distortion in rutile-related mixed metal ternary uranium oxides: a structural, spectroscopic, and theoretical investigation, Inorg. Chem., 2021, 60, 2246–2260.
54. K. O. Kvashnina, P. M. Kowalski, S. M. Butorin, G. Leinders, J. Pakarinen, R. Bès, H. Li and M. Verwerft, Trends in the valence band electronic structures of mixed uranium oxides, Chem. Commun., 2018, 54, 9757–9760.
55. P. M. Kowalski, Z. He and O. Cheong, Electrode and electrolyte materials from atomistic simulations: properties of Li_xFePO₄ electrode and zircon-based ionic conductors, Front. Energy Res., 2021, 9, 653542.
56. Y.-Y. Ting and P. M. Kowalski, Refined DFT+U method for computation of layered oxide cathode materials, Electrochim. Acta, 2023, 443, 141912.
57. K. Jakob and H. Oberhofer, Master's thesis: Self-Consistency in the Hubbard-Corrected DFT+U Method, Faculty of Chemistry, Technical University of Munich, 2021.
58. W. B. Begna, G. S. Gurmesa and C. A. Geffe, Ortho-atomic projector assisted DFT+U study of room temperature Ferro- and antiferromagnetic Mn-doped TiO₂ diluted magnetic semiconductor, Mater. Res. Express, 2022, 9, 076102.
59. B. Dorado, B. Amadon, M. Freyss and M. Bertolus, DFT + U calculations of the ground state and metastable states of uranium dioxide, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 235125.
60. J. P. Allen and G. W. Watson, Occupation matrix control of d- and f-electron localisations using DFT+U, Phys. Chem. Chem. Phys., 2014, 16, 21016–21031.
61. R. Ouyang, S. Curtarolo, E. Ahmetcik, M. Scheffler and L. M. Ghiringhelli, SISSO: A compressed-sensing method for identifying the best low-dimensional descriptor in an immensity of offered candidates, Phys. Rev. Mater., 2018, 2, 083802.
62. T. A. R. Purcell, M. Scheffler, C. Carbogno and L. M. Ghiringhelli, SISSO++: A C++ Implementation of the Sure-Independence Screening and Sparsifying Operator Approach, J. Open Source Softw., 2022, 7, 3960.
63. T. A. R. Purcell, M. Scheffler and L. M. Ghiringhelli, Recent advances in the SISSO method and their implementation in the SISSO++ code, J. Chem. Phys., 2023, 159, 114110.
64. Y. Zhang, J. W. Furness, B. Xiao and J. Sun, Subtlety of TiO₂ phase stability: Reliability of the density functional theory predictions and persistence of the self-interaction error, J. Chem. Phys., 2019, 150, 014105.
65. T. Arlt, M. Bermejo, M. A. Blanco, L. Gerward, J. Z. Jiang, J. Staun Olsen and J. M. Recio, High-pressure polymorphs of anatase TiO₂, Phys. Rev. B: Condens. Matter Mater. Phys., 2000, 61, 14414–14419.
66. M. Setvin, C. Franchini, X. Hao, M. Schmid, A. Janotti, M. Kaltak, C. G. Van de Walle, G. Kresse and U. Diebold, Direct View at Excess Electrons in TiO₂ Rutile and Anatase, Phys. Rev. Lett., 2014, 113, 086402.
67. A. Chaudhari, A. J. Logsdail and A. Folli, Polymorph-Induced Reducibility and Electron Trapping Energetics of Nb and W Dopants in TiO₂, J. Phys. Chem. C, 2025, 129, 15453–15461.
68. F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot and E. Duchesnay, Scikit-learn: Machine Learning in Python, J. Mach. Learn. Res., 2011, 12, 2825–2830.
69. The GPyOpt authors, GPyOpt: A Bayesian Optimization framework in Python, http://github.com/SheffieldML/GPyOpt, 2016, Accessed June 2024.
70. D. R. Jones, M. Schonlau and W. J. Welch, Efficient Global Optimization of Expensive Black-Box Functions, J. Glob. Optim., 1998, 13, 455–492.
71. M. D. Mckay, R. J. Beckman and W. J. Conover, A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output From a Computer Code, Technometrics, 2000, 42, 55–61.
72. S. Dalton and D. Richardson, pyDOE: The experimental design package for Python, https://pythonhosted.org/pyDOE/, 2014, Accessed June 2024.
73. R. Ramakrishnan, P. O. Dral, M. Rupp and O. A. von Lilienfeld, Big Data Meets Quantum Chemistry Approximations: The -Machine Learning Approach, J. Chem. Theory Comput., 2015, 11, 2087–2096.
74. D.-H. Seo, A. Urban and G. Ceder, Calibrating transition-metal energy levels and oxygen bands in first-principles calculations: Accurate prediction of redox potentials and charge transfer in lithium transition-metal oxides, Phys. Rev. B: Condens. Matter Mater. Phys., 2015, 92, 115118.
75. P. Mahalanobis, On the Generalised Distance in Statistics, Sankhya, Ser. A, 1936, 80, 1–7.
76. V. I. Anisimov, F. Aryasetiawan and A. I. Lichtenstein, First-principles calculations of the electronic structure and spectra of strongly correlated systems: the LDA+U method, J. Phys.: Condens. Matter, 1997, 9, 767–808.
77. B. Himmetoglu, A. Floris, S. de Gironcoli and M. Cococcioni, Hubbard-corrected DFT energy functionals: The LDA+U description of correlated systems, Int. J. Quantum Chem., 2014, 114, 14–49.
78. T. Schäfer, N. Daelman and N. López, Cerium Oxides without U: The Role of Many-Electron Correlation, J. Phys. Chem. Lett., 2021, 12, 6277–6283.
79. R. song Li, D. qiang Xin, S. qi Huang, Z. jian Wang, L. Huang and X. hua Zhou, A full potential all-electron calculation on electronic structure of NiO, Chin. J. Phys., 2018, 56, 2829–2836.
80. G. Gebreyesus, L. Bastonero, M. Kotiuga, N. Marzari and I. Timrov, Understanding the role of Hubbard corrections in the rhombohedral phase of BaTiO₃, Phys. Rev. B, 2023, 108, 235171.
81. T.-c. Liu, D. Gaines, H. Kim, A. Salgado-Casanova, S. B. Torrisi and C. Wolverton, Anomalous reversal of stability in Mo-containing oxides: A difficult case exhibiting sensitivity to DFT + U and distortion, Phys. Rev. Mater., 2025, 9, 055402.
82. S. Phoka, P. Laokul, E. Swatsitang, V. Promarak, S. Seraphin and S. Maensiri, Synthesis, structural and optical properties of CeO₂ nanoparticles synthesized by a simple polyvinyl pyrrolidone (PVP) solution route, Mater. Chem. Phys., 2009, 115, 423–428.
83. F. B. Baker, E. J. Huber, C. E. Holley and N. Krikorian, Enthalpies of formation of cerium dioxide, cerium sesquicarbide, and cerium dicarbide, J. Chem. Thermodyn., 1971, 3, 77–83.
84. M. Kick, C. Grosu, M. Schuderer, C. Scheurer and H. Oberhofer, Mobile Small Polarons Qualitatively Explain Conductivity in Lithium Titanium Oxide Battery Electrodes, J. Phys. Chem. Lett., 2020, 11, 2535–2540.
85. M. Kick, C. Scheurer and H. Oberhofer, Formation and stability of small polarons at the lithium-terminated Li₄Ti₅O₁₂ (LTO) (111) surface, J. Chem. Phys., 2020, 153, 144701.
86. J. Aarons, M. Sarwar, D. Thompsett and C.-K. Skylaris, Perspective: Methods for large-scale density functional calculations on metallic systems, J. Chem. Phys., 2016, 145, 220901.
87. M. Arrigoni and G. K. H. Madsen, A comparative first-principles investigation on the defect chemistry of TiO₂ anatase, J. Chem. Phys., 2020, 152, 044110.
88. T. D. Pham and N. A. Deskins, Efficient Method for Modeling Polarons Using Electronic Structure Methods, J. Chem. Theory Comput., 2020, 16, 5264–5278.
89. R. Zhang, A. Chutia, A. A. Sokol, D. Chadwick and C. R. A. Catlow, A computational investigation of the adsorption of small copper clusters on the CeO₂(110) surface, Phys. Chem. Chem. Phys., 2021, 23, 19329–19342.
90. G. P. Kerker, Efficient iteration scheme for self-consistent pseudopotential calculations, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 23, 3082–3084.
91. K. Warda, E. Macke, I. Timrov, L. C. Ciacchi and P. M. Kowalski, Getting the manifold right: The crucial role of orbital resolution in DFT+ U for mixed df electron compounds, arXiv, 2025, preprint, arXiv:2508.16435,  DOI:10.48550/arXiv.2508.16435.
92. C. Spreafico and J. VandeVondele, The nature of excess electrons in anatase and rutile from hybrid DFT and RPA, Phys. Chem. Chem. Phys., 2014, 16, 26144–26152.
93. C. Bigi, Z. Tang, G. M. Pierantozzi, P. Orgiani, P. K. Das, J. Fujii, I. Vobornik, T. Pincelli, A. Troglia, T.-L. Lee, R. Ciancio, G. Drazic, A. Verdini, A. Regoutz, P. D. C. King, D. Biswas, G. Rossi, G. Panaccione and A. Selloni, Distinct behavior of localized and delocalized carriers in anatase TiO₂ (001) during reaction with O₂, Phys. Rev. Mater., 2020, 4, 025801.
94. P. Nandi, S. Shin, H. Park, Y. In, U. Amornkitbamrung, H. J. Jeong, S. J. Kwon and H. Shin, Large and Small Polarons in Highly Efficient and Stable Organic-Inorganic Lead Halide Perovskite Solar Cells: A Review, Sol. RRL, 2024, 8, 2400364.
95. P. Gono, J. Wiktor, F. Ambrosio and A. Pasquarello, Surface Polarons Reducing Overpotentials in the Oxygen Evolution Reaction, ACS Catal., 2018, 8, 5847–5851.
96. Y. Huang, Y. Zhu, H. Fu, M. Ou, C. Hu, S. Yu, Z. Hu, C.-T. Chen, G. Jiang, H. Gu, H. Lin, W. Luo and Y. Huang, Mg-Pillared LiCoO₂: Towards Stable Cycling at 4.6V, Angew. Chem., Int. Ed., 2021, 60, 4682–4688.
97. J. Xia, N. Zhang, D. Yi, F. Lu, Y. Yang, X. Wang and Y. Wang, Stabilizing 4.6 V LiCoO₂ via Er and Mg Trace Doping at Li-Site and Co-Site Respectively, Small, 2024, 20, 2311578.
98. M. Mladenov, R. Stoyanova, E. Zhecheva and S. Vassilev, Effect of Mg doping and MgO-surface modification on the cycling stability of LiCoO₂ electrodes, Electrochem. Commun., 2001, 3, 410–416.
99. S. Levasseur, M. Ménétrier and C. Delmas, On the Dual Effect of Mg Doping in LiCoO₂ and Li₁+δCoO₂: Structural, Electronic Properties, and 7Li MAS NMR Studies, Chem. Mater., 2002, 14, 3584–3590.
100. S. Shi, C. Ouyang, M. Lei and W. Tang, Effect of Mg-doping on the structural and electronic properties of LiCoO₂: A first-principles investigation, J. Power Sources, 2007, 171, 908–912.
101. X. G. Xu, C. Li, J. X. Li, U. Kolb, F. Wu and G. Chen, Electronic Structure of Li(Co, Mg)O₂ Studied by Electron Energy-Loss Spectrometry and First-Principles Calculation, J. Phys. Chem. B, 2003, 107, 11648–11651.
102. H. Tukamoto and A. R. West, Electronic Conductivity of LiCoO₂ and Its Enhancement by Magnesium Doping, J. Electrochem. Soc., 1997, 144, 3164.
103. W. Luo, X. Li and J. R. Dahn, Synthesis and Characterization of Mg Substituted LiCoO₂, J. Electrochem. Soc., 2010, 157, A782.
104. F. Zhou, M. Cococcioni, C. A. Marianetti, D. Morgan and G. Ceder, First-principles prediction of redox potentials in transition-metal compounds with LDA + U, Phys. Rev. B: Condens. Matter Mater. Phys., 2004, 70, 235121.
105. K. Hoang and M. D. Johannes, Defect chemistry in layered transition-metal oxides from screened hybrid density functional calculations, J. Mater. Chem. A, 2014, 2, 5224–5235.
106. L. Dahéron, H. Martinez, R. Dedryvère, I. Baraille, M. Ménétrier, C. Denage, C. Delmas and D. Gonbeau, Surface Properties of LiCoO₂ Investigated by XPS Analyses and Theoretical Calculations, J. Phys. Chem. C, 2009, 113, 5843–5852.
107. R. Sathiyamoorthi, P. Shakkthivel, R. Gangadharan and T. Vasudevan, Layered LiCo_1-xMg_xO₂ (x = 0.0, 0.1, 0.2, 0.3 and 0.5) cathode materials for lithium-ion rechargeable batteries, Ionics, 2007, 13, 25–33.
108. N. E. Kirchner-Hall, W. Zhao, Y. Xiong, I. Timrov and I. Dabo, Extensive benchmarking of DFT+U calculations for predicting band gaps, Appl. Sci., 2021, 11, 2395.
109. I. Timrov, F. Aquilante, M. Cococcioni and N. Marzari, Accurate Electronic Properties and Intercalation Voltages of Olivine-Type Li-Ion Cathode Materials from Extended Hubbard Functionals, Phys. Rev. X, 2022, 1, 033003.
110. Y.-C. Wang, Z.-H. Chen and H. Jiang, The local projection in the density functional theory plus U approach: A critical assessment, J. Chem. Phys., 2016, 144, 144106.
111. S. Bhowmik, A. J. Medford and P. Suryanarayana, Real-space Hubbard-corrected density functional theory, arXiv, 2025, preprint, arXiv:2507.23612,  DOI:10.48550/arXiv.2507.23612.

---