Pseudopotential and structural effects on electronic properties in transition metal oxides using Wannier functions

Abstract

How do pseudopotential choices affect the electronic transport properties in first-principles calculations? How does the atomic structure affect the computed properties? In this work, we present a comparison of six different pseudopotential libraries and their impact on computed material properties, including band gap, high-frequency permittivity, and electrical conductivity, using a set of transition metal oxides. We also highlight that Wannier functions can be used to calculate electronic transport properties of materials. Our workflow demonstrates that Wannier functions provide a reliable, accurate method to handle complex systems using any pseudopotential library. We validate and compare our results with existing theoretical works and experimental values from the literature. We find that the effect of various pseudopotentials on the calculated band gap, electrical conductivity, and high-frequency permittivity is negligible. Therefore, Wannier functions can be used to benchmark calculations of band gaps and electrical conductivities of metallic and semiconducting, as well as magnetic TMOs. We also find that the effect of the choice of the atomic structure (i.e., experimental vs. relaxed/optimized) on the electrical conductivity is not negligible. The experimental structure generally gives a more accurate result in comparison with experimental values of intrinsic conductivity than the relaxed atomic structure.

---

1 Introduction

Transition metal oxides (TMOs) have unique electronic properties and complex structures. In energy storage, TMOs serve as electrode materials, particularly as cathodes in lithium-ion batteries, due to their tunable conductivity.¹ In electronics, TMOs are active layer materials in sensors, transistors, and memory devices, due to the tunability of their properties and optical transparency.² In electrocatalysis³ and photocatalysis,⁴ TMOs are used as electrolytes and photocatalysts due to the strong covalent interaction between the transition metal ‘d’ orbitals and oxygen ‘p’ orbitals. This orbital interaction determines key catalytic properties, such as the reactivity of oxide surfaces and the adsorption of reactants, by facilitating charge transfer. While there are many practical applications of TMOs, ab initio modeling of electronic transport properties in TMOs remains challenging, due to a limited understanding of electronic behavior within these complex structures.⁵

Several recent theoretical studies have focused on calculating the electrical properties of TMOs. Borlido et al.⁶ reported on the effect of pseudopotentials (PPs) compared to the all-electron approach on band gap properties of solids. The all-electron approach was shown to agree very well with existing experimental band gap results. The absolute error in band gaps with respect to the experimental values can in some cases increase dramatically when PPs are used with a different exchange-correlation (XC) functional than that used for their construction. Kraus et al.⁷ showed that inconsistent PP – XC-functional combinations do not affect the prediction of the lattice constants in first-principles calculations. This tells us that the combination of PP and XC-functional might be critical, and this possibility will be assessed in the current work.

Petousis et al.⁸ reported that density functional perturbation theory (DFPT⁹) can be used to find the high-frequency permittivity of a wide range of inorganic materials (88 structures) with a band gap greater than 0.1 eV. They presented a workflow to calculate the electronic and ionic contributions to the permittivity of the materials. They stated that GGA is, in general, more accurate than LDA for finding dielectric permittivity.

Ricci et al.¹⁰ reported that electrical conductivity can be calculated from the band structure of the material using smoothed Fourier interpolation of the bands (as implemented in the BoltzTrap¹¹ code) for a wide range of inorganic materials. However, the direct Fourier interpolation method is computationally costly, as a dense k-point mesh has to be used to obtain convergence. Instead of the direct band interpolation method, an alternative method is the Wannier-based approach, which uses maximally localized Wannier functions (MLWFs)¹² to interpolate band structures (as implemented in the BoltzWann¹³ code). Since the BoltzTrap code uses direct Fourier interpolation of the bands to calculate electrical conductivity, it requires a denser k-point grid in comparison to the BoltzWann code. A comparison of the computational cost of BoltzTrap and BoltzWann, including different k-mesh density, is included in the SI, Section 2.3.4.

Neither of these high-throughput (HT) works^8,10 evaluated the impact of different pseudopotential libraries on the calculated electrical properties. Both of the works used relaxed structures instead of the available experimental data. Since one of the purposes of this work is to compare theoretical data with existing experimental values from literature, we will compare the effect of the choice of atomic structure on the materials' properties, such as electrical conductivity.

Finally, Youssef et al.¹⁴ also showed that Wannier functions can accurately capture field-dependent polarization in binary oxides. They used the Berry phase approach with MLWFs to show the effect of a high electric field on defect formation energy.

In the current work, we focus on three topics: (i) the use of Wannierization as a method for band structure interpolation as applied to the calculation of electronic properties of TMOs, (ii) the effect of atomic structure (experimental or relaxed) on the electronic properties, and (iii) an assessment of the impact of PP selection by comparing results obtained with the following six PP libraries commonly used in the literature:

• Optimized Norm Conserving Vanderbilt PPs (SG15-NC¹⁵),

• Ultra-Soft PPs from Rutgers (GBRV-US¹⁶),

• Projected Augmented Wave PPs of Dal Corso (PS-PAW¹⁷),

• Ultra-Soft PPs of Dal Corso (PS-US¹⁷),

• Standard Solid-State Pseudopotential (SSSP¹⁸).

• PseudoDojo PPs (Ps-Dojo) from van Setten¹⁹

In particular, we focus on the calculation of the electronic transport properties, including the electronic band gap, high-frequency dielectric permittivity, and electrical conductivity. Our approach systematically uses MLWFs for interpolation of electronic bands, even for complex systems such as TMOs.

The TMOs systems were selected because of their importance in applications such as energy storage, electronics, and catalysis, as well as to evaluate the methodology for different polymorphs of the same formula (TiO₂-R and TiO₂-A), various oxidation states (MoO₂ and MoO₃), effects of layered structures and tunable band gap (V₂O₅ and V₂O₃), effect of magnetism in Cr₂O₃, and ZnO, with an experimental band gap of ∼3.4 eV which is notoriously difficult to capture by traditional density functional theory (DFT). We validate our methodology using data from existing computational datasets, including data from the Materials Project (MP²⁰), and also using available experimental data. Finally, we present a study of the impact of PP selection on the electronic transport properties of selected TMOs, using the local density (LDA) as well as a generalised gradient (GGA) density functional approximation (DFAs). Additionally, we investigated the Hubbard + U correction with the PBE functional on electronic properties of the studied materials.

We find that our approach using MLWFs enables the calculation of high-frequency permittivity in TMOs, regardless of the choice of PPs. Another finding of this work is that the electrical conductivity of materials calculated using the BoltzWann code implemented with MLWFs is in good agreement with the results obtained from the BoltzTrap code implemented with direct interpolation by Ricci et al.¹⁰ Finally, we showed that the choice of PPs has a negligible impact on the calculated band gap, electrical conductivity, and high-frequency permittivity, while the atomic structure plays an important role in predicted electrical conductivity.

2 Methods

All DFT calculations were performed using Quantum ESPRESSO version 7.2 (QE)²¹ and Wannier90 version 3.1.0 (W90).²² All calculations were carried out using the Perdew–Zunger (PZ) LDA²³ or the Perdew–Burke–Ernzerhof (PBE) GGA,²⁴ and Hubbard corrections with +U values from AFLOW.²⁵ Spin-unpolarized calculations were used, except for the magnetic Cr₂O₃, where spin-polarization was allowed. BoltzWann version 3.1.0 (post-W90)¹³ was used to calculate the electrical conductivity of materials. The PBE + U_AFLOW calculations were not possible with the SG15-NC PP library, since the SG15-NC library does not have the required orbital information for the application of Hubbard corrections. However, the PS-Dojo library is included, representing norm-conserving PPs.

