ACS Publications. Most Trusted. Most Cited. Most Read
My Activity

Figure 1Loading Img
RETURN TO ISSUEPREVC: Physical Properti...C: Physical Properties of Materials and InterfacesNEXT

Properties of Oxygen Vacancy and Hydrogen Interstitial Defects in Strontium Titanate: DFT + Ud,p Calculations

  • Szymon Winczewski*
    Szymon Winczewski
    Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, Poland
    *Email: [email protected]
  • Jacek Dziedzic
    Jacek Dziedzic
    Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, Poland
    School of Chemistry, University of Southampton, Highfield, SouthamptonSO17 1BJ, U.K.
  • Tadeusz Miruszewski
    Tadeusz Miruszewski
    Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, Poland
  • Jarosław Rybicki
    Jarosław Rybicki
    Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, Poland
    TASK Computer Centre, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, Poland
  • , and 
  • Maria Gazda
    Maria Gazda
    Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, Poland
    More by Maria Gazda
Cite this: J. Phys. Chem. C 2022, 126, 43, 18439–18465
Publication Date (Web):October 24, 2022

Copyright © 2022 The Authors. Published by American Chemical Society. This publication is licensed under

CC-BY 4.0.
  • Open Access

Article Views





PDF (8 MB)
Supporting Info (1)»


This work presents extensive theoretical studies focused on the mixed ion-electron transport in cubic strontium titanate (STO). A new approach to the description of this difficult system was developed within the framework of linear-scaling Kohn–Sham density functional theory, as realized in the ONETEP program. The description we present is free of any empirical parameters and relies on the Hubbard U and Hund’s J corrections applied to both Ti and O atoms. The proposed methodology was validated by considering perfect STO. Its calculated properties were found to be in close agreement with experiments and calculations at higher levels of theory. The validated approach was subsequently employed to study the oxygen vacancy (VO) and the hydrogen interstitial (IH), using very large supercells (625 ± 1 atoms). The relaxed configurations of defects were obtained through fastidious energy minimization and later analyzed from a number of perspectives. The calculated defect formation energies and charge transition levels (CTLs) were found to be in close agreement with the experiment. With the exception of the charge-neutral VO, all considered defects were found to introduce shallow states, located down to 0.2 eV below the conduction band. Our calculations revealed a large 1 eV difference in the thermodynamic and optical CTLs of the neutral VO, explaining the inconsistencies observed─till now─between conduction and spectroscopic measurements. The influence of defects on the bonding characteristics and the crystalline structure of STO was quantified, revealing that both VO and IH defects lead to a significant polarization and strong tilting of the TiO6 octahedra.

This publication is licensed under

CC-BY 4.0.
  • cc licence
  • by licence

1. Introduction

Jump To

1.1. Strontium Titanate

Strontium titanate (SrTiO3, STO) is a prototype perovskite material. At high temperatures, it crystallizes in a cubic structure (space group Pmm, see Figure 1), while at low temperatures, it adopts a tetragonal structure (space group I4/mcm) due to an antiferrodistortive phase transition which occurs at 105 K. (1−3)

Figure 1

Figure 1. Unit cell of the cubic SrTiO3 (panel a) and the 5 × 5 × 5 supercell used in the calculations (panel b). Panel c presents the initial positioning of the hydrogen interstitial with respect to the closest TiO6 octahedra. All visualizations presented in this work were prepared using OVITO. (156)

In the last few decades, significant effort has been dedicated to understanding the structural, (4−9) mechanical, (10−12) optical, (13−24) and transport (25−32) properties of STO. The influence of various dopants on these properties has also been extensively studied. (28,33−35) This high interest in STO is driven by the variety of its potential applications, which span fields such as electronics, optoelectronics, and photocatalysis. (36)
In recent years, STO has also become the subject of an increased interest in electrochemistry, as it has been demonstrated to be a good base material for the fabrication of the so-called triple-conducting oxides (TCOs). (34,37,38) These materials are capable of conducting oxygen (O2–) and hydrogen (H+) ions, and electrons or holes, and are therefore considered promising for use as novel cathodes in proton-conducting fuel cells. (39,40) One of the examples of promising STO-based TCO is SrTixFe1–xO3−δ, which was found to be electrochemically stable, offering high ionic conductivity and good ion surface exchange rates. (35,41−43)
In the last few years, significant advances were made in the experimental studies focusing on the synthesis, testing, and application of various TCOs. (37−40,44) However, their better understanding still requires considerable theoretical work because even the base materials are still poorly understood in the context of mixed electron-ion transport. This knowledge gap mainly originates from the difficulties encountered in theoretical modeling studies. (45) The related issues can be well illustrated using the example of STO.

1.2. Difficulties in Theoretical Description

Understanding the mixed transport properties of STO requires a suitable description of all three constituent types of transport. However, even a description of the static properties of charge carriers constitutes a challenging task, as the theoretical methods used to date have been either incapable of describing the electronic structure correctly or have been too computationally expensive for studies concerning defects.
The local and semi-local density functional theory (DFT) methods that have been used typically in this context are known to poorly describe strong electron correlation effects, which leads to a poor description of the band structure with significantly underestimated band gaps. For STO, DFT calculations employing the local density approximation (LDA) and the generalized gradient approximation (GGA) for exchange–correlation (XC) interactions typically yield an indirect band gap of 1.6–2.0 eV, (46−52) underestimating the experimental value (3.25 eV (31)) by as much as 40–50%.
Methods based on higher levels of theory (e.g., DFT using hybrid functionals (46,47,49,53−57) or various GW approaches (58−61)) perform considerably better in describing the band structure (band gaps very close to the experiment were reported in refs (49,58,59,and61)). However, because of their extremely high computational cost, they cannot, as of today, be fully applied to defects.
The investigation of defect properties needs to begin with structure relaxation, which is a highly demanding task, even when local and semi-local DFT methods are used. It cannot be performed with the higher level of theory applied throughout the entire minimization with the currently available computational resources. The apparent exceptions found in the literature (62) report on the defect properties calculated at a higher level of theory but for structures optimized at a lower level.
One of the possible remedies to this accuracy-cost problem is the application of DFT + U approaches. These methods are characterized by a computational cost comparable to semi-local DFT, while providing significant improvement in the quality of the description, which is often comparable to DFT with hybrid functionals. The last two decades have witnessed considerable advances in the theory of DFT + U. (63−65) As a consequence, DFT + U has gained significant popularity, becoming a standard tool for treating difficult (strongly correlated) systems containing transition-metal atoms. (66)
Within DFT + U, the standard DFT energy functional is augmented by a Hubbard term, whereby additional Coulomb repulsion is applied between the strongly localized (typically d and f) electrons belonging to the same Hubbard atom. The magnitude of this interaction is proportional to the occupation numbers and is controlled by the Hubbard U parameter. In so doing, DFT + U addresses the self-interaction error, a well-known deficiency of standard DFT.
DFT + U can be extended to DFT + U + J by additionally correcting the interactions between unlike spins and adding an appropriate exchange interaction term in the corrective functional. The magnitude of thus introduced Hund’s exchange term is controlled by another parameter, typically denoted with J. The +J correction has been demonstrated to be crucial in describing materials with strong spin–orbit coupling. (63,64)
The application of DFT + U + J requires a prior determination of the correction parameters, U and J. Their choice can be motivated on theoretical grounds using the linear response approach, (67) its recent reformulations, (68,69) or other approaches. (70,71)
Although theoretical foundations for finding U and J exist, in many cases these parameters are determined empirically by fitting to experimental properties, like the band gap width. (66) Such an approach often leads to unphysical (i.e., very large) U values, (72−76) often resulting in a poor description of the properties which were not taken into account when choosing U. More often than not, the values of U and J are transferred between various systems containing the same corrected element, (77) without regard for their dependence on the material. It should be emphasized that U and J are also sensitive to the methodology, in the sense that they depend on the details of the theoretical approach used, as the exact form of the corrective functional, the form of the projectors used to calculate the occupancies for the correction, and the XC functional used.
Because of the above, the (U, J) values used in the literature often differ significantly between different studies concerning the same material. (66) For example, two different (U, J) sets were used for correcting Ti 3d electrons in STO (all given in eV): (5.0, 0.64) proposed by Pavarini et al. (78) (used in refs (48,77,79–87)) and (3.2, 0.90) proposed by Solovyev et al. (88) (used in refs (48,80,and89)). We believe that neither of the above sets is, in fact, well justified, as they were transferred from studies concerning LaTiO3 (Solovyev’s set, Pavarini’s set) and YTiO3 (Pavarini’s set). There are also other choices made by individual authors. (62,75,90−92)
This inability to describe the underlying physics in a justified manner is reflected in the published results. The existing studies, which employed DFT + U (and +J) methods and focused on oxygen vacancies (VO) in STO, have provided conflicting reports concerning the location and nature of defect levels introduced by the neutral VO. DFT + U studies reported levels that are shallow or located within the conduction band (CB), while higher-level calculations typically predicted deeper locations (0.7 eV below the CB). This problem has been recently extensively addressed by Ricca et al., (93) who questioned many of the previous results obtained for VO in STO with DFT + U approaches. These failures can be attributed to the fact that standard DFT + U (which, when applied to STO, corrects only Ti 3d electrons through on-site terms only) is still unable to faithfully reproduce the band structure of perfect STO, offering only a minor improvement over semi-local DFT and yielding band gaps (typically 2.0–2.4 eV), which are still in poor agreement with the experiment, among other deficiencies.
The last few years have provided new developments in the corrective approaches. The proposed extensions (e.g., methods such as DFT + U + V, (94) which also corrects off-site interactions) are now emerging as a valuable tool for studying STO and similar challenging systems. Ricca et al. (93) very recently demonstrated that DFT + U + V is capable of reproducing STO’s band gap width with an accuracy comparable to that of hybrid DFT calculations. Crucially, this was done using U and V parameters determined on theoretical grounds.
Recent work by O’Regan et al. (95) has shown that a similar improvement can also be achieved within the DFT + U + J framework. They demonstrated that the application of the +U +J corrections also to the 2p electrons of the oxygen atoms (in addition to Ti 3d electrons) considerably improves the predicted fundamental band gap of TiO2, which constituted a long-standing challenge for DFT calculations. The band gaps calculated agreed almost perfectly with the experiment (differing by only 0.01 and 0.03 eV, for rutile and anatase, respectively). Reference (95) also achieved this using well-justified methods of choosing the corrective parameters. The U and J parameters were determined from the theory using the so-called minimum-tracking linear-response (MT-LR) method. This approach─also proposed recently by O’Regan and co-workers (69) ─can be considered as an analogue of the conventional linear response approach of Cococcioni and de Gironcoli, (67) with their procedure being adapted to the case of direct-minimization codes, like ONETEP (96,97) used in refs (69and95).

1.3. Difficulties in Modeling of Defects

The accurate theoretical description is not the sole difficulty. By default, the calculations on point defects employ periodic boundary conditions (PBCs) and the so-called supercell approach. The use of PBCs means that the simulated defect is in fact a periodic array of defects. Particular care must be taken to guarantee that the calculated properties are not affected by spurious defect–defect interactions (mechanical, electrostatic, or magnetic).
This problem can be overcome─at least in principle─by extrapolating to infinity the results obtained for a series of finite systems, thus reaching the dilute defect limit. However, the available computational power constraints the system sizes tractable in practice. These typically do not exceed 320 atoms (corresponding to a supercell containing 4 × 4 × 4 replicas of the five atom unit cell) for perovskite-type materials, for semi-local DFT methods or methods with a similar cost.
Additionally, not all system sizes are equivalent. A study carried out by Zhang et al. (86) has shown that there are better and worse system sizes n × n × n. This conclusion has been reached based on the calculated defect formation energies (VO in STO was considered). Energies obtained for even n’s turned out to be systematically higher than those obtained for odd n’s. Reference (86) explained this effect referring to the relaxation of the crystal structure, which occurred after the introduction of the defect. The possible relaxation was limited in supercells with even values of n, which restricted the rotation of the TiO6 octahedra.
Studies concerning defects often consider their various charge states in order to predict which state is energetically preferred at the given conditions specified by the Fermi level. (98−103) By doing so, other quantities of interest─like thermodynamic and optical charge transition levels (CTLs)─can also be accessed.
In the above context, charge-neutral defects do not pose significant difficulty. In contrast, charged defects─which most often constitute the main interest─require additional treatment, that is, the calculation of an image interaction correction (IIC), with the aim of removing the spurious contributions to the energy which originate from the electrostatic interactions of the defect charge density with its periodic images (an artifact of PBCs).
Another correction must be applied in calculations on charged defects performed in PBCs. In such calculations, a uniform neutralizing background charge, the so-called jellium, must be introduced to make the Ewald summation convergent. The introduction of jellium means that the electrostatic potentials of charged and neutral defects cannot be directly compared (the introduction of jellium corresponds to setting the average potential to zero). Therefore, before evaluating the formation energy of the charged defect, a so-called potential alignment (PA) must be performed to correctly account for the energy required to exchange the electrons in the process of defect charging.
A number of approaches to computing both corrections (IIC and PA) have been proposed over the last three decades. (98,99,104−117) More often than not, these corrections are either omitted (81,83) or are treated without the diligence they deserve (the underlying assumptions are not verified), even though they might give contributions to the energy which are significant in the context of mixed electron–ion transport. This requires verification.
Comparing theoretical calculations with the experiment also constitutes a challenge, and there are several reasons for that.
Theoretical calculations are typically carried out at 0 K. Therefore, comparing the properties measured in the experiment requires extrapolating results obtained for finite temperatures down to 0 K, which may introduce discrepancies. For example, the experimental determination of defect formation energies can be done by measuring the temperature dependence of conductivity and applying the methods of defect chemistry, which rely on a number of assumptions. (29)
Differences between the conditions used in experiments and those assumed in calculations also manifest themselves through differences in the chemical potentials. Typically, in theoretical calculations, chemical potentials are computed, and as such they are known precisely. However, they are affected by the errors originating from the limitations of the chosen theoretical description. On the other hand, in the experiment chemical potentials are difficult to precisely control and measure. Many properties of defects, for example, their formation energies, strongly depend on chemical potentials, which, in fact, specify thermodynamic conditions. An appropriate comparison between calculations and experiments requires ensuring that the same conditions are considered, that is, that both use the same sets of chemical potentials.
Some possibilities of comparison are also provided─at least in principle─by the measurements of optical properties. The locations of defect levels, that is, the electron states introduced by defects are accessible from calculations. However, also here, considerable difficulties arise. Again, the related problems can be well illustrated on the example of STO.
With the exception of the charge-neutral VO, the considered defects (oxygen vacancies in the 2+ and 1+ charge states, hydrogen interstitial) are problematic for optical observation because they do not introduce any states, or only introduce states which are very close to the CB (≈0.1 eV below). This proximity makes it challenging to experimentally obtain reliable information, which could be used to validate the theory.
The charge-neutral VO introduces states, which might be seen as deep (≈1.1 eV below the CB (14,15)). However, some optical experiments measured shallow locations (≈0.1 eV below the CB (16)). The conduction measurements of the ionization energies also predicted a shallow location. (25,26) This significant discrepancy is the reason why the charge-neutral VO attracted considerable theoretical interest. (53,56,57,62,77,81−83,86,87,93,118−132)
However, the calculations carried out to date could not explain the discrepancies seen in the experimental characterization of the charge-neutral VO. Various theoretical studies gave locations of defect states which differed significantly and disagreed with spectroscopic and conduction measurements. This problem was discussed in detail in ref (93), to which we refer the interested reader. We will return to this problem later in this article, as it constitutes an essential aspect of this work. The methodology developed and applied in this work was partially validated by considering this problematic case of the charge-neutral VO. We believe that a careful validation of the proposed approach should always constitute an indispensable test of its adequacy and quality.
Although there are many theoretical reports focusing on the oxygen vacancies in STO, the properties of hydrogen interstitials were not a subject of any broader study. The existing reports have serious limitations: using small supercells (72,73,133−135) (typically 2 × 2 × 2 or 3 × 3 × 3) and/or employing inadequate descriptions (local and semi-local DFT, or DFT + U with unphysically large U = 8.5 eV, (72,73,134,136) or even simpler methods (137−139)). Furthermore, to the best of our knowledge, no single study focused on all three types of charge carriers (electrons, oxygen and hydrogen ions), studying them on the same theoretical footing within the context of mixed electron-ion transport. That hinders the assessment of their role in STO.

1.4. Aims and Outline of the Work

Building on the above observations and recent theoretical developments, we decided to fastidiously investigate the oxygen vacancy and hydrogen interstitial defects in STO. Our focus on these defects is motivated by the long-term scope of the research undertaken, a theoretical characterization of various STO-based (and similar) TCOs, with a particular interest in their mixed electron-ion transport. From this perspective, this work paves the way for future studies, constituting a proof of concept focused on developing and validating the underlying methodology (presented in Section 2).
In particular, we studied VO and IH defects in STO using the DFT + Ud,p approach, recently proposed by O’Regan et al. (95) After validating the methodology (on perfect STO, Section 3), we considered a number of charge states of these defects, performed their geometry relaxation and characterized their various properties. Here, particular care was taken to correctly account for all essential theoretical and computational ingredients discussed in the above introduction. The results of this extensive study are presented in Section 4, where we report on formation energies, electronic structures, defect levels, CTLs, structural relaxations, and other properties. We conclude with a discussion and summary in Section 5, outlining the directions for future works.

2. Methodology

Jump To

2.1. ONETEP Approach

All the calculations reported in this article were performed using ONETEP (96,97) (order-N electronic total energy package), a linear-scaling approach for quantum-mechanical calculations based on DFT.
Within ONETEP, the cubic scaling limitation of conventional DFT is overcome by working directly with the density matrix
which is constructed from localized orbitals ϕα and ϕβ centered on atoms. The orbitals are nonorthogonal-generalized Wannier functions (NGWFs), which can be chosen to be real-valued. They are expressed in an underlying basis of periodic cardinal sine (psinc) functions defined in real space on a regular Cartesian grid. The grid spacing dpsinc controls the quality of the basis set and can be related to the kinetic energy cutoff Ecut of traditional plane-wave (PW) DFT. Another parameter that affects the quality of the computations is the localization radius of the NGWFs (Rϕ), beyond which they strictly vanish.
The linear-scaling property of the ONETEP approach is a consequence of the so-called nearsightedness of electronic matter. (140,141) This phenomenon makes the density kernel Kαβ short-ranged in systems with a large band gap. Therefore, it can be truncated for |RαRβ|>RK, leading to a highly sparse matrix representation, when large systems are studied. The parameter RK is known as the density kernel cutoff.
The sparsity of Kαβ and the localization of {ϕα} are exploited to obtain a linear-scaling formalism. For this purpose, ONETEP utilizes sparse matrix algebra and the so-called FFT box technique. (142) Here, as we only want to outline the most critical aspects of ONETEP approach briefly, we omit a detailed discussion of the underlying formalism, referring the interested reader to refs (96,97,143–147).
Within ONETEP, the ground state is found as in standard DFT, that is, through the variational principle. The minimization of the system’s total energy is achieved using a two-level loop. The inner loop minimizes the energy with respect to the density kernel elements with the NGWFs fixed, while the outer loop optimizes the psinc coefficients of the NGWFs. The minimization continues until the convergence with respect to both quantities─that is, Kαβ and {ϕα}─is achieved.
In ONETEP─similarly to other approaches to DFT─the quality of the obtained results also depends on the basis set used to represent the orbitals. In this respect, the ONETEP approach is similar to the linear combination of atomic orbitals (LCAO) approach, as it also uses atom-centered basis functions, specifically the NGWFs. A different number of NGWFs can be assigned to different atoms. The type (i.e., symmetry) and the spatial extent of NGWFs can also be controlled, making it possible to work with basis sets of different qualities (like split-valence basis sets) and characteristics (like bases containing diffuse functions).
What strongly distinguishes the ONETEP approach from LCAO approaches is the fact that during the calculation the shapes of the NGWFs change, as they are optimized in situ. Consequently, the basis functions adapt to the local chemical environment experienced by the atom on which they are centered. This approach─as it facilitates expressing the sought Kohn–Sham eigenlevels─significantly improves the description of various systems, which within ONETEP can be accurately modeled with near-complete basis set accuracy even with minimal basis sets. The use of a minimal basis also reduces the computational cost, facilitating studying large systems, which in the case of ONETEP may contain hundreds of thousands of atoms, depending on the system type and chosen methodology (see ref (96) and works cited therein).

