A computational protocol is developed for efficient studies of partially reduced redox-active oxides using the self-consistent charge density functional tight-binding method. The protocol is demonstrated for ceria, which is a prototypical reducible oxide material. The underlying idea is to achieve a consistent (and harmonized) set of Slater–Koster (SK) tables with connected repulsive potentials that enable switching on and off the in-valence description of the Ce 4*f* states without serious loss of accuracy in structure and energetics. The implicit treatment of the Ce 4*f* states, with the use of *f*-in-core SK-tables, is found to lead to a significant decrease in computational time. More importantly, it allows for explicit control of the oxidation states of individual Ce atoms. This makes it possible to “freeze” the electronic configuration, thereby allowing the exploration of the energetics for various meta-stable configurations. We anticipate that the outlined strategy can help to shed light on the interplay between the size, shape, and redox activity for nanoceria and other related materials.

## I. INTRODUCTION

The Density Functional Theory (DFT) is one of the pillars of modern computational materials science. The success of the DFT method rests on its favorable scaling behavior and its (often) respectable accuracy. However, the method is still too computationally demanding to address the complexity of many materials under realistic conditions. Therefore, the method often serves as a reference method in modern multi-scale simulation protocols.

One example of a challenging—and interesting—group of materials is the reducible oxides with their strongly correlated electronic *d* and *f* states.^{1,2} Accurate treatment of these systems requires many-body theory (higher-order quantum mechanical methods) or, within DFT, hybrid (non-local exchange) functionals, together with efficient minimization algorithms to ensure proper convergence to the correct electronic ground state. These challenges stimulate research activities with the aim of developing new methods or finding schemes and protocols for electronic structure calculations with improved computational efficiency.

Among the reducible oxides, cerium dioxide (CeO_{2}, ceria) has become something of a prototypical material, with its signature Ce 4*f* states being populated in the reduced state. The interest in ceria is not only fundamental but also technical with a multitude of useful applications in fields ranging from catalysis^{3} to solid-oxide fuel cells,^{4} and nanomedicine.^{5,6} Most of these applications utilize the ability of Ce to relatively easily toggle between the Ce^{3+} and Ce^{4+} oxidation states, thereby creating spatially localized “redox hotspots.”^{7}

The redox activity of ceria is often correlated with its ability to create and heal oxygen vacancies. In fact, one of the most important applications of ceria is as an oxygen storage (buffer) material in the three-way catalyst.^{8} However, other reaction mechanisms, without the inclusion of explicit oxygen vacancies, have also been shown to be important. One example is the enhanced low-temperature redox activity found for very small (d < 5 nm) ceria nanoparticles (NPs).^{9} The observed low-temperature redox activity could effectively be explained through a mechanism where the redox active centers are intrinsic nanospecific motifs, rather than being accompanied by explicit oxygen vacancies.^{10}

Despite different mechanistic pictures, i.e., with and without the formation of explicit oxygen vacancies, these examples highlight the Ce 4*f* state as a key parameter (descriptor) for the stability and reactivity of ceria. Of particular importance is the eigenvalue position of the occupied Ce 4*f* state relative to the valence band maximum (VBM), which mainly consists of O 2*p* electronic states. An accurate description of the Ce 4*f* electronic state is central in the approximate electronic structure computational protocol presented in this paper.

It is well known that the effects of the strongly correlated Ce 4*f* states are difficult to capture using semi-local DFT due to the so-called self-interaction error. This is also an inherent drawback in many methods derived from DFT, such as the self-consistent charge density functional-based tight-binding (SCC-DFTB) method.^{11} However, the SCC-DFTB method is significantly less computationally demanding than DFT, and the method has been shown to be relatively accurate in the description of correlated electronic states when used in combination with a Hubbard + U approach.^{12} In this work, we develop an approach (a protocol) aiming at an efficient prediction of structural, energetic, and electronic properties of reduced ceria systems within the SCC-DFTB methodology. The underlying idea is to achieve a consistent (and harmonized) methodology that enables switching on and off the explicit (i.e., in-valence) description of the Ce 4*f* states. This allows us in advance to decide where to place the Ce^{3+} ions and perform structural relaxations for meta-stable electronic configurations. Then, by switching to the *f*-in-valence description, the protocol allows for efficient exploration of the occupied *f* states of many structural configurations in a controlled manner. For large systems, say, reasonably sized NPs, it could allow us to explore and predict the reactivity and stability of different Ce^{3+} patterns originating from oxygen vacancies or as a consequence of the shape of the particle as discussed in Ref. 12.