2.1 Electronic band localization

PZ and PBE, being LDA and GGA functionals, have some known limitations, including incorrect description of non-covalent interactions as well as poor electron localization due to the self-interaction error. Therefore, it is challenging to capture localized electron–electron interactions using PZ or PBE in complex systems such as TMOs, which may exhibit strong electronic correlation.⁵ One approach to overcome this limitation is by transforming the delocalized electronic states into a localized orbital representation and then computing the electronic properties of TMOs using such localized descriptions. The result of the localization procedure strongly depends on the quality of the input wavefunction, the choice of the bands that should be localized, as well as the method by which they are localized. One option is to minimize the spread of these orbitals, i.e., computing maximally localized Wannier functions (MLWFs).

A common way of constructing MLWFs is by projecting the input wavefunction onto a set of atomic or hydrogenic orbitals (e.g., s, p, d, f) or their common hybridization (e.g., sp, sp², sp³d, etc.). These then form an initial guess for the localization procedure. If suitable projections are chosen, the localization converges quickly. However, this approach requires manual specification of the projections. Generally, some “chemical intuition” is required to ensure the initial guess is appropriate. There are two common ways of automating the construction of MLWFs: the selected columns of the density matrix (SCDM)²⁶ method, and the projectability-disentangled Wannier functions (PDWFs)²⁷ method. In SCDM, an optimal set of basis vectors is generated directly from the density matrix. SCDM uses all selected bands (i.e., both valence and conduction bands) to construct MLWFs, so that the free electron contribution in the system can be screened. This is useful, especially for materials that have no band gap, i.e., materials with a metallic character. In PDWFs, for each Bloch state (i.e., the delocalized electronic wavefunctions from QE), PDWFs resemble a predefined set of physically intuitive local atomic orbitals, which are taken from the PPs used in the calculation. This is useful, especially for materials with a band gap. If the calculation is concentrated around the Fermi level, where the electronic transport happens, and if the initial guess for localized orbitals is chosen correctly, then the spread of the MLWFs will be smaller, and the calculation will be less costly. As long as the spread is small, the calculation gives a more accurate result, because Wannierization localizes well.

For the metallic systems, we employed the SCDM method in conjunction with automatic projection for the initial guesses for Wannier functions. This method enables the initialization of Wannierization without specifying frozen or disentanglement windows. We set sufficient iterations to have a localized spread of Wannier functions. The number of Wannier functions is selected based on the number of bands, keeping the metallic character of the system while generating localized orbitals to calculate the electronic conductivity of the system.

For insulators as well as magnetic cases (E_G > 1 eV), the PDWFs method was carried out manually for the selected materials, according to the workflow of Qiao et al.²⁷ Automated projection (e.g. using an AiiDA workflow) will be carried out in future work focusing on benchmarking of DFAs. In this method, several terms have to be defined, including disentanglement windows, the number of Wannierized bands, and orbital projections, to construct the initial projectors. The inner (frozen) window is set to 2 eV above the conduction band minimum to reproduce states around the Fermi energy. The outer (disentanglement) window is set to the maximum and minimum eigenvalues of the system. The number of Wannierized bands is the number of projectable atomic orbitals in the selection of the disentanglement window. The orbital projections can be selected by checking the density of states representation of the bands.

2.2 Choice of pseudopotentials

Pseudopotentials (PPs) are used to model the atomic nuclei and core electrons of the studied system. The use of PPs is a simplification that reduces computational demands while maintaining accuracy in predicting electronic structure.²⁸ PPs replace core states with an effective potential, so core states act like a point charge. This simplification makes the wavefunctions represented by a small number of plane waves, enabling efficient electronic-structure calculations.

Six different PP libraries were used in this work to investigate their effect on the electrical properties of the studied materials:

• The norm-conserving pseudopotentials (NC PPs, represented by SG15-NC) are constructed to explicitly treat valence states. The inside of the core states is replaced by using pseudo wavefunctions, while keeping the total charge (the norm) consistent with the all-electron wavefunctions. The outside of the core states in NC PPs is the same as in all-electron bases. This means the NC PPs can reproduce all-electron-like wavefunctions around the nucleus, achieving transferability, i.e., the potential is similarly suitable for different chemical environments. However, NC PPs need a higher number of plane waves and a higher cut-off energy; therefore, the computational time is not as fast as other PP libraries.

• The ultra-soft pseudopotentials (US PPs, represented by GBRV-US and PS-US) relax (or “soften”) the norm-conservation constraint of the NC PPs. As a consequence, they require a lower cut-off energy. They are designed to be faster but maintain accuracy by including augmented charges localized in the core region.

• The projector-augmented wave pseudopotentials (PAW PPs, represented by PS-PAW) are based on the all-electron PAW method, which explicitly represents the interactions between the nucleus, core, and valence electrons. The US PPs and PAW PPs can both provide better transferability at lower plane-wave cut-offs than NC PPs. PAW is a mixed approach, combining ideas from US PPs²⁹ and the linearized augmented-plane-wave method (LAPW³⁰). Notably, in calculations within the Materials Project (MP), PAW-based PPs are used as the standard PP library.³¹ There is no conceptual difference between PS-PAW, which is used in the current work, and the PAW PPs used in the MP calculations (carried out using VASP). The only practical difference is in the valence states: in the MP calculations used as a reference in this work, the _pv variants are used (for Ti, 10 explicit electrons in the 3p, 3d, and 4s shells), while the PS-PAW PPs are more analogous to the _sv variants available in VASP (for Ti, 12 explicit electrons, including also the semi-core 3s shell), see ref. 32.

• The SSSP library is a set of PPs optimized specifically for solid-state systems (i.e., bulk structures). This library enables the use of a high-performance PP for each element, based on a rigorous comparison of PPs of various construction. It includes PPs from the libraries mentioned earlier, and can be used for applications in high-throughput calculations, such as benchmarking. There are two variants of the library: precision and efficiency. The difference between these libraries comes from the benchmarking procedures, such as the equation of states (Δ-factor) or band structure deviations. The efficiency variant is more suitable for benchmarking, as it is less computationally costly than the precision variant, due to lower cut-off energies for most of the elements. In this work, we have used the SSSP-efficiency library version 1.3.0. For the materials studied, the difference between SSSP-efficiency and SSSP-precision is only in the value of the wavefunction and density cut-off.

• The PseudoDojo library (v0.5) is a set of NC PPs designed to provide high accuracy across a wide range of materials, including elemental solids and their oxides. Unlike SSSP, which is a composite library, PseudoDojo provides a consistent family of NC PPs generated using a unified methodology. The library is built using tight plane-wave cut-off energy tests, while reproducing formation energies across different chemical compounds.¹⁹ The PseudoDojo PP library was constructed by testing with both the PBE and the PZ DFAs.

2.3 Materials properties

2.3.1 Band gaps. The band gap (E_G) is the energy difference between the valence band maximum and the conduction band minimum. As we use MLWFs to calculate the electronic properties of TMOs, they must be able to reproduce the structures of both valence and conduction bands. Therefore, for materials with a lower band gap (<1 eV), the SCDM method was used to capture the delocalized electronic behavior in metallic systems. For materials with a higher band gap, we use fixed occupations to obtain the wavefunction, and the projection is performed using all atomic orbitals listed in the respective PP files.