2.2. Theoretical Description

In this work, we used the ONETEP approach to study defects (oxygen vacancy and hydrogen interstitial) in STO. We used the DFT method within the GGA with the Perdew–Burke–Ernzerhof (PBE) (148) XC functional, which was augmented with the Hubbard +U and Hund’s +J corrections using the DFT + Ud,p approach. (95) We motivate our choice of the PBE functional with the fact that in the calculations carried out, we also had to consider H2 and O2 molecules (used as the reference states in the calculation of defect formation energies). The versatility of PBE is the reason why we preferred it over GGA functionals that were specifically designed to model solids (like PBEsol (149)).
The innermost core electrons of Sr, Ti, and O were treated with norm-conserving pseudopotentials obtained from the Rappe library. (150,151) For Ti atoms, these sets do not include 3s or 3p semicore states, which were treated explicitly. Therefore, we used 13 NGWFs for each Ti atom, with the initial electron configuration being that of a neutral atom (3s23p63d24s24p0). The same number of 13 NGWFs was used for Sr atoms (initial electron configuration 4s24p64d05s25p0). In both cases (Sr and Ti), the additional p-type NGWFs served as diffuse functions and were found by numerical experiments to be necessary for obtaining a well-converged ground state. For all O atoms, we used 4 NGWFs (2s22p4 initial electron configuration). The hydrogen atom was described with one NGWF function.
For all NGWFs, we used the same localization radius Rϕ = 12 a0, where a0 denotes the Bohr radius. The density kernel was not truncated. The cutoff energy was taken as Ecut = 1090 eV, meaning the NGWFs were expanded with a dpsinc = 0.43 a0 spacing. The above values ensured well-converged results, as was concluded from the numerical tests performed for the perfect system.
The use of a relatively high energy cutoff is motivated by the fact that the performed calculations mainly comprised geometry optimizations, which require accurate forces for stable convergence. For the same reason, the stopping criterion for NGWF optimization was set to a more stringent value of 8 × 10–7Eha0–3/2 (Eh denotes the Hartree energy) and the density kernel was optimized on each update of the NGWFs and also in trial steps taken during the conjugate-gradient-based optimization of NGWFs. Denser than default working grids were also used for expressing the charge density and potential (with 340 grid points in each coordinate direction) to reduce numerical noise. Γ-Point sampling was used for reciprocal space integration. In the modeled supercells, this corresponds to sampling on a Γ-centered regular grid containing 5 × 5 × 5 (or 7 × 7 × 7, see later) points in the first Brillouin zone (FBZ).
In our studies, we also considered systems with an odd number of electrons (oxygen vacancy with a q = 1 + formal charge and neutral hydrogen interstitial) and systems with locally broken spin symmetry (singlet and triplet states of the neutral oxygen vacancy and a negatively charged hydrogen interstitial). Therefore, all the calculations were performed in the spin-polarized scheme, treating both spin channels independently. We note that within ONETEP spin-polarization is introduced through the dependence of the density kernel Kαβ on spin, but the basis functions (NGWFs) are spin-independent. The employed spin-polarization was also a consequence of the fact that we applied +J corrections.
What distinguishes our approach from other DFT + U + J studies on STO is that we applied +U and +J corrections to both Ti (3d subspace) and O (2p subspace) atoms, following the guidelines provided in refs (95,101,152,and153). Using the +U +J corrections requires choosing the so-called Hubbard projectors, which directly enter the definition of Hubbard and Hund’s terms (therefore, affecting the effectiveness of the +U +J corrections), as they are used to calculate occupancies of the strongly localized orbitals being corrected. This is done by projecting the calculated Kohn–Sham states (molecular orbitals or bands) onto a basis of localized atomic-like states. In this work, we used Hubbard projectors corresponding to the atomic NGWFs (obtained from the pseudoatomic solver), according to a method proposed by O’Regan et al. (144,145) and implemented in ONETEP. Here, we underline the particular importance of using localized Wannier-like projectors, which─according to many authors─is the most appropriate choice. (69,95,154,155) Our procedure did not include projector self-consistency. (144)
Although there are some sets of U and J parameters for correcting Ti in STO available in the literature, to the best of our knowledge, there are no parameters suitable for correcting oxygen 2p electrons in STO. Here, similarly to refs (64and69), we emphasize that the considered parameters are known to be case-specific, in the sense that they should not be transferred between different systems and, what is also important, between different methodologies. For these reasons, at the very beginning of our studies, we determined all the required parameters (i.e., UTi, JTi, UO and JO) from scratch for STO and the chosen methodology.

2.3. Calculations Performed

All the calculations were performed at zero temperature. A supercell consisting of 5 × 5 × 5 repetitions of the 5-atom cubic unit cell (Figure 1) was used, with N = 625 atoms in the perfect system. Where clearly indicated, a 7 × 7 × 7 supercell (with N = 1715 atoms) was used instead. PBCs were used in all directions. The lattice constant was set to the experimental value of a = 3.91 Å. (37)
Our simulation protocol consisted of the following steps. At the very beginning, we determined the U and J parameters for Ti and O atoms. This was done by studying the response of the perfect system using the MT-LR method (see Supporting Information, Section S1).
The obtained parameters were subsequently validated by comparing the results obtained with the chosen methodology with the results obtained with different methodologies and with experimental data. Our validation included comparing measures describing the band structure (band gap width Eg, total and partial densities of states) and measures characterizing the structure and energetics of the perfect system. To this end, we have calculated the equation of state (EOS) and, based on it, determined the equilibrium lattice constant a, the bulk modulus B0, and its pressure derivative B0 = ∂B0/∂p (both for zero pressure).
With U and J parameters at hand, we moved to systems containing defects, studying the oxygen vacancy with a formal charge of q = 0, 1+, and 2+ (denoted with VO0, VO1+, and VO2+) and the hydrogen interstitial with a formal charge of 1–, 0 and 1+ (IH1–, IH0, and IH1+, respectively). For the neutral oxygen vacancy, we considered two spin states, with two “excess” electrons (i.e., electrons remaining after the removal of the neutral O atom and located on the two neighboring Ti atoms) having initially either antiparallel (singlet state, VO0,S=0) or parallel (triplet, VO0,S=2) spins. We also considered a non-magnetic singlet state, with two excess electrons located on the same Ti atom and having antiparallel spins. However, we found this state to be highly unpreferable. Singlet and triplet spin states were also considered for IH1– (denoted as IH1–,S=0 and IH1–,S=2).
For all considered defects, we performed geometry relaxation, which started from the perfect crystal with one O atom removed (VOq) or one H atom added (IHq). In the initial structure, the hydrogen atom was positioned on the TiO2 plane, with the O–H bond length and the Ti–O–H valence angle taken as dO–H = 1 Å and θTi–O–H = 90° (Figure 1c).
To check the influence of the initial conditions (and numerical noise) on the obtained minima, we have also performed additional (and independent) geometry optimization, which started from structures with slightly distorted atomic positions (all atoms were displaced by 10–3 Å in random directions). We have considered one additional structure for each charge/spin state of the VOq defect and three additional structures for each state of the IHq defect. In total, we have considered 4 × 2 (for VO) + 4 × 4 (for IH) different structures, performing 24 separate geometry optimizations. For this purpose, we employed a combination of two geometry optimization techniques, which were used alternately in a loop.
In the beginning, the optimization was carried out in the delocalized internal coordinates (157) (DI), according to a method proposed in ref (158). Through numerical experiments, we found this method to be more efficient in the early stages of geometry relaxation.
The structure obtained from the DI minimization was later subjected to further relaxation, carried out in Cartesian coordinates with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm. (159−161) In these two initial stages of relaxation, we used the following convergence criteria ΔE < 10–7Eh/ion (change in the system’s total energy between two successive iterations), Fmax < 2 × 10–4Eh/a0 (the maximum force acting on any ion), and ΔRmax < 5 × 10–4a0 (maximum displacement of any ion between two successive iterations), requiring all of them to be met simultaneously at convergence.
This two-stage process─consisting of repeated DI- and BFGS-based minimizations─was continued until no significant drop (less than 3 meV) in the system’s total energy was observed throughout four consecutive minimizations (DI-, BFGS-, DI-, and BFGS-based). During this stage, more stringent (twice as strict as before) convergence criteria were used in the individual DI/BFGS minimizations. The criterion for NGWF convergence was also tightened to 2 × 10–7Eha0–3/2, to avoid “empty” iterations (corresponding to updates of the geometry with no significant update of the electronic structure).
Here, we note that the simulation cell was kept constant, as changing it would change the grid spacing dpsinc, resulting in varying Ecut and, in turn, the basis set. However, we performed a complete relaxation of the atomic coordinates, that is, without constraining the positions of any (even distant) atoms. The performed complete relaxation allowed us to observe effects such as structural polarization and octahedral tiltings and to account for their influence on the defect’s properties.

2.4. Formation Energies

The formation energies of defects were calculated using the well-known Zhang and Northup formula, (162) according to which the formation energy of a defect D in charge state q is given as (99,114)
Here, ED,q and EP represent the total energies of the defective and perfect systems, respectively. With ni, we denote the number of removed (ni < 0) or added (ni > 0) i-type atoms, μi refers to their chemical potential in the chosen reference state, and the third term in eq 2 represents contributions to ΔH that are due to the exchange of the atoms with the environment. With ϵVBM and ϵF, we denote the energy level of the valence band maximum (VBM) of the perfect system and the location of the Fermi level (with respect to the VBM), respectively. The corresponding fourth term represents the energy change originating from the exchange of electrons (charging of a defect).
The two remaining terms, EIIC and EPA, represent corrections. These corrections are applied a posteriori in calculations on charged defects in PBCs to remove artificial contributions to ED,q, which arise from spurious interactions between the considered charged defect, its periodic images, and the neutralizing background charge. The first of these (EIIC) is the so-called IIC, while the latter (EPA, typically written as qΔvPA) is the PA term. Both of these can be calculated using various approaches of different sophistication. (98,99,104−117)
When determining the IIC and PA corrections, we used a very recent method proposed by Durrant et al. (114) This method (for details, see Supporting Information, Section S3) is in many aspects similar to the well-known Freysoldt–Neugebauer–Van de Walle (FNV) procedure, (105) as it also relies on the comparison of the charge densities (determination of IIC) and electrostatic potentials (determination of PA) of perfect and defected (both neutral and charged) supercells. However, it does not assume any functional form of the charge density associated with the defect’s charge (like a Gaussian for the FNV procedure) and therefore, it is better suited to situations where the obtained defect is difficult to describe analytically (we found this to be the case in our calculations).

2.5. Charge Transition Levels

The calculated formation energies ΔH were used to determine CTLs (see Figure 2 for their graphical interpretation). This was done based on the framework proposed by Lany and Zunger. (99)

Figure 2

Figure 2. Thermodynamic (a) and optical (b,c) CTLs. In the thermodynamic transition, the increase of the Fermi level ϵF above the corresponding TTL (ϵth(q/q′), see eq 3) causes a qq′ = q – 1 charge transition. This process corresponds to a thermally induced excitation of an electron (obtained from the host material’s valence band) to a state introduced by the defect. The opposite process may also occur during the lowering of ϵF. In the optical transition, the absorption (or emission) of a photon causes an excitation (or recombination) of the electron, which passes from the defect’s level to the CB of the host material (or vice versa). Both optical processes result in a change of the defect’s charge: qq′ = q + 1 (for the absorption) or qq′ = q – 1 (for the emission). The energy of the absorbed (or emitted) photon defines the location of the corresponding OTL (ϵth(q/q′), see eq 4). In the thermodynamic transition, both states of the defect (initial and final) correspond to the relaxed configurations (D, q and D, q′, respectively), while in the optical transitions, only the initial state is relaxed (D, q).

The thermodynamic transition levels (TTLs), ϵth(q/q′), were calculated from the equation
The TTL ϵth(q/q′), which in the above expression is calculated relative to VBM of the perfect system, describes a thermally induced charge transition qq′ (electron excitation or trapping). ϵth(q/q′) corresponds to the Fermi level at which the formation energies of the defect in charge states q and q′ become equal. Here, because thermal excitations are considered, the defect is allowed to relax in both charge states q and q′. Therefore, both formation energies from eq 3 are evaluated for the relaxed configurations. Based on ϵth(q/q′), the thermal ionization energy of the defect can be determined. For example, the distance of ϵth(1+/0) from CBM (and of ϵth(0/1−) from VBM) measures the thermal ionization energy of simple donors (and acceptors).
We have also calculated the optical transition levels (OTLs), ϵop(q/q′), which measure the optical absorption (or emission) energy related to a qq′ excitation (or recombination). The ϵop(q/q′) levels (relative to CBM, with ϵF = Eg) were calculated from the equation
Because optical (vertical) transitions are considered (which occur over extremely short time scales), both formation energies are evaluated for the atomic configuration corresponding to the initial charge state q, according to the Franck–Condon principle.
The expression for ϵop(q/q′) can be written in the following form
Here, ϵCBM stands for the location of the conduction band minimum (CBM), while ED,q|D,q and ED,q represent the total energies of two supercells containing the defect with a charge state q′ and q, respectively. In the interest of brevity, we omitted electrostatic corrections EIIC and EPA in eq 5. They were, however, calculated (using Durrant’s procedure, see Supporting Information, Section S3) and added to the two total energies present in eq 5.
Based on Janak’s theorem, (163) the energy difference in eq 5 can be approximated with: (100,102,103)
which holds for q′ = q – 1 (recombination), and with
which is valid for q′ = q + 1 (excitation). Here, ϵiKS(q,N) denotes the i-th single-particle Kohn–Sham eigenvalue of an N-electron system with charge q. The notation |D,q serves as a reminder that the eigenvalue must be evaluated for the atomic configuration corresponding to charge state q. The eigenvalues, which appear in eqs 6 and 7, correspond to the highest occupied (ϵN+1KS(q,N+1) and ϵNKS(q,N)) and lowest unoccupied (ϵN+1KS(q,N) and ϵNKS(q,N1)) molecular orbitals. Similar─although more sophisticated─expressions can be derived for transitions (recombinations/excitations), which correspond to the removal/addition of two electrons (with q′ = q ± 2). The notation used in eqs 6 and 7 does not distinguish spin channels. It is noted that the relevant (majority- or minority-spin) eigenvalues should be considered.
The approximations defined by eqs 6 and 7 can be used to calculate ϵop(q/q′) based on the Kohn–Sham eigenvalues ϵiKS. It must be emphasized that─because of the spurious defect–defect interactions─the eigenvalues calculated for charged systems must be also corrected. This can be done using the formula (103,111)
which defines the correction that must be added to the eigenvalue corresponding to the defect state ϵdKS. Here, the electrostatic corrections (EIIC and EPA) were also calculated using Durrant’s procedure (see Supporting Information, Section S3).
In the literature there are many reports (see e.g., refs (100and103)) which employed the approximations defined by eqs 6 and 7 in the determination of OTLs ϵop(q/q′). More often than not, the accuracy of these approximations was not tested (for a discussion of their correctness see ref (164)). In this work, we verified their validity practically by comparing the approximated OTLs with those calculated using the exact expression given by eq 5 (for results of this comparison, see Supporting Information, Section S4).

3. Methodology Validation─Perfect STO

Jump To

3.1. Validation Outline

From the MT-LR method, we obtained the following values of the U and J parameters (for details, see Supporting Information, Section S1): UTi = 4.111 eV, JTi = 0.374 eV, UO = 9.097 eV, and JO = 0.935 eV. To validate these and the chosen approach, we calculated the electronic ground state of perfect STO at the experimental density (a = 3.91 Å). We also calculated its EOS by finding the ground state for eight other lattice constants. Here, to maintain the same accuracy, we had to restrict ourselves to system sizes that differ by an integer number of psinc grid spacings dpsinc. (165−167)
To directly demonstrate how the chosen approach (employing +U +J corrections for both Ti and O atoms, denoted as DFT + Ud,p) improves the description of STO, similar calculations were also performed for two different methodologies. Within these, (i) the +U +J corrections were not applied (neither for Ti nor for O) or (ii) were applied only to Ti atoms (with the U and J parameters for Ti taken as in DFT + Ud,p). These two approaches will be denoted as DFT and DFT + Ud. The first one corresponds to a standard DFT calculation (GGA-based, without any corrections). The second one corresponds to a typical approach from the literature, in which only Ti atoms are corrected. (48,62,72−77,79−87,89−92) In the DFT and DFT + Ud calculations, all other (i.e., not related to +U +J corrections) parameters were the same as in the DFT + Ud,p calculations to permit a direct comparison. The EOS calculations (Section 3.2) used a 5 × 5 × 5 supercell. The band structure calculations (Section 3.3) used a 7 × 7 × 7 supercell.

3.2. Equation of State

In Table 1, we present structural and mechanical properties of STO obtained with the three considered approaches, calculated by fitting the obtained E(a) dependencies (the total energy vs lattice constant, see Figure 3) with the Birch–Murnaghan EOS.

Figure 3

Figure 3. EOS of the perfect STO, calculated with the three considered approaches. Points correspond to the computational results, while solid lines represent the obtained fits (of the Birch–Murnaghan form). Blue dashed lines represent the adjusted experimental EOSs. To increase readability, the curves corresponding to DFT and DFT + Ud calculations were offset by + 0.1 and +0.2 eV, respectively.

Table 1. Properties of the Perfect STO as Predicted by Various Approaches Considered in This Work, Compared with Other Theoretical Calculations, and Experimental Data (for Room Temperature)a,b
methoda (Å)B0 (GPa)B0EgI(eV)EgD(eV)
this work     
 DFT (PBE)3.966161.9 (0.2)4.37 (0.02)1.661.96
 DFT + Ud3.986154.2 (0.3)4.35 (0.02)1.952.32
 DFT + Ud,p3.957172.1 (0.2)4.25 (0.02)2.923.20
other computational studies     
(a) DFT (all with PBE)     
 GTO3.93a, 3.96b171a,b 1.8b, 1.99a2.1b, 2.35a
 LAPW3.94c,d, 3.95e167e, 169c, 175d4.61c1.63d, 1.97c 
 PW3.94b,f,g, 3.95h168f,r, 169b, 171h 1.8bf, 1.85h2.1b, 2.2f,r
 hybrid3.94i172i 1.74i 
