Intrinsic Material Properties Dictating Oxygen Vacancy Formation ...
Intrinsic Material Properties Dictating Oxygen Vacancy Formation Energetics in Metal Oxides Ann M. Deml, Aaron M. Holder, Ryan P. OβHayre, Charles B. Musgrave, and Vladan StevanoviΔ Supplementary Information Contents ................................................................................................................................. Pg - Methods .................................................................................................................................... 2 - Table S1: List of 45 oxides used for fitting the πΈπ model and their calculated properties. ...... 4 - Figure S1: Plot of calculated πΈπ for 45 oxides versus a linear combination of Ξπ»π and πΈπ . ..... 6 - Figure S2: Plot of calculated πΈπ for 45 oxides versus a linear combination of 3 terms ............ 6 - Descriptions of candidate basis descriptors examined for the πΈπ model................................ 7 - Table S2: List of 18 oxides used for validation of the πΈπ model and their calculated properties. ................................................................................................................................. 8 - Table S3: Predicted πΈπ of ~1800 oxide materials. o Main group metal oxides .............................................................................................. 9 o Nonmagnetic transition metal oxides......................................................................... 14 o Magnetic transition metal oxides................................................................................ 19 - References ............................................................................................................................... 25
Deml S1
Methods Direct πΈπ calculations. Our study focuses on the relationship between the intrinsic properties of oxides and the formation energetics of charge neutral oxygen vacancies. Admittedly, depending on the position of the Fermi energy of the host material, the oxygen vacancy can assume a thermodynamically more favorable charged state. However, the energy to form the neutral ππ is relevant (i) as the upper limit for πΈπ regardless of whether the defect is charged or neutral and (ii) for n-type oxides in which the transition level between the charged and neutral ππ states resides inside the band gap. Standard supercell methods1β3 were used to compute the EV of a single, charge neutral ππ as ππππππ‘
ππππππ‘ πΈπ‘ππ‘
πΈπ = πΈπ‘ππ‘
βππ π‘ β πΈπ‘ππ‘ + ππ ,
(S1)
βππ π‘ πΈπ‘ππ‘
where and are the total energies of a supercell with and without the ππ , respectively, and ππ corresponds to the oxygen chemical potential characterizing the reservoir of oxygen atoms. Supercells with ππ were created by removal of a neutral oxygen atom followed by self-consistent optimization of the electron density. In the language of KrΓΆger-Vink notation, both the 2+ charged ππ and the charge compensating electrons are included in the supercell. Due to the explicit charge neutrality of the supercell, finite size corrections such as the charge-image and band-alignment corrections1 are not needed. Likewise, band gap corrections are unnecessary because all 45 oxides exhibit βdeepβ oxygen vacancies, i.e. the electrons that were previously participating in bonds with the removed oxygen atom occupy states well within the band gap. ππππππ‘
βππ π‘ To calculate πΈπ‘ππ‘ and πΈπ‘ππ‘ , supercells with 40-80 atoms were generated by replicating DFTrelaxed bulk unit cells taken from the Inorganic Crystal Structure Database (ICSD).4 Different supercell sizes result from differences in the unit cell size and symmetry. Distances between periodic ππ in neighboring supercells were on the order of 8-10 Γ with ππ concentrations ranging from 2-6 %. For the defect supercell calculations, only atomic positions were optimized while cell volumes and shapes were fixed at their bulk values. All non-equivalent O sites were sampled. Energy differences between different O sites were found to be small (β€0.1 eV); therefore, the energies were averaged. We find that the πΈπ resulting from this scheme are converged to within 0.1 eV with respect to supercell size (e.g. compared to 160 atoms in the case of BaTiO3). In all calculations of πΈπ the PBE exchange-correlation functional5 is employed together with the addition of the onsite Hubbard term (PBE+U)6 as implemented within the projector augmented wave (PAW) method7 in the VASP code.8 A Monkhorst-Pack k-point sampling9 was applied with all total energies converged to within 3 meV/atom with respect to the number of kpoints. We chose an energy cutoff of 340 eV corresponding to a value 20 % greater than the highest cutoff energy suggested by the employed pseudopotentials (282 eV for oxygen). In this work, we consider a broad range of main group and transition metal oxides. A constant Hubbard correction of U=3 eV was applied to d-orbitals of all transition metals except Cu and Ag for which U=5 eV was used consistent with the numerical setup and findings from Ref. 10. Spin degrees of freedom were treated explicitly. For systems containing transition metals, we enumerated all possible magnetic configurations on a primitive unit cell and used the lowest energy configurations in the subsequent defect calculations. In the case of magnetic rock salt structures we used the (111) antiferromagnetic superlattice, known from experiments to be the magnetic ground state.11
Deml S2
The set of computational values described above produces accurate defect energetics1 as well as accurate oxide enthalpies of formation using the Fitted Elemental Reference Energies (FERE) approach.10 FERE elemental energies (chemical potentials) correspond to standard state elemental phases; therefore, our calculated πΈπ using the FERE oxygen chemical potential (πππΉπΈπ πΈ ) correspond to the standard state conditions of gaseous oxygen but can easily be adjusted for other T and ππ2 conditions. Electronic structure of the host. To accurately model the electronic structure properties of the bulk (host) systems, in particular the band gap (πΈπ ), spin polarized many-body GW calculations12 within the PAW implementation of the VASP code were performed. The structures were relaxed with DFT+U, as detailed above, to obtain the initial DFT+U structures, eigenenergies, and wave functions prior to the quasiparticle energy (QPE) calculations in GW. The lattice parameters of compounds containing Ge and Sn were scaled to experimental values due to the βsoftβ nature of these systems. Occupation independent on-site potentials were used for the 3d states of transition metals Ti, Cr, Mn, Fe, Co, and Ni for improved agreement with experimental band gap energies according to the method developed by Lany.13 No on-site potentials are available for elements possessing f electrons (e.g. La, Ce); therefore, we consider only those without f electrons. The same k-point sampling and energy cutoff used for DFT+U calculations were applied. The total number of bands was taken as 60 times the number of atoms in the unit cell. As shown in Figure 3, πΈπ calculated using GW exhibit close agreement with experimental values while πΈπ calculated using DFT+U generally underestimate experimental values, consistent with numerous previous reports.13β15 All gaps correspond to minimum (fundamental) gaps. Analysis tools. A stepwise linear regression approach implemented in JMP16 was used to investigate possible models for πΈπ and to select the most statistically-significant subset of candidate descriptors using the corresponding p-values. The resultant πΈπ descriptors include, in the order of decreasing significance: i) the oxide enthalpy of formation (Ξπ»π ), ii) the band gap (πΈπ ), iii) the average atomic electronegativity difference between O and its first-nearest neighbors (β©Ξπβͺ), and iv) the O-2p band center relative to the valence band maximum (πΈππ΅π β πΈπ2π ). The cumulative inclusion of these contributions results in R2 values of 0.71, 0.85, 0.92, and 0.97, respectively, while inclusion of additional terms increases R2 by
Deml S1
Methods Direct πΈπ calculations. Our study focuses on the relationship between the intrinsic properties of oxides and the formation energetics of charge neutral oxygen vacancies. Admittedly, depending on the position of the Fermi energy of the host material, the oxygen vacancy can assume a thermodynamically more favorable charged state. However, the energy to form the neutral ππ is relevant (i) as the upper limit for πΈπ regardless of whether the defect is charged or neutral and (ii) for n-type oxides in which the transition level between the charged and neutral ππ states resides inside the band gap. Standard supercell methods1β3 were used to compute the EV of a single, charge neutral ππ as ππππππ‘
ππππππ‘ πΈπ‘ππ‘
πΈπ = πΈπ‘ππ‘
βππ π‘ β πΈπ‘ππ‘ + ππ ,
(S1)
βππ π‘ πΈπ‘ππ‘
where and are the total energies of a supercell with and without the ππ , respectively, and ππ corresponds to the oxygen chemical potential characterizing the reservoir of oxygen atoms. Supercells with ππ were created by removal of a neutral oxygen atom followed by self-consistent optimization of the electron density. In the language of KrΓΆger-Vink notation, both the 2+ charged ππ and the charge compensating electrons are included in the supercell. Due to the explicit charge neutrality of the supercell, finite size corrections such as the charge-image and band-alignment corrections1 are not needed. Likewise, band gap corrections are unnecessary because all 45 oxides exhibit βdeepβ oxygen vacancies, i.e. the electrons that were previously participating in bonds with the removed oxygen atom occupy states well within the band gap. ππππππ‘
βππ π‘ To calculate πΈπ‘ππ‘ and πΈπ‘ππ‘ , supercells with 40-80 atoms were generated by replicating DFTrelaxed bulk unit cells taken from the Inorganic Crystal Structure Database (ICSD).4 Different supercell sizes result from differences in the unit cell size and symmetry. Distances between periodic ππ in neighboring supercells were on the order of 8-10 Γ with ππ concentrations ranging from 2-6 %. For the defect supercell calculations, only atomic positions were optimized while cell volumes and shapes were fixed at their bulk values. All non-equivalent O sites were sampled. Energy differences between different O sites were found to be small (β€0.1 eV); therefore, the energies were averaged. We find that the πΈπ resulting from this scheme are converged to within 0.1 eV with respect to supercell size (e.g. compared to 160 atoms in the case of BaTiO3). In all calculations of πΈπ the PBE exchange-correlation functional5 is employed together with the addition of the onsite Hubbard term (PBE+U)6 as implemented within the projector augmented wave (PAW) method7 in the VASP code.8 A Monkhorst-Pack k-point sampling9 was applied with all total energies converged to within 3 meV/atom with respect to the number of kpoints. We chose an energy cutoff of 340 eV corresponding to a value 20 % greater than the highest cutoff energy suggested by the employed pseudopotentials (282 eV for oxygen). In this work, we consider a broad range of main group and transition metal oxides. A constant Hubbard correction of U=3 eV was applied to d-orbitals of all transition metals except Cu and Ag for which U=5 eV was used consistent with the numerical setup and findings from Ref. 10. Spin degrees of freedom were treated explicitly. For systems containing transition metals, we enumerated all possible magnetic configurations on a primitive unit cell and used the lowest energy configurations in the subsequent defect calculations. In the case of magnetic rock salt structures we used the (111) antiferromagnetic superlattice, known from experiments to be the magnetic ground state.11
Deml S2
The set of computational values described above produces accurate defect energetics1 as well as accurate oxide enthalpies of formation using the Fitted Elemental Reference Energies (FERE) approach.10 FERE elemental energies (chemical potentials) correspond to standard state elemental phases; therefore, our calculated πΈπ using the FERE oxygen chemical potential (πππΉπΈπ πΈ ) correspond to the standard state conditions of gaseous oxygen but can easily be adjusted for other T and ππ2 conditions. Electronic structure of the host. To accurately model the electronic structure properties of the bulk (host) systems, in particular the band gap (πΈπ ), spin polarized many-body GW calculations12 within the PAW implementation of the VASP code were performed. The structures were relaxed with DFT+U, as detailed above, to obtain the initial DFT+U structures, eigenenergies, and wave functions prior to the quasiparticle energy (QPE) calculations in GW. The lattice parameters of compounds containing Ge and Sn were scaled to experimental values due to the βsoftβ nature of these systems. Occupation independent on-site potentials were used for the 3d states of transition metals Ti, Cr, Mn, Fe, Co, and Ni for improved agreement with experimental band gap energies according to the method developed by Lany.13 No on-site potentials are available for elements possessing f electrons (e.g. La, Ce); therefore, we consider only those without f electrons. The same k-point sampling and energy cutoff used for DFT+U calculations were applied. The total number of bands was taken as 60 times the number of atoms in the unit cell. As shown in Figure 3, πΈπ calculated using GW exhibit close agreement with experimental values while πΈπ calculated using DFT+U generally underestimate experimental values, consistent with numerous previous reports.13β15 All gaps correspond to minimum (fundamental) gaps. Analysis tools. A stepwise linear regression approach implemented in JMP16 was used to investigate possible models for πΈπ and to select the most statistically-significant subset of candidate descriptors using the corresponding p-values. The resultant πΈπ descriptors include, in the order of decreasing significance: i) the oxide enthalpy of formation (Ξπ»π ), ii) the band gap (πΈπ ), iii) the average atomic electronegativity difference between O and its first-nearest neighbors (β©Ξπβͺ), and iv) the O-2p band center relative to the valence band maximum (πΈππ΅π β πΈπ2π ). The cumulative inclusion of these contributions results in R2 values of 0.71, 0.85, 0.92, and 0.97, respectively, while inclusion of additional terms increases R2 by