The cut-off energies are taken directly from recommended values in the PP files. The highest cut-off energy among the elements in the TMO is used. For example, in the SSSP-efficiency library, the kinetic energy cut-off for oxygen wavefunctions is 50 Ry, and the kinetic energy cut-off for charge density is 400 Ry. The rest of the cut-off energies for each material can be found in the SI, Table S1. The k-point mesh for the system is determined by performing a convergence test for Brillouin zone sampling. The convergence criterion is based on the change of the total energy of the system with mesh density. When the total energy per atom is converged within 1 meV per atom between meshes, the smaller k-point mesh is selected. This ensures that the k-point selection will give reliable results while maintaining the speed of the calculation.

2.3.2 High-frequency permittivities. The dielectric permittivity is a frequency-dependent quantity that consists of electronic and ionic contributions. However, in this work, we focus only on the electronic contribution, i.e., the high-frequency part (ε^∞). The high-frequency permittivity can be calculated from the macroscopic polarization P of the material as a response to an external field (ε) or as the second derivative of the total energy (E) with respect to :³³where α and β represent the indices of the coordinates of the system, and δ_αβ is the Kronecker delta, which equals 1 when the indices are identical (α = β), otherwise 0.

The high-frequency dielectric permittivity can be computed by applying a homogeneous electric field within DFT calculations with periodic boundary conditions. The dipole moment can also be calculated using the modern theory of polarization,³⁴ where polarization is described as a Berry phase of the occupied electronic wavefunctions, see e.g. Spaldin.³⁵ One can derive eqn (1) using the change in the dipole moment Δp as a function of the change in the electric field strength: (ref. 36)

where Ω is the unit cell volume. Additionally, at the high frequency limit, the change in the total dipole moment Δp can be replaced by the electronic component, Δp_el.

An alternative approach based on MLWFs has been proposed by Vanderbilt and Resta.³⁴ Each MLWF represents a localized electron pair. Since Wannier functions are real-space representations of the constructed wavefunctions, this approach should be equivalent to the Berry phase formalism.³⁴ The sum over the Wannier centers of the occupied bands gives the electronic contribution to dipole moments,

where τ_i represents the position vector of the ith Wannier center, and e is the elementary charge. After finding Δp_el at two different field strengths ε, the high-frequency permittivity can be calculated using eqn (2).

2.3.3 Electrical conductivities. In solids, electrical conductivity arises from the transport of charge carriers such as electrons, holes, or ions. In TMOs, three common conduction mechanisms describing charge transport can be expected: intrinsic electronic conduction, a small or large polaron conduction mechanism, and variable-range hopping (VRH). Intrinsic conduction is governed by delocalized carriers moving in bands from the valence band to the conduction band at high temperatures.³⁷ VRH occurs by carrier transport between localized states in a disordered energy landscape at low temperatures, as described by Mott.³⁸ The polaron model occurs at mid-range temperatures through electron-phonon interaction via hopping. A polaron can be defined as a quasiparticle consisting of an electron that is polarized around the displaced ions, leading to self-trapping. The polaron model can be broadly categorized into two different mechanisms, namely small and large polarons. A small polaron forms when the short-range electron–lattice interaction is dominant, resulting in the self-trapped state localizing to a single site. An example of small polaron conduction is TiO₂-R.³⁹ In contrast, a large polaron forms when long-range electron–lattice interaction is dominant, resulting in the self-trapped state localizing to more than a single site.⁴⁰ In TiO₂-A, the conduction mechanism at low temperature is driven by large-polaron conduction.⁴¹ However, in this work, using the MLWFs approach within the constant relaxation time approximation, we focus only on the intrinsic electronic transport regime.

The Boltzmann transport equation can be used to calculate transport properties, describing how a particle moves in a solid.⁴² Using first-principles methods such as DFT, the electronic transport properties can be calculated from the electronic structure of a material using Boltzmann transport theory. The derivation from the Boltzmann equation into the first-principles approach is based on parametrized band structure (e.g., WFs) and scattering mechanism (e.g., relaxation time). For a detailed explanation, see Scheidemantel et al.⁴³ This can be done using the direct band interpolation method, as implemented in the BoltzTrap code, which was used by Ricci et al..¹⁰ However, this technique requires the use of a highly dense k-mesh, which is costly. An alternative approach to the direct interpolation method is to use Wannier functions, as implemented in the BoltzWann code. This approach uses MLWFs to interpolate the band structure of the material and, with a scattering mechanism, to calculate the electrical conductivity. It has the advantage of using a less dense k-mesh while maintaining the accuracy of the calculation in comparison to the direct interpolation method. For a computational time and convergence comparison between the two approaches, see SI, Section 2.3.4.

In this work, the BoltzWann code using MLWFs is used to calculate electrical conductivity properties of selected materials. In our scope of work, the thermal conductivity is not investigated. The indirect interpolation using MLWFs, as implemented in BoltzWann, is as follows:

where µ represents the chemical potential (which is set to be around Fermi energy), α and β are the Cartesian indices, f(E, µ, T) is the Fermi–Dirac distribution function, and Σ_αβ(E) is the transport distribution function (TDF). The Fermi–Dirac function represents the number of charge carriers, while the TDF represents their mobility. As such, the TDF includes an explicit sum of the group band velocities v_α,β(n, k) over all quantum numbers of the system, the energy of the bands E_n,k and the relaxation time τ_nk, as shown in eqn (5).

In eqn (5), δ is the Dirac delta (δ = 1 when E = E_n,k, otherwise 0). τ_nk is, on the other hand, a critical parameter, affecting Σ_αβ(E) and therefore σ_αβ(µ, T) with a linear dependency. This parameter has to be specified in the input file before the calculation. The effects of the relaxation time on the electron-phonon interaction were investigated for semiconductors and compared to experimental results for anisotropic bulk structures by Allen et al.⁴⁴ In previous works,¹⁰ a constant relaxation time of 10 fs was used. Here, we used the same relaxation time for the sake of comparability. An investigation of the effect of the relaxation time on selected systems will be carried out in future work.

The electrical conductivity σ_αβ(µ, T) depends on the chemical potential and temperature. In a first-principles calculation, the chemical potential must be specified. Its value varies with doping or temperature. To obtain the conductivity at a given temperature, we set the chemical potential to a range that covers both the valence band maximum (VBM) and conduction band minimum (CBM). For metals, the temperature-dependent Fermi level values are found using Fermi–Dirac smearing via QE for each temperature. For insulators, the Fermi level is approximately in the middle of the gap and shifts with increasing temperature (see ref. 45, Section 28). In the current work, no doped systems were considered.

2.4 Method validation

We first validated our method before the calculations were done, see SI, Section 1. The band structures obtained after Wannierization versus the bands obtained directly from QE are shown in Fig. S1–S8, showing good agreement. The high-frequency permittivity values obtained using the Wannier charge center method were compared with the Berry phase method, see Fig. S9–S11, showing that our result closely follows the reference data. We validated our workflow for electrical conductivity calculation by comparing with results from MP, see Fig. S12–S14, showing a good agreement given the difference in methodology.