(b) DFT + Ud (all with PBE and PW)     
Ueff = 4 eV3.97j164j4.35j≈2.0j 
U = 5 eV, J = 0.64 eV3.90g, 3.93k, 3.98l  2.39l, 2.68k 
(c) hybrid DFT     
 B3LYP, GTO3.94a, 3.95m179m, 180a 3.57a, 4.38m3.89a
 B3PW, GTO3.89n, 3.90o,p, 3.92b, 3.93q190b 3.4b, 3.6n, 3.63o, 3.64p3.7b, 3.96p
 HSE06, PW3.90f,r192f,r 3.01r, 3.07f,r3.47f,r
 PBE0, GTO3.90b201b 3.8b4.2b
 PBE0, PW3.90b193b 4.0b4.4b
(d) various GW approaches (all based on LDA)     
 QSGW   3.32s 
G0W0   3.40t3.76t
 QSGW80 + SO   3.56u 
 QSGW + LP + 0.8Σ   3.2v3.7v
experimental data3.905w,x178.8 (4.6)y4.31 (0.1)y3.25z3.75z

Values given in brackets denote errors, which for the calculated lattice constants were lower than 10–3 Å.


References: a, (46) b, (49) c, (51) d, (50) e, (168) f, (47) g, (86) h, (48) i, (52) j, (62) k, (83) l, (84) m, (53) n, (56) o, (55) p, (54) q, (125) r, (57) s, (58) t, (59) u, (60) v, (61) w, (9) x, (37) y, (12) z. (31)

Standard DFT significantly underestimates B0, by 9.4% when compared to the experimental result. (12) Interestingly, the application of the +U +J corrections for Ti atoms only (DFT + Ud) does not improve the bulk modulus and even makes it worse (the calculated B0 is by 13.8% lower than the experimental value). In contrast to these two, the result obtained for DFT + Ud,p is lower by only 3.8%. For this approach, the calculated pressure derivative of the bulk modulus (B0) was also found as very close to the experiment being─in fact─within the experimental uncertainty (the same holds for DFT and DFT + Ud calculations). The calculated lattice parameter a differs from the experimental result significantly (by 0.06–0.08 Å), but for DFT + Ud,p, the agreement is─still─the best.
The results obtained from the standard DFT calculations deviate slightly from analogous calculations reported in the literature but are comparable. The calculated lattice constant (3.97 Å) and bulk modulus (162 GPa) are higher by 0.02–0.03 Å and lower by 5–9 GPa than those obtained from PW calculations employing the same XC functional (PBE). Our DFT + Ud results (a = 3.99 Å, B0 = 154 GPa, and B0 = 4.35) are similar to those presented in ref (62) (3.97 Å, 164 GPa, and 4.35, respectively), where a similar Ueff = UJ = 4 eV was employed (3.74 eV in our case). The results obtained from DFT + Ud,p cannot be compared with the literature because no studies to date employed similar methodology.
In Figure 3, we present three calculated EOSs and compare them with the corresponding adjusted experimental EOSs, which were obtained by using (in the Birch–Murnaghan EOS) the experimental values of B0 (178.8 GPa) and B0 (4.31). Their adjustment consisted in using the calculated (for the given approach) equilibrium lattice constant a.
Even a cursory glance reveals an excellent agreement obtained for DFT + Ud,p. The calculated EOS coincides with the (adjusted) experimental curve almost in the entire range (notable deviations are observed only at high compressions). For the two other approaches, the agreement is visibly much worse.
The presented comparison clearly demonstrates that the chosen methodology (DFT + Ud,p) significantly improves the description of energetics of STO, with good accuracy obtained in a wide range of strains. The chosen methodology correctly describes mechanical properties of STO both in the linear (reproduced B0) and non-linear (correct B0) regimes. This gives us confidence that the lattice distortions, which appear around defects (and the related energetic effects) will also be correctly described by DFT + Ud,p.

3.3. Electronic Structure

The examination of the calculated band gaps (Table 1) clearly shows that the application of DFT + Ud,p also significantly improves the description of the band structure around the Fermi level. Standard DFT calculations significantly underestimate the experimental indirect band gap (by more than 1.5 eV) and DFT + Ud improves the band gap only slightly (by 0.3 eV). The application of the +U +J corrections also to O atoms yields a band gap, which agrees well with the experiment. Both calculated band gaps (indirect R → Γ, EgI = 2.92 eV; direct Γ → Γ, EgD = 3.20 eV) compare favorably with the experiment (3.25 and 3.75 eV, (31) respectively) and with calculations based on a higher level of theory, like those employing hybrid functionals (46,47,49,53−57,125) or various GW approaches (58−61) (see Table 1). The calculated valence band width is also in excellent agreement with experimental results of 5.5–6.5 eV: (169,170) 5.5 eV was obtained using DFT + Ud,p, while standard DFT and DFT + Ud calculations predicted 4.5 and 4.6 eV, respectively. DFT + Ud,p also offers an improved description of the dispersion of the valence and conduction bands. To illustrate this, in Figure 4, we present the Kohn–Sham single-particle energies evaluated at high-symmetry points of the FBZ. Visual inspection reveals that the energies obtained from DFT + Ud,p match quasi-particle energies obtained from G0W0 calculations. (59)

Figure 4

Figure 4. Electronic structure of the perfect STO. Single-particle energies of the nine highest valence bands and the lowest CB are represented with horizontal line segments. Four panels correspond to the X, Γ, M, and R points of the FBZ. Each panel presents results obtained with DFT (black, left), DFT + Ud (red), and DFT + Ud,p (green), and also results of G0W0 calculations (59) (blue, right), which serve as a reference. The zero of energy (dashed gray line) was set to the Fermi energy EF.

To explain this success of DFT + Ud,p, we use Figure 5 to present the local density of states (l-DOS) calculated for all three considered approaches. A comparison of characteristics obtained for DFT and DFT + Ud shows that the valence band edge essentially remains unaffected when correcting the Ti 3d subshell only. The additional application of corrections to the O 2p subshell lowers its energy. This re-establishes the possibility of hybridization with Ti 3d (see Figure 5), resulting in a significant band gap opening with DFT + Ud,p. Similar observations were made in ref (95).

Figure 5

Figure 5. Total and l-DOSs of perfect STO as predicted by various approaches considered in this work. To allow direct comparison, all panels use the mid-gap energy of the DFT calculation for 0 eV (dashed gray line). Gaussian broadening of 0.1 eV was used.

Figure 6 presents the total DOS calculated with DFT + Ud,p and its decomposition into local angular-momentum projected densities (lp-DOS). (171) The low-lying three narrow maxima (located at −55.7, −32.6, and −29.6 eV) correspond to Ti 3s, Ti 3p, and Sr 4s semicore states. Two broad groups of features located at energies between −16 and −12 eV correspond to bands formed from O 2s and Sr 4s levels. The valence band, which starts at −5.5 eV, is mainly composed of O 2p states, with a small contribution of Ti 3d states. The first two groups of conduction states (located between 3 and 8 eV) have a Ti 3d–t2g and Ti 3d–eg character. Another group of features─originating from the Sr 4d levels─is found above 8 eV.

Figure 6

Figure 6. Total DOS of perfect STO and its decomposition into local angular momentum-projected densities (red─s, green─p, and blue─d). The zero of the energy corresponds to the Fermi energy EF. Two vertical dashed lines indicate the VBM and CBM.

The above results are consistent with G0W0 calculations, (59) which─similarly─located the Ti 3p state at −32.5 eV below EF and yielded a similar position of the O 2p valence band center. Some quantitative differences are─however─visible in the locations of other bands, which in our picture are positioned at higher energies (see Table 2). The same conclusions can be drawn from a comparison with recent X-ray photoelectron spectroscopy (XPS) data, (172) which also provided locations similar to ours. Interestingly, while the position of the Sr 4s state was not identified in ref (172), our DFT + Ud,p results suggest that this level may lie slightly (3 eV) above the Ti 3p level. This could explain the asymmetry of the peak assigned to Ti 3p in the referenced photoelectron spectroscopy experiment.
Table 2. Electronic Structure of Perfect STOa
 Ti 3sTi 3pSr 4sO 2sSr 4pO 2p
DFT + Ud,p–
G0W0 0016.018.129.8
XPS–24.80 15.718.3≈31.7

For DFT + Ud,p (this work) and G0W0 (ref (59)), weighted centers of various features are presented (calculated from the lp-DOS). All values are given in eV and are measured relative to the location of Ti 3p state (the sharpest feature) to allow comparison with the XPS experiment of ref (172).

The presented validation clearly shows that DFT + Ud,p significantly improves the description of STO. The calculated measures characterizing its mechanical properties and the electronic structure agree with recent experiments. They are also consistent with recent theoretical calculations based on a higher level of theory, and this can be obtained within DFT + Ud,p for a fraction of their cost. The well-reproduced energetics and good reproduction of valence band-related features (character, width, dispersion, forbidden gaps) give us confidence that the defect-related properties (such as formation energies, transition levels, location of defect levels) will also be correctly described by DFT + Ud,p.
One more issue deserves comment here, namely our choice not to apply the +U +J corrections to Sr 4d states. These states are located high in energy, mostly form the upper CBs, and contribute only marginally to the VB and lower CBs (see Figures 5 and 6). Therefore, it can be expected that the Sr 4d states (and─consequently─their correcting) will have negligible influence on the investigated ground-state properties. We verified this by performing DFT + Ud,p calculation in which the Sr 4d states were also corrected (we used USr = 0.934 eV and JSr = 0.152 eV, which were obtained from the MT-LR method). These calculations showed that the additional correcting of the Sr atoms does not appreciably influence the results. For example, the calculated direct band gap width increased by only 4 meV.

4. Results─Defective Systems

Jump To

4.1. Structural Relaxation

We begin by discussing the structural relaxations, which proved to be computationally very demanding. In the most difficult cases, we found it necessary to perform up to 15 cycles of DI + BFGS optimizations, with up to 1000–1600 updates of the geometry required to reach the assumed convergence criterion (energy change below 3 meV/2 cycles). One relaxation took between 6 and 9 months of continuous calculations on 16 nodes (equipped with two 14 core processors each). The total wall time of all carried relaxations exceeded 70 M CPU-hours.
The difficulty originates from the very high dimensionality of the considered problem (1869 ± 3 degrees of freedom). However, to a large degree, it is also a consequence of the problem’s specificity, that is, the fact that the introduction of a defect leads to significant rotations and deformations of TiO6 octahedra, and to notable displacements of the nearest Sr ions. These effects feed back into one another, propagating through the crystal lattice, and significantly affecting the system’s energy.
To better illustrate the above problem, we use Figures 7 and 8 to present the evolution of various measures during the relaxation of a VO0,S=0 defect. The data corresponds to the structure whose minimization started from perfect (undistorted) positions. For VO0,S=0, this initial structure led to a lower final energy (however, this is not always the case, see later discussion).

Figure 7

Figure 7. Structural relaxation of the system containing a neutral oxygen vacancy VO0,S=0. Panel (a) presents the evolution of the total energy (relative to the converged value). Note that in the right part (i.e., for energy evaluations >500), the vertical axis was adjusted to better show detail. Two insets present the evolution of forces and displacements (the latter relative to the previous iteration). Panel (b) presents the evolution of the components of the total energy. These are also measured with respect to their converged values. Blue and red identify DI- and BFGS-based minimizations.

Figure 8

Figure 8. Structural relaxation of the neutral oxygen vacancy VO0,S=0 (cont’d). Configurations obtained after the first minimization cycle (a–c) and at convergence (d–f). All panels present a cross-sectional view through the supercell: the central TiO2 (001) plane containing the defect (located in the center) is shown, and (for a,b and d,e) SrO (001) plane located above. Panels (b,e) present atomic displacements, measured relative to the converged structure (b) or the initial structure (e). Panels (c,f) visualize two electronic states introduced by the defect, using their single-particle wavefunctions ψdKS. Here, isosurfaces |ψdKS|2=102e/Å3 are shown.

The analysis of the total energy, forces, and displacements (Figure 7a) shows that the employed strategy overcame many local barriers. As expected, the energy drop was the highest in the first minimization cycle (2.8 eV). However, the subsequent minimizations reduced the energy further (by approx. 0.2 eV), yielding a more suitable configuration for calculating CTLs or studying ion transport. For such transport, the energy barriers are of the order of tenths of eV and─therefore─very fine minima are required to obtain reliable (not underestimated) barrier heights.
The validity of the carried out relaxation can also be supported with an analysis of the total energy decomposition (Figure 7b), which shows that all contributions to the total energy converged to their final values before the relaxation was complete.
We note that the main features of the electronic states introduced by the defects, such as their spatial character (described by the corresponding single-particle wavefunction, ψdKS) and position within the gap (given by the corresponding single-particle Kohn–Sham eigenvalue, ϵdKS), are already well reproduced after the first minimization cycle. For example, at this point, the two ψdKS states of VO0,S=0 were located ϵCBM–ϵdKS = 0.714 eV below the CB edge, which is very close to the converged result (difference of <2 meV). After the first DI + BFGS cycle, the single-particle wavefunctions ψdKS were also very similar to those at convergence (see Figure 8).
However, a single-cycle relaxation does not provide an accurate picture of the defect’s influence on the crystal structure, with the atomic positions undergoing significant changes during later cycles. For example, for VO0,S=0, the positions obtained after the first cycle differed by as much as 0.4 Å (and by 0.12 Å on average) from those at convergence (see Figure 8). This clearly shows that multiple DI + BFGS cycles are required to obtain a reliable picture of the defect.
The relaxation proved to be comparatively difficult for hydrogen interstitials, although the opposite could be expected due to the smaller size of IH. This can be explained by the fact that in the considered range of distances (limited by half of the supercell side, ≈10 Å), the lattice distortions introduced by IH are comparable to those associated with the presence of VO.
We found the equilibrium structure obtained from the relaxation to depend─to some extent─on the initial configuration. However, as it concerns the defects’ static (i.e., not related to their migration) properties, the differences between the minima corresponding to the same defect state were negligible. Tables 3 (data for VO) and 4 (IH) summarize all the relaxations carried out. The differences in the calculated total energies of the corresponding minima are small (below 0.05–0.15 eV). The locations of the defect levels are also comparable (differences in ϵdKS are below 0.01 eV). The similarity of the spin density distributions shows that the determined ground states are also similar.
Table 3. Minima Identified from the Relaxations Carried outa
defect state   distances (Å)
/ energy (eV)defect level position (eV)integral of |ρ – ρ| (e)ΔdTi–TiΔdSr–Sr
VO2+d+404.640n/a0+0.221+0.144, +0.155
 p+0.039n/a0+0.240+0.132, +0.170
VO1+p+415.791–0.0991.88+0.209+0.145, +0.155
 d+0.155–0.0981.88+0.212+0.143, +0.160
VO0,S=0p+427.045–0.716, −0.7162.39+0.046+0.134, +0.208
 d+0.046–0.719, −0.7182.45+0.045+0.167, +0.186
VO0,S=2d+427.044–0.724, −0.7212.47+0.046+0.165, +0.186
 p+0.015–0.719, −0.7142.45+0.046+0.134, +0.208

Data for various charge/spin states (given in the first column) of the VO defect is presented. The second column indicates if the relaxation started from perfect (p) or distorted (d) positions. The total energies of the identified lowest energy minima (first rows) are measured relative to the perfect system (ED,qEP is given). For other─higher-lying─minima, the energy difference to the corresponding lowest energy minimum is given. Positions of the defect levels are given as ϵdKS–ϵCBM, that is, they are measured with respect to the CBM of the perfect system. The fifth column gives the integral of |ρ(r) – ρ(r)| (spin up and spin down densities) over the entire supercell. The last two columns give the changes in the Ti–Ti and Sr–Sr distances for Ti and Sr atoms closest to the defect. For ΔdSr–Sr, the minimum and maximum is given.

Table 4. Minima Identified from the Relaxations Carried Out (cont’d)a
defect state   distances (Å)angles (degrees)
/ energy (eV)defect level position (eV)integral of |ρ – ρ| (e)ΔdTi–TiΔdSr–SrdO–HθTi–O–HθSr–O–H
IH1-,S=0d–5.288–0.115, −0.1153.67+0.282+0.3340.96377.646.3
 d+0.014–0.108, −0.1043.66+0.279+0.3250.96277.946.2
 d+0.092–0.115, −0.1153.66+0.288+0.3260.95786.244.3
 p+0.105–0.107, −0.1053.65+0.280+0.3230.95686.544.4
IH1-,S=2d–5.298–0.116, −0.1103.75+0.282+0.3340.96377.546.2
 p+0.009–0.107, −0.0993.75+0.280+0.3240.96178.545.9
 d+0.021–0.106, −0.0993.75+0.279+0.3260.96278.046.2
 d+0.045–0.096, −0.0893.73+0.282+0.3310.96278.146.0

Data for the IH defect. Columns 1–6 have the same meaning as in Table 3. Column 7 presents the change in the distance between two Sr atoms closest to the H atom. The last three columns describe the position of the H atom.

For the oxygen vacancy, the measures which characterize lattice distortions that occurred in the very vicinity of the defect (the outward displacements of two closest Ti atoms and four closest Sr atoms) are also comparable between the corresponding minima. The positions of the H atoms are also similar for all minima corresponding to IH0, as well as for IH1+.
For IH1–,S=0, considerable differences are, however, seen in the calculated θTi–O–H valence angle, which is either ≈78 or ≈86°, for the two lowest and two highest energy minima, respectively. Because of this difference, these two pairs of minima can be classified as two distinct motifs. They also differ (by 0.005–0.007 Å) in the calculated dO–H bond length. Similar motifs (with θTi–O–H = 90 and 76°) were identified by T-Thienprasert et al. (136) for the IH1+ defect.
The absence of the low θTi–O–H angle motif for the q = 0 and 1+ charge states of IH does not originate from differences in the used initial configurations. The relaxations performed for each charge state of IH used the same initial atomic positions. The underlying reasons are difficult to identify.
In general, we could not ascertain any pattern that would make one initial structure more suitable for the relaxation of VO. For VO0,S=0 and VO1+, relaxations which started from perfect (undistorted) positions led to a lower final energy. For the remaining states introducing initial distortions proved to be beneficial. This was also the case for IH, for which all identified minimum energy structures corresponded to distorted initial guesses.

4.2. Formation Energies

Table 5 presents the calculated formation energies and the obtained IIC and PA corrections, which were calculated (108,173−175) using Durrant’s procedure (114) (for details, see Supporting Information, Section S3).
Table 5. Calculated Formation Energies ΔH (Uncertainty in Brackets) and the Contributions of the IIC and PA Correctionsa
 formation energycorrections
defect stateΔH (eV)EIIC (eV)EPA (eV)
VO2++1.45 (0.12)[0, 0.009][−0.062, +0.139]
VO1++4.36 (0.11)[0, 0.002][−0.058, +0.017]
VO0,S=0+7.32 (0.02)00
VO0,S=2+7.31 (0.01)00
IH1+–3.56 (0.13)[0, 0.002][−0.050, +0.063]
IH0–0.54 (0.07)00
IH1–,S=0+2.40 (0.11)[0, 0.002][−0.093, +0.068]
IH1–,S=2+2.34 (0.07)[0, 0.002][−0.085, +0.058]

Here, we do not distinguish different minima of the same defect. Instead, the results obtained for the various minima were used to determine the presented uncertainties. These uncertainties also account for the uncertainties of the EIIC and EPA corrections, for which bounds are given. Note that both corrections are zero by definition for the neutral defects.