In this work, the switching mentioned above is achieved by using a “*f*-in-core” approach, i.e., by treating the correlated Ce 4*f* electronic states either as inactive core electrons or as active valence electrons in our Slater–Koster (SK) tables. The approach is similar to the *f*-in-core pseudopotentials for actinides and lanthanides developed by Weigand *et al.*^{13} To allow for consistency between an explicit and implicit description of the Ce 4*f* electrons, we make use of the recently developed Curvature-Constrained-Splines (CCS) method^{14,15} for efficient generation of the so-called repulsive potentials that constitute an important part of the DFTB energy expression. The CCS method ensures that a harmonized description is obtained when both in-core Ce^{3+} and Ce^{4+} and in-valence Ce^{3+} and Ce^{4+} ions are treated in one and the same system. The reduced system in focus in this work is bulk ceria with a single oxygen vacancy. This rather simple system allows us to effectively demonstrate and validate the protocol in depth.

The outline of this paper is as follows. Section II gives a short summary of the SCC-DFTB method (relating it to DFT) and provides the computational details. Section III discusses the computational protocol in detail by going through the various building blocks. In Sec. IV, the results are reported and discussed. Section V concludes this paper.

## II. THEORY AND COMPUTATIONAL DETAILS

The method in focus in this study is SCC-DFTB (see, e.g., Refs. 11, 16, and 17 and references therein). The current section will review the main ingredients of the SCC-DFTB method. Additionally, some features of the CCS regression approach will be discussed as it is used to generate repulsive potentials in this work.

### A. SCC-DFTB formalism

^{17}):

In the equation above, *f*_{α} is the occupation of a single-particle state *ψ*_{α} and *n* [or more correctly *n*(**r**), where (**r**) has been left out for space-saving reasons] is the electron density. The first term within the bracket is the kinetic energy of non-interacting electrons in the system, and the second (*V*_{ext}) is the electron–nuclei interaction. The first term outside the bracket is the Hartree (or electron–electron Coulomb) energy contribution formed by the Hartree potential V_{H} acting on the electron density, the second term is the exchange–correlation energy term, and the third one is the nuclei–nuclei repulsive interaction. Among these terms, $Excn$ is the one that includes many unknowns and, in general, determines the accuracy of the DFT functional.

*δn*, around a reference electronic density,

*n*

_{0}, which is usually taken to be equal to the superposition of atomic densities. The SCC-DFTB energy expression equivalent to Eq. (1) is

The first term, *E*_{BS}, called the band structure term, corresponds to the bracketed expression in Eq. (2). It is equal to the sum of the occupied orbital energies of the reference electron density *n*_{0}. The second term, the so-called Coulomb term *E*_{Coul}, consists of the first two terms outside the brackets in Eq. (2). This term handles charge fluctuations from the reference density. Finally, the repulsive term *E*_{rep} is the part that contains all other interactions that were not yet taken into account. It corresponds to the last four terms in Eq. (2) where the ion–ion interactions E_{II} constitute the dominant part.

In practice, *E*_{rep}[*n*_{0}] is approximated and expressed as a sum of pairwise two-body terms (*V*_{rep}). This term is conceptually the most problematic one, and the parametrization procedure (discussed further in Sec. III D) is often not trivial.

For strongly correlated electrons, in general, such as for 4*f* electrons in Ce, the SCC-DFTB method suffers from the same problems as its parent DFT method. This is the tendency to delocalize such states due to the well-known self-interaction problem, resulting in an incorrect description. In the literature, there are numerous studies dedicated to benchmarking and exploring different DFT approximations and practices to describe ceria systems.^{18–22} Herein, we will use the Hubbard + U corrected DFT as a reference owing to its wide use in the community. This implies that a Hubbard + U correction is needed in the SCC-DFTB calculations. In this work, we have used the fully localized limit (FLL) form of the Hubbard correction, as implemented by Hourahine *et al.*^{23}

### B. Computational details

** DFT calculations**. The electronic structure reference calculations presented in this work were performed using DFT in the implementation with plane waves and pseudopotentials. The exchange–correlation energy was described by the PBE density functional developed by Perdew, Burke, and Ernzerhof.

^{24}To treat the highly correlated Ce 4

*f*electrons, we used the rotationally invariant Hubbard correction of Dudarev

*et al.*

^{25}using a value of 5 eV for Ce 4

*f*states, which has previously been deemed appropriate in Refs. 19 and 22 and other previous studies.

^{18}