After successfully validating our MLWF-based framework, we can now confidently calculate the material properties and evaluate the effect of different PP libraries and structural effects. We note that in this work, PZ, PBE, and PBE + U_AFLOW DFAs are used for all calculations; benchmarking of other DFAs will be carried out in further work.

3 Results

As discussed in the introduction, the effect of the inconsistency between the DFA used to parametrize the PPs and that used for the calculations is to be assessed. This comparison aims to show whether an inconsistency between PP parametrization and the XC functional used can affect materials properties such as band gap and conductivity. In the inconsistency calculations, we use the PBE DFA in combination with various pslibrary (PS) PPs, which were parametrized using the BP, PBEsol, and PBE functionals, and are provided in both PAW and ultrasoft (US) forms. The results for ZnO and V₂O₅ are shown in Section 2.1 of the SI. For the band structures of the two materials (Fig. S15–S18), the use of inconsistent PPs can lead to ∼0.1 eV difference in the band gaps, consistent with previous work.⁸ In the extreme case of the already underpredicted band gap of ZnO, the use of the PS-US PPs with PZ-parametrization can cause a 0.3 eV difference (i.e., a 40% deviation). For the V₂O₅ case, the effect is reduced (to below 7%) as the band gap is larger at ∼1.6 eV. The results are analogous for the electrical conductivity (Fig. S19–S22). There is no significant deviation for the V₂O₅ case, with the activation energies around 0.01 eV across functionals. On the other hand, there is a non-negligible effect for ZnO, especially when the PZ-parametization is used (a 0.06 eV or a 30% deviation in the activation energy from PBE and BP results). Our current results are therefore consistent with the work of Petousis et al.⁸ the effect of the XC functional – PP parametization inconsistency cannot be neglected for materials with a band gap below ∼1 eV; in most other cases the effect is negligible.

The effect of the atomic structure on the computed material properties, such as electrical conductivity, also has to be assessed. This was investigated for insulating materials, including the magnetic Cr₂O₃, using different functionals (PZ, PBE, PBE + U_AFLOW) using Ps-Dojo PPs by comparing relaxed/optimized computational structure from MP (solid lines) versus the experimental structure from the Crystallography Open Database (dashed lines), see the Fig. S23–S28. The experimental structure provides a better agreement with the literature values of activation energy of intrinsic electrical conductivity for all materials except ZnO, discussed below in Section 3.3. In the case of MoO₃, shown in Fig. S24, the improved agreement when the experimental structure is used is quite striking.

3.1 Band gaps

Fig. 1 shows the effect of various PPs on the computed band gaps. The band gap values are shown only for the semiconducting materials in our set, since metals do not have any band gap near the Fermi level. The available experimental values of E_G are indicated by the dashed lines in Fig. 1. For all calculations, the computed band gaps are underestimated compared to the experimental values. This underestimation is expected since the PZ and PBE DFAs are well-known for predicting lower band gaps than experimental values; as expected, the use of AFLOW + U values improves the agreement with experiment somewhat.

Fig. 1 Band gap comparison for Cr₂O₃ across different PPs with various functionals: (a) PZ, (b) PBE, (c) PBE + U_AFLOW. Experimental data (dashed line) from Abdullah et al.⁴⁶

The Fig. 1 shows the typical, systematic dependence of E_G of Cr₂O₃ on the DFA, while PP dependence remains relatively minor. For PZ (Fig. 1(a)), the band gaps are underestimated across all PPs (∼0.2–0.4 eV), which is expected due to the well-known over-delocalization error of LDA. Changing the functional to PBE (Fig. 1(b)) provides the typical improvement of GGA over LDA across all PPs (∼0.8–1.2 eV), but still far from the experimental value (∼3.2 eV (ref. 46)). In contrast to PBE, for Cr₂O₃, the Hubbard correction using AFLOW + U values significantly improves the band gap prediction with respect to the experimental value. Since Cr₂O₃ is an antiferromagnetic oxide, all calculations are performed within a spin-polarized framework. Cr₂O₃ serves as an important stress test for our MLWFs-based framework, as it extends the material scope to include magnetic cases where spin polarization and localized d states play a crucial role. The results for the rest of the materials can be found in the SI, Fig. S29–S34.

The influence of different PP libraries on the band gap is minor. For example, MoO₃ (see Fig. S29) shows nearly the same values of E_G among all PPs for PZ, PBE, and PBE + U_AFLOW. The largest relative deviation is obtained for the ZnO case when PZ is used (see Fig. S33(a)), with around 40% deviation between PS-PAW and SG15-NC or PS-US, as the predicted band gap is very low (0.515 eV on average). With PBE, the deviation between PPs drops to 14% as the E_G values increase. Note that the Hubbard correction improves the band gap prediction in ZnO only slightly, to 1.25 eV on average (i.e., still a 2 eV underprediction), and decreases the deviation between PPs to 7%. Finally, despite the PS-PAW and PS-US libraries being generated from the same atomic data,¹⁷ the band gaps obtained using these two libraries can differ. No consistent trend can be observed between the different PPs in this set of materials. However, the band gap values are largely robust concerning the choice of PPs for the selected materials.

3.2 High-frequency permittivity

The structural effect on predicting high-frequency permittivity was investigated in SI, Section 1.3, for different atomic structures, such as conventional and primitive non-relaxed cells, as well as primitive relaxed cells. Fig. S10 and S11 show that perfect agreement between the Wannier function and Berry phase methods can be obtained when one uses a relaxed primitive cell to predict high-frequency permittivity. The reason for that can be explained by poorer convergence of Wannierization at structures further from the energy minimum, see Tables S4 and S5.

In the following, we use primitive cells obtained from MP, i.e., not relaxed cells, which leads to a discrepancy between the Wannier function and Berry phase results. The aim of this study is not to optimize structural parameters but to compare the performance of the different PP libraries and/or functionals under identical, experimentally relevant conditions. Additionally, a geometry optimization using a hybrid functional may be too costly for routine work.

Fig. 2 shows the variation of ε^∞ for the two most interesting cases, TiO₂-R and V₂O₅; the results for the rest of the materials can be found in the SI, Fig. S34–S38. The colors indicate the direction of the applied electric field as well as the lattice vector (a, b, c). As mentioned in SI, Section S1.3, the electric field is always applied parallel to the studied axis, i.e., the Kronecker delta δ_αβ is always 1. All the calculations were done by using a finite external field ( = 0.001 Ry a.u).

Fig. 2 High-frequency permittivities (ε^∞) computed from MLWF centers (solid bars) and the Berry phase (hatched bars) for TiO₂-R (left) and V₂O₅ (right). Six pseudopotential libraries are compared. The squares above the PS-PAW data are taken from ref. 8. Colors represent the lattice vectors of the materials: a (blue), b (orange), and c (green).

For TiO₂-R, with a tetragonal lattice, all of the computed ε^∞ values along the a- and b-axis are identical due to the isotropic ab plane, while the c-axis values are different, but only slightly (see Fig. 2 (left)). The effect of PPs on Berry phase results is negligible, with the a and b components in the 6.95–7.10 range and the c component in the 8.56–8.59 range across all six PPs. For this nearly isotropic case, the differences between the MLWF-based results (solid bars) and Berry phase data (hatched bars) are minor. The agreement between the two methods is within 1%, which shows that MLWFs can reproduce Berry phase polarization very reliably. With both methods, SG15-NC gives slightly smaller values for both lattice directions by 0.1 compared to the rest.