The IIC corrections proved to be negligible: for both of VO and IH defects, they are in the order of 2–10 meV. The main reason for that is the very high dielectric constant of STO (we assumed 330 (27)). A similar observation─about the practical insignificance of IICs for STO─was made in ref (93).
The obtained PA corrections were proved to be larger than IICs. However, they are still small (typically below ±0.1 eV), and their neglect would not change the obtained picture of the defects’ relative stability. The negligibility of the PA corrections can be explained by referring to the details of the calculations we carried out. The large size of the supercells causes the jellium to have a low charge density and, therefore, it does not appreciably shift the average electrostatic potential.
The formation energies presented in Table 5 correspond to the Fermi level ϵF set to the VBM, and, crucially, oxygen- and hydrogen-rich limits (the H2 and O2 molecules were taken as the reference states, for more details, see Supporting Information, Section S2).
Our DFT + Ud,p calculations reveal that the oxygen vacancy prefers the q = 2+ charge state, which has the lowest formation energy (1.45 eV for the oxygen-rich conditions). This result is in line with experimental observations. The formation energies of VO1+ and VO0 are higher by ca. 2.9 and 5.9 eV, respectively, than the formation energy of VO2+. Keeping the uncertainties in mind, the singlet and triplet states of the neutral VO are characterized by the same formation energy, which equals 7.3 eV.
The hydrogen interstitial strongly prefers the 1+ charge state, which has the lowest formation energy (≈−3.6 eV). Its negative value informs that under the considered conditions (hydrogen-rich conditions, Fermi level set to the VBM), hydrogen will be naturally incorporated into STO in the form of protons, and this also agrees with experimental observations. Although the formation energy of IH0 was also negative (−0.5 eV), its significantly higher value (by ≈3 eV) suggests that in STO, the formation of neutral hydrogen defects will not be as favorable. The formation of negatively charged hydrogen defects is strongly unfavorable (their formation energy is higher by 6 eV than that of IH1+).
We now turn to comparing our results against those of other computational studies. There are many literature reports, which considered the formation of VO using various theoretical approaches. LDA–DFT calculations (120) yielded significantly smaller formation energies: 4.5 eV (VO0), 2.4 eV (VO1+), and 0.3 eV (VO2+). Results in line with our calculations were obtained by Carrasco et al., (122) with ΔH = 7.65 eV for VO0 based on GGA–DFT calculations with a PW91 functional. DFT + U calculations predicted a lower formation energy of VO0 (5.3, (62) 6.1, (86) 6.7, (83) and 6.8 eV (77)), but a similar─to ours─formation energy of VO2+ (1.2 eV (83)). Formation energies obtained from DFT calculations with hybrid functionals are also lower than those reported here. Reference (53) reported 6.8 eV for VO0 (obtained with B3LYP), while Mitra et al. (57) obtained 5.9 eV (VO0), 3.2 eV (VO1+), and 0.7 eV (VO2+) using the HSE functional. A comparison with other existing studies (56,123,125) is difficult, as they used a different reference chemical potential μO.
The above comparison highlights significant discrepancies among the existing reports and our work. However, we found a reasonable agreement with a very recent work of Ricca et al., (93) who used a similar method (DFT + U + V, with self-consistently adjusted U and V parameters). Reference (93) reported ΔH = 6.4 eV for VO0 (7.3 eV in this work), 3.4 eV for VO1+ (here, 4.4 eV), and 0.6 eV for VO2+ (1.4 eV in the present study). While the corresponding formation energies differ (in this study, they are higher by ≈1 eV), the differences between the formation energies of various charge states are similar. We found that the formation energy of a neutral oxygen vacancy is by 2.9 and 5.9 eV higher than that of a single and double charged vacancy, respectively. Reference (93) reported similar differences: 3.0 and 5.8 eV, respectively. This similarity suggests that the apparent shift of 1 eV in the formation energies may originate from the choice of the reference chemical potential μO. Reference (93) also used μO=EO2/2; however, this energy was obtained directly from DFT calculations and not from an adjustment procedure, as in the present work (see Supporting Information, Section S2). This may introduce a difference in the chemical potential μO and, consequently, a uniform shift in the calculated formation energies. We note that ref (93) found that the same formation energy characterizes the singlet and triplet states of the neutral oxygen vacancy, and this agrees with our result.
The formation energy of VO2+ can be compared with experimental data. ΔH (see eq 2) is defined through the following reaction (written using Kröger–Vink notation (176))
that is, a process in which the formation of a doubly ionized oxygen vacancy (denoted with VO••, in our notation VO2+) and gaseous oxygen reduces the hole (h) concentration. This process occurs in the region of p-type conductivity (ϵF close to VBM). Its activation energy (the so-called oxidation enthalpy) was measured by Moos and Härdtl (29) to be ΔHOx ≈ 1 eV (see Section 4.7 therein). A similar value was measured by Choi and Tuller, (177) who obtained ΔHOx = 1.29 eV. Our result, ΔH(q = 2+,ϵF = 0) = (1.45 ± 0.12) eV, agrees very well with both experimental measurements.
Similarly, the calculated ΔH(q = 2+,ϵFEg) can also be compared with the experiment, as it corresponds to the activation energy (the so-called reduction enthalpy, ΔHRed) of the following reaction
This reaction, in which two excess electrons are emitted into the CB of the host material, is seen in the region of n-type conductivity (for high ϵF). Moos and Härdtl (29) measured ΔHRed = 6.1 eV, which compares well with our calculated ΔH(q = 2+,ϵFEg) ≈ 7.3 eV (see Figure 9 for ΔH vs ϵF dependencies). The observed discrepancy of ≈1 eV might originate from the difference in the Fermi level ϵF (its exact value was not specified in ref (29)).

Figure 9

Figure 9. Dependence of the formation energy ΔH on the Fermi level ϵF. Two panels present the result obtained for the oxygen vacancy (a) and hydrogen interstitial (b). Individual lines correspond to various considered charge states q. The error bars drawn at ϵF = 0, 1.25, and 2.5 eV depict the uncertainties of ΔH. Here, no distinction is made between singlet and triplet states of VO0, which gave similar results. The presented data encompasses both cases. The same holds for IH1–.

The fact that the calculated formation energies ΔH are higher than the experimental results should be expected. This can be explained by accounting for the limitations of the modeling technique used. The application of the electrostatic corrections means the obtained ΔH is unaffected by spurious electrostatic interactions. However, the calculated total energies still suffer from other finite-size effects. The application of the supercell approach limits the possible long-range relaxation of the structure around the defect, making the obtained total energies higher than those which would be obtained for larger systems (for demonstration, see, e.g., refs (86and120)) and, consequently, in the dilute defect limit. Thus, we expect the formation energies calculated for finite systems to be higher than those measured in the experiment, as we found in this work for VO2+.
Many of the previous theoretical studies obtained formation energies that were lower than those measured experimentally. In the context of the above discussion, this suggests that the theoretical approach used in these studies must suffer from one or more deficiencies, as it would give even more underestimated formation energies in the dilute defect limit.
Our results obtained for the hydrogen interstitial disagree with those presented by Varley et al., (178) who obtained ΔH = −2.00 eV (for IH1+), 1.59 eV (for IH0), and 5.55 eV (for IH1–), based on DFT calculations with hybrid functionals. Their formation energies are higher by ≈1.6 eV (IH1+), 2.1 eV (IH0), and 3.2 eV (IH1–) than the formation energies obtained in this work. This discrepancy can be explained by the observation that ref (178) used a smaller supercell (135 atoms), which strongly limits relaxation. A comparison with other studies which considered IH (72,73,134,136) is difficult, as the formation energies were either not calculated or their presentation omits important details (like the chosen reference state).

4.3. Charge Transition Levels

Figure 9 presents the dependence of the formation energy ΔH on the position of the Fermi level ϵF. For VO, the q = 2+ charge state is preferred in almost the entire range of ϵF. VO2+ thermodynamically transitions to states with lower nominal charge q (1+ and 0) when the Fermi level is ≈0.1 eV below the CBM (location of the (2+/1+) and (2+/0) TTLs). Keeping in mind the uncertainties, we cannot say if VO2+ transitions directly into VO0 (this would suggest a negative-U behavior, a similar observation was made in ref (120)), or through a VO1+ charge state. All three ΔH versus ϵF lines presented in Figure 9 cross at ϵF = 2.95 eV.
For the Fermi level close to the CBM, the formation energies of all three charge states of VO are very similar. This suggests that neutral, single-charged, and double-charged oxygen vacancies may coexist under the corresponding conditions.
The band gap Eg = 3.06 eV, as shown in Figure 9, corresponds to the result obtained for a 5 × 5 × 5 supercell (the same system size was used for defect calculations). Calculations performed for a larger 7 × 7 × 7 supercell gave a slightly lower band gap Eg = 2.92 eV (see Table 1, the main reason for this is the better sampling of reciprocal space), which is below the location of the transition levels observed for VO. This would mean that for VO in cubic STO, there are no TTLs in the band gap. A similar observation was made by Choi et al. (82) based on DFT + U calculations.
Figure 9b shows that the hydrogen interstitial prefers a q = 1+ charge state over almost the entire range of ϵF. At ϵF = 2.97 eV, it thermodynamically transitions to IH1–. In this case, it also cannot be concluded if a direct transition takes place, or if the intermediate q = 0 state is involved.
Similarly to VO, all hydrogen interstitials have comparable formation energies for ϵFEg. This suggests that they might coexist under the corresponding conditions.
Table 6 presents the calculated TTLs and OTLs (see eqs 3 and 4 for their definition), which are also schematically illustrated in Figure 10. VO introduces both shallow and deep OTLs, which lie 0.08 eV ((1+/2+) level) and 1.0–1.1 eV ((0/1+) and (0/2+) levels) below the CBM when the excitations are considered. The OTLs related to the excitation of the hydrogen interstitial are shallow, located ≈0.1 eV below the CB edge.

Figure 10

Figure 10. Thermodynamic ϵth(q/q′) and optical ϵop(q/q′) CTLs. The (q/q′) labels describe the transition type. For absorption, the levels correspond to the opposite transition (q′/q). The Kohn–Sham eigenvalues ϵdKS+ΔϵdKS (corrected values) are also presented. Black error bars represent uncertainties. For VO0 and IH1– no distinction is made between singlet and triplet states, for which results were similar.

Table 6. Thermodynamic and Optical Transition Levelsa
 charge transition level (eV)
transition typethermodynamicabsorptionemission
oxygen vacancy
(2+/1+)–0.11 (0.20)–0.08 (0.08)–0.09 (0.10)
(1+/0, S = 0)–0.11 (0.13)–0.98 (0.06)–0.02 (0.04)
(1+/0, S = 2)–0.12 (0.12)–0.98 (0.06)–0.02 (0.04)
(2+/0, S = 0)–0.11 (0.05)–1.09 (0.09)–0.07 (0.09)
(2+/0, S = 2)–0.12 (0.05)–1.09 (0.09)–0.08 (0.08)
hydrogen interstitial
(1+/0)–0.04 (0.20)–0.11 (0.06)–0.04 (0.06)
(0/1 −, S = 0)–0.12 (0.18)–0.07 (0.08)–0.02 (0.07)
(0/1–, S = 2)–0.18 (0.14)–0.08 (0.07)–0.03 (0.06)
(1+/1–, S = 0)–0.08 (0.12)–0.08 (0.07)–0.03 (0.06)
(1+/1–, S = 2)–0.11 (0.10)–0.07 (0.07)–0.04 (0.06)

Locations of all levels are measured with respect to CBM. A negative value informs that the level is located below CB. The value given in brackets is the uncertainty. Note that the OTL presented in column 3 corresponds to the opposite transition q′ → q.

For the neutral oxygen vacancy, significant differences are seen between OTLs related to excitation and those corresponding to recombination. While the (0/1+) and (0/2+) excitations of VO0 produce deep features (1.0–1.1 eV below CBM), the recombination of electrons from the CB to the states introduced by VO0 results in features that are very close to the CB edge (less than 0.1 eV below CBM, for all recombinations to VO0). This observation explains the significant differences in the experimental studies concerning oxygen vacancies in STO.
Absorption techniques reported deep locations of defect levels associated with VO. Ultraviolet photoemission spectroscopy measurements identified a feature located between 1.0 (14) and 1.2 eV (15) below the CB, which agrees almost perfectly with our calculated (0/1+) and (0/2+) transition levels (1.0–1.1 eV below CBM). On the other hand, spectroscopy based on emission analysis reported shallow locations. For example, Kan et al., (16) based on photoluminescence (PL) measurements, observed a state which is located slightly (0.1 eV) below the CB. This agrees very well with our calculated (2+/1+), (1+/0), and (2+/0) OTLs. We note that recent PL studies also reported another level (located ≈0.4 eV below the CB). However, this feature was attributed (16,18) to exciton recombination not considered in this work.
The calculated TTLs also agree well with conduction measurements of donor ionization energies, which all predicted a very shallow location (down to 0.16 eV below the CBM (25,26)) of the TTLs of VO.
Based on the above comparison, we conclude that the DFT + Ud,p results obtained for VO are in very close agreement with the experiment. A lack of experimental data prevents comparing the OTLs obtained for the hydrogen interstitials. However, the calculated TTLs are qualitatively consistent with recent experimental studies, (32,34) which observed hydride conduction in STO subjected to reducing conditions and high temperatures. These experiments suggested that the (1+/1−) TTL of IH is located very close to the CB edge, which agrees very well with our result.
For IH, the obtained values of ϵdKS (≈0.1 eV below CBM) correctly reflect the positioning of the TTLs and OTLs. For VO1+, the calculated ϵdKS also agrees reasonably with the obtained TTLs and OTLs. However, the ϵdKS eigenvalues obtained for VO0 (≈0.75 eV below CBM) differ strongly from the calculated TTLs (0.1 eV below CBM) and also do not reflect the positioning of the corresponding OTLs (1.0–1.1 eV below CBM for excitations of VO0 and 0.1 eV for recombinations to VO0). The obtained results show that─in general─the ϵdKS Kohn–Sham eigenvalues are not accurate measures of the real locations of the defects’ levels. Their correct estimation requires the calculation of the CTLs.

4.4. Defect Levels and Electronic Structure

The following analysis will be based on Kohn–Sham single-particle levels (ψiKS wavefunctions) and the corresponding single-particle energies (ϵiKS eigenvalues). Although both measures─as argued in the previous section─have some limitations, they are still helpful in the qualitative description of how the considered defects influence the electronic structure of STO. The presented analysis will also show why the Kohn–Sham eigenvalues ϵdKS failed in describing the VO0 defect.
Figure 11a visualizes the two levels introduced in the forbidden gap by VO0,S=0. Each of these is localized on one neighboring Ti atom (Ti1 or Ti2), being mainly composed of its 3d–t2g level (with a small addition of the O 2p levels of the nearby O atoms, see the DOS in Figure 12f). Because of this trapping, the Ti1 and Ti2 atoms adopt a 3+ oxidation state. This change is visible in the DOS and manifests itself in a 1–2.5 eV increase of the Kohn–Sham eigenvalues corresponding to the 3s and 3p levels of the Ti1 and Ti2 atoms (see Figure 12, panels b and c).

Figure 11

Figure 11. Levels introduced in the forbidden gap by the VO (a–c) and IH (d–f) defects. All panels show a cross-sectional view through the supercell (the central TiO2 plane is shown, the defect is located in the center) and visualize isosurfaces corresponding to an electronic density of |ψdKS|2=102e/Å3. In the interest of readability, planes containing lobes of ψdKS were framed with color borders. For the triplet states (VO0,S=2 and IH1–,S=2), the two spin majority levels are presented separately.

Figure 12

Figure 12. Influence of the VO and IH defects on the electronic structure of STO. Panel (a) (bottom) presents the total DOS of perfect STO and serves as a guide. It indicates the position of regions shown in panels (b–f), which present features introduced by VO0,S=0 (b,c,f) and IH1+ (d,e) defects in the low-energy parts (b–e) and within the forbidden gap (f). Panels (b–f) present the local (atom-decomposed) DOS. All panels use the same reference energy (EF of perfect STO).

Similar electronic behavior is observed for VO0,S=2, which also introduces two levels with the t2g character in the forbidden gap. In this case, the decrease of the oxidation state is also visible in the lower energy parts of the DOS (not shown here). It also manifests itself in the outward displacements of the Ti1 and Ti2 atoms, which for VO0,S=2 (and also VO0,S=0) are by ≈0.2 Å lower than for the other charge states of VO (see Table 3).
The levels introduced by VO0,S=2 (see Figure 11b) differ in symmetry from those introduced by VO0,S=0. For the triplet state, each level is composed of two orthogonal t2g levels of the Ti1 and Ti2 atoms, such that the lobes of each level are located around two perpendicular planes. Two levels of VO0,S=2 differ in the sign of their single-particle wavefunctions.
The above analysis explains why the Kohn–Sham eigenvalues ϵdKS are inadequate for predicting the CTLs of VO0. During the charge transition from (and to) the neutral VO, a change in the oxidation state of the Ti1 and Ti2 atoms manifests. This change is reflected in a considerable change of the Kohn–Sham eigenvalues of the deep levels (Ti 3s and 3p). Because the ϵdKS eigenvalues obtained for VO0 do not account for this change in the low-energy eigenvalues, they cannot correctly describe the total energy change in the transitions to (and from) VO0. From the above analysis, a general conclusion can be drawn: single (i.e., corresponding to one charge state) Kohn–Sham eigenvalues ϵdKS are insufficient for predicting the locations of CTLs, which involve significant, that is, chemical, changes in the electronic structure.
In contrast to VO0, VO1+ introduces a shallow state, which is delocalized. This level is formed from the 3d–t2g levels of multiple Ti atoms belonging to the same TiO2 (010) plane (see Figure 11c). The participation of the 2p levels of the O atoms belonging to the same TiO2 plane is also visible. For VO1+, there is no change in the oxidation state of the Ti1 and Ti2 atoms, as evidenced from the DOS (not shown here) and an analysis of the Ti1–Ti2 distance.
Similarly to VO1+, the levels introduced by IH0 and IH1– are shallow, delocalized, and also formed from the 3d–t2g and 2p levels of the Ti and O atoms belonging to the same TiO2 plane (see Figure 11, panels d–f).
VO2+ does not introduce any levels in the forbidden gap and no other new features in the DOS. A similar situation (no levels within the forbidden gap) is observed for IH1+, which, however, manifests its presence in the lower energy parts of the DOS.
For all IH defects, the addition of the H atom leads to new features appearing in the DOS below the O 2s and 2p valence bands (see Figure 12, panels d and e). These features are also composed of the H 1s level and should be attributed to the formation of an OH bond (and accompanying energy shifts). These new levels lie 2.31 eV below the O 2s band and 1.88 and 0.33 eV below the O 2p band. These locations are─within the obtained uncertainty ±0.04 eV─independent of the charge state q. The given locations correspond to the Kohn–Sham eigenvalues and therefore should be considered only as estimates. A more accurate evaluation of the real positioning of these levels (energies which might be observed in the optical spectroscopy and which correspond to excitations from the identified OH levels) is not possible within the applied approach, requiring methods suitable for studying excited states.
The analysis of these (O 2s band- and O 2p band-derived, OH bond-related) states shows that they are predominantly formed by the states provided by the bonded O and H atoms. However, to a large degree they are also composed of the states provided by the nearest Ti atoms (see Figure 13). This reveals the importance of the Ti sublattices for protonic transport. With this in mind, future research could focus on verifying how the OH bonding could be affected by doping STO with other cations. Possibilities include enhancing proton transport by weakening the OH bond via electron drain to a more electronegative dopant.