All calculations were performed using the Vienna *Ab initio* Simulation Package (VASP)^{26–29} using projector augmented wave (PAW) pseudopotentials.^{30} Cerium was described with 12 valence electrons (5*s*^{2}5*p*^{6}6*s*^{2}5*d*^{1}4*f*^{1}) and oxygen with six valence electrons (2*s*^{2}2*p*^{4}). The VASP code by default only accounts for the spherical contribution to the exchange–correlation energy gradient inside the PAW spheres. Aspherical contributions are essential for accurate total energies and band structures for elements with 4*f* electrons. They were accounted for in our calculations.

The cutoff energy for the plane wave basis was set to 400 eV, and the calculations were performed with Γ-point sampling of the Brillouin zone. The maximum forces on any atom were converged to 0.02 eV/Å.

** DFTB calculations**. All SCC-DFTB calculations were performed using the DFTB+ code.

^{31}The geometry optimizations were performed until the maximum force on each atom was lower than 10

^{−4}a.u. The self-consistent charge convergence criteria (SCC tolerance) that we used was either 10

^{−6}(all Ce atoms treated with

*f*-in-core) or 10

^{−3}(

*f*states described as valence on some Ce atoms).

In addition, we employed a Methfessel–Paxton Brillouin zone sampling method, which produces a Fermi-like electron distribution with much lower electron entropy. This was done to ensure proper electron occupancies. The Hubbard correction was also introduced in our DFTB calculations. We used the fully localized limit (FLL) form of the Hubbard correction.^{23}

## III. BUILDING THE WORKFLOW TOWARD CE IN-CORE/IN-VALENCE HARMONIZATION

### A. Overview

This subsection presents an overview of the scheme that we have used to generate a new and consistent set of SK-tables for CeO_{2} involving an in-core description of the empty and occupied 4*f* states. Our reference model structure for the occupied case will be a single oxygen vacancy in bulk ceria with different Ce^{3+} localization patterns. Figure 1 shows a representation of the CeO_{2} bulk structure alongside an illustration depicting the local environment around an oxygen vacancy. With a single oxygen vacancy in the structure, two nearby Ce ions are reduced from Ce^{4+} to Ce^{3+}. These tend to occupy two of the Ce positions indicated in Fig. 1(b). Their exact location will be discussed in more detail below.

The purpose is to design a protocol that allows us to treat the strongly localized *f* states associated with the Ce^{3+} and Ce^{4+} cations either implicitly, using an in-core approach, or explicitly, as regular valence electronic state. To accomplish this, we follow these steps:

Step 1 (Sec. III B). Generate the SK-tables and the corresponding SCC-related parameters for the E

_{BS}and E_{Coul}energy terms in Eq. (3) for three types of Ce atoms: Ce_{f-in-valence}, $Cef\u2212in\u2212core4+$, and $Cef\u2212in\u2212core3+$.Step 2 (Sec. III C). Optimize the U-parameter for the

*f*-states (through linear fitting) in search of a good description of its occupied 4*f*wavefunction and its corresponding eigenvalue (its location in the bandgap).Step 3 (Sec. III D). Fit the repulsive energy contributions E

_{rep}[*n*_{0}] (as part of the total DFTB total energy) against a training set of key ceria structures and their accompanying energies at the DFT level using the robust CCS method.^{14}

Thus, first, we compute the necessary integrals and parameters (as described in Sec. II) to generate SK-tables that offer a good description in terms of the electronic density of states (DOS). In the second step, we optimize the U parameter to ensure that we obtain a good description of the occupied Ce 4*f* states, both in terms of spatial distribution and energy. In the last step, we fit the repulsive potential to ensure that we reproduce the reference structure and key energetic properties.

### B. Step 1: Generating *f*-in-core SK-tables

The SCC-DFTB formalism requires an SK-table to be generated for each unique atom pair. The “connectivity table” in Table I lists the pair types relevant to the present study and gives references to the origin of the respective SK-table (i.e., to the entailing parametrizations, for example, of the compression radii of the atomic wave functions, and more). For consistency with previous calculations, we use the established mio-1.0 set of Ref. 11 for O–O and the Ce–O parameters of Ref. 12, where the 4*f* electrons are treated in-valence. To account for the two new, *f*-in-core, charge states of Ce, we have generated seven new SK-tables, in addition to the previous three.^{32}

. | O . | Ce . | $Cef\u2212in\u2212core3+$ . | $Cef\u2212in\u2212core4+$ . |
---|---|---|---|---|

O | mio-1.0^{11} | Reference 12 | ✓ | ✓ |

Ce | Reference 12 | Reference 12 | ✓ | ✓ |