With both methods, the calculated permittivities are lower than the MP reference values, which are 7.55 for the a and b components, and 9.20 for the c component. The reason for the deviation may be that the reference values are calculated using a relaxed primitive cell as well as a k-point mesh that is almost five times denser than in our work (ρ_k = 4.13a₀⁻³ for MP, 0.83a₀⁻³ in our work).

By contrast, V₂O₅ has an orthorhombic structure with three different axes and a significant anisotropy. As shown in Fig. 2 (right), with the MLWF approach, the a-axis permittivity ranges from about 6.37 (SG15-NC) to 6.60 (PS-PAW), the b-axis changes between 6.45 (SG15-NC) and 6.66 (GBRV-US), and the c-axis is between 3.8 and 4.0 across all PPs. The Berry phase calculations yield lower values, narrower spreads between PPs, but the same trends as the MLWF results.

The largest deviation for V₂O₅ between the MLWF and the Berry phase methods appears along the a-axis, with around a 10% difference. The deviation is moderate in the b-axis (3%) and negligible in the c-axis. The structure used for V₂O₅ in our work is a primitive cell based on the experimental structure, with a volume of 175.71 ų. By contrast, the relaxed cell used to obtain the reference MP values (a = 5.36, b = 5.90, c = 3.40) has a volume of 196.78 ų. This means the experimental structure is far from the PBE equilibrium structure, likely due to van der Waals effects. If we used the relaxed cell from the MP calculation, there would be no deviation between the reference data, the MLWF, and the Berry phase results, see Fig. S11 in the SI. It is worth noting that the MP values are calculated with almost five times denser k-point mesh density than our current work (3.41 a₀⁻³ and 0.71 a₀⁻³, respectively).

In conclusion, the Berry phase values for high-frequency dielectric permittivity can be reproduced using the Wannier function approach. Importantly, the two approaches differ in the representation of the electronic polarization: the Berry phase approach evaluates the phase of occupied states in reciprocal space, while the Wannier function method expresses the polarization using localized orbitals as charge centers in real space. Nevertheless, they give similar results for high-frequency dielectric permittivity using experimental structures, and identical results with relaxed primitive cells. The other main finding is that the choice of the PPs does not play an important role in high-frequency dielectric permittivity calculations.

3.3 Electrical conductivity

After validation of the MLWF-based band structures and obtaining a good agreement with the available literature data (see SI, Section S1.4),¹⁰ we calculated the electrical conductivity for the whole dataset, including materials with a metallic character (MoO₂ and V₂O₃) and the magnetic Cr₂O₃ treated with spin-polarization. The temperature ranges used in the calculations are determined from the available experimental studies, since the purpose of this work is to compare with experimental values. In all cases, the experimental data were obtained from the existing experimental literature,^47–54 focusing on the high temperature range, where the intrinsic electronic conduction mechanism ought to be dominant. Among of the reference literature works, TiO₂-R, MoO₃, and MoO₂ are from single crystal samples, while the rest are polycrystalline (ZnO), annealed films (V₂O₅), nanoporous (TiO₂-A), or powder (Cr₂O₃) samples. Note that our calculations of electrical conductivity are performed using constant relaxation time τ_nk = 10 fs.

The results of the conductivity for molybdenum oxides (MoO₂ and MoO₃) are shown in Fig. 3, while the results using PZ and PBE + U_AFLOW can be found in the SI, Fig. S40 and S41. The effect of the PPs library on the conductivity prediction of both molybdenum oxides is negligible, see Fig. S40 and S41. This means that our workflow can predict the electrical conductivity for insulators (such as MoO₃) and for metallic cases (such as MoO₂). A notable difference, however, arises from the choice of atomic structure, as shown in Fig. 3 (left). For MoO₃, the optimized structure significantly overestimates the activation energy (∼1.20 eV) compared to the experimental value (∼0.60 eV), while the experimental structure (dashed lines) estimates more accurately (∼0.42 eV). On the other hand, MoO₂ shows almost perfect agreement with experimental values, regardless of the structure used, even though the temperature is extremely low, see Fig. 3 (right). As aforementioned, metallic materials do not have a band gap around the Fermi level; therefore, their electrons are free to move from the valence to the conduction band without thermal activation. Thus, the activation energy for metals is expected to be around zero. Note that in both cases, the experimental conductivity values are obtained using single crystals.

Fig. 3 Electrical conductivity (σ) with PBE and PS-Dojo PPs as a function of temperature for MoO₃ (left) and MoO₂ (right). Calculations using optimized atomic structure (solid lines) compared to experimental atomic structure (dashed lines) and experimental data (black dots) from ref. 47 for MoO₃, and ref. 53 for MoO₂.

The results of the conductivity for TiO₂-A and TiO₂-R are presented in Fig. 4. The effect of PPs on the electrical conductivity of TiO₂-A and TiO₂-R is negligible for the three DFAs used, see Fig. S42 and S43. TiO₂-A is known to exhibit a large polaron formation at low temperature, see ref. 41, 56 and 57. The phase diagram of polarons in TiO₂-A clearly shows that above 240 K the intrinsic conduction mechanism starts dominating over the polaron mechanism of conductivity, see Fig. 9 in ref. 41. Therefore, we calculated the electrical conductivity of TiO₂-A above this threshold, in the range from 400 K to 800 K, where experimental data is also available.^51,55 Fig. 4 compares the calculated and experimental datasets, showing different activation energies. This indicates a strong sensitivity to sample-dependent effects such as microstructure, defect concentration, and experimental conditions. However, the deviation between the two experimental datasets (∼0.12 eV) is much lower than the deviation between the optimized atomic structure and experimental atomic structure (∼0.30 eV).

Fig. 4 Electrical conductivity (σ) with PBE and PS-Dojo PPs as a function of temperature for TiO₂-A (left) and TiO₂-R (right). Legend as in Fig. 3, with experimental data from Ref. 55 and 51 for TiO₂-A, and ref. 50 for TiO₂-R.

On the other hand, TiO₂-R is known to exhibit a small polaron formation at low temperature.³⁹ Above 400 K, the conduction mechanism is governed by band transport.³⁹ We calculated the electrical conductivity values in the range of 850 K to 1250 K, as shown in Fig. 4 (right). The effect of the choice of structure is also here not negligible. The results from the experimental atomic structure calculations give more accurate conductivity behavior in this temperature range (∼1.44 eV) when compared to the experimental value (∼1.48 eV) obtained from single crystals, than when using the optimized structure (∼1.06 eV). This is similar to the MoO₃ case. Finally, an improvement in the band gap prediction does not correlate with an improvement in the conductivity prediction for either of the titanium oxides, see Fig. S42 and S43 in the SI.

The results of the conductivity for ZnO and Cr₂O₃ are shown in Fig. 5. For ZnO, despite the significantly underestimated band gap with PBE (see Fig. S33), the conductivity is captured well, with reasonable agreement in the absolute value of σ (within an order of magnitude), and excellent agreement in the activation energy of intrinsic conductivity E_a(σ). The conductivity does not depend only on band gap energies (i.e., the number of charge carriers), but also on band velocities (i.e., charge carrier mobility), see eqn (5). Given the underprediction of the ZnO band gap and associated overprediction in the number of carriers, the charge carrier mobility here must be strongly underpredicted, leading to conductivity values consistent with experiment due to fortuitous error cancellation. This is confirmed by the use of PBE + U_AFLOW, see Fig S46, which improves the band gap energy, but reduces the agreement in E_a(σ). Additionally, ZnO exhibits an opposite sensitivity to the choice of the atomic structure compared to the other materials, as discussed in Section 3. When the experimental atomic structure is used for ZnO, a much worse agreement with experimental conductivity values is obtained, again due to the breaking of the above error cancellation. ZnO is different from the rest of the materials set, because it has a fully occupied 3d closed shell. Normally, TMOs have partially filled d-states near the Fermi level. The interaction between O 2p and Zn 3d is strong and sensitive to the bond length, which is a structural parameter.⁵⁸ Therefore, small changes in the Zn–O distance can affect p–d repulsion. We can conclude that ZnO is very sensitive to the choice of PP, DFA, inclusion of Hubbard corrections, and the atomic structure. On the other hand, the conductivity results of Cr₂O₃ shown in Fig. 5 (right) confirm that our methodology can calculate the temperature dependence of the conductivity values using spin polarization, with a reasonable agreement obtained as long as the experimental atomic structure is used. Even here is the effect of the PP choice negligible, see Fig. S47.

Fig. 5 Electrical conductivity (σ) with PBE and PS-Dojo PPs as a function of temperature for ZnO (left) and Cr₂O₃ (right). Legend as in Fig. 3, with experimental data from ref. 48 for ZnO, and ref. 54 for Cr₂O₃.

The results of the conductivity for V₂O₅ and V₂O₃ are presented in Fig. 6. The effect of PPs on the electrical conductivity of V₂O₅ and V₂O₃ is negligible, see Fig. S44 and S45. One can easily recognize in Fig. 6 that the prediction of the conductivity for V₂O₅ is in poor agreement with the existing experimental literature values. This can be explained by any of the following reasons: first, the temperature range studied may still be in the polaron regime, see Fig. 5 from Wernbacher et al.⁵⁹ Therefore, the polaron hopping mechanism is contributing to the conduction. Since our model can only predict intrinsic conduction, the prediction does not show good agreement. Second, the experimental conductivity data⁵² is from an annealed film deposited on glass substrates, which is not a single crystal sample. Furthermore, V₂O₅ is known to form oxygen vacancies readily, even more so than titanium oxides.⁶⁰ The effects of defects are also not modelled using our methodology. Finally, V₂O₅ is well-known as a layered structure, therefore the use of a dispersion-corrected DFA may be necessary to capture these effects. It is worth remembering that materials that have the highest anisotropy in our dataset also showed the highest deviation between the Berry phase approach and the MLWF approach in the high-frequency dielectric permittivity when using the experimental structures, see Fig. 2 and S38. However, the calculated activation energy of conductivity is around half of the experimental band gap,⁶¹ which corroborates that our methodology predicts a reasonable intrinsic conductivity. V₂O₅ is therefore a good candidate to investigate further, considering electron-phonon coupling effects. On the other hand, the results for the metallic V₂O₃ are shown in Fig. 6 (right). The absolute values of σ are predicted to within a factor of 1.5 and the choice of atomic structure plays no effect; the trends are therefore similar to the metallic MoO₂ (see Fig. 3 (right)), as discussed above.

Fig. 6 Electrical conductivity (σ) with PBE and PS-Dojo PPs as a function of temperature for V₂O₅ (left) and V₂O₃ (right). Legend as in Fig. 3, with experimental data from ref. 52 for V₂O₅, and ref. 49 for V₂O₃.

The choice of PP library therefore makes no significant difference in σ or E_a(σ). However, in general, the inclusion of a Hubbard correction using the tabulated AFLOW + U values also does not lead to a systematic improvement in the predicted E_a(σ). This is a common trend for all materials in the dataset. The correlation between band-gap accuracy and intrinsic electrical conductivity within the constant-relaxation-time framework is weak. On the other hand, for semiconducting materials where experimental data from single crystals is available (MoO₃, TiO₂-R), the use of an experimental atomic structure is necessary to obtain a good agreement with the activation energy of the intrinsic conductivity. Accurate treatment of materials such as V₂O₅, ZnO, and Cr₂O₃ may therefore require the use of an experimental structure with a more appropriate DFT approach, or a treatment beyond the constant relaxation time τ_nk, in addition to comparing with single-crystal experimental data.

4 Conclusion

In this work, we tested the effect of the PPs and atomic structure on the electronic structure and electronic transport properties of materials using PZ, PBE, and PBE + U DFAs, including the band gap, high-frequency dielectric permittivity, and electrical conductivity. We have also assessed the feasibility of using maximally localized Wannier functions to calculate these properties. We found that the effect of PPs on the results of band gap, electrical conductivity calculations, and high-frequency dielectric permittivity is generally negligible. However, for perfect consistency between the MLWF and Berry phase approaches, the MLWFs have to be well converged, the occupied bands have to be properly disentangled, the primitive cells have to be used, and the structures should be relaxed. Moreover, we were able to calculate the temperature response of the conductivity of TMOs with both insulating and metallic character, as well as magnetic cases, allowing for material screening using a Wannier function-based workflow. We showed that the choice of atomic structure has a non-negligible impact on the electrical conductivity by comparing results computed with optimized/relaxed and experimental structures. Finally, we investigated the effect of the inconsistency between the PP parametrization and the XC-functional used, with significant impact on materials with a low calculated band gap (<1 eV), e.g., ZnO. Our findings pave the road for benchmarking of density functional approximations on a wider range of materials in our future work.

Conflicts of interest

There are no conflicts to declare.

Data availability

Additional supplemental information including all QE and W90 input and output files as well as postprocessing scripts is available at DOI: https://doi.org/10.5281/zenodo.15630691.

Supplementary information (SI): methodological validation and additional results for the effect of PP-XC inconsistency, atomic structure, and PP choices on the material properties such as band gap, high-frequency permittivity, and electrical conductivity as well as the computational cost comparison. See DOI: https://doi.org/10.1039/d5ra09637e.

Acknowledgements

This work was funded by the Deutsche Forschungsgemeinschaft (DFG) under project number 490703766. The calculations were performed using the computational resources of the Institute of Mathematics at Technische Universität Berlin. The authors would like to thank G. Pizzi and No. Paulish from PSI for valuable discussions and support.

References

1. P. Mallick, S. Moharana, L. Biswal and S. K. Satpathy, in Emerging Nanodielectric Materials for Energy Storage, ed. S. Moharana, D. H. Gregory and R. N. Mahaling, Springer International Publishing, Cham, 2024, pp. 331–355, DOI:  DOI:10.1007/978-3-031-40938-7_12.
2. Y. Jeon, D. Lee and H. Yoo, Recent Advances in Metal-Oxide Thin-Film Transistors: Flexible/Stretchable Devices, Integrated Circuits, Biosensors, and Neuromorphic Applications, Coatings, 2022, 12, 204,  DOI:10.3390/coatings12020204.
3. Y. Dang, T. Wu, H. Tan, J. Wang, C. Cui, P. Kerns, W. Zhao, L. Posada, L. Wen and S. L. Suib, Partially reduced Ru/RuO₂ composites as efficient and pH-universal electrocatalysts for hydrogen evolution, Energy Environ. Sci., 2021, 14, 5433–5443,  10.1039/D1EE02380B.
4. S. Lin, H. Huang, T. Ma and Y. Zhang, Photocatalytic Oxygen Evolution from Water Splitting, Adv. Sci., 2021, 8, 2002458,  DOI:10.1002/advs.202002458.
5. C. J. Cramer and D. G. Truhlar, Density functional theory for transition metals and transition metal chemistry, Phys. Chem. Chem. Phys., 2009, 11, 10757,  10.1039/b907148b.
6. P. Borlido, J. Doumont, F. Tran, M. A. L. Marques and S. Botti, Validation of Pseudopotential Calculations for the Electronic Band Gap of Solids, J. Chem. Theory Comput., 2020, 16, 3620–3627,  DOI:10.1021/acs.jctc.0c00214.
7. P. Kraus, P. Raiteri and J. D. Gale, Computational workflows for perovskites: case study for lanthanide manganites, Phys. Chem. Chem. Phys., 2023, 25, 14799–14811,  10.1039/D3CP00041A.
8. I. Petousis, W. Chen, G. Hautier, T. Graf, T. D. Schladt, K. A. Persson and F. B. Prinz, Benchmarking density functional perturbation theory to enable high-throughput screening of materials for dielectric constant and refractive index, Phys. Rev. B, 2016, 93, 115151,  DOI:10.1103/PhysRevB.93.115151.
9. S. Baroni, S. De Gironcoli, A. Dal Corso and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys., 2001, 73, 515–562,  DOI:10.1103/RevModPhys.73.515.
10. F. Ricci, W. Chen, U. Aydemir, G. J. Snyder, G.-M. Rignanese, A. Jain and G. Hautier, An ab initio electronic transport database for inorganic materials, Sci. Data, 2017, 4, 170085,  DOI:10.1038/sdata.2017.85.
11. G. K. Madsen and D. J. Singh, BoltzTraP. A code for calculating band-structure dependent quantities, Comput. Phys. Commun., 2006, 175, 67–71,  DOI:10.1016/j.cpc.2006.03.007.
12. N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza and D. Vanderbilt, Maximally localized Wannier functions: Theory and applications, Rev. Mod. Phys., 2012, 84, 1419–1475,  DOI:10.1103/RevModPhys.84.1419.
13. G. Pizzi, D. Volja, B. Kozinsky, M. Fornari and N. Marzari, BoltzWann: a code for the evaluation of thermoelectric and electronic transport properties with a maximally-localized Wannier functions basis, Comput. Phys. Commun., 2014, 185, 422–429,  DOI:10.1016/j.cpc.2013.09.015.
14. M. Youssef, K. J. Van Vliet and B. Yildiz, Polarizing Oxygen Vacancies in Insulating Metal Oxides under a High Electric Field, Phys. Rev. Lett., 2017, 119, 126002,  DOI:10.1103/PhysRevLett.119.126002.
15. M. Schlipf and F. Gygi, Optimization algorithm for the generation of ONCV pseudopotentials, Comput. Phys. Commun., 2015, 196, 36–44,  DOI:10.1016/j.cpc.2015.05.011.
16. K. F. Garrity, J. W. Bennett, K. M. Rabe and D. Vanderbilt, Pseudopotentials for high-throughput DFT calculations, Comput. Mater. Sci., 2014, 81, 446–452,  DOI:10.1016/j.commatsci.2013.08.053.
17. A. Dal Corso, Pseudopotentials periodic table: from H to Pu, Comput. Mater. Sci., 2014, 95, 337–350,  DOI:10.1016/j.commatsci.2014.07.043.
18. G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet and N. Marzari, Precision and efficiency in solid-state pseudopotential calculations, npj Comput. Mater., 2018, 4, 72,  DOI:10.1038/s41524-018-0127-2.
19. M. Van Setten, M. Giantomassi, E. Bousquet, M. Verstraete, D. Hamann, X. Gonze and G.-M. Rignanese, The PseudoDojo: Training and grading a 85 element optimized norm-conserving pseudopotential table, Comput. Phys. Commun., 2018, 226, 39–54,  DOI:10.1016/j.cpc.2018.01.012.
20. A. Jain, S. P. Ong, G. Hautier, W. Chen, W. D. Richards, S. Dacek, S. Cholia, D. Gunter, D. Skinner, G. Ceder and K. A. Persson, Commentary: The Materials Project: A materials genome approach to accelerating materials innovation, APL Mater., 2013, 1, 011002,  DOI:10.1063/1.4812323.
21. P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. De Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari and R. M. Wentzcovitch, QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials, J. Phys.: Condens. Matter, 2009, 21, 395502,  DOI:10.1088/0953-8984/21/39/395502.
22. G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi and J. R. Yates, Wannier90 as a community code: new features and applications, J. Phys.: Condens. Matter, 2020, 32, 165902,  DOI:10.1088/1361-648X/ab51ff.
23. J. P. Perdew and A. Zunger, Self-interaction correction to density-functional approximations for many-electron systems, Phys. Rev. B, 1981, 23, 5048–5079,  DOI:10.1103/PhysRevB.23.5048.
24. J. P. Perdew, K. Burke and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett., 1996, 77, 3865–3868,  DOI:10.1103/PhysRevLett.77.3865.
25. A. Jain, G. Hautier, C. J. Moore, S. Ping Ong, C. C. Fischer, T. Mueller, K. A. Persson and G. Ceder, A high-throughput infrastructure for density functional theory calculations, Comput. Mater. Sci., 2011, 50, 2295–2310,  DOI:10.1016/j.commatsci.2011.02.023.
26. A. Damle and L. Lin, Disentanglement via Entanglement: A Unified Method for Wannier Localization, Multiscale Model. Simul., 2018, 16, 1392–1410,  DOI:10.1137/17M1129696.
27. J. Qiao, G. Pizzi and N. Marzari, Projectability disentanglement for accurate and automated electronic-structure Hamiltonians, npj Comput. Mater., 2023, 9, 208,  DOI:10.1038/s41524-023-01146-w.
28. P. Schwerdtfeger, The Pseudopotential Approximation in Electronic Structure Theory, ChemPhysChem, 2011, 12, 3143–3155,  DOI:10.1002/cphc.201100387.
29. D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue formalism, Phys. Rev. B, 1990, 41, 7892–7895,  DOI:10.1103/PhysRevB.41.7892.
30. O. K. Andersen, Linear methods in band theory, Phys. Rev. B, 1975, 12, 3060–3083,  DOI:10.1103/PhysRevB.12.3060.
31. G. Kresse and D. Joubert, From ultrasoft pseudopotentials to the projector augmented-wave method, Phys. Rev. B, 1999, 59, 1758–1775,  DOI:10.1103/PhysRevB.59.1758.
32. M.-T. Huebsch and M. Wolloch, VASP-WiKi, English, Wikipedia, 2025.
33. X. Wu, D. Vanderbilt and D. R. Hamann, Systematic treatment of displacements, strains, and electric fields in density-functional perturbation theory, Phys. Rev. B, 2005, 72, 035105,  DOI:10.1103/PhysRevB.72.035105.
34. R. Resta and D. Vanderbilt, in Physics of Ferroelectrics, Springer Berlin Heidelberg, Berlin, Heidelberg, 2007, vol. 105, pp. 31–68,  DOI:10.1007/978-3-540-34591-6_2.
35. N. A. Spaldin, A beginner's guide to the modern theory of polarization, J. Solid State Chem., 2012, 195, 2–10,  DOI:10.1016/j.jssc.2012.05.010.
36. P. Umari and A. Pasquarello, Ab initio Molecular Dynamics in a Finite Homogeneous Electric Field, Phys. Rev. Lett., 2002, 89, 157602,  DOI:10.1103/PhysRevLett.89.157602.
37. M. G. Hutchins, O. Abu-Alkhair, M. M. El-Nahass and K. Abdel-Hady, Electrical conduction mechanisms in thermally evaporated tungsten trioxide (WO₃) thin films, J. Phys.: Condens. Matter, 2006, 18, 9987–9997,  DOI:10.1088/0953-8984/18/44/001.
38. N. Mott, Conduction in glasses containing transition metal ions, J. Non-Cryst. Solids, 1968, 1, 1–17,  DOI:10.1016/0022-3093(68)90002-1.
39. V. N. Bogomolov and D. N. Mirlin, Optical Absorption by Polarons in Rutile (TiO₂) Single Crystals, Phys. Status Solidi B, 1968, 27, 443–453,  DOI:10.1002/pssb.19680270144.
40. D. Emin, Optical properties of large and small polarons and bipolarons, Phys. Rev. B, 1993, 48, 13691–13702,  DOI:10.1103/PhysRevB.48.13691.
41. C. Franchini, M. Reticcioli, M. Setvin and U. Diebold, Polarons in materials, Nat. Rev. Mater., 2021, 6, 560–586,  DOI:10.1038/s41578-021-00289-w.
42. B. Nag, in Electron Transport in Compound Semiconductors, ed. B. Nag, Springer Berlin Heidelberg, Berlin, Heidelberg, 1980, pp. 129–170,  DOI:10.1007/978-3-642-81416-7_7.
43. T. J. Scheidemantel, C. Ambrosch-Draxl, T. Thonhauser, J. V. Badding and J. O. Sofo, Transport coefficients from first-principles calculations, Phys. Rev. B, 2003, 68, 125210,  DOI:10.1103/PhysRevB.68.125210.
44. P. B. Allen, W. E. Pickett and H. Krakauer, Anisotropic normal-state transport properties predicted and analyzed for high-T c oxide superconductors, Phys. Rev. B, 1988, 37, 7482–7490,  DOI:10.1103/PhysRevB.37.7482.
45. N. Ashcroft and N. Mermin, Solid State Physics, Saunders College Publishing, United State of America, English Edition, 1976.
46. M. M. Abdullah, F. M. Rajab and S. M. Al-Abbas, Structural and optical characterization of Cr2O3 nanostructures: Evaluation of its dielectric properties, AIP Adv., 2014, 4, 027121,  DOI:10.1063/1.4867012.
47. S. K. Deb, Physical properties of a transition metal oxide: optical and photoelectric properties of single crystal and thin film molybdenum trioxide, Proc. R. Soc. London, A, 1997, 304, 211–231,  DOI:10.1098/rspa.1968.0082.
48. K. H. Kim, S. H. Shim, K. B. Shim, K. Niihara and J. Hojo, Microstructural and Thermoelectric Characteristics of Zinc Oxide-Based Thermoelectric Materials Fabricated Using a Spark Plasma Sintering Process, J. Am. Ceram. Soc., 2005, 88, 628–632,  DOI:10.1111/j.1551-2916.2005.00131.x.
49. S. A. Carter, T. F. Rosenbaum, P. Metcalf, J. M. Honig and J. Spalek, Mass enhancement and magnetic order at the Mott-Hubbard transition, Phys. Rev. B, 1993, 48, 16841–16844,  DOI:10.1103/PhysRevB.48.16841.
50. E. H. Greener, F. J. Barone and W. M. Hirthe, Electrical Conductivity of Single and Polycrystalline Near-Stoichiometric Rutile in the Range 600° to 1400 °C, J. Am. Ceram. Soc., 1965, 48, 623–627,  DOI:10.1111/j.1151-2916.1965.tb14692.x.
51. T. Dittrich, J. Weidmann, F. Koch, I. Uhlendorf and I. Lauermann, Temperature- and oxygen partial pressure-dependent electrical conductivity in nanoporous rutile and anatase, Appl. Phys. Lett., 1999, 75, 3980–3982,  DOI:10.1063/1.125513.
52. Y. Liu, Q. Chen, X. Du, X. Liu and P. Li, Effects of substrate on the structure and properties of V₂O₅ thin films prepared by the sol-gel method, AIP Adv., 2019, 9, 045028,  DOI:10.1063/1.5095718.
53. L. Ben-Dor and Y. Shimony, Crystal structure, magnetic susceptibility and electrical conductivity of pure and NiO-doped MoO₂ and WO₂, Mater. Res. Bull., 1974, 9, 837–844,  DOI:10.1016/0025-5408(74)90120-2.
54. H. Nagai, T. Fujikawa and K.-i. Shoji, Electrical Conductivity of Cr₂O₃ Doped with La₂O₃, Y₂O₃ and NiO*, Trans. Japan Inst. Met., 1983, 24, 581–588,  DOI:10.2320/matertrans1960.24.581.
55. P. Knauth and H. L. Tuller, Electrical and defect thermodynamic properties of nanocrystalline titanium dioxide, J. Appl. Phys., 1999, 85, 897–902,  DOI:10.1063/1.369208.
56. F. Liu, Z. Yang, Y. Luo, S. Guo, C. Zhang, S. Choo, X. Wang, A. Mkhoyan, M. Bernardi and B. Jalan, “Understanding Polaronic Transport in Anatase TiO₂ Films by Combining Precise Synthesis and First-Principles Many-Body Theory”, 2025.
57. Y.-L. Han, Z.-Z. Luo, C.-J. Li, S.-C. Shen, G.-L. Qu, C.-M. Xiong, R.-F. Dou, L. He and J.-C. Nie, Carrier-mediated Kondo effect and Hall mobility by electrolyte gating in slightly doped anatase TiO₂ films, Phys. Rev. B, 2014, 90, 205107,  DOI:10.1103/PhysRevB.90.205107.
58. L. Y. Lim, S. Lany, Y. J. Chang, E. Rotenberg, A. Zunger and M. F. Toney, Angle-resolved photoemission and quasiparticle calculation of ZnO: the need for d band shift in oxide semiconductors, Phys. Rev. B, 2012, 86, 235113,  DOI:10.1103/PhysRevB.86.235113.
59. A. M. Wernbacher, M. Eichelbaum, T. Risse, S. Cap, A. Trunschke and R. Schlögl, Operando electrical conductivity and complex permittivity study on vanadia oxidation catalysts, J. Phys. Chem. C, 2019, 123, 8005–8017,  DOI:10.1021/acs.jpcc.8b07417.
60. Z. Zhang and V. E. Henrich, Surface electronic structure of V₂O₅(001): defect states and chemisorption, Surf. Sci., 1994, 321, 133–144,  DOI:10.1016/0039-6028(94)90034-5.
61. C. Kittel, Introduction to Solid State Physics, Wiley, Hoboken, NJ, vol. 8, 2005.

---