Toward an efficient f-in-core/f-in-valence switchable description for DFTB calculations of Ce 4f states in ceria

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₂, ceria) has become something of a prototypical material, with its signature Ce 4f 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³ to solid-oxide fuel cells,⁴ and nanomedicine.^5,6 Most of these applications utilize the ability of Ce to relatively easily toggle between the Ce³⁺ and Ce⁴⁺ oxidation states, thereby creating spatially localized “redox hotspots.”⁷ 

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.⁸ 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).⁹ 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.¹⁰ 

Despite different mechanistic pictures, i.e., with and without the formation of explicit oxygen vacancies, these examples highlight the Ce 4f state as a key parameter (descriptor) for the stability and reactivity of ceria. Of particular importance is the eigenvalue position of the occupied Ce 4f state relative to the valence band maximum (VBM), which mainly consists of O 2p electronic states. An accurate description of the Ce 4f 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 4f 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.¹¹ 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.¹² 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 4f states. This allows us in advance to decide where to place the Ce³⁺ 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³⁺ 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 4f 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.¹³ To allow for consistency between an explicit and implicit description of the Ce 4f 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³⁺ and Ce⁴⁺ and in-valence Ce³⁺ and Ce⁴⁺ 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.

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.

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, $Exc n$ is the one that includes many unknowns and, in general, determines the accuracy of the DFT functional.

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₀. 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₀] 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 4f 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.²³ 

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.²⁴ To treat the highly correlated Ce 4f electrons, we used the rotationally invariant Hubbard correction of Dudarev et al.²⁵ using a value of 5 eV for Ce 4f states, which has previously been deemed appropriate in Refs. 19 and 22 and other previous studies.¹⁸ 

All calculations were performed using the Vienna Ab initio Simulation Package (VASP)^26–29 using projector augmented wave (PAW) pseudopotentials.³⁰ Cerium was described with 12 valence electrons (5s²5p⁶6s²5d¹4f¹) and oxygen with six valence electrons (2s²2p⁴). 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 4f 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.³¹ The geometry optimizations were performed until the maximum force on each atom was lower than 10⁻⁴ a.u. The self-consistent charge convergence criteria (SCC tolerance) that we used was either 10⁻⁶ (all Ce atoms treated with f-in-core) or 10⁻³ (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.²³ 

This subsection presents an overview of the scheme that we have used to generate a new and consistent set of SK-tables for CeO₂ involving an in-core description of the empty and occupied 4f states. Our reference model structure for the occupied case will be a single oxygen vacancy in bulk ceria with different Ce³⁺ localization patterns. Figure 1 shows a representation of the CeO₂ 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⁴⁺ to Ce³⁺. 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³⁺ and Ce⁴⁺ 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−in−core4+$⁠, and $Cef−in−core3+$⁠.

• 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 4f wavefunction and its corresponding eigenvalue (its location in the bandgap).

• Step 3 (Sec. III D). Fit the repulsive energy contributions E_rep[n₀] (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.¹⁴ 

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 4f 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.

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 4f 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.³² 

The two new Ce_f-in-core types correspond to the one with a valence of 4 but with the empty 4f electronic levels treated as core and the one with a valence of 3, again treating the 4f 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 4f states here, the atomic references for $Cef−in−core4+$ and $Cef−in−core3+$ are [Xe]4f⁰5d²6s² and [Xe]4f¹5d¹6s², 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 4f 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₂ system is not appreciably affected by the in-core $Cef−in−core4+$ model. Figure 2 shows the resulting electronic density of states (DOS) for stoichiometric bulk ceria where all Ce ions are described with the $Cef−in−core4+$ 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 4f states). We conclude that, overall, the electronic structure is not much affected by locking the unoccupied f orbitals into the core.

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 4f-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³⁺ 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³⁺ 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 4f state in a similar place as PBE+U (U = 5 eV).²² These 4f 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 4f 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 4f 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).²² 

We end this subsection by commenting on the fact that the PBE+U (and SCC-DFTB+U) eigenvalues of the occupied Ce 4f band are underestimated when compared to the rather large span of experimentally measured values for the location of the occupied Ce 4f 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.

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 4f 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−in−core4+$–O pairs.

• Isotropic bulk scans for bulk CeO₂ 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₂(111) surface slab (11 structures).

• Optimized slab geometries for the CeO₂(111), CeO₂(110), and CeO₂(100) surfaces as well as a reconstructed CeO₂(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₂ 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₂, as shown in Fig. 5. All phonon spectra were calculated using phonopy,³⁴ and the reference structure used in these calculations was a 2 × 2 × 2 supercell of the optimized PBE+U structure.

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.

In this section, we again (as in Sec. III B) focus on the Ce³⁺ 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³⁺ 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⁴⁺ ions with the f-in-core SK-files, and for the sites where we chose to place the Ce³⁺ 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.

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³⁺ 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₂ bulk di-vacancies and all unique configurations of Ce³⁺ 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 4f-levels in the CeO₂ bandgap.

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 4f 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 4f 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 4f 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 4f 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₂ and inspire the generation of efficient DFTB models for related systems.

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.

The authors have no conflicts to disclose.

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).

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