Figure 13

Figure 13. Influence of the IH1+ defect on the electronic structure of STO. All panels present a cross-sectional view through the supercell, visualizing atoms from the central TiO2 plane (with the H atom in the center). Panel (a) serves as a guide, showing the location of the region visualized in panels (b–d). These panels present the Kohn–Sham single-particle wavefunctions ψiKS (|ψiKS|2=102e/Å3 isosurfaces are shown) of deep lying states, found below O 2p (b,c) and O 2s (d) bands. Only spin up ψiKS are presented.

The picture obtained in this work for VO0 disagrees (in many aspects) with previous theoretical studies. Their predictions are schematically summarized in Figure 14 (scenarios a–d), which also presents the findings of this work (Figure 14e).

Figure 14

Figure 14. Defect levels introduced by the neutral oxygen vacancy VO0 in the forbidden gap, as predicted in the literature (a–d). Scenario (e) corresponds to the picture obtained in this work, in which the Kohn–Sham eigenvalues (ϵdKS) still─as in previous works─disagree with the experiment (f). However, the calculated CTLs agree well with the available experimental data and explain the differences between the spectroscopic and conduction measurements done for VO0. The presented schematic was partly (scenarios a–d) redrawn after ref (93) (Figure 2 therein): C. Ricca, I. Timrov, M. Cococcioni, N. Marzari, U. Aschauer, Self-consistent DFT + U + V study of oxygen vacancies in SrTiO3, Phys. Rev. Res. 2, 023313, 2020, DOI:; licensed under the Creative Commons Attribution 4.0 International license. (a) References (80and119–121), (b) refs (53,56,57,122,123,125,and126), (c) refs (57,80–82,93,and130), (d) refs (81,82,and93), (f) refs (14–16,25,and26).

Standard DFT calculations and DFT + U studies, which used low U values, (80,119−121) predicted that VO0 introduces two delocalized t2g levels, which are positioned close to or in the CB (scenario a). This would mean that the q = 0 charge state is unstable. This failure can be attributed to a significant underestimation of the STO band gap. Hybrid DFT calculations (53,56,57,122,123,125,126) predicted deep levels for VO0 but found that they have a Ti 3d–eg character. DFT + U studies that used larger U values (yielding a reasonable band gap) predicted deep localized levels with an eg character (57,80−82,130) (scenario c, obtained for the singlet state) or positioned one of the defect levels in the CB (81,82) (scenario d, obtained for the triplet state). Our results also disagree with those obtained recently based on DFT + U + V calculations by Ricca et al. (93) They also found that VO0,S=0 introduces two deep localized levels but with an eg character (scenario c). For VO0,S=2, two distinct levels were identified in ref (93) : one deep, localized level (with an eg character) and one delocalized level (with t2g character) located within the CB (scenario d). This also disagrees with our result.
To the best of our knowledge, none of the theoretical studies carried out to date reported defect levels with characteristics similar to that found in this work. We do not set out to identify the reason or reasons for the observed discrepancies (like differences in the supercell size, accounting for structural relaxation, the type and the quality of the basis set used, the XC functional chosen, the details related to the Hubbard + U and Hund’s + J corrections, and others), only noting that some efforts toward establishing these were undertaken in ref (93). Instead, we emphasize that the results obtained in this work are supported by the agreement of the calculated CTLs with those measured experimentally (14−16,25,26) (compare scenario e─our work─with scenario f─experiment).

4.5. Atomic Charges

We assessed the influence of the VO and IH defects on the bonding of STO using natural population analysis (NPA). (179,180) The NPA atomic charges calculated for the perfect system are +1.75 for Sr, +2.12 for Ti, and −1.31 and −1.25 for equatorial and apical O atoms, respectively. The atomic charge of Sr is very close to its oxidation state (2+), showing the ionic character of its bonding. In contrast to this, the atomic charges on Ti and O differ strongly from their oxidation states (4+ and 2–, respectively), revealing a degree of covalency in their bonding.
The atomic charges obtained for the perfect system were used as reference charges in the analysis of charge redistribution in the defective systems. For this purpose, we calculated the acquired charges, that is, the differences
of the NPA atomic charges corresponding to the defective (D, q) and perfect (P) systems, with i indexing atoms. All the studied defective systems were analyzed in this manner.
Figure 15 shows the dependence of the acquired charge ΔqiNPA on the distance from the defect rid = |Rird| (Ri denotes the position of atom i and rd represents the position of the defect). Results for two example systems are presented: the lowest energy structures of VO2+ and IH1+. Polarization of the crystal structure is apparent. The introduction of either defect causes the Ti and O atoms to acquire additional charge─mean values (averaged over all Ti/O atoms) are given: +0.10 for Ti and −0.034 for O (for VO2+) and +0.14 for Ti and −0.050 for O (for IH1+). The Sr atoms, in contrast, acquire only negligible charge.

Figure 15

Figure 15. NPA carried out for the VO2+ (a) and IH1+ (b) defects, and all considered systems (c). Panels (a,b) show how the charges acquired by individual atoms (ΔqiNPA) depend on their separation from the defect (rid). Here, each point represents a single atom, dashed lines represent mean acquired charges of Sr, Ti, and O atoms, and the shaded regions depict the corresponding standard deviations. Panel (c) summarizes results obtained for all identified minima (the ordering of minima is the same as in Tables 3 and 4 if read from left to right). Here, the points represent the mean acquired charges of Sr, Ti, and O atoms, and the error bars depict standard deviations [corresponding to the regions highlighted in panels (a,b)].

Although the observed polarization is modest (for Ti and O atoms the mean acquired charges constitute 5–7% and 2–4% of their reference charges, respectively), it extends over the entire system and─in the studied systems─does not diminish appreciably even at the largest considered distances rid (limited by the size of the employed supercell). Our calculations do not set out to determine how localized this effect is.
The ΔqiNPA(rid) characteristics obtained for the other considered systems (other charge states and other minima) are qualitatively similar to those presented. However, some quantitative differences were identified. These are shown in Figure 15c, which summarizes all the population analyses we carried out.
The polarization identified above occurs for all systems and is manifested through the charging of the Ti and O atoms but not Sr atoms. The fact that the polarization is observed for the completely ionized (VO2+ and IH1+) defects means that it is not related to the redistribution of electrons from the defect levels. It is also not a consequence of a non-zero defect charge because it is also observed for VO0 and IH0. These two observations lead us to attributing the observed polarization to the removal/addition of the defect’s atom.
As seen from Figure 15c, the magnitude of the polarization varies between systems and can even be different for systems containing the same type of defect (compare results obtained for various minima of VO1+, VO0,S=0, IH1+, and IH0). This is illustrated in Tables 7 and 8, where we present the atomic charges of the atoms close to the defect. These are the two nearby Ti atoms (Ti1 and Ti2), the four closest Sr atoms, and the 10 closest O atoms (neighbors of Ti1 and Ti2), as well as (only for IH defects) the O and H atoms forming the OH group.
Table 7. Atomic Charges for Atoms Close to the VO Defect. Data Corresponding to Various Charge/Spin States (First Column) and Various Identified Minima (Consecutive Rows) Is Presenteda
 acquired charge ΔqiNPA
defect stateTi1Ti2O (10 atoms)Sr (4 atoms)
VO2++0.089+0.243–0.069 (0.024)–0.010 (0.002)
 +0.092+0.253–0.067 (0.032)–0.003 (0.002)
VO1++0.163+0.288–0.088 (0.029)–0.005 (0.004)
 +0.012+0.132–0.034 (0.023)–0.021 (0.002)
VO0,S=0–0.077+0.001–0.145 (0.012)+0.012 (0.001)
 –0.187–0.065–0.099 (0.026)–0.008 (0.002)
VO0,S=2–0.125–0.056–0.115 (0.022)–0.015 (0.002)
 –0.130–0.023–0.127 (0.018)–0.008 (0.002)

For the O and Sr atoms, the mean value (averaged over the specified number of atoms) and its standard deviation (in brackets) are given. They specify the charge acquired by a single Sr or O atom. The ordering of the minima is the same as in Table 3.

Table 8. Atomic Charges for Atoms Close to the IH Defecta
 acquired charge ΔqiNPAatomic charge qiNPA
defect stateTi1Ti2O (10 atoms)Sr (4 atoms)OH
IH1++0.197+0.202–0.059 (0.018)+0.013 (0.011)–1.152+0.340
 +0.142+0.141–0.044 (0.014)+0.006 (0.012)–1.207+0.403
 +0.061+0.056–0.008 (0.015)+0.006 (0.012)–1.209+0.427
 +0.175+0.181–0.050 (0.016)+0.019 (0.009)–1.218+0.408
IH0+0.071+0.066–0.008 (0.016)+0.002 (0.006)–1.197+0.419
 +0.057+0.075–0.014 (0.017)+0.004 (0.008)–1.230+0.437
 +0.037+0.036–0.010 (0.019)+0.002 (0.005)–1.207+0.429
 +0.053+0.058–0.009 (0.017)–0.001 (0.013)–1.199+0.416
IH1–,S=0+0.076+0.126–0.020 (0.011)–0.005 (0.010)–1.132+0.343
 +0.083+0.053–0.014 (0.009)–0.008 (0.011)–1.140+0.362
 +0.015+0.015+0.001 (0.019)–0.004 (0.009)–1.197+0.417
 +0.012+0.027–0.002 (0.015)–0.006 (0.009)–1.217+0.437
IH1–,S=2+0.016+0.029–0.002 (0.017)–0.007 (0.010)–1.132+0.366
 +0.048+0.084–0.016 (0.013)–0.005 (0.010)–1.192+0.393
 +0.003+0.008+0.001 (0.017)–0.006 (0.009)–1.188+0.409
 +0.070+0.104–0.020 (0.013)–0.006 (0.010)–1.155+0.375

The meaning of columns 1–5 is the same as in Table 7. The last two columns present the atomic charges of the O and H atoms from the OH group. The ordering of the minima is the same as in Table 4.

The most significant differences are observed for VO1+, VO0,S=0, and IH1+ defects, where the charge acquired by Ti1 and Ti2 atoms differs (between the identified minima) by as much as 0.15 (VO1+), 0.11 (VO0,S=0), and 0.14 (IH1+). They are also visible in the charges acquired by the nearest O atoms, where the total (summed) atomic charge may differ even by 0.54 (VO1+), 0.46 (VO0,S=0), and 0.51 (IH1+). The strongly differing minima identified here correspond to those which visibly differ in the mean acquired charges (cf. Figure 15c). This reveals the connection between the polarization observed in the very vicinity of the defect and the entire structure.
The charges of the H and O atoms from the OH group (see Table 8) are also─to some degree─sensitive to the atomic positions (they differ slightly between the identified minima). By itself, the hydrogen atom carries a charge, which ranges from +0.34 to +0.44. The remaining charge of the IH defect can be either localized or delocalized.
The localized charge is also accommodated on the nearby O atom (from the OH group), whose atomic charge increases by 0.08–0.16 after introducing IH. The delocalized charge is distributed in the entire supercell, as either a polarization charge (see Figure 15) or in the form of delocalized defect levels (for IH0 and IH1–, see Figure 11).
Population analysis reveals that the IH defect carries a localized +0.5 e charge (here, ΔqiNPA of the O atom and qiNPA of the H atom were summed), which seems to be independent of the formal charge q of the IH defect. A similar conclusion was made by Bork et al., (134) who, by combining DFT calculations with Bader charge analysis, found the same (+0.5) charge for IH1+ and IH0.

4.6. Crystalline Structure

The influence of the defects on the crystalline structure was quantified in terms of interatomic distances, with the aid of the pair correlation function (PCF), defined as
Here, N(r) represents the number of pairs of atoms (of any type) that are located within the distance [r, r + Δr), while V and N stand for the volume and the total number of atoms in the system, respectively. The PCF relates the local atomic number density at a distance r from any atom to the average atomic number density of the entire system, N/V.
In Figure 16, we present the total PCFs of all considered defects and decompositions of PCFs obtained for VO2+ and IH1+ into partial functions, corresponding to various pairs of atom types.

Figure 16

Figure 16. Influence of the VO and IH defects on the short- and medium-range order of STO. Panels (a,b) compare the total PCFs g(r) obtained for various charge/spin states of VO (a) and IH (b). Panels (c,d) present the decompositions of the total g(r) obtained for VO2+ and IH1+ into partial PCFs, corresponding to various pairs of atom types. The vertical dashed lines correspond to distances in the perfect system (described at the top). Note that all panels use vertical offsets. Their values can be directly read from the labeled tick-marks.

The similarity of all the presented PCFs shows that in the considered short and medium ranges, the VO and IH defects have a similar effect on the crystalline structure of STO. The degree of this influence is comparable among various charge/spin states.
As already mentioned in Section 4.1 (see also Figure 8), the introduction of defects mainly leads to rotations and distortions of the TiO6 octahedra. Changes of another character, like the displacements of the Sr and Ti ions, are only visible in the very vicinity of the defect. This picture allows explaining the obtained PCFs.
The rotations of the TiO6 octahedra cause a significant change in the O–O, Ti–O, and Sr–O distances. These changes are well visible in all PCFs and manifest through the broadening of the corresponding peaks: second (O–O and Sr–O distances), fourth (O–O distances), fifth (Ti–O distances), and sixth (O–O and Sr–O distances). In contrast, no significant broadening (or shift) of the first peak is observed. This shows that the shortest (intra-octahedral) Ti–O distance remained constant, on average, in defect-containing systems.
The sharp fourth peak (corresponding to the shortest Sr–Sr and Ti–Ti distances) indicates no significant changes in either the Sr or the Ti sublattice. Their relative displacement is also not observed (sharp third peak, Sr–Ti distances).
Figure 17 reveals that in the defect containing systems the valence angles θSr–Sr–Sr and θTi–Ti–Ti (formed by a triple of the nearest Sr- or Ti-atoms) are close to 90 and 180°. This shows that the Sr and Ti sublattices remained regular. The θTi–Sr–Ti angles also are similar to those seen in the perfect system, demonstrating no relative displacement of the Sr and Ti sublattices (except for the very vicinity of the defect).

Figure 17

Figure 17. Influence of the VO and IH defects on the short- and medium-range order of STO (cont’d). The two panels present distributions of various valence angles (described within the plots) for VO2+ (a) and IH1+ (b) defects. Distributions were normalized to values given in the round brackets. The vertical dashed lines correspond to angles in the perfect system. Note that all panels use vertical offsets. Their values can be directly read from the labeled tick-marks.

The distributions of all valence angles incorporating O atoms display significant broadening, up to ±10°. This, again, shows the rotations (changes visible in the θO–Sr–O angles) and distortions (θO–Ti–O and θO–O–O angles) of the TiO6 octahedra.
There is one more valence angle, which is especially important in the context of proton conductivity, namely, the θTi–O–Ti angle. As shown in Figure 18e, this angle (spanned by two Ti atoms belonging to two neighboring octahedra, and the O atom being their common neighbor) can be used to measure the relative tilting of the TiO6 octahedra, with the tilt angle defined as
Its importance originates from the fact that the octahedral tilting is suspected to reduce the activation energy for proton migration, enhancing proton conductivity. (181,182)

Figure 18

Figure 18. Defect-induced octahedral tilting in STO. Panels (a,b) present distributions of the tilt angle θtilt calculated for all considered defects (labeled in the plot). Panels (c,d) show the dependence of the tilt angle on the distance from the defect for two exemplary systems (minimum energy structures of VO2+ and IH1+). Here, each point represents a single O atom, presenting its distance from the defect and the value of the tilt angle θtilt centered on it. Panel (e) illustrates the meaning of θtilt.

Figure 18 demonstrates that the introduction of VO and IH defects alike results in a well pronounced octahedral tilting of ≈7–8° on average (a similar mean was found for all systems), which may reach even 20–30° close to the defect. This tilting effect does not decay with the separation from the defect. However, it can be suspected that the limited size of the simulated systems (≈20 Å supercells), as well as the use of PBCs, may have some influence on this non-vanishing behavior. The similarity of the θtilt(rid) and ΔqiNPA(rid) (see Figure 15) dependencies shows that the octahedral tilting and the observed polarization of the TiO6 octahedra are related.

5. Conclusions

Jump To

5.1. Summary

This work investigated oxygen vacancy and hydrogen interstitial defects in cubic STO. A new approach to describing this difficult system was developed within the framework of linear-scaling Kohn–Sham DFT. The band gap problem was overcome by applying Hubbard + U and Hund’s + J corrections to both Ti (3d subshell) and O (2p subshell) atoms, following the DFT + Ud,p method recently proposed by O’Regan et al. (95)
The U and J parameters were determined on theoretical grounds, using the MT-LR approach (69) and later validated for perfect STO. DFT + Ud,p significantly improved the description of STO’s electronic structure close to the Fermi energy, as well as its mechanical properties, both of which are essential for the accurate modeling of defects.
The validated methodology was applied to the study of oxygen vacancy (VOq) and hydrogen interstitial (IHq) defects. Various charge and spin states of these defects were considered and modeled using large supercells, containing N = 625 ± 1 atoms. The relaxed configurations of defects were obtained through fastidious energy minimization, which consisted of a series of continued minimizations, carried out alternately in the delocalized internal and Cartesian coordinates.
Reproducibility tests showed that the potential energy surface of STO contains multiple local minima, which differ in the atomic positions and the observed internal polarization, but have comparable energies and correspond to similar static (i.e., not related to the ion transport) properties of defects.
For charged defects, the calculated total energies were corrected using a method proposed recently by Durrant et al. (114) It was found that for STO (which has a high dielectric constant) and for large supercells, the image interaction and PA corrections are small, not exceeding 0.01 and 0.15 eV, respectively.
The obtained relaxed configurations of defects were analyzed from a number of perspectives. The calculated formation energies showed that the oxygen vacancy and hydrogen interstitial strongly prefer a fully ionized state, with a 2+ and 1+ formal charge, respectively. The predicted formation energy of VO2+ was in close agreement with experimental oxidation and reduction enthalpies.
The locations of defect levels were determined using the CTL formalism. (99) With the exception of the neutral VO, all considered defects introduced states that would be seen as shallow in both conductivity and optical measurements being located down to 0.2 eV below the CB. The neutral oxygen vacancy introduced shallow (thermodynamic, optical for the emission) and deep (optical for the absorption) CTLs. This observation explains the inconsistencies seen in experimental characterizations of VO0.
It was revealed that, in general, the locations of defect levels cannot be reliably predicted based on the Kohn–Sham eigenvalues corresponding to a single charge state, as these were unable to reflect the positioning of the defect levels of VO0. The underlying reason was the change in the oxidation state of two nearby Ti atoms.
The VO1+, IH0, and IH1- defects introduced shallow and delocalized states in the forbidden gap, which were composed of the 3d–t2g and 2p states of the Ti and O atoms belonging to the same TiO2 plane. The deep-lying states of the VO0 defect were composed of the 3d–t2g states of two nearby Ti atoms, which is in contrast to previous theoretical reports that predicted a Ti 3d–eg character.
Both VO and IH defects manifested their presence also in the lower parts of the DOS. The neutral oxygen vacancy considerably (by 1.0–2.5 eV) shifted the 3s and 3p levels of two nearby Ti atoms. All hydrogen interstitials introduced states located below the O 2s and O 2p bands, which can be attributed to the formation of the OH bond.
Analysis revealed that VO and IH defects have a similar influence on the bonding and crystalline structure of STO, causing polarization and tilting of the TiO6 octahedra. Both effects occurred due to the defect’s atom removal/addition. Our calculations did not answer how localized these effects are. The Sr sublattice was found to be insensitive to the presence of the VO and IH defects.

5.2. Discussion and Future Work

This work demonstrated the DFT + Ud,p approach to be successful for describing STO, continuing along the direction set by O’Regan et al. (95) (and also other authors (101,152,153)), who demonstrated its adequacy for describing other difficult transition-metal oxides, like TiO2, CuO, and CeO2.
The calculations we carried out yield a picture that is highly consistent with the experiment. We attribute this success to: (i) an adequate formal description (theoretically motivated, free of empirical parameters, and validated for the perfect system), (ii) high quality of the calculations (near-complete basis set accuracy, accounting for the semi-core states, and a tailored procedure for carrying out geometry relaxations), and (iii) physically sound calculations and analysis (large supercells, careful accounting for structural relaxations, diligent treatment of electrostatic corrections, and explicit consideration of various charge states and calculation of the CTLs).
The adopted methodology can be applied almost straightforwardly to other perovskite materials. Here, we underline the importance of the MT-LR method, which provides a unified framework for finding the essential (correction defining) U and J parameters for various elements, even those that display magnetization. (69)
The MT-LR method relies on the perturbation of individual sites and, therefore, allows determining the U and J parameters individually for each atom chosen for correcting. This opens the possibility of studying dopants (like Fe substitutes at the Ti sites in SrTi1–xFexO3−δ) and better accounting for chemical changes, which may occur due to the introduction of defects or their migration. The MT-LR approach also enables the application of the +U and +J corrections in a self-consistent fashion (similarly as it was done in ref (93)). We did not account for the potential site dependence of the U and J parameters in the present study. Therefore, in future work, it would be prudent to check, for example, how the predicted CTLs of VO0 change if the changes in the oxidation state of two nearby Ti atoms are explicitly taken into account (e.g., by reevaluating their U and J parameters).
This work used the ONETEP program, which provides a formalism for linear-scaling DFT calculations. However, linear scaling was not exploited here, as we did not (for now) truncate the density kernel as we wanted to ensure that the results are not affected by additional approximations. In the next step, it will be essential to verify if the quality of the results can be maintained in the presence of truncation (the large band gap of STO suggests such possibility), also checking how other numerical parameters (like the localization radius of the NGWFs) influence the theoretical predictions. In this respect, the obtained results serve as a trustful checkpoint, allowing to test what further approximations can be made to reduce the computational cost. Its significant reduction will enable studying larger systems and efficiently investigating the ionic transport through transition state searches and even Born–Oppenheimer molecular dynamics simulations.
The obtained results suggest that defected STO (and, perhaps, other perovskites) might be characterized by a potential energy surface with multiple local minima, all with a similar total energy but with markedly different spatial distributions of the charge. These differences─which manifest in the observed polarization─may prove to be significant in the context of ion transport. It should be expected that the paths for O and H ion migration (and the corresponding energy barriers) depend on the atomic charges seen in the immediate vicinity of the defect and are strongly felt by the migrating ion. Our results show that these charges may be different for the minima of the same defect type. This reveals the necessity to verify how reproducible the minimum energy paths are (and the corresponding barrier heights) calculated based on the theory using, for example, the nudged elastic band method (see e.g., refs (72), (73), (134), (183), and (184)).
The obtained picture of how the defects influence the structure and bonding of STO also requires a discussion within the context of protonic transport. According to the commonly accepted view, in perovskites, it occurs through the so-called Grotthuss mechanism, (185−187) which consists of two basic events: proton rotations (around the O atoms that they are bonded to) and proton hopping (to another─neighboring─O atom). It is also believed that proton rotations are energetically less expensive than proton hopping.
Our calculations suggest that the IH defect causes a notable polarization and tilting of the TiO6 octahedra. Let us─within this context─consider a scenario, in which the migrating proton moves from one equilibrium position to another (the exact type of the migration event is not relevant here). The associated polarization and octahedral tilting can be expected not to follow the proton immediately. Thus, after the migration event (rotation or hopping), the crystal lattice has locally higher energy, temporarily displaying non-equilibrium polarization and tilting. We note that the above reasoning can be motivated based on purely inertial and configurational grounds: rotations and distortions of many octahedra consisting of much heavier (than the proton) O ions are required for the crystal lattice to reach its new equilibrium state.
With the above in mind, in perovskite materials, protonic transport might be characterized by more than two energy scales (barrier heights for the rotation and jump) because, in the discussed picture, the knowledge of the lattice excitation energies is also required. This suggests that theoretical characterizations of the ion transport, which considered only the migration between two minima of the same energy, might be incomplete.
In this work, we do not attempt to answer the questions which naturally arise from the above discussion, like: “what are the heights of the energy barriers and energies of the related lattice excitations?”, “is the proton rotation/jump to its previous position more likely?”, “what are the related probabilities?”, “how long does it take for the crystal lattice to reach new equilibrium after the proton migration event?”, “what influence does the temperature have on all of these?”. Here, we only draw the scientific community’s attention, noting that the related problems are the topic of our ongoing research.

Supporting Information

Jump To

The Supporting Information is available free of charge at

  • Detailed information related to the choice of U and J parameters, choice of the reference chemical potentials, calculation of the electrostatic corrections, and determination of the optical transition levels (PDF)

Terms & Conditions

Most electronic Supporting Information files are available without a subscription to ACS Web Editions. Such files may be downloaded by article for research use (if there is a public use license linked to the relevant article, that license may permit other uses). Permission may be obtained from ACS for other uses through requests via the RightsLink permission system:

Author Information

Jump To

  • Corresponding Author
  • Authors
    • Jacek Dziedzic - Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, PolandSchool of Chemistry, University of Southampton, Highfield, SouthamptonSO17 1BJ, U.K.Orcid
    • Tadeusz Miruszewski - Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, Poland
    • Jarosław Rybicki - Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, PolandTASK Computer Centre, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, Poland
    • Maria Gazda - Faculty of Applied Physics and Mathematics, Gdansk University of Technology, Narutowicza 11/12, Gdańsk80-233, PolandOrcid
  • Notes
    The authors declare no competing financial interest.


Jump To

This work was supported by the National Science Centre in Poland (grant number UMO-2016/23/B/ST5/02137). The authors gratefully acknowledge the PSNC Poznan Supercomputing and Networking Center (Poznan, Poland) and the TASK Academic Computer Centre (Gdansk, Poland) for providing computer time and facilities. This research was supported by the PL-Grid Infrastructure. S.W. gratefully acknowledges the help of and valuable discussions with D. D. O’Regan.


Jump To

This article references 187 other publications.

  1. 1
    Cowley, R. A.; Buyers, W. J. L.; Dolling, G. Relationship of normal modes of vibration of strontium titanate and its antiferroelectric phase transition at 110 °K. Solid State Commun. 1969, 7, 181184,  DOI: 10.1016/0038-1098(69)90720-0
  2. 2
    Shirane, G.; Yamada, Y. Lattice-dynamical study of the 110 °K phase transition in SrTiO3. Phys. Rev. 1969, 177, 858863,  DOI: 10.1103/physrev.177.858
  3. 3
    Cowley, R. A. The phase transition of strontium titanate. Philos. Trans.: Math., Phys. Eng. Sci. 1996, 354, 27992814
  4. 4
    Lytle, F. W. X-ray diffractometry of low-temperature phase transformations in strontium titanate. J. Appl. Phys. 1964, 35, 22122215,  DOI: 10.1063/1.1702820
  5. 5
    Wyckoff, R. W. G. Crystal Structures─Volume 2: Inorganic Compounds RXn, RnMX2, RnMX3; John Wiley & Sons: New York, London, Sydney, 1964.
  6. 6
    Unoki, H.; Sakudo, T. Electron spin resonance of Fe3+ in SrTiO3 with special reference to the 110 °K phase transition. J. Phys. Soc. Jpn. 1967, 23, 546552,  DOI: 10.1143/jpsj.23.546
  7. 7
    Abramov, Y. A.; Tsirelson, V. G.; Zavodnik, V. E.; Ivanov, S. A.; Brown, I. D. The chemical bond and atomic displacements in SrTiO3 from X-ray diffraction analysis. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 1995, 51, 942951,  DOI: 10.1107/s0108768195003752
  8. 8
    Cao, L.; Sozontov, E.; Zegenhagen, J. Cubic to tetragonal phase transition of SrTiO3 under epitaxial stress: An X-ray backscattering study. Phys. Status Solidi A 2000, 181, 387404,  DOI: 10.1002/1521-396x(200010)181:2<387::aid-pssa387>;2-5
  9. 9
    Schmidbauer, M.; Kwasniewski, A.; Schwarzkopf, J. High-precision absolute lattice parameter determination of SrTiO3, DyScO3 and NdGaO3 single crystals. Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater. 2012, 68, 814,  DOI: 10.1107/s0108768111046738
  10. 10
    Edwards, L. R.; Lynch, R. W. The high pressure compressibility and Grüneisen parameter of strontium titanate. J. Phys. Chem. Solids 1970, 31, 573574,  DOI: 10.1016/0022-3697(70)90098-3
  11. 11
    Beattie, A. G.; Samara, G. A. Pressure dependence of the elastic constants of SrTiO3. J. Appl. Phys. 1971, 42, 23762381,  DOI: 10.1063/1.1660551
  12. 12
    Fischer, G. J.; Wang, Z.; Karato, S. Elasticity of CaTiO3, SrTiO3 and BaTiO3 perovskites up to 3.0 GPa: The effect of crystallographic structure. Phys. Chem. Miner. 1993, 20, 97103,  DOI: 10.1007/bf00207202
  13. 13
    Capizzi, M.; Frova, A. Optical gap of strontium titanate (deviation from Urbach tail behavior). Phys. Rev. Lett. 1970, 25, 12981302,  DOI: 10.1103/physrevlett.25.1298
  14. 14
    Henrich, V. E.; Dresselhaus, G.; Zeiger, H. J. Surface defects and the electronic structure of SrTiO3 surfaces. Phys. Rev. B: Solid State 1978, 17, 49084921,  DOI: 10.1103/physrevb.17.4908
  15. 15
    Courths, R. Ultraviolet photoelectron spectroscopy (UPS) and LEED studies of BaTiO3 (001) and SrTiO3 (100) surfaces. Phys. Status Solidi B 1980, 100, 135148,  DOI: 10.1002/pssb.2221000114
  16. 16
    Kan, D.; Terashima, T.; Kanda, R.; Masuno, A.; Tanaka, K.; Chu, S.; Kan, H.; Ishizumi, A.; Kanemitsu, Y.; Shimakawa, Y. Blue-light emission at room temperature from Ar+-irradiated SrTiO3. Nat. Mater. 2005, 4, 816819,  DOI: 10.1038/nmat1498
  17. 17
    Dejneka, A.; Tyunina, M.; Narkilahti, J.; Levoska, J.; Chvostova, D.; Jastrabik, L.; Trepakov, V. A. Tensile strain induced changes in the optical spectra of SrTiO3 epitaxial thin films. Phys. Solid State 2010, 52, 20822089,  DOI: 10.1134/s1063783410100124
  18. 18
    Ravichandran, J.; Siemons, W.; Scullin, M. L.; Mukerjee, S.; Huijben, M.; Moore, J. E.; Majumdar, A.; Ramesh, R. Tuning the electronic effective mass in double-doped SrTiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 035101,  DOI: 10.1103/physrevb.83.035101
  19. 19
    Tarun, M. C.; McCluskey, M. D. Infrared absorption of hydrogen-related defects in strontium titanate. J. Appl. Phys. 2011, 109, 063706,  DOI: 10.1063/1.3561867
  20. 20
    Rice, W. D.; Ambwani, P.; Bombeck, M.; Thompson, J. D.; Haugstad, G.; Leighton, C.; Crooker, S. A. Persistent optically induced magnetism in oxygen-deficient strontium titanate. Nat. Mater. 2014, 13, 481487,  DOI: 10.1038/nmat3914
  21. 21
    Li, Y.; Lei, Y.; Shen, B. G.; Sun, J. R. Visible-light-accelerated oxygen vacancy migration in strontium titanate. Sci. Rep. 2015, 5, 14576,  DOI: 10.1038/srep14576
  22. 22
    Crespillo, M. L.; Graham, J. T.; Agulló-López, F.; Zhang, Y.; Weber, W. J. Isolated oxygen vacancies in strontium titanate shine red: Optical identification of Ti3+ polarons. Appl. Mater. Today 2018, 12, 131137,  DOI: 10.1016/j.apmt.2018.04.006
  23. 23
    Crespillo, M. L.; Graham, J. T.; Agulló-López, F.; Zhang, Y.; Weber, W. J. Recent advances on carrier and exciton self-trapping in strontium titanate: Understanding the luminescence emissions. Crystals 2019, 9, 95,  DOI: 10.3390/cryst9020095
  24. 24
    Siebenhofer, M.; Viernstein, A.; Morgenbesser, M.; Fleig, J.; Kubicek, M. Photoinduced electronic and ionic effects in strontium titanate. Mater. Adv. 2021, 2, 75837619,  DOI: 10.1039/d1ma00906k
  25. 25
    Tufte, O. N.; Chapman, P. W. Electron mobility in semiconducting strontium titanate. Phys. Rev. 1967, 155, 796802,  DOI: 10.1103/physrev.155.796
  26. 26
    Lee, C.; Yahia, J.; Brebner, J. L. Electronic conduction in slightly reduced strontium titanate at low temperatures. Phys. Rev. B: Solid State 1971, 3, 25252533,  DOI: 10.1103/physrevb.3.2525
  27. 27
    Neville, R. C.; Hoeneisen, B.; Mead, C. A. Permittivity of strontium titanate. J. Appl. Phys. 1972, 43, 21242131,  DOI: 10.1063/1.1661463
  28. 28
    Maier, J.; Schwitzgebel, G.; Hagemann, H.-J. Electrochemical investigations of conductivity and chemical diffusion in pure and doped cubic SrTiO3 and BaTiO3. J. Solid State Chem. 1985, 58, 113,  DOI: 10.1016/0022-4596(85)90264-6
  29. 29
    Moos, R.; Hardtl, K. H. Defect chemistry of donor-doped and undoped strontium titanate ceramics between 1000° and 1400 °C. J. Am. Ceram. Soc. 1997, 80, 25492562,  DOI: 10.1111/j.1151-2916.1997.tb03157.x
  30. 30
    Raevski, I. P.; Maksimov, S. M.; Fisenko, A. V.; Prosandeyev, S. A.; Osipenko, I. A.; Tarasenko, P. F. Study of intrinsic point defects in oxides of the perovskite family: II. Experiment. J. Phys.: Condens. Matter 1998, 10, 80158032,  DOI: 10.1088/0953-8984/10/36/012
  31. 31
    van Benthem, K.; Elsässer, C.; French, R. H. Bulk electronic structure of SrTiO3: Experiment and theory. J. Appl. Phys. 2001, 90, 61566164,  DOI: 10.1063/1.1415766
  32. 32
    Wideroe, M.; Waser, R.; Norby, T. Transport of hydrogen species in a single crystal SrTiO3. Solid State Ionics 2006, 177, 14691476,  DOI: 10.1016/j.ssi.2006.07.029
  33. 33
    Kreuer, K. D. Aspects of the formation and mobility of protonic charge carriers and the stability of perovskite-type oxides. Solid State Ionics 1999, 125, 285302,  DOI: 10.1016/s0167-2738(99)00188-5
  34. 34
    Steinsvik, S.; Larring, Y.; Norby, T. Hydrogen ion conduction in iron-substituted strontium titanate, SrTi1–xFexO3–x/2 (0 ≤ x ≤ 0.8). Solid State Ionics 2001, 143, 103116,  DOI: 10.1016/s0167-2738(01)00838-4
  35. 35
    Merkle, R.; Maier, J. How is oxygen incorporated into oxides? A comprehensive kinetic study of a simple solid-state reaction with SrTiO3 as a model material. Angew. Chem., Int. Ed. 2008, 47, 38743894,  DOI: 10.1002/anie.200700987
  36. 36
    Phoon, B. L.; Lai, C. W.; Juan, J. C.; Show, P.-L.; Chen, W.-H. A review of synthesis and morphology of SrTiO3 for energy and other applications. Int. J. Energy Res. 2019, 43, 51515174,  DOI: 10.1002/er.4505
  37. 37
    Miruszewski, T.; Dzierzgowski, K.; Winiarz, P.; Wachowski, S.; Mielewczyk-Gryń, A.; Gazda, M. Structural properties and water uptake of SrTi1–xFexO3–x/2−δ. Materials 2020, 13, 965,  DOI: 10.3390/ma13040965
  38. 38
    Miruszewski, T.; Dzierzgowski, K.; Winiarz, P.; Jaworski, D.; Wiciak-Pawłowska, K.; Skubida, W.; Wachowski, S.; Mielewczyk-Gryń, A.; Gazda, M. Structure and transport properties of triple-conducting BaxSr1–xTi1–yFeyO3−δ oxides. RSC Adv. 2021, 11, 1957019578,  DOI: 10.1039/d0ra10048j
  39. 39
    Papac, M.; Stevanović, V.; Zakutayev, A.; O’Hayre, R. Triple ionic–electronic conducting oxides for next-generation electrochemical devices. Nat. Mater. 2021, 20, 301313,  DOI: 10.1038/s41563-020-00854-8
  40. 40
    Samat, A. A.; Darus, M.; Osman, N.; Baharuddin, N. A.; Anwar, M. A short review on triple conducting oxide cathode materials for proton conducting solid oxide fuel cell. AIP Conf. Proc. 2021, 2339, 020233,  DOI: 10.1063/5.0044224
  41. 41
    Rothschild, A.; Menesklou, W.; Tuller, H. L.; Ivers-Tiffée, E. Electronic structure, defect chemistry, and transport properties of SrTi1–xFexO3–y solid solutions. Chem. Mater. 2006, 18, 36513659,  DOI: 10.1021/cm052803x
  42. 42
    Jung, W.; Tuller, H. L. Impedance study of SrTi1–xFexO3−δ (x=0.05 to 0.80) mixed ionic-electronic conducting model cathode. Solid State Ionics 2009, 180, 843847,  DOI: 10.1016/j.ssi.2009.02.008
  43. 43
    Mroziński, A.; Molin, S.; Jasiński, P. Effect of sintering temperature on electrochemical performance of porous SrTi1–xFexO3−δ (x = 0.35, 0.5, 0.7) oxygen electrodes for solid oxide cells. J. Solid State Electrochem. 2020, 24, 873882,  DOI: 10.1007/s10008-020-04534-0
  44. 44
    Winiarz, P.; Mielewczyk-Gryń, A.; Wachowski, S.; Jasiński, P.; Witkowska, A.; Gazda, M. Structural and electrical properties of titanium-doped yttrium niobate. J. Alloys Compd. 2018, 767, 11861195,  DOI: 10.1016/j.jallcom.2018.07.134
  45. 45
    Freysoldt, C.; Grabowski, B.; Hickel, T.; Neugebauer, J.; Kresse, G.; Janotti, A.; Van de Walle, C. G. First-principles calculations for point defects in solids. Rev. Mod. Phys. 2014, 86, 253305,  DOI: 10.1103/revmodphys.86.253
  46. 46
    Piskunov, S.; Heifets, E.; Eglitis, R. I.; Borstel, G. Bulk properties and electronic structure of SrTiO3, BaTiO3, PbTiO3 perovskites: An ab initio HF/DFT study. Comput. Mater. Sci. 2004, 29, 165178,  DOI: 10.1016/j.commatsci.2003.08.036
  47. 47
    Wahl, R.; Vogtenhuber, D.; Kresse, G. SrTiO3 and BaTiO3 revisited using the projector augmented wave method: Performance of hybrid and semilocal functionals. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 104116,  DOI: 10.1103/physrevb.78.104116
  48. 48
    Kinaci, A.; Sevik, C.; Çağin, T. Electronic transport properties of SrTiO3 and its alloys: Sr1–xLaxTiO3 and SrTi1–xMxO3 (M = Nb, Ta). Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 155114,  DOI: 10.1103/PhysRevB.82.155114
  49. 49
    Evarestov, R. A.; Blokhin, E.; Gryaznov, D.; Kotomin, E. A.; Maier, J. Phonon calculations in cubic and tetragonal phases of SrTiO3: A comparative LCAO and plane-wave study. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 134108,  DOI: 10.1103/physrevb.83.134108
  50. 50
    Sakhya, A. P.; Maibam, J.; Saha, S.; Chanda, S.; Dutta, A.; Sharma, B. I.; Thapa, R.; Sinha, T. P. Electronic structure and elastic properties of ATiO3 (A = Ba, Sr, Ca) perovskites: A first principles study. Indian J. Pure Appl. Phys. 2015, 53, 102109
  51. 51
    Mubarak, A. A. The first-principle study of the electronic, optical and thermoelectric properties of XTiO3 (X = Ca, Sr and Ba) compounds. Int. J. Mod. Phys. B 2016, 30, 1650141,  DOI: 10.1142/s0217979216501411
  52. 52
    Holmström, E.; Spijker, P.; Foster, A. S. The interface of SrTiO3 and H2O from density functional theory molecular dynamics. Proc. R. Soc. A 2016, 472, 20160293,  DOI: 10.1098/rspa.2016.0293
  53. 53
    Ricci, D.; Bano, G.; Pacchioni, G.; Illas, F. Electronic structure of a neutral oxygen vacancy in SrTiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 68, 224105,  DOI: 10.1103/physrevb.68.224105
  54. 54
    Eglitis, R. I.; Piskunov, S.; Heifets, E.; Kotomin, E. A.; Borstel, G. Ab initio study of the SrTiO3, BaTiO3 and PbTiO3 (001) surfaces. Ceram. Int. 2004, 30, 19891992,  DOI: 10.1016/j.ceramint.2003.12.176
  55. 55
    Heifets, E.; Piskunov, S.; Kotomin, E. A.; Zhukovskii, Y. F.; Ellis, D. E. Electronic structure and thermodynamic stability of double-layered SrTiO3(001) surfaces: Ab initio simulations. Phys. Rev. B: Condens. Matter Mater. Phys. 2007, 75, 115417,  DOI: 10.1103/physrevb.75.115417
  56. 56
    Alexandrov, V. E.; Kotomin, E. A.; Maier, J.; Evarestov, R. A. First-principles study of bulk and surface oxygen vacancies in SrTiO3 crystal. Eur. Phys. J. B 2009, 72, 5357,  DOI: 10.1140/epjb/e2009-00339-4
  57. 57
    Mitra, C.; Lin, C.; Robertson, J.; Demkov, A. A. Electronic structure of oxygen vacancies in SrTiO3 and LaAlO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 155105,  DOI: 10.1103/physrevb.86.155105
  58. 58
    Hamann, D. R.; Vanderbilt, D. Maximally localized Wannier functions for GW quasiparticles. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 045109,  DOI: 10.1103/physrevb.79.045109
  59. 59
    Sponza, L.; Véniard, V.; Sottile, F.; Giorgetti, C.; Reining, L. Role of localized electrons in electron-hole interaction: The case of SrTiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 235102,  DOI: 10.1103/physrevb.87.235102
  60. 60
    Deguchi, D.; Sato, K.; Kino, H.; Kotani, T. Accurate energy bands calculated by the hybrid quasiparticle self-consistent GW method implemented in the ecalj package. Jpn. J. Appl. Phys. 2016, 55, 051201,  DOI: 10.7567/jjap.55.051201
  61. 61
    Bhandari, C.; van Schilfgaarde, M.; Kotani, T.; Lambrecht, W. R. L. All-electron quasiparticle self-consistent GW band structures for SrTiO3 including lattice polarization corrections in different phases. Phys. Rev. Mater. 2018, 2, 013807,  DOI: 10.1103/physrevmaterials.2.013807
  62. 62
    Triggiani, L.; Muñoz-García, A. B.; Agostiano, A.; Pavone, M. Promoting oxygen vacancy formation and p-type conductivity in SrTiO3 via alkali metal doping: A first principles study. Phys. Chem. Chem. Phys. 2016, 18, 2895128959,  DOI: 10.1039/c6cp05089a
  63. 63
    Himmetoglu, B.; Floris, A.; de Gironcoli, S.; Cococcioni, M. Hubbard-corrected DFT energy functionals: The LDA+U description of correlated systems. Int. J. Quantum Chem. 2014, 114, 1449,  DOI: 10.1002/qua.24521
  64. 64
    Tolba, S. A.; Gameel, K. M.; Ali, B. A.; Almossalami, H. A.; Allam, N. K. Density Functional Calculations; Yang, G., Ed.; IntechOpen: London, 2018; Chapter 1.
  65. 65
    Harun, K.; Salleh, N. A.; Deghfel, B.; Yaakob, M. K.; Mohamad, A. A. DFT + U calculations for electronic, structural, and optical properties of ZnO wurtzite structure: A review. Results Phys. 2020, 16, 102829,  DOI: 10.1016/j.rinp.2019.102829
  66. 66
    Capdevila-Cortada, M.; Łodziana, Z.; López, N. Performance of DFT+U approaches in the study of catalytic materials. ACS Catal. 2016, 6, 83708379,  DOI: 10.1021/acscatal.6b01907
  67. 67
    Cococcioni, M.; de Gironcoli, S. 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,  DOI: 10.1103/physrevb.71.035105
  68. 68
    Timrov, I.; Marzari, N.; Cococcioni, M. Hubbard parameters from density-functional perturbation theory. Phys. Rev. B 2018, 98, 085127,  DOI: 10.1103/physrevb.98.085127
  69. 69
    Linscott, E. B.; Cole, D. J.; Payne, M. C.; O’Regan, D. D. Role of spin in the calculation of Hubbard U and Hund’s J parameters from first principles. Phys. Rev. B 2018, 98, 235157,  DOI: 10.1103/physrevb.98.235157
  70. 70
    Aryasetiawan, F.; Imada, M.; Georges, A.; Kotliar, G.; Biermann, S.; Lichtenstein, A. I. Frequency-dependent local interactions and low-energy effective models from electronic structure calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 70, 195104,  DOI: 10.1103/physrevb.70.195104
  71. 71
    Aryasetiawan, F.; Karlsson, K.; Jepsen, O.; Schönberger, U. Calculations of Hubbard U from first-principles. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 74, 125106,  DOI: 10.1103/physrevb.74.125106
  72. 72
    Bork, N.; Bonanos, N.; Rossmeisl, J.; Vegge, T. Ab initio charge analysis of pure and hydrogenated perovskites. J. Appl. Phys. 2011, 109, 033702,  DOI: 10.1063/1.3536484
  73. 73
    Bork, N.; Bonanos, N.; Rossmeisl, J.; Vegge, T. Thermodynamic and kinetic properties of hydrogen defect pairs in SrTiO3 from density functional theory. Phys. Chem. Chem. Phys. 2011, 13, 1525615263,  DOI: 10.1039/c1cp20406h
  74. 74
    Lee, J.; Demkov, A. A. Charge origin and localization at the n-type SrTiO3/LaAlO3 interface. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 193104,  DOI: 10.1103/physrevb.78.193104
  75. 75
    Kim, M. S.; Park, C. H. GW calculation of the electronic structure of cubic perovskite SrTiO3. J. Korean Phys. Soc. 2010, 56, 490493,  DOI: 10.3938/jkps.56.362
  76. 76
    Sikam, P.; Moontragoon, P.; Sararat, C.; Karaphun, A.; Swatsitang, E.; Pinitsoontorn, S.; Thongbai, P. DFT calculation and experimental study on structural, optical and magnetic properties of Co-doped SrTiO3. Appl. Surf. Sci. 2018, 446, 92113,  DOI: 10.1016/j.apsusc.2018.02.161
  77. 77
    Brown, J. J.; Ke, Z.; Geng, W.; Page, A. J. Oxygen vacancy defect migration in titanate perovskite surfaces: Effect of the A-site cations. J. Phys. Chem. C 2018, 122, 1459014597,  DOI: 10.1021/acs.jpcc.8b03322
  78. 78
    Pavarini, E.; Biermann, S.; Poteryaev, A.; Lichtenstein, A. I.; Georges, A.; Andersen, O. K. Mott transition and suppression of orbital fluctuations in orthorhombic 3d1 perovskites. Phys. Rev. Lett. 2004, 92, 176403,  DOI: 10.1103/physrevlett.92.176403
  79. 79
    Okamoto, S.; Millis, A. J.; Spaldin, N. A. Lattice relaxation in oxide heterostructures: LaTiO3/SrTiO3 superlattices. Phys. Rev. Lett. 2006, 97, 056802,  DOI: 10.1103/physrevlett.97.056802
  80. 80
    Cuong, D. D.; Lee, B.; Choi, K. M.; Ahn, H.-S.; Han, S.; Lee, J. Oxygen vacancy clustering and electron localization in oxygen-deficient SrTiO3: LDA + U study. Phys. Rev. Lett. 2007, 98, 115503,  DOI: 10.1103/physrevlett.98.115503
  81. 81
    Hou, Z.; Terakura, K. Defect states induced by oxygen vacancies in cubic SrTiO3: First-principles calculations. J. Phys. Soc. Jpn. 2010, 79, 114704,  DOI: 10.1143/jpsj.79.114704
  82. 82
    Choi, M.; Oba, F.; Kumagai, Y.; Tanaka, I. Anti-ferrodistortive-like oxygen-octahedron rotation induced by the oxygen vacancy in cubic SrTiO3. Adv. Mater. 2013, 25, 8690,  DOI: 10.1002/adma.201203580
  83. 83
    Liu, B.; Cooper, V. R.; Xu, H.; Xiao, H.; Zhang, Y.; Weber, W. J. Composition dependent intrinsic defect structures in SrTiO3. Phys. Chem. Chem. Phys. 2014, 16, 1559015596,  DOI: 10.1039/c4cp01510j
  84. 84
    Gogoi, P. K.; Trevisanutto, P. E.; Yang, M.; Santoso, I.; Asmara, T. C.; Terentjevs, A.; Della Sala, F.; Breese, M. B. H.; Venkatesan, T.; Feng, Y. P. Optical conductivity renormalization of graphene on SrTiO3 due to resonant excitonic effects mediated by Ti 3d orbitals. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 035424,  DOI: 10.1103/physrevb.91.035424
  85. 85
    Jeschke, H. O.; Shen, J.; Valentí, R. Localized versus itinerant states created by multiple oxygen vacancies in SrTiO3. New J. Phys. 2015, 17, 023034,  DOI: 10.1088/1367-2630/17/2/023034
  86. 86
    Zhang, L.; Liu, B.; Zhuang, H.; Kent, P. R. C.; Cooper, V. R.; Ganesh, P.; Xu, H. Oxygen vacancy diffusion in bulk SrTiO3 from density functional theory calculations. Comput. Mater. Sci. 2016, 118, 309315,  DOI: 10.1016/j.commatsci.2016.02.041
  87. 87
    Brovko, O. O.; Tosatti, E. Controlling the magnetism of oxygen surface vacancies in SrTiO3 through charging. Phys. Rev. Mater. 2017, 1, 044405,  DOI: 10.1103/physrevmaterials.1.044405
  88. 88
    Solovyev, I.; Hamada, N.; Terakura, K. t2g versus all 3d localization in LaMO3 perovskites (M=Ti–Cu): First-principles study. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 53, 71587170,  DOI: 10.1103/physrevb.53.7158
  89. 89
    Xie, Y.; Cao, H.-Y.; Zhou, Y.; Chen, S.; Xiang, H.; Gong, X.-G. Oxygen vacancy induced flat phonon mode at FeSe/SrTiO3 interface. Sci. Rep. 2015, 5, 10011,  DOI: 10.1038/srep10011
  90. 90
    Kim, Y. S.; Kim, J.; Moon, S. J.; Choi, W. S.; Chang, Y. J.; Yoon, J.-G.; Yu, J.; Chung, J.-S.; Noh, T. W. Localized electronic states induced by defects and possible origin of ferroelectricity in strontium titanate thin films. Appl. Phys. Lett. 2009, 94, 202906,  DOI: 10.1063/1.3139767
  91. 91
    Carlotto, S.; Natile, M. M.; Glisenti, A.; Vittadini, A. Electronic structure of SrTi1–xMxO3−δ (M=Co, Ni, Cu) perovskite-type doped-titanate crystals by DFT and DFT+U calculations. Chem. Phys. Lett. 2013, 588, 102108,  DOI: 10.1016/j.cplett.2013.10.020
  92. 92
    Maldonado, F.; Maza, L.; Stashans, A. Electronic properties of Cr-, B-doped and codoped SrTiO3. J. Phys. Chem. Solids 2017, 100, 18,  DOI: 10.1016/j.jpcs.2016.09.003
  93. 93
    Ricca, C.; Timrov, I.; Cococcioni, M.; Marzari, N.; Aschauer, U. Self-consistent DFT + U + V study of oxygen vacancies in SrTiO3. Phys. Rev. Res. 2020, 2, 023313,  DOI: 10.1103/physrevresearch.2.023313
  94. 94
    Campo, V. L., Jr.; Cococcioni, M. Extended DFT + U + V method with on-site and inter-site electronic interactions. J. Phys.: Condens. Matter 2010, 22, 055602
  95. 95
    Orhan, O. K.; O’Regan, D. D. First-principles Hubbard U and Hund’s J corrected approximate density functional theory predicts an accurate fundamental gap in rutile and anatase TiO2. Phys. Rev. B 2020, 101, 245137,  DOI: 10.1103/physrevb.101.245137
  96. 96
    Prentice, J. C. A.; Aarons, J.; Womack, J. C.; Allen, A. E. A.; Andrinopoulos, L.; Anton, L.; Bell, R. A.; Bhandari, A.; Bramley, G. A.; Charlton, R. J. The ONETEP linear-scaling density functional theory program. J. Chem. Phys. 2020, 152, 174111,  DOI: 10.1063/5.0004445
  97. 97
    Skylaris, C.-K.; Haynes, P. D.; Mostofi, A. A.; Payne, M. C. Introducing ONETEP: Linear-scaling density functional simulations on parallel computers. J. Chem. Phys. 2005, 122, 084119,  DOI: 10.1063/1.1839852
  98. 98
    Lany, S.; Zunger, A. Light- and bias-induced metastabilities in Cu(In,Ga)Se2 based solar cells caused by the (VSe-VCu) vacancy complex. J. Appl. Phys. 2006, 100, 113725,  DOI: 10.1063/1.2388256
  99. 99
    Lany, S.; Zunger, A. Assessment of correction methods for the band-gap problem and for finite-size effects in supercell defect calculations: Case studies for ZnO and GaAs. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 235104,  DOI: 10.1103/physrevb.78.235104
  100. 100
    Gallino, F.; Pacchioni, G.; Di Valentin, C. Transition levels of defect centers in ZnO by hybrid functionals and localized basis set approach. J. Chem. Phys. 2010, 133, 144512,  DOI: 10.1063/1.3491271
  101. 101
    Mattioli, G.; Alippi, P.; Filippone, F.; Caminiti, R.; Amore Bonapasta, A. Deep versus shallow behavior of intrinsic defects in rutile and anatase TiO2 polymorphs. J. Phys. Chem. C 2010, 114, 2169421704,  DOI: 10.1021/jp1041316
  102. 102
    Chakrabarty, A.; Patterson, C. H. Transition levels of defects in ZnO: Total energy and Janak’s theorem methods. J. Chem. Phys. 2012, 137, 054709,  DOI: 10.1063/1.4739316
  103. 103
    Gerosa, M.; Bottani, C. E.; Caramella, L.; Onida, G.; Di Valentin, C.; Pacchioni, G. Defect calculations in semiconductors through a dielectric-dependent hybrid DFT functional: The case of oxygen vacancies in metal oxides. J. Chem. Phys. 2015, 143, 134702,  DOI: 10.1063/1.4931805
  104. 104
    Makov, G.; Payne, M. C. Periodic boundary conditions in ab initio calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 1995, 51, 40144022,  DOI: 10.1103/physrevb.51.4014
  105. 105
    Freysoldt, C.; Neugebauer, J.; Van de Walle, C. G. Fully ab initio finite-size corrections for charged-defect supercell calculations. Phys. Rev. Lett. 2009, 102, 016402,  DOI: 10.1103/PhysRevLett.102.016402
  106. 106
    Hine, N. D. M.; Frensch, K.; Foulkes, W. M. C.; Finnis, M. W. Supercell size scaling of density functional theory formation energies of charged defects. Phys. Rev. B: Condens. Matter Mater. Phys. 2009, 79, 024112,  DOI: 10.1103/physrevb.79.024112
  107. 107
    Corsetti, F.; Mostofi, A. A. System-size convergence of point defect properties: The case of the silicon vacancy. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 035209,  DOI: 10.1103/physrevb.84.035209
  108. 108
    Hine, N. D. M.; Dziedzic, J.; Haynes, P. D.; Skylaris, C.-K. Electrostatic interactions in finite systems treated with periodic boundary conditions: Application to linear-scaling density functional theory. J. Chem. Phys. 2011, 135, 204103,  DOI: 10.1063/1.3662863
  109. 109
    Taylor, S. E.; Bruneval, F. Understanding and correcting the spurious interactions in charged supercells. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 075155,  DOI: 10.1103/physrevb.84.075155
  110. 110
    Komsa, H.-P.; Rantala, T. T.; Pasquarello, A. Finite-size supercell correction schemes for charged defect calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 045112,  DOI: 10.1103/physrevb.86.045112
  111. 111
    Chen, W.; Pasquarello, A. Correspondence of defect energy levels in hybrid density functional theory and many-body perturbation theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 115104,  DOI: 10.1103/physrevb.88.115104
  112. 112
    Murphy, S. T.; Hine, N. D. M. Anisotropic charge screening and supercell size convergence of defect formation energies. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 094111,  DOI: 10.1103/physrevb.87.094111
  113. 113
    Kumagai, Y.; Oba, F. Electrostatics-based finite-size corrections for first-principles point defect calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 195205,  DOI: 10.1103/physrevb.89.195205
  114. 114
    Durrant, T. R.; Murphy, S. T.; Watkins, M. B.; Shluger, A. L. Relation between image charge and potential alignment corrections for charged defects in periodic boundary conditions. J. Chem. Phys. 2018, 149, 024103,  DOI: 10.1063/1.5029818
  115. 115
    Naik, M. H.; Jain, M. CoFFEE: corrections for formation energy and eigenvalues for charged defect simulations. Comput. Phys. Commun. 2018, 226, 114126,  DOI: 10.1016/j.cpc.2018.01.011
  116. 116
    Tsunoda, N.; Kumagai, Y.; Takahashi, A.; Oba, F. Electrically benign defect behavior in zinc tin nitride revealed from first principles. Phys. Rev. Appl. 2018, 10, 011001,  DOI: 10.1103/physrevapplied.10.011001
  117. 117
    Chagas da Silva, M.; Lorke, M.; Aradi, B.; Farzalipour Tabriz, M.; Frauenheim, T.; Rubio, A.; Rocca, D.; Deák, P. Self-consistent potential correction for charged periodic systems. Phys. Rev. Lett. 2021, 126, 076401,  DOI: 10.1103/PhysRevLett.126.076401
  118. 118
    Astala, R.; Bristowe, P. D. Ab initio study of the oxygen vacancy in SrTiO3. Modell. Simul. Mater. Sci. Eng. 2001, 9, 415422,  DOI: 10.1088/0965-0393/9/5/306
  119. 119
    Tanaka, T.; Matsunaga, K.; Ikuhara, Y.; Yamamoto, T. First-principles study on structures and energetics of intrinsic vacancies in SrTiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2003, 68, 205213,  DOI: 10.1103/physrevb.68.205213
  120. 120
    Buban, J. P.; Iddir, H.; Öğüt, S. Structural and electronic properties of oxygen vacancies in cubic and antiferrodistortive phases of SrTiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 69, 180102,  DOI: 10.1103/physrevb.69.180102
  121. 121
    Luo, W.; Duan, W.; Louie, S. G.; Cohen, M. L. Structural and electronic properties of n-doped and p-doped SrTiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2004, 70, 214109,  DOI: 10.1103/physrevb.70.214109
  122. 122
    Carrasco, J.; Illas, F.; Lopez, N.; Kotomin, E. A.; Zhukovskii, Y. F.; Evarestov, R. A.; Mastrikov, Y. A.; Piskunov, S.; Maier, J. First-principles calculations of the atomic and electronic structure of F centers in the bulk and on the (001) surface of SrTiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, 73, 064106,  DOI: 10.1103/physrevb.73.064106
  123. 123
    Evarestov, R. A.; Kotomin, E. A.; Zhukovskii, Y. F. DFT study of a single F center in cubic SrTiO3 perovskite. Int. J. Quantum Chem. 2006, 106, 21732183,  DOI: 10.1002/qua.20855
  124. 124
    Zhukovskii, Y. F.; Kotomin, E. A.; Evarestov, R. A.; Ellis, D. E. Periodic models in quantum chemical simulations of F centers in crystalline metal oxides. Int. J. Quantum Chem. 2007, 107, 29562985,  DOI: 10.1002/qua.21483
  125. 125
    Zhukovskii, Y. F.; Kotomin, E. A.; Piskunov, S.; Ellis, D. E. A comparative ab initio study of bulk and surface oxygen vacancies in PbTiO3, PbZrO3 and SrTiO3 perovskites. Solid State Commun. 2009, 149, 13591362,  DOI: 10.1016/j.ssc.2009.05.023
  126. 126
    Zhukovskii, Y. F.; Kotomin, E. A.; Piskunov, S.; Mastrikov, Y. A.; Ellis, D. E. The effect of oxygen vacancies on the atomic and electronic structure of cubic ABO3 perovskite bulk and the (001) surface: Ab initio calculations. Ferroelectrics 2009, 379, 191198,  DOI: 10.1080/00150190902852240
  127. 127
    Djermouni, M.; Zaoui, A.; Kacimi, S.; Bouhafs, B. Vacancy defects in strontium titanate: Ab initio calculation. Comput. Mater. Sci. 2010, 49, 904909,  DOI: 10.1016/j.commatsci.2010.06.045
  128. 128
    Ertekin, E.; Srinivasan, V.; Ravichandran, J.; Rossen, P. B.; Siemons, W.; Majumdar, A.; Ramesh, R.; Grossman, J. C. Interplay between intrinsic defects, doping, and free carrier concentration in SrTiO3 thin films. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 85, 195460,  DOI: 10.1103/physrevb.85.195460
  129. 129
    Evarestov, R.; Blokhin, E.; Gryaznov, D.; Kotomin, E. A.; Merkle, R.; Maier, J. Jahn-Teller effect in the phonon properties of defective SrTiO3 from first principles. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 85, 174303,  DOI: 10.1103/physrevb.85.174303
  130. 130
    Lin, C.; Mitra, C.; Demkov, A. A. Orbital ordering under reduced symmetry in transition metal perovskites: Oxygen vacancy in SrTiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 86, 161102,  DOI: 10.1103/physrevb.86.161102
  131. 131
    Al-Hamadany, R.; Goss, J. P.; Briddon, P. R.; Mojarad, S. A.; Al-Hadidi, M.; O’Neill, A. G.; Rayson, M. J. Oxygen vacancy migration in compressively strained SrTiO3. J. Appl. Phys. 2013, 113, 024108,  DOI: 10.1063/1.4775397
  132. 132
    Souto-Casares, J.; Spaldin, N. A.; Ederer, C. Oxygen vacancies in strontium titanate: A DFT + DMFT study. Phys. Rev. Res. 2021, 3, 023027,  DOI: 10.1103/physrevresearch.3.023027
  133. 133
    Xiong, K.; Robertson, J.; Clark, S. J. Behavior of hydrogen in wide band gap oxides. J. Appl. Phys. 2007, 102, 083710,  DOI: 10.1063/1.2798910
  134. 134
    Bork, N.; Bonanos, N.; Rossmeisl, J.; Vegge, T. Simple descriptors for proton-conducting perovskites from density functional theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 014103,  DOI: 10.1103/physrevb.82.014103
  135. 135
    Baker, J. N.; Bowes, P. C.; Irving, D. L. Hydrogen solubility in donor-doped SrTiO3 from first principles. Appl. Phys. Lett. 2018, 113, 132904,  DOI: 10.1063/1.5047793
  136. 136
    T-Thienprasert, J.; Fongkaew, I.; Singh, D. J.; Du, M.-H.; Limpijumnong, S. Identification of hydrogen defects in SrTiO3 by first-principles local vibration mode calculations. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 85, 125205,  DOI: 10.1103/physrevb.85.125205
  137. 137
    Münch, W.; Kreuer, K.-D.; Seifertli, G.; Majer, J. A quantum molecular dynamics study of proton diffusion in SrTiO3 and CaTiO3. Solid State Ionics 1999, 125, 3945,  DOI: 10.1016/s0167-2738(99)00156-3
  138. 138
    Münch, W.; Kreuer, K.-D.; Seifert, G.; Maier, J. Proton diffusion in perovskites: Comparison between BaCeO3, BaZrO3, SrTiO3, and CaTiO3 using quantum molecular dynamics. Solid State Ionics 2000, 136-137, 183189,  DOI: 10.1016/s0167-2738(00)00304-0
  139. 139
    Villamagua, L.; Barreto, R.; Prócel, L. M.; Stashans, A. Hydrogen impurity in SrTiO3: Structure, electronic properties and migration. Phys. Scr. 2007, 75, 374378,  DOI: 10.1088/0031-8949/75/3/024
  140. 140
    Kohn, W. Density functional and density matrix method scaling linearly with the number of atoms. Phys. Rev. Lett. 1996, 76, 31683171,  DOI: 10.1103/physrevlett.76.3168
  141. 141
    Prodan, E.; Kohn, W. Nearsightedness of electronic matter. Proc. Natl. Acad. Sci. U.S.A. 2005, 102, 1163511638,  DOI: 10.1073/pnas.0505436102
  142. 142
    Skylaris, C.-K.; Mostofi, A. A.; Haynes, P. D.; Pickard, C. J.; Payne, M. C. Accurate kinetic energy evaluation in electronic structure calculations with localized functions on real space grids. Comput. Phys. Commun. 2001, 140, 315322,  DOI: 10.1016/s0010-4655(01)00248-x
  143. 143
    O’Regan, D. D.; Hine, N. D. M.; Payne, M. C.; Mostofi, A. A. Linear-scaling DFT + U with full local orbital optimization. Phys. Rev. B: Condens. Matter Mater. Phys. 2012, 85, 085107
  144. 144
    O’Regan, D. D.; Hine, N. D. M.; Payne, M. C.; Mostofi, A. A. Projector self-consistent DFT + U using nonorthogonal generalized Wannier functions. Phys. Rev. B: Condens. Matter Mater. Phys. 2010, 82, 081102
  145. 145
    O’Regan, D. D.; Payne, M. C.; Mostofi, A. A. Subspace representations in ab initio methods for strongly correlated systems. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 245124
  146. 146
    Hine, N. D. M.; Robinson, M.; Haynes, P. D.; Skylaris, C.-K.; Payne, M. C.; Mostofi, A. A. Accurate ionic forces and geometry optimization in linear-scaling density-functional theory with local orbitals. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 83, 195102,  DOI: 10.1103/physrevb.83.195102
  147. 147
    Ruiz-Serrano, Á.; Hine, N. D. M.; Skylaris, C.-K. Pulay forces from localized orbitals optimized in situ using a psinc basis set. J. Chem. Phys. 2012, 136, 234101,  DOI: 10.1063/1.4728026
  148. 148
    Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized gradient approximation made simple. Phys. Rev. Lett. 1996, 77, 38653868,  DOI: 10.1103/physrevlett.77.3865
  149. 149
    Perdew, J. P.; Ruzsinszky, A.; Csonka, G. I.; Vydrov, O. A.; Scuseria, G. E.; Constantin, L. A.; Zhou, X.; Burke, K. Restoring the density-gradient expansion for exchange in solids and surfaces. Phys. Rev. Lett. 2008, 100, 136406,  DOI: 10.1103/physrevlett.100.136406
  150. 150
    The Rappe group. PBE GGA PSP database. (accessed Feb 3, 2022).
  151. 151
    Rappe, A. M.; Rabe, K. M.; Kaxiras, E.; Joannopoulos, J. D. Optimized pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 41, 12271230,  DOI: 10.1103/physrevb.41.1227
  152. 152
    Himmetoglu, B.; Wentzcovitch, R. M.; Cococcioni, M. First-principles study of electronic and structural properties of CuO. Phys. Rev. B: Condens. Matter Mater. Phys. 2011, 84, 115108,  DOI: 10.1103/physrevb.84.115108
  153. 153
    Plata, J. J.; Márquez, A. M.; Sanz, J. F. Communication: Improving the density functional theory+U description of CeO2 by including the contribution of the O 2p electrons. J. Chem. Phys. 2012, 136, 041101,  DOI: 10.1063/1.3678309
  154. 154
    Marzari, N.; Mostofi, A. A.; Yates, J. R.; Souza, I.; Vanderbilt, D. Maximally localized Wannier functions: Theory and applications. Rev. Mod. Phys. 2012, 84, 14191475,  DOI: 10.1103/revmodphys.84.1419
  155. 155
    Sathe, P.; Harper, F.; Roy, R. Compactly supported Wannier functions and strictly local projectors. J. Phys. A: Math. Theor. 2021, 54, 335302,  DOI: 10.1088/1751-8121/ac1167
  156. 156
    Stukowski, A. Visualization and analysis of atomistic simulation data with OVITO–the open visualization tool. Modell. Simul. Mater. Sci. Eng. 2009, 18, 015012,  DOI: 10.1088/0965-0393/18/1/015012
  157. 157
    Baker, J.; Kessi, A.; Delley, B. The generation and use of delocalized internal coordinates in geometry optimization. J. Chem. Phys. 1996, 105, 192212,  DOI: 10.1063/1.471864
  158. 158
    Andzelm, J.; King-Smith, R. D.; Fitzgerald, G. Geometry optimization of solids using delocalized internal coordinates. Chem. Phys. Lett. 2001, 335, 321326,  DOI: 10.1016/s0009-2614(01)00030-6
  159. 159
    Fletcher, R. Practical Methods of Optimization, 2nd ed.; John Wiley & Sons: New York, USA, 1987.
  160. 160
    Pfrommer, B. G.; Côté, M.; Louie, S. G.; Cohen, M. L. Relaxation of crystals with the quasi-Newton method. J. Comput. Phys. 1997, 131, 233240,  DOI: 10.1006/jcph.1996.5612
  161. 161
    Packwood, D.; Kermode, J.; Mones, L.; Bernstein, N.; Woolley, J.; Gould, N.; Ortner, C.; Csányi, G. A universal preconditioner for simulating condensed phase materials. J. Chem. Phys. 2016, 144, 164109,  DOI: 10.1063/1.4947024
  162. 162
    Zhang, S. B.; Northrup, J. E. Chemical potential dependence of defect formation energies in GaAs: Application to Ga self-diffusion. Phys. Rev. Lett. 1991, 67, 23392342,  DOI: 10.1103/physrevlett.67.2339
  163. 163
    Janak, J. F. Proof that ∂E∂ni=ϵ in density-functional theory. Phys. Rev. B: Condens. Matter Mater. Phys. 1978, 18, 71657168,  DOI: 10.1103/physrevb.18.7165
  164. 164
    Sanna, S.; Frauenheim, T.; Gerstmann, U. Validity of the Slater-Janak transition-state model within the LDA + U approach. Phys. Rev. B: Condens. Matter Mater. Phys. 2008, 78, 085201,  DOI: 10.1103/physrevb.78.085201
  165. 165
    Skylaris, C.-K.; Haynes, P. D.; Mostofi, A. A.; Payne, M. C. Using ONETEP for accurate and efficient density functional calculations. J. Phys.: Condens. Matter 2005, 17, 57575769,  DOI: 10.1088/0953-8984/17/37/012
  166. 166
    Skylaris, C.-K.; Haynes, P. D. Achieving plane wave accuracy in linear-scaling density functional theory applied to periodic systems: A case study on crystalline silicon. J. Chem. Phys. 2007, 127, 164712,  DOI: 10.1063/1.2796168
  167. 167
    Ruiz-Serrano, Á.; Skylaris, C.-K. A variational method for density functional theory calculations on metallic systems with thousands of atoms. J. Chem. Phys. 2013, 139, 054107,  DOI: 10.1063/1.4817001
  168. 168
    Tinte, S.; Stachiotti, M. G.; Rodriguez, C. O.; Novikov, D. L.; Christensen, N. E. Applications of the generalized gradient approximation to ferroelectric perovskites. Phys. Rev. B: Condens. Matter Mater. Phys. 1998, 58, 1195911963,  DOI: 10.1103/physrevb.58.11959
  169. 169
    Nakamatsu, H.; Adachi, H.; Ikeda, S. Electronic structure of the valence band for perovskite-type titanium double oxides studied by XPS and DV-Xα cluster calculations. J. Electron Spectrosc. Relat. Phenom. 1981, 24, 149159,  DOI: 10.1016/0368-2048(81)80002-3
  170. 170
    Kowalczyk, S. P.; McFeely, F. R.; Ley, L.; Gritsyna, V. T.; Shirley, D. A. The electronic structure of SrTiO3 and some simple related oxides (MgO, Al2O3, SrO, TiO2. Solid State Commun. 1977, 23, 161169,  DOI: 10.1016/0038-1098(77)90101-6
  171. 171
    Hine, N. D. M.; Avraam, P. W.; Tangney, P.; Haynes, P. D. Linear-scaling density functional theory simulations of polar semiconductor nanorods. J. Phys.: Conf. Ser. 2012, 367, 012002,  DOI: 10.1088/1742-6596/367/1/012002
  172. 172
    Nono Tchiomo, A. P.; Babu-Geetha, G.; Carleschi, E.; Ngabonziza, P.; Doyle, B. P. Surface characterization of clean SrTiO3(100) substrates by x-ray photoelectron spectroscopy. Surf. Sci. Spectra 2018, 25, 024001,  DOI: 10.1116/1.5041734
  173. 173
    Fuka, V. PoisFFT – A free parallel fast Poisson solver. Appl. Math. Comput. 2015, 267, 356364,  DOI: 10.1016/j.amc.2015.03.011
  174. 174
    Budiardja, R. D.; Cardall, C. Y. Parallel FFT-based Poisson solver for isolated three-dimensional systems. Comput. Phys. Commun. 2011, 182, 22652275,  DOI: 10.1016/j.cpc.2011.05.014
  175. 175
    Womack, J. C.; Anton, L.; Dziedzic, J.; Hasnip, P. J.; Probert, M. I. J.; Skylaris, C.-K. DL_MG: A Parallel Multigrid Poisson and Poisson-Boltzmann Solver for Electronic Structure Calculations in Vacuum and Solution. J. Chem. Theory Comput. 2018, 14, 14121432,  DOI: 10.1021/acs.jctc.7b01274
  176. 176
    Kröger, F. A.; Vink, H. J. Solid State Physics; Seitz, F., Turnbull, D., Eds.; Academic Press: New York, USA, 1956; Vol. 3, pp 307435.
  177. 177
    Choi, G. M.; Tuller, H. L. Defect structure and electrical properties of single-crystal. Ba0.03Sr0.97TiO3. J. Am. Ceram. Soc. 1988, 71, 201205,  DOI: 10.1111/j.1151-2916.1988.tb05848.x
  178. 178
    Varley, J. B.; Janotti, A.; Van de Walle, C. G. Hydrogenated vacancies and hidden hydrogen in SrTiO3. Phys. Rev. B: Condens. Matter Mater. Phys. 2014, 89, 075202,  DOI: 10.1103/physrevb.89.075202
  179. 179
    Lee, L. P.; Cole, D. J.; Payne, M. C.; Skylaris, C.-K. Natural bond orbital analysis in the ONETEP code: Applications to large protein systems. J. Comput. Chem. 2013, 34, 429444,  DOI: 10.1002/jcc.23150
  180. 180
    Reed, A. E.; Weinstock, R. B.; Weinhold, F. Natural population analysis. J. Chem. Phys. 1985, 83, 735746,  DOI: 10.1063/1.449486
  181. 181
    Gomez, M. A.; Griffin, M. A.; Jindal, S.; Rule, K. D.; Cooper, V. R. The effect of octahedral tilting on proton binding sites and transition states in pseudo-cubic perovskite oxides. J. Chem. Phys. 2005, 123, 094703,  DOI: 10.1063/1.2035099
  182. 182
    Eriksson Andersson, A. K.; Selbach, S. M.; Grande, T.; Knee, C. S. Thermal evolution of the crystal structure of proton conducting BaCe0.8Y0.2O3−δ from high-resolution neutron diffraction in dry and humid atmosphere. Dalton Trans. 2015, 44, 1083410846,  DOI: 10.1039/c4dt03948c
  183. 183
    Sheppard, D.; Terrell, R.; Henkelman, G. Optimization methods for finding minimum energy paths. J. Chem. Phys. 2008, 128, 134106,  DOI: 10.1063/1.2841941
  184. 184
    Zhou, Z.; Chu, D.; Cazorla, C. Ab initio description of oxygen vacancies in epitaxially strained SrTiO3 at finite temperatures. Sci. Rep. 2021, 11, 11499,  DOI: 10.1038/s41598-021-91018-4