$Cefin\u2212core3+$ | ✓ | ✓ | ✓ | ✓ |

$Cefin\u2212core4+$ | ✓ | ✓ | ✓ | ✓ |

The two new Ce_{f-in-core} types correspond to the one with a valence of 4 but with the empty 4*f* electronic levels treated as core and the one with a valence of 3, again treating the 4*f* electronic levels as core, but now with an occupation of 1 electron. The *l*-dependent compression radii used herein are the same as the one used in Ref. 12. However, as we do not explicitly treat the Ce 4*f* states here, the atomic references for $Cef\u2212in\u2212core4+$ and $Cef\u2212in\u2212core3+$ are [Xe]4*f*^{0}5*d*^{2}6*s*^{2} and [Xe]4*f*^{1}5*d*^{1}6*s*^{2}, respectively. The SK-tables involving oxygen are kept the same as in Ref. 12.

The primary aim of our paper is to explore to what extent it is feasible to give an in-core description of the empty and/or occupied Ce 4*f* states in partially reduced ceria within the SCC-DFTB framework. However, here we first want to ascertain that the electronic structure of stoichiometric bulk CeO_{2} system is not appreciably affected by the in-core $Cef\u2212in\u2212core4+$ model. Figure 2 shows the resulting electronic density of states (DOS) for stoichiometric bulk ceria where all Ce ions are described with the $Cef\u2212in\u2212core4+$ model. The structure was fixed at an optimized PBE+U (U = 5 eV) structure (with a cell parameter of 5.50 Å). The corresponding DFT-computed DOS is shown for reference (excluding the Ce 4*f* states). We conclude that, overall, the electronic structure is not much affected by locking the unoccupied *f* orbitals into the core.

### C. Step 2: Parametrizing electron localization (Hubbard U)

Next, we aim to find the U-value that best reproduces the DOS for reduced ceria compared to data from the PBE+U reference method. First, we investigate the effect the U-value has on the electronic eigenvalue of the occupied 4*f*-states of ceria within the SCC-DFTB framework. We also look at the electronic spin density to ensure that the electrons localize on the designated Ce ions. Partially reduced ceria, CeO_{2−x}, with one oxygen vacancy and two different localization patterns for the two excess *f* electrons will be used as test cases. Both correspond to removing a single oxygen atom from a 4 × 4 × 4 super-cell. The structures differ with respect to where the two associated Ce^{3+} ions are located [cf. Fig. 1(b)]. In the nearest-neighbor (NN+NN) configuration, the two reduced Ce ions are both located on the inner rim around the vacancy of Fig. 1(b). In the next-nearest-neighbor (NNN+NNN) configuration, the two Ce^{3+} ions are both located on the outer rim.

We performed single-point SCC-DFTB+U calculations with varying U values using the NN+NN and NNN+NNN structures from the PBE+U (U = 5 eV) optimizations in Ref. 22. In the following, we will often refer to NN+NN and NNN+NNN by just “NN” and “NNN,” respectively. For each localization pattern, a SCC-DFTB+U scan was performed in search of a U value that would position the occupied 4*f* state in a similar place as PBE+U (U = 5 eV).^{22} These 4*f* positions are 1.22 and 1.29 eV above the VBM for the NN and NNN cases, respectively. These values are marked with red and blue horizontal dashed lines in Fig. 3(a). With SCC-DFTB+U, the position of the occupied 4*f* energy eigenvalues in the bandgap is seen to display a linear dependence with U (the solid tilted lines in the figure; in fact, a linear behavior was obtained for PBE+U as well in Ref. 22).

Moreover, Fig. 3 shows that it is possible to find a U-value within the SCC-DFTB+U framework that mimics the PBE+U placement of the occupied *f*-eigenvalue in the bandgap. The optimal U-values calculated from the intercept between the solid and dashed lines in Fig. 3 are 1.95 eV for NN and 2.03 eV for NNN. We will use the average of these two U-values, U = 1.99 eV, which gives an O_{2p,VBM} → $Ce4focc$ bandgap of 1.21 eV for NN and 1.25 eV for NNN. This is in good agreement with the 4*f* location in the PBE+U case. Figure 3(b) shows the corresponding spin density isosurface for the NNN case computed using the SCC-DFTB+U method with an optimal U-value. Similar localization behavior is obtained for all tested U-values. The spatial distribution of the spin density seen in Fig. 3(b) is very similar to the one obtained using PBE+U (U = 5 eV).^{22}

We end this subsection by commenting on the fact that the PBE+U (and SCC-DFTB+U) eigenvalues of the occupied Ce 4*f* band are underestimated when compared to the rather large span of experimentally measured values for the location of the occupied Ce 4*f* states in the bandgap. This is interesting to know, but our objective is to parameterize the SCC-DFTB method to mimic a selected reference space, here being the PBE+U (U = 5 eV), and not to mimic experiments.

### D. Step 3: Generating the repulsive potentials

The electronic properties fitted in the previous steps have all been done on static structures geometry optimized at the PBE+U level. To ensure that the good electronic description is preserved also for unknown structures, we need make certain that the structure is as close to the reference PBE+U result as possible. In the SCC-DFTB framework, this is done by the fitting of the so called repulsive potential [see Eq. (3) above]. In this work, we use the Curvature Constrained Splines (CCS) method as developed in Ref. 14 to generate the repulsive potentials. As the name suggests, the method makes use of cubic splines with special constraints and can be utilized to obtain general pair potentials or repulsive potentials as in this work. In the following, we will provide a brief description of the repulsive potential generation. For a full tutorial of the procedure, we refer to the repository at https://github.com/Teoroo-CMC/Ceria_DFTB.git, which contains the data and tools used to produce our repulsive potentials.

The shifting of the 4*f* levels to the core will result in a discrepancy in the band structure energy contribution between the new and the old SK-tables. To ensure that the new SK-tables give the same relative energies and structural relaxations, we need to parameterize the repulsive energy contribution for each atom type combination in Table I.

We use the CCS fitting method to parameterize the repulsive potential using a training set of PBE+U energies and their accompanying structures. More precisely, the reference energy for each structure in the training set is calculated as *E*^{tot,PBE+U} − $ECoulDFTB$ − $EBSDFTB$ derived from Eq. (3). Our procedure is described in detail in Ref. 15 with application to Si. In our experience, the quality of the training set used in the parametrization has a large influence on the accuracy and transferability of the repulsive potential. In this work, we used a training set comprised of the following stoichiometric structures to obtain the repulsive potential for the Ce–O and $Cef\u2212in\u2212core4+$–O pairs.

Isotropic bulk scans for bulk CeO

_{2}in the fluorite structure using a 4 × 4 × 4 supercell, including random displacements (±0.1 Å along each xyz-direction) of internal coordinates (59 structures).Isotropic scans for a three atomic layer thin CeO

_{2}(111) surface slab (11 structures).Optimized slab geometries for the CeO

_{2}(111), CeO_{2}(110), and CeO_{2}(100) surfaces as well as a reconstructed CeO_{2}(110) surface (four structures).

In our previous efforts on the generation of repulsive potentials for oxide materials,^{12,33} we have found that it is essential to add structures with a variety of coordination numbers in the training-set. In this work, such structures are represented by the surface slabs. Note that we have omitted Ce metal and O_{2} from the training-set and only parameterize the Ce–O repulsive potentials. All Ce–Ce repulsive potentials in this work are set to zero. All structures in the training-set were weighted equally, and version 0.11.2 of the CCS package was used to obtain the repulsive potential (available at https://github.com/Teoroo-CMC/CCS). We used a cutoff radius of 10 bohrs and 40 knots in the cubic spline.

The quality of the repulsive potential can be appreciated from Figs. 4(a) and 4(b) for the in-core and in-valence cases, respectively. These figures show scatter plots between the target repulsive energy and the energy obtained from the CCS for all structures in the training-set. Both fittings are quite good, but the one involving the in-valence description is slightly better. This finding is also supported by the comparison of the calculated phonon spectra of bulk CeO_{2}, as shown in Fig. 5. All phonon spectra were calculated using phonopy,^{34} and the reference structure used in these calculations was a 2 × 2 × 2 supercell of the optimized PBE+U structure.

^{3+}using an in-core approach, we also need a repulsive potential describing the Ce

^{3+}–O interactions. In principle, we could obtain this contribution with the CCS regression approach using a training set with a sufficient number of structurally diverse reduced structures. However, we found that it is sufficient to adopt a simpler approach where the Ce

^{3+}–O repulsive potential, $VrepCe3+\u2212O$, is generated from the Ce

^{4+}–O repulsive potential, $VrepCe4+\u2212O$, by scaling according to the following expression:

*α*is a scalar.

The scaling parameter *α* was optimized to achieve a minimal relaxation around the NN vacancy structure. *α* was adjusted in a range from 0.8 to 1.5 in steps of 0.01, and we found that an *α* value of 1.30 gave the smallest relaxation energy, cf. Fig. 4(c). The starting point, and reference point, in this relaxation was the optimized PBE+U geometry. The degree of structural distortion caused by the vacancy will be discussed in Sec. IV, as well as other quality measures.

## IV. RESULTS OF THE HARMONIZATION TESTS

In this section, we again (as in Sec. III B) focus on the Ce^{3+} species in bulk ceria with one oxygen vacancy but now from a different perspective. By “harmonization” (cf. the heading), we refer to the construction of a Ce description that is consistent in such a way that going from an *f*-in-valence description to *f*-in-core will not cause drastic changes in computed properties. Consistency in terms of structure, electronic band gaps, and relative stabilities of different Ce^{3+} localization patterns.

To test the harmonization of the different SK-files (SK-file = SK-tables + repulsive), i.e., how interchangeable they are, we performed a detailed assessment of electronic, structural, and energetic properties of the NN and NNN electron localizations around one oxygen vacancy in our 4 × 4 × 4 ceria bulk model. We used either the *f*-in-core SK-files for all Ce atoms or employed a “*mixed*” *approach*, treating all Ce^{4+} ions with the *f*-in-core SK-files, and for the sites where we chose to place the Ce^{3+} ions, the *f*-in-valence SK-files were used.

Ideally, we would like to be able to determine accurate structures (“geometries”) and energies at the least demanding computational level. This would here correspond to the SCC-DFTB in-core approach. It would also be very useful if the electronic structure could be determined from a single-point calculation with an in-valence description only at those sites where we have *f*-electrons. As before, our reference is the PBE+U geometry, energy, and electronic structure. Therefore, in order to examine how close to the ideal PBE+U results our SCC-DFTB+U calculations come, we compare five different workflows (“cases”) as summarized in the left column of Table II and elaborated below. Case 4 would refer to our ideal “least computationally demanding” case.

Computational details . | $Egap(O2p\u2192Ce4focc)$ (NN,NNN) . | $d\u0304Ce3+\u2212O2\u2212$ (NN,NNN) . | ΔE (V_{NN} − V_{NNN})
. |
---|---|---|---|

PBE+U Ref. 22 | (1.22, 1.29) | (2.426, 2.473) | 0.07 |

HSE06’ Ref. 22 | (2.62, 2.65) | (−, −) | 0.11 |

PBE (in-core) | ⋯ | (2.429, 2.485) | 0.17 |

DFTB, case 1: (Ce^{4+} in-core, Ce^{3+} in-val, SP@GGA-geom) | (1.19, 1.31) | ⋯ | 0.12 |

DFTB, case 2: (Ce^{4+} in-core, Ce^{3+} in-core, SP@GGA-geom) | ⋯ | ⋯ | 0.15 |

DFTB, case 3: (Ce^{4+} in-core, Ce^{3+} in-core, geom opt) | ⋯ | (2.423, 2.477) | 0.15 |

DFTB, case 4: (Ce^{4+} in-core, Ce^{3+} in-val, SP@Case 3-geom) | (1.15, 1.28) | ⋯ | 0.20 |

DFTB, case 5: (Ce^{4+} in-core, Ce^{3+} in-val, geom opt) | (1.08, 1.07) | (2.427, 2.498) | 0.12 |

Computational details . | $Egap(O2p\u2192Ce4focc)$ (NN,NNN) . | $d\u0304Ce3+\u2212O2\u2212$ (NN,NNN) . | ΔE (V_{NN} − V_{NNN})
. |
---|---|---|---|

PBE+U Ref. 22 | (1.22, 1.29) | (2.426, 2.473) | 0.07 |

HSE06’ Ref. 22 | (2.62, 2.65) | (−, −) | 0.11 |

PBE (in-core) | ⋯ | (2.429, 2.485) | 0.17 |

DFTB, case 1: (Ce^{4+} in-core, Ce^{3+} in-val, SP@GGA-geom) | (1.19, 1.31) | ⋯ | 0.12 |

DFTB, case 2: (Ce^{4+} in-core, Ce^{3+} in-core, SP@GGA-geom) | ⋯ | ⋯ | 0.15 |

DFTB, case 3: (Ce^{4+} in-core, Ce^{3+} in-core, geom opt) | ⋯ | (2.423, 2.477) | 0.15 |

DFTB, case 4: (Ce^{4+} in-core, Ce^{3+} in-val, SP@Case 3-geom) | (1.15, 1.28) | ⋯ | 0.20 |

DFTB, case 5: (Ce^{4+} in-core, Ce^{3+} in-val, geom opt) | (1.08, 1.07) | (2.427, 2.498) | 0.12 |

*data column*in Table II reports the occupied

*f*level energy in terms of the O

_{2p}→ $Ce4focc$ bandgap for each of the two Ce

^{3+}ions in the NN and NNN localization cases. The first

*data rows*in Table II report these bandgaps calculated with our reference DFT method (PBE+U). The HSE06 values are included to provide perspective, namely, to highlight that the variation between the five Cases, i.e., when the in-core/in-valence

*f*description is being varied, should be considered in relation to the (well-known) notion that varying the DFT functional can lead to very large shifts of the bandgap values.

**Case 1**refers to single-point SCC-DFTB calculations on the DFT geometries with the*f*-in-core description for the Ce^{4+}ions and*f*-in-valence description for the Ce^{3+}ions (i.e., the “mixed” approach, like in Cases 4 and 5 below). We compare the quality of the resulting electronic structures and of ΔE (V_{NN}− V_{NNN}) with the DFT results. The SCC-DFTB+U*f*-level energies compare very well with the DFT results and ΔE (V_{NN}− V_{NNN}) is small as it is with DFT, and it has the right sign compared to DFT.**Case 2**makes use of single-point calculations at the DFT-optimized geometries of case 1, but now with SCC-DFTB calculations using the*f*-in-core description for both the Ce^{3+}and Ce^{4+}ions. Since the structure is fixed and the*f*electrons are in-core, our only performance probe here is the quality of the energy description for ΔE (V_{NN}− V_{NNN}). The PBE+U value is small. The SCC-DFTB computed value is of the right sign but slightly overestimated.**Case 3**has the same setup as case 2 in terms of the*f*electron description (*f*-in-core for all Ce), but now the calculations also allow for structural relaxation. A measure of the relaxation is provided in the*second data column*of Table II) in terms of the average distance, $d\u0304Ce3+\u2212O2\u2212$, between Ce^{3+}and all oxygen ions within their first coordination shells. Figure 6 illustrates the relaxation in more detail. In general, the Ce–O distances around the Ce^{3+}ions are longer than the equilibrium Ce^{4+}–O in the stoichiometric bulk. However, in the DFT V_{NN}geometry, two of Ce^{3+}–O distances, related to oxygen ions located on the opposite side of vacancy, are significantly shorter than the equilibrium bulk value.This somewhat unexpected behavior could perhaps be understood in terms of a combination of two interactions that can counterbalance the expected outward relaxation of the Ce

^{3+}–O bond, namely, (i) the repulsion between the vacancy and Ce^{3+}ion, which would push the Ce^{3+}ion closer to the oxygen ion on the opposite side of the vacancy, and (ii) the attraction between the vacancy and the oxygen io,n which would pull the oxygen closer to the Ce^{3+}ion.This effect is smaller in the case 3 geometry where the “short” Ce

^{3+}–O bond is only very slightly smaller than the corresponding equilibrium value bulk value.We conclude that the PBE+U distortion pattern is quite well reproduced by the scc-DFTB in-core description while maintaining a reasonable ΔE (V

_{NN}− V_{NNN}) value (cf. Table II).

**Case 4**corresponds to a single-point calculation using the geometries obtained in Case 3 but now with a*f*-in-valence description for the Ce^{3+}ions. Here, we can compare the quality of the electronic structure, including the possible effects caused by changes in geometry, with the results obtained when the reference DFT structures were used (case 1). We find that the calculated occupied Ce 4*f*eigenvalues in our SCC-DFTB approach change by less than 0.05 eV (and remain in good agreement with the DFT reference) when based on the structures that were fully relaxed using the in-core description for the Ce ions. The relative stability of V_{NN}and V_{NNN}is bit worse.**Case 5**makes use of geometry relaxation with the*f*-in-core description for Ce^{4+}ions and*f*-in-valence description for Ce^{3+}ions. The aim of this calculation is to compare the relaxation around Ce^{3+}ions with the in-core description as well as with the reference DFT geometries and the possible effect of the new structure on the bandgaps and ΔE due to the new geometry. The consequences on the in-valence relaxation are seen in Fig. 6. The essential features of the PBE+U results are not captured as well as in the in-core model (case 3), which is somewhat counter-intuitive.Moreover, comparing case 4 with case 5, we find that geometry relaxation at the in-valence level induces a significant down-shift of the occupied Ce 4

*f*eigenvalues for the NNN localization pattern. The eigenvalues shift from 1.31 eV above the valance band using the DFT geometry (case 1) to 1.07 eV in case 5.

In summary, so far, the results suggest that our “mixed” in-core/in-valance SCC-DFTB approach well captures the electronic structure of the DFT reference. This is true both in terms of the whole DOS of the pristine bulk (excluding the *f*-band) and the position of the *f*-level eigenvalues when using an in-valance description for Ce^{3+} ions surrounding a single oxygen vacancy, provided the structure is the same or similar to the DFT reference. This condition is met using case 4, which relies on a geometry optimization using an in-core only description but with a single point calculation using the “mixed” approach. Further optimizing the geometry with the “mixed” approach (case 5) leads to a worsened description of the structure compared to the DFT reference. This is particularly seen for the V_{NNN} configuration. Concerning the energetics, the computational protocol qualitatively describes the relative stabilities of V_{NN} and V_{NNN}. However, the relative energies of these structures are overestimated in all five cases compared to the PBE+U reference.

As a further test of our approach, we considered four types of CeO_{2} bulk di-vacancies and all unique configurations of Ce^{3+} ions at NN positions around them as follows:

A di-vacancy along the [100] direction (five configurations).

A di-vacancy along the [110] direction (11 configurations).

A di-vacancy along the [111] direction with a Ce ion right between the vacancies (five configurations). We label these configurations [111]Ce.

A di-vacancy along the [111] direction without a Ce ion between the vacancies (ten configurations).

Altogether, data using the case 3 and case 4 protocols for 31 different structures have been computed. This particular choice of protocol corresponds to the ideal situation where energies and *f*-levels are both computed in the most efficient computational way without having to resort to any prior knowledge of the structure from a reference calculation at the PBE+U level. The results are given in Fig. 7 for the relative energies computed according to case 3 and the *f*-levels computed according to case 4 for all di-vacancy structures. In Figs. 6 and 7, a comparison is made to the corresponding PBE+U results. It is clear that our “mixed” in-core/in-valance DFTB approach manages to place the stability of the various di-vacancies groups [100], [110], [111]Ce, and [111] in the correct order. Furthermore, it also gives a good prediction of the corresponding localized Ce 4*f*-levels in the CeO_{2} bandgap.

## V. CONCLUSION

A computational protocol for a more efficient calculation of structural, energetic, and electronic properties of partially reduced metal-oxide systems within the SCC-DFTB+U method has been developed. The protocol is demonstrated for ceria, which is famous for its redox chemistry originating from its energetically accessible, but at the same time theoretically problematic, correlated 4*f* electronic states.

The protocol is built around the generation of consistent (and harmonized) SK-tables with corresponding repulsive potentials (here constructed using the curvature constrained splines method) and aims to enable an efficient and fairly seamless switching between an explicit and implicit description of the correlated electronic Ce 4*f* states. The switch implies that the problematic occupied *f* states are treated as core states to the largest possible extent without a serious loss of quality in the description. The latter was tested and assessed in a five-step process (“cases 1–5”) in Sec. IV. The explicit in-valence description of the Ce 4*f* states is only used when necessary, for example, when electronic properties are needed or for studying the actual redox reaction (charge transfer) in detail.

In summary, the implicit treatment of the Ce 4*f* states, with the use of *f*-in-core potentials, leads to a significant decrease in computational time (faster self-consistent charge convergence of the SCC-DFTB method). In addition, our computational protocol allows for explicit control of the oxidation states of individual Ce atoms, which makes it possible to “freeze” the electronic configuration of desired sites and explore the redox properties of others.

We anticipate that the outlined computational protocol can be used to shed light on the interplay between *f*-electrons and more complex defect structures (e.g., large vacancy clusters) in reduced CeO_{2} and inspire the generation of efficient DFTB models for related systems.

## ACKNOWLEDGMENTS

We acknowledge the support from the Swedish Research Council (VR), the Center for Interdisciplinary Mathematics (CIM) at Uppsala University, and the Swedish National Strategic e-Science programme (eSSENCE). The calculations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at UPPMAX and NSC.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Bojana Kocmaruk**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal). **Akshay Krishna Ammothum Kandy**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Kersti Hermansson**: Conceptualization (equal); Funding acquisition (equal); Writing – original draft (equal); Writing – review & editing (equal). **Jolla Kullgren**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Project administration (equal); Writing – original draft (equal); Writing – review & editing (equal). **Peter Broqvist**: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in https://github.com/Teoroo-CMC/Ceria_DFTB.git.

## REFERENCES

*U*instead of Stoner I

*U*for electron localization and structure at oxygen vacancies in ceria

*Ab initio*molecular dynamics for liquid metals

*Ab initio*molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium

*ab initio*total-energy calculations using a plane-wave basis set

*Ab initio*molecular dynamics with full wave functions