Electronic interaction $U_{pp}$ on oxygen $p$ orbitals in oxides: exemples of UO$_2$ and TiO$_2$, and role of correlated orbitals

We carry out a detailed study of the role of electronic interaction on $p$ oxygen orbitals in Mott insulators oxides and charge transfer oxide using UO$_2$ and TiO$_2$ as prototypical systems. First, we calculate values of effective interactions $U_{ff}$, $U_{pp}$ and $U_{fp}$ in UO$_2$ and $U_{dd}$, $U_{pp}$ and $U_{dp}$ in TiO$_2$. Second, we analyse the role of electronic interactions $U_{pp}$ on $p$ orbitals of oxygen on spectral and structural properties. Finally, we show that this role depends strongly on the definition of correlated orbitals and that using Wannier functions lead to more physical results for spectral and structural properties.


I. INTRODUCTION
Since its development the Density Functional Theory (DFT) [1,2] has been applied successfully to a wide variety of systems [3]. For some elements however, the Local Density Approximation (LDA) [4] or the Generalized Gradient Approximation (GGA) [5] have a tendency to overdelocalize electrons, due both to the self interaction error and to the approximate description of interaction effects. This error is particularly visible on systems containing localized d or f orbitals, such as transition metal oxides or lanthanides and actinides compound. To address this problem, and more generally to improve these functionals or the limitations of DFT concerning excited states, several methods have been designed, such as e.g. hybrid functionals, the GW method [6], DFT+U [7][8][9], DFT+DMFT [10,11].
More specifically, the description of systems that contain atoms with spatially localized atomic orbitals (such as d or f ) requires a dedicated treatment, because the large electronic interaction between electrons inside these orbitals needs to be taken into account. It can be handled at the static mean field level (in DFT+U ) or in dynamical mean field theory (DMFT).
However, these two schemes rely on two parameters, the direct interaction U and the exchange interaction J. These parameters can be adjusted such that e.g band gap or cell parameter are in agreement with experiment. But in this case, such calculation is no more first principles. A more fruitful solution is to try to calculate U and J from first principles. The main used methods to obtain U from first principles are the constrained DFT (cDFT) [7], the linear response method [12,13], method using Hartree Fock orbitals [14,15] and the constrained Random Phase Approximation (cRPA) [16]. These methods have been successful to describe a wide range of systems including pure metals and oxides. In these systems, most of the calculations took into account only the Coulomb interaction on d or f orbitals of the metallic atom (e.g. Ti, U, Fe). However, it has been shown also that other orbitals (e.g. p orbitals of oxygen), are important, as well as intershell interactions and have an impact on the electronic structure and/or structural parameters of e.g. TiO 2 [17,18], transition metal oxides [19][20][21][22],ZrO 2 [23], actinides and lanthanides oxides [24][25][26], cerium [25], lanthanides compounds [27] and high T c cuprates [28]. Such idea has also been applied to design a scheme for high throughput [15,29], as effective interactions can also be used within generalized DFT+U schemes that can be seen as simplified but fast hybrid functionals. However, only a few work focused on an ab-initio calculation of these Coulomb interaction [15,25,26,[29][30][31]. Moreover, a detailed understanding of the effect of these interactions on electronic or structural properties is needed: What is the effect of U pp on spectral functions? The inclusion of U pp often leads to a surprising decrease of atomic volume (see eg. [17,23,24,32]): what is its physical interpretation? Such basic and fundamental questions are important because the U pp interaction is a building block of recent proposed scheme for high throughput calculations [15,29,32].
Whichever group of orbitals may be chosen as correlated states, the DFT+U and DFT+DMFT need a precise definition of correlated orbitals. It has been emphasized recently that the calculation of U itself, the spectral function as well as structural properties can be impacted by this choice [33][34][35][36][37][38]. The impact of this choice on p states was only briefly discussed [29].
In this paper, we carry out a detailed study of interaction effects on oxygen atoms in oxides. We first compute interaction in the full dp (U dd , U dp and U pp ) or df states (U f f , U f p and U pp ) using cRPA implementation [33,39] with the ABINIT code [40,41]. Then we study the relative impact of interactions on spectral and structural properties. We propose an explanation for the observed unusual behavior over volume observed here and in other studies. Lastly, we emphasize the key role of the definition of correlated orbitals, namely atomic orbitals versus Wannier functions: such definition impacts quantitatively the results of the calculations as it has previously been discussed [33][34][35][36][37][38]. However, in our case, and especially for structural properties, results are opposite. We focus on two prototypical systems, UO 2 , a Mott Hubbard insulator and TiO 2 , a charge transfer insulator, both containing strongly-correlated orbital (d or f ) and oxygen p orbital close to the Fermi level.
The section II presents the DFT+U and cRPA methodologies, definition of correlated orbitals as well as computational details. Section III and IV present results on UO 2 and TiO 2 . The conclusion is in section V.

II. METHODS AND COMPUTATIONAL DETAILS
A. DFT+U

Expression of energy in DFT+U
The standard expression for DFT+U total energy is where E DFT is the energy of the system in DFT and E U is the energy due to the DFT+U correction. E U can be split in two terms, with E ee , the mean field electron-electron interaction in DFT+U and E dc the double counting correction (see below). We use the rotationally invariant expression of E ee [42] 1 where m 1 , m 2 , m 3 , m 4 are indexes of real spherical harmonics of angular momentum l, |m 1 is a generic correlated orbital, m 1 m 2 |V ee |m 3 m 4 is an element of the electron-electron interaction matrix V ee 2 and σ is the spin. n σ m1,m2 is the element m 1 , m 2 of the occupation matrix for correlated orbital, calculated as: with Ψ σ ν,k the Kohn-Sham wave function and f σ ν,k the occupation factor, for band ν, k-point k and spin σ. 1 Summations over correlated atoms and species are implied in this section 2 It is expressed using Slater integrals F k and Gaunt coefficient m 1 |m|m 2 of real spherical harmonics as 4π k=0,2,4,6 Considering a basis where the occupation matrix is diagonal (n m4,m2 = n m2 δ m4,m2 and n m3,m1 = n m1 δ m3,m1 ), we can simplify E ee expression to: The second part in E U is E dc . It can take various forms [7,43], on this study, we will focus on the Full Localized Limit (FLL) formulation. The role of this double counting correction is to cancel the local LDA electronelectron interaction. It is written with N σ being the total number of electrons for the considered orbital (the trace of the occupation matrix for spin σ) and N = σ N σ .
In this work, we use the PAW [44] implementation [45,46] of DFT+U in ABINIT [40,41]. In the next subsection we discuss the choice of correlated orbitals.

Choice of correlated orbitals in DFT+U
In this subsection, we precise several possibilities to define correlated orbitals |m 1 .
Indeed, the occupation matrix n σ m1,m2 contains the major information about localization, hybridization and orbital polarization or anisotropy in the system under study. Precising its detailed calculation is thus central: It is computed using equation 3. Inside this equation, the term |m 1 is central [33][34][35][36][37][38].
We can separate the local orbital in a radial and angular part, yielding with |Y l,m1 , a spherical harmonic accounting for the angular part, and |φ the radial part, that can take different formulation (and depends on l).
In this work, we use 2 different ways to define correlated orbitals namely, atomic orbitals and Wannier orbitals.
a. Atomic orbitals. The simplest formulation for |φ is to use an atomic local orbital |φ 0 3 . Ref. [36] showed that in the PAW scheme and when atomic orbitals are truncated, the best way to describe atomic p orbitals and especially the atomic limit, is to use a re-normalized scheme for truncated orbitals (the discussion is summarized in appendix B). In the notation of [36], this scheme is called o C . If not precised explicitly in the text, this scheme will be used in this paper as the "atomic" one. Note that, as discussed in [36], unrenormalized truncated orbitals (this is the scheme o B detailed in appendix B) commonly used in PAW [46,48] must not be used for orbitals that are not fully localized in the PAW sphere, such as p orbitals.
b. Wannier orbitals Wannier Functions can be seen as local orbitals adapted to the solid (this scheme will be noted o W ). In this work we use Projected Localized Orbital Wannier functions [33,47] defined in two steps. In the first step, we built functions |χ lσ km that are not normalized and not orthogonal, by projecting atomic orbitals over Kohn-Sham wave-functions as: with |Ψ σ kν being Kohn-Sham wave function and |χ l km an atomic-like orbital. As the sum is limited to a subset of bands (W),χ lσ km are not orthonormal. In a second step, we thus proceed to an orthonormalization using the overlap matrix [47]. After orthonormalization, we obtain a set of Wannier functions|W lσ km . The choice of the subset of bands W is guided by correlation effects we want to consider in our system [47]. 3. DFT+U with more than one correlated orbital At this step it is important to stress the fact that it is possible to consider more than one orbital in DFT+U . If one neglects the inter-orbital interaction [25], we have for the total energy, considering a system with +U term, on both f orbitals and p orbitals: B. The cRPA method to compute effective Coulomb interactions Constrained Random Phase Approximation (cRPA) intends to compute effective interactions, by carefully separating screening effect [16,31,49,50]: in order to avoid double counting of screening effects, effective interactions between correlated electrons should not contain screening effects arising from correlated electrons.
In this work, we have generalized the implementation of Ref. [33] in ABINIT to the calculation of effective interactions among different correlated orbitals. This implementation has been used in Ref. [26].
In this section, we present the method and notations we will be using in the rest of the paper.
1. From the non-interacting polarisability to effective interactions In a system of interacting electrons, with a Coulomb interaction v = 1 |r−r | , each interaction is screened by the rest of the system. Using linear response theory, it is possible to compute this screening in RPA using non interacting polarisability of the system χ 0 (ω) and v. χ 0 (ω) is obtained from first order perturbation as electron-hole excitations (see e.g. [51]). The main idea of cRPA is to constrained a dynamical polarisability calculation, by splitting this term as with χ correl 0 is the part due to excitations among correlated orbitals and χ r 0 the part due to the rest of the electronic excitations. Such splitting enables to define χ r 0 . The cRPA polarisability is thus with σ k,ν the energy for band index ν at k-point k,with spin σ and f σ k,ν the factor of occupation. w(k, k , ν, ν , σ) is a weight function, that permits to exclude transition among the correlated bands. Different definitions of this term are detailed in the next section II B 2 for entangled and non-entangled bands.
The polarisability is linked to the susceptibility, in matrix notation for the position variables, by ε r (ω) = 1 − vχ r 0 (ω). We can define the dynamically screened Coulomb interaction matrix using this Finally, we compute the U scalar using this interaction matrix We can also define J as [39,43,52] 2. Practical calculation of χ r 0 We detail here the calculation of χ r 0 in eq.10. To constrained electronic transition in χ r 0 to transitions other that the one among correlated bands, we use the term w(k, k , ν, ν , σ) in equation 10. First we have to define the subset C, which is the subset of bands considered correlated. In the simplest cases C = W, but it can be useful to take different definitions [33]. If correlated bands are completely isolated from the other ones, then one have [16]: when (νk) and (ν k ) are both correlated bands and w = 1 otherwise. This weighting scheme is called model (a). If correlated bands are hybridized with non correlated bands, then the previous scheme cannot be used: Due to the partial character of some bands, it is mandatory to consider some non-integer weight. To achieve such goal, we follow what has been proposed by Şaşıoglu et al. [50], and use a weight function proportional to correlated orbitals weight on each Kohn Sham functions (see also [33,53]). This weighting function takes the following form and is called model (b).

Notations and models for cRPA calculations
The cRPA scheme is based on those three parameters : • The bands used to construct Wannier functions W • The bands used to remove transition C • The model to remove transition, either (a) or (b) We use a simplified notation for cRPA calculation that reads as C − W (a/b). We note a model using the main character of the band constituting the subset. It may be useful to take a large C, using model (b), in this case, we note ext (for extended) as the name of bands (see Ref. 33). For example, in UO 2 , if we take W as containing exclusively f bands, choose C = W, and remove transition using model (a), the results will be noted f − f (a) (see also [33]). Table I gives the model used in cRPA calculations in UO 2 in this work. Bands 7 to 19 covers the 6 p-like bands and the 7 f -like bands.

C. Calculations parameters
All the calculation are done using the ABINIT code, in the PAW formalism. Spectral calculations are converged to a precision of 0.1 eV on the gap, structural parameter calculations are converged to 0.1Å 3 on the volume and cRPA calculations of effective interactions are converged to 0.1 eV.
For all the spectral and structural calculations, we have used the occupation matrix control method [46,[54][55][56] to find the antiferromagnetic ground state of UO 2 in DFT+U and we use the LDA approximation for DFT (a comparison to GGA and ferromagnetic case is given as a supplementary material). We have used a plane wave cut-off of 20 Ha for the wave function and 35 Ha for the compensation quantity of PAW, and a 64 k-points grid.
For the cRPA calculations, we used the ferromagnetic GGA due to the prohibitive cost of calculation for the anti-ferromagnetic state. In the supplementary, we show that FM LDA and FM GGA give very similar results. The atomic data are the same as for the other calculations. The plane wave cut-off are set to, 25 Ha for the wave function, 50 Ha for the PAW compensation quantity, 7 Ha for the polarization and 35 Ha for the effective interaction. We have used 100 bands and a 64 k-points grid.

TiO2
All the calculations are done in LDA. As the p band is full and the d band is empty, the ground state is found, without the need to use the occupation matrix control method. The PAW atomic data consider 3s, 3p, 4s, and 3d as valence electrons for titanium and 2s and 2p for oxygen.
The cut-off are set to 30 Ha for the wave function and 60 Ha for the PAW compensation quantity. An 8 × 8 × 8 k-points grid is used, except for the HSE06 calculation with a 6 × 6 × 6 grid.
For the cRPA calculations, we used cut-off of 90 Ha for the plane wave, 100 Ha for the PAW compensation quantity, 5 Ha for the polarization and 30 Ha for the effective interaction. We used a 64 k-points grid and 70 bands.

III. ROLE OF U ff AND Upp IN UO2
This section focus on UO 2 which is is a prototypical Mott Hubbard insulator. It orders antiferromagnetically under 30K and exhibits a 2 eV gap [57]. It is also important because it is fuel for nuclear reactor. Its first principles description requires to go beyond LDA or GGA approximations and include explicitly the electronic interactions using DFT+U [55,58] or DFT+DMFT. In this work, we use DFT+U to show the impact of pp interactions on spectral and structural properties.

A. Spectral properties
We first remind the role of U ff on density of states of UO 2 before studying the role of U pp .

Spectral properties with DFT+U ff
We emphasize that UO 2 is particularly interesting to investigate the effect of U pp because O p-like bands and lower U f -like Hubbard bands are separated in energy so that the relative impact of U ff and U pp can be disentangled. This separation is clearly seen in Fig. 1. This figure compares (a) the experimental spectral function [57], and our calculation of the (b) LDA spectra, (c) DFT+U spectra using U ff =5 eV and (d) HSE 06 functional spectra. The experimental photo-emission spectra shows that UO 2 is insulating. The conduction band is composed of a broad band centered at 5 eV, it can be interpreted as a upper f -like Hubbard band. The valence band exhibits a peak at -1 eV -which can be seen as a lower Hubbard band -and a broad band localized at -5 eV which can be interpreted as a p-like band. In the DFT-LDA band structure (graph (b)) we recover the f -like band at the Fermi level, thus the system is described as metallic in contradiction to the experimental photo-emission spectra (graph (a)) [59].
DFT+U with U = 4.5 eV (graph (c)) recovers an insulating density of state, by splitting the 5f band in two parts, creating a gap of 2.0 eV in good agreement with experiment [59]. Nevertheless, the position of the p band is at -3 eV, too high in comparison to experiment.
Lastly, using the hybrid functional HSE06 improves on the position of p like band, but do not correctly describes its width, in agreement with previous work [60].
Thus the comparison of DFT+U ff results with experiment or hybrid functional calculations highlights the wrong placement of the O-2p-like band in DFT+U ff .
In order to improve this description, we include electronic interaction U pp on the 2p orbitals. Such correction is expected to lower the position of the p-like band. It is discussed in the following subsection.

Spectral properties with DFT+U ff + Upp
The goal of this section is to investigate the role of U pp on spectral properties.
a. Results with atomic orbitals Fig. 2 (a) first reproduces as a reference the density of states obtained in DFT+U ff with U = 4.5 eV.
Then we use DFT+U ff +U pp with U pp = 5.0 eV. The result of this calculation is reproduced in graph (b). We first compare these fully converged DFT+U pp and DFT+U pp +U ff calculations (graph (a) and (b)). As the +U potential for p orbitals is −U pp (n p −0.5), and as p orbitals are filled (n p = 1), we could naively expect a shift of U /2 for the p-band. U pp =5 eV should thus shift the p band by 2.5 eV. 4 However, we observe a surprisingly small shift in graph (b) with respect to graph (a).
In order to investigate this behavior, we report on the graph (c) of Fig.2, the density of states without self-consistency. This is similar to apply only once the DFT+U ff +U pp potential to the DFT+U ff density and then diagonalize the Hamiltonian without recomputing the density. The only difference between the graph (b) and the graph (c) in Fig.2 is the density and thus the number of electrons that appears in the DFT+U potential: In graph (b) the number of p electrons used to compute the Hamiltonian is the fully converged DFT+U ff +U pp number of electrons, whereas in graph (c), it is the DFT+U ff number of electrons. In this case, we observe that the shift of the p-like band is much larger.
We investigate the origin of these differences in the next subsection.
b. Discussion Before studying the DFT+U potential, we first focus on the DFT+U energy. The DFT+U energy due to the 4 Because of hybridization of Kohn Sham states, the true shift for Kohn Sham states Ψ σ kν should be smaller and equal to m Ψ σ kν |χ l m Upp(0.5 − np) χ l m |Ψ σ kν with | χ l m |Ψ σ kν | 2 close but lower than 1. p electrons contains a Coulomb interaction term and a double counting term and is expressed as where U pp is the strength of the effective interaction between p electrons, n p is the number of electrons in a p orbital, and N p the total number of p electrons. This energy can be derived with respect to n p to obtain the potential associated to the p orbitals, We obtain for the U pp potential.
In a similar way, we can also express a total Hartree energy, restricted to p − p interactions with V pp being the strength of the unscreened Coulomb interaction between two p electrons (on the same atom) 5 .
We can rewrite the Hartree energy by considering that even at the Hartree level, the Coulomb interaction is screened by the other electrons of the system.   Thus the Coulomb interactions can be expressed with U pp the effective Coulomb interactions. We can also rewrite this expression by using that p1,p2 n p1 n p2 = p1 n p1 p2 n p2 = N 2 p . The effective Hartree energy, limited to p − p interactions is thus Again, this energy can be derived with respect to n p to obtain the potential associated to the p orbitals, We obtain for the Hartree potential. The equation for the Hartree potential (Eq.22) and +U potential (Eq.18) thus depends only on the number of occupations n p and interaction U pp . We can thus use this expression to understand graph (b) and (c) of Fig.2. For the effective interactions, we use re-normalized cRPA values for U pp (U pp = U cRPA pp − U cRPA pf = 5 eV). The numbers of electrons obtained from DFT+U ff and DFT+U ff +U pp are gathered in Table II. They are coherent with a decrease of hybridization when U pp is added to the Hamiltonian: N p increases and N f decreases. From these values, it is possible to compute the values of +U and Hartree potentials using Eqs. 18 and 22.
We can now use these expressions to understand Fig.2. The physical picture briefly discussed in the caption of Tab.II is the following: When U pp is added to the system without self-consistency on the charge (graph (c) of Fig. 2 and first column of Tab. II), the only change in the Hamiltonian is V U = U pp ( 1 2 − n p ) −1.4 eV so that the p is shifted by 1.4 eV to low energy. Self-consistency over charge (graph (b) of Fig. 2 and column two of Tab. II) induces an increase of p electron number [20]. Such increase mainly increases the Coulomb repulsion among p electrons that arises from Hartree energy. Such increase thus pushes the p-like band to high energy, nearly canceling the shift due to V U . It thus explains that the effect of U pp is rather small at self-consistency as observed previously [17].
c. Results with Wannier orbitals The Fig. 3 show results using Wannier orbitals as correlated orbitals and compare (a) DFT+U ff , (b) Fully converged DFT+U ff +U pp and (c) DFT+U ff +U pp with the DFT+U ff density and no self-consistency (as in (c) of Fig. 2).
Let's now discuss calculations using atomic and Wannier orbitals. Graph (a) shows that DFT+U ff density of states are similar except that the p-like band width is larger. This larger bandwidth can be attributed to the larger hybridization between the p-like band and the lower f Hubbard band using Wannier orbitals.
For DFT+U ff +U pp , the non self-consistent calculations (graph (c)) are similar also except the width of the p-like band. For self-consistent DFT+U ff +U pp (graph (b)), there is an important difference. Whereas with the atomic orbitals, the p bands are weakly shifted with respect to DFT+U ff , with the Wannier orbitals, the shift is much larger (1 eV). The difference can be explained using the number of electrons in the Wannier calculation (see Tab. III). Indeed the change of the number of elec-   trons as U pp increases is much smaller in the Wannier case than in the atomic case and as a consequence the p-like band stays low in energy as discussed in Tab.III. Concerning the final results, the density of states of the fully converged calculation is in rather good agreement with HSE06 calculation and experiment concerning the position of the bands.
B. Calculation of U ff , Upp and U f p in UO2 In this section, we use our recent cRPA implementation which allows calculations of effective interactions among multiple orbitals in systems with entangled bands to compute U ff , U pp and U f p in UO 2 .
The goal is first to understand the role of the initial band structure on calculated values and to determine values that could be used to compute structural and electronic properties. The results are gathered on Tab. III B 1. cRPA calculations with the f p−f p model lead to large values of U , whereas with the f p − ext model, values of U ff and U f p are weaker. It can be explained with the same argument that in Ref. [33] for the difference between f − f p (a) and f − f p (b) calculations: In the f p − ext, there is a residual weight at the Fermi level, which is neither of f nor p character. It contribute a lot to the screening of the interaction. We can also note that in contrast the value of U pp increases a bit comparing to U f p or U pp . This is because the bare value of interactions increased importantly in the f p − ext model for these orbitals that are more delocalized. Note that for GGA+U with U ff =2 eV, the difference between f p − f p model and f p − ext model is negligible because the system is insulating and thus there is no residual contribution at the Fermi level. Using LDA instead of GGA leads to similar results, as reported in the supplementary.

Role of U in
ff on computed U values.
As can be seen on Tab. III B 1, when U in ff increases from 0 eV to 6.0 eV, values of U ff , U f p , and U pp increases. This is simply due to the fact that the screening is lower as U ff and thus the gap increases.

Role of U in
pp on computed U values.
We now compare the values of U , with respect to the input values of U pp . The effect of U pp is to increase the partial gap between f and p bands as discussed in section III A 2. As a consequence, some transitions are occuring at a larger transition and thus the screening is reduced. Comparing values computed for U pp ranging from 0 to 10 eV, we indeed obtain increasing values of all interactions.

Hybrid functionals
Using the HSE06 hybrid functional, the p band is even lower in energy. As a consequence, the screening is even lower, and thus values of interactions are found larger.

C. Effect of U on structural parameter
In this section, we investigate the role of U pp on structural properties. In DFT+U ff , results are nearly independent of correlated orbitals whereas in DFT+U ff +Upp, atomic and Wannier orbitals lead to opposite variation of volume with respect to DFT+U ff . As discussed in the text, this difference comes from the variation of the number of p electrons as a function of volume, which is more physical using Wannier orbitals.
The figure 4 report the calculated values for the equilibrium volume in different cases.
a. LDA. The LDA calculation underestimates the equilibrium volume for UO 2 , with a value of 147.4Å 3 where the experimental one is 163.7Å 3 .
b. LDA+U ff . As shown previously in a large number of studies (e.g [58]), adding a U ff in the calculation increase the value of the equilibrium volume, and we get 159.9Å 3 . Using Wannier orbitals, the volume is very close (159.7Å 3 ).
c. LDA+U ff +U pp . The application of U pp leads to drastically different results depending on the choice of correlated orbitals. If we choose the atomic basis, we find a decrease of the volume when U pp is applied: the volume decreases from 159.9Å 3 to 154.7Å 3 . This is in agreement with several studies on oxides [17,22,24,32]. In Ref. [32], for example, DFT+U dd +U pp calculations (called ACBN0 in this context) for oxides gives smaller volume than DFT calculations.
In contrast, when using the Wannier basis, we observe a much larger equilibrium volume in comparison to the LDA+U ff case, with a value of 164.1Å 3 .
We interpret this difference in the next subsection.

Influence of the DFT+U correlated orbitals on equilibrium volume
Understanding this difference of behavior according to correlated orbitals requires to go back to the expression of total energy. The DFT+U total energy can be written as: In order to understand the origin of the change of volume in various DFT+U calculations, we have plotted E DFT+U , E DFT and E U on Fig. 5 for four different calculations: DFT+U ff using atomic orbitals or Wannier orbitals and DFT+U ff +U pp using atomic orbitals or Wannier orbitals. In those graphs, the different energies have been shifted in order to make the comparison of the variation easier. The blue curve is the total energy of the system E DFT+U , as a function of the volume. The minimum of this function gives the equilibrium energy and the equilibrium volume. The green curve represent the DFT energy, calculated with the DFT+U density E DFT [n DFT+U ]. The orange graph represents E U that depends explicitly on U ff and U pp terms.
In DFT+U ff , the energy variation are similar (see graphs (a) and (b) of Fig. 5) for the calculations using either atomic or Wannier orbitals. This is expected as the equilibrium volumes are similar as shown previously on Fig. 4. In DFT+U ff +U pp , the energy variations are different (graphs (c) and (d)), coherently with the fact that the equilibrium volume are different. The analysis of these graphs shows that a large contribution to the change of volume comes from E U . Indeed, by comparing the evolution of this term, we observe that for atomic orbitals, it decreases with the volume, whereas for Wannier orbitals, it increases with the volume. The E U contains two contributions coming respectively from U pp and U ff . The U ff term in E U has a similar variation as seen in (a) and (b). Thus the difference in the variation of E U in (c) and (d) comes mainly from the U pp term in E U . As U pp is fixed, this difference is linked to the evolution of occupations as a function of volume which are plotted in Fig. 6. As expected, the p occupation of each orbital have a different variation when volume increases in the two cases: it decreases for atomic orbitals and increases for Wannier orbitals. The variation of the number of p electrons in the calculation using Wannier orbitals is expected on a physical basis: When volume increases, oxygen p orbitals and uranium f orbitals are less hybridized so that p like bands contains a lower contribution of f orbitals, and thus occupations of p orbitals increases. So we can expect the Wannier calculation to give more physical results. The variation of n p in the calculation using atomic orbitals appears nonphysical. It might be related to the fact that atomic orbitals are not adequate to describe accurately these ionic systems as a function of volume 6 .
Let's understand now how these variations of n p impact the energy versus volume curve.

Mechanism of structure modification
The link between structural properties and occupations is simply the relation between E U and occupations which is (with J = 0) [12,59]: where n i,σ is the number of electrons with spin σ on the orbital i. This graph is plotted as function of n i in Fig.7 and Fig.8. Before explaining the role of U pp and the impact of the variation of the number of p electrons on equilibrium volume, we first focus on U ff whose impact on equilibrium volume is well established.  Fig. 6). As discussed in the text, the Wannier description is coherent with the decrease of hybridization as the volume increases: it leads to a decrease of energy as a function of volume and thus an increase of the equilibrium volume.
a. U ff In UO 2 , there is 14 f orbitals and 2 f electrons. In DFT, the 2 electrons are shared between all the orbitals taking into account the crystal field (and/or spin orbit coupling). In DFT+U , applying U ff to the f orbital split orbitals in two parts in order to lower the DFT+U energy [55]: 2 orbitals getting nearly full and below the Fermi level (See Fig. 1) and the 12 other getting nearly empty [55].
The Fig. 7 illustrates the change of energy due to the change of occupations when the volume increases in DFT+U ff . First, when the volume increases, hybridization lowers. Thus, the 12 f orbitals corresponding to empty bands will get more empty (n i closer to 0): their occupation is the black line on Fig.7, and this occupation change according to arrows when the volume increases. On the other side, the two f orbitals corresponding to filled bands will get more full. By looking on Fig. 7, we can understand that it will lead to a lower DFT+U energy as volume increases, so that the equilibrium volume will be greater in agreement with Fig. 4. In other words, if numbers of electron are closer to 0 or 1, then the self-interaction correction decreases.
b. U pp We now discuss the evolution of energy for the E Upp term. p orbitals are different from f orbitals: p-like bands are all below Fermi level and so all fully occupied. Because of the hybridization p orbitals are however not completely filled.
The Fig. 8 illustrates the DFT+U energy variation due to occupation of p orbitals. As discussed above the number of p electrons variation upon increase of volume is different for the atomic and Wannier calculations. Straightforwardly, the increase/decrease of the number of electron in DFT+U using Wannier/atomic orbitals induces a decrease/increase of the energy. As a consequence, using the atomic orbitals, the equilibrium volume lowers when U pp increases. Using Wannier orbitals, we have the opposite (more physical and expected) behavior: when interaction increases, volume increases. Such results thus raises doubts on the physical results obtained with U pp applied on atomic orbitals.
In conclusion of this section on UO 2 , we have shown that the effect of U pp can be understood and its role depends crucially on the choice of correlated orbitals.

IV. ROLE OF U dd AND Upp IN TIO2
Here, we study the role and calculation of U pp on the rutile phase of TiO 2 . Contrary to UO 2 , which is a Mott insulator, TiO 2 is a prototypical charge transfer insulator with empty d orbitals. It has been the subject of several studies using DFT+U dd or DFT+U dd +U pp (eg. [17,18,22,29,56,61,62]).

A. Density of states of TiO2
In this section, we compare and analyse successively the respective role of U dd and U pp on the spectral properties of TiO 2 .

Density of states calculated in LDA+U dd
The figure 9 compares the densities of states calculated in LDA, LDA+U dd and HSE 06 functional. As for several charge transfer insulators with filled shells, DFT/LDA describes the system as an insulator, but underestimates the band gap (1.9 eV instead of 3.0 eV experimentally [63]) (see e.g. [64]). The effect of U dd on the density of states is straightforward. The empty d-like band is shifted upward and the band gap increases (see Fig. 9). For U dd =4.5 eV (graph (b) on Fig.9), LDA+U dd gives a gap of 2.7 eV. The last calculation (graph (c)) uses the self-consistent values of U dd = 8.1 eV calculated in the dp − dp model (see tab.V) and gives a band gap of 3.5 eV. We recover thus the main results of previous work [17,18]. Concerning HSE06, we find a gap of 3.2 eV in agreement with previous work [64].

Density of states calculated in LDA+U dd +Upp
We now include the interaction on both titanium d orbitals and oxygen p orbitals . The figure 10 compare the densities of states in LDA+U dd and LDA+U dd +U pp . As for UO 2 , in order to disentangled physical effects, we show self-consistent calculations (b) and non selfconsistent calculations (c). We show that adding U pp increase the band gap to 5.3 eV in the first self-consistent field iteration (c). At convergence of the self-consistency, we find a band gap value of 4.1 eV.
As it has already been discussed [15,62], using DFT+U on e.g. TiO 2 affects the p band width. In order to decorrelate the effect of band width from the band shift, we represent on the density of states the barycenter of the p band, which should not be affected by a band width change. The study of the barycenter gives us a shift of -1.27 eV when U pp is included in a non self-consistent way from the DFT+U dd calculations (comparison of (c) and (a)). Self-consistency gives a shift of 1.00 eV (from the comparison of (b) and (c)). We carry out the same kind of study as for UO 2 , by calculating numerically the band shift using the change in the occupation numbers. The results are reported in the table VI.
As for UO 2 , the second and third columns of the table correspond to the numbers of electrons and corresponding potentials used to produce non-self consistent density of states (graph (c)) and self-consistent one (graph (b)) of Fig. 10. In the non self-consistent DOS (graph (c)), the p band is approximately shifted by V Upp with respect to the case where U pp =0 eV. Here V Upp = -1.59 eV is of the same order of magnitude as the shift of the barycenter of the p band (1.27 eV). 7 At convergence we see that even if V Upp is decreased (-0.21 eV), the increasing of the Hartree potential leads to a total increase of the potential of the p band (+1.06 eV) with respect to the non self-consistent case. Fig. 10 shows indeed that the barycenter of p band is shifted by 1.00 eV (from graph (c) to graph (b)), not far from the calculated value. In summary, U pp has a small effect on the density of states. We underline that the effect would be even smaller if we would not re-normalize the atomic wave-functions used to define correlated orbitals.

Density of states calculated in LDA+U dd +Upp using Wannier functions
Again as for UO 2 , we present here DFT+U dd +U pp calculations to study the dependence with correlated orbitals. The Fig. 11 represents TiO 2 total density of states similar to Fig. 10 but computed using Wannier orbitals. In Fig. 11, we see, by comparing graph (a), (b) to (c), that in comparison to atomic-like orbitals calculations (10), the calculation using Wannier functions gives slightly larger shift of the p-like band : The shift is -2.41 eV for the p-band between graph (a) to (c) of using Wannier functions occupations. It arises from this analysis that number of p electrons are larger, which explains the large shift, but its variation is smaller, which explains the weak self-consistency effect. In practice, the observed self-consistency effect is not so weak (about 1.5 eV) and maybe the large variation of N d could explain it. However and globally, comparing DFT+U dd (a) and DFT+U dd +U pp (b) we see that the global effect of U pp is slightly larger using Wannier orbitals.

B. cRPA calculations for TiO2
In this section, we use our cRPA implementation [26,33] to compute the U interactions parameters in TiO 2 (see Tab.V).
In these calculations, it appears that Ti-d-like and Op-like bands are separated from the rest of the bands in the system, thus model (a) and (b) as discussed in section dp − dp 0 0 11.7 9.8 3.9 LDA dp − dp 8.1 5.3 11.6 8.8 3.7  [65]. We carry out a selfconsistent calculations in the sense of Ref. [26,33,66] and we find a value of 3 eV, very close to the non selfconsistent value.
In the dp − dp model, values of U are larger because screening is weaker (see e.g. [33]). Self-consistency effect is weak the maximum difference being 1.0 eV on U pp . As for UO 2 , the final values we used for U dd and U pp in the dp−dp model are re-normalized by U dp (see [25]). Values we retain are U dd = 8.1 eV and U pp = 5.3 eV.

C. Structural parameters
We compare on Fig 12 equilibrium  First the LDA-DFT underestimates largely the volume with a value of 60.48Å 3 . Using U dd , the equilibrium volume increases. Using atomic like correlated orbitals leads to a larger increase than using Wannier correlated orbitals.
We then carried out the study using DFT+U dd +U pp , using U dd =8.1 eV and U pp =5.3 eV. The application of U pp leads to the same effect as for UO 2 : With atomic orbitals, volume decreases (down to a value of 60.70Å 3 ) whereas with Wannier orbitals, volume increase (up to a value of 63.74Å 3 ), as expected from a physical point of view. Besides, the Wannier calculation is somewhat closer to experiment.
In appendix, we decompose the energy as we have done in UO 2 . The conclusions are similar. Also the number of p electrons in Wannier and atomic orbitals show similar trends as in UO 2 . So the relative role of Wannier orbitals and atomic orbitals concerning structural properties is similar in TiO 2 and UO 2 .

V. CONCLUSION
We conducted here a detailed study on the role of electronic interactions in p orbitals of oxides using UO 2 and TiO 2 as prototypical Hubbard and charge transfer insulators. We used our new cRPA implementation (allowing   multi orbitals interaction calculations) to obtain effective interactions among f /d and p orbitals. Using obtained U values and DFT+U , we investigate the effect of U pp on spectral and structural properties and discuss its physical origin in terms of electrons numbers. For structural properties, we find a reduction of the volume when U pp is added, which is a counter-intuitive result, in agreement with previous studies [17,23,24,32], when correlated orbitals are atomic orbitals. We show that using Wannier orbitals as correlated orbitals restore expected results, mainly because the number of electrons and its evolution as a function of volume are more physical. Such results put light on the physics brought by the U pp interaction in oxides. It is especially important because several scheme use U pp as a way to mimic hybrid functionals in the goal of performing high throughput calculations [15,29,32]. Our implementation is in this line, as our combination of DFT+U pp and cRPA using a common Wannier basis gives results in rather good agreement with experiment or HSE06 hybrid functionals. More calculations are needed to confirm this view.
The results are very similar to UO 2 . First using only U dd = 8.3 eV, we found that E U decreases while the volume increases, using either atomic or Wannier orbitals.
Using U dd = 8.1 eV and U pp 5.3 eV, the effects are also similar to what we observed in UO 2 . For an increase of volume, E U slightly increases using atomic orbitals, and largely decreases using Wannier orbitals. It can also be explained using the occupations of the p orbitals. Results are shown in Fig. 14. This figure clearly show that for a volume increase, atomic like orbitals leads to reduced occupation of the p orbitals, contrary to Wannier orbitals. As discussed in section III C 3, this mechanism leads to a reduction of the equilibrium volume for atomic orbitals.
Where χ l,m is an atomic-like orbital. In this appendix, we sum-up the different ways to define correlated orbitals |φ t 0 in PAW as discussed in Ref. [36].
• In first one (noted o B in Ref. [36]) we re-normalize the truncated radial part and so |φ t 0 are replaced by As well explained in Appendix C of [36], applying this formula straightforward leads to nonphysical occupations, with a maximum occupation of around 4.5 electrons for p orbitals (depending on PAW radius).
• In the second one (noted o C ), we use the renormalized wave-function as in o B , and then renormalize again the scalar product of the KS wavefunction and the correlated orbitals define in o B . This leads to a scalar product which is: As discussed in [36], such formulation recovers the atomic limit.
• In the third one (noted o W ), that we introduce here, we drop the local projection with atomic orbitals and use Projected Localized Wannier functions, detailed in part II B.
In the supplementary material, we compare calculations using o B and o C . The solution used in [36] is to re-normalized the occupation in order to have a maximum occupation of 6 electrons for a p orbital. In the case of d or f orbital, the re-normalization of occupation has a small impact as they are well localized. To highlight this point, we compare in the supplementary, 'the density of states for DFT+U f f and DFT+U f f + U pp calculated using o B and o C .