Charge self-consistency in density functional theory + dynamical mean field theory: k-space reoccupation and orbital order

We study effects of charge self-consistency within the combination of density functional theory (DFT; Wien2k) with dynamical mean field theory (DMFT; w2dynamics) in a basis of maximally localized Wannier orbitals. Using the example of two cuprates, we demonstrate that even if there is only a single Wannier orbital with fixed filling, a noteworthy charge redistribution can occur. This effect stems from a reoccupation of the Wannier orbital in k-space when going from the single, metallic DFT band to the split, insulating Hubbard bands of DMFT. We analyze another charge self-consistency effect beyond moving charge from one site to another: the correlation-enhanced orbital polarization in a freestanding layer of SrVO3.


I. INTRODUCTION
Density functional theory (DFT) 1,2 is highly successful in predicting various material properties such as crystal structures, ionization energies, electrical, magnetic and vibrational properties. Indeed DFT is the de facto standard for calculating materials' physical properties. But even the best approximations for the DFT exchangecorrelation functional fail to describe one class of materials, known as strongly correlated systems. In these materials, the interaction between electrons is insufficiently screened to be amenable to the available functionals. One might add a static Coulomb correction within the socalled DFT+U formalism. 3 This often yields an improved description, in particular of strongly correlated insulators, but it has its own limitations: DFT+U is essentially a Hartree-Fock-like treatment with a single-Slaterdeterminant ground state. In this situation, the energy cost of the Coulomb interaction can only be avoided by symmetry breaking, which is hence largely overestimated.
Dynamic, albeit local, correlations can be taken into account by dynamical mean field theory (DMFT) [12][13][14] which has been merged with DFT for realistic calculations of correlated materials. [15][16][17][18] Here electrons can stay on or leave lattice sites dynamically so as to greatly suppress double occupation and the cost of the Coulomb interaction, even in a paramagnetic phase without any symmetry breaking. If one has a three dimensional material at elevated temperatures, say room temperature, and if there is no magnetic or other phase transition close-by, 19 these local DMFT correlations prevail. Already the first applications showed that DFT+DMFT well describes transition metals 16 , their oxides 20 , and felectron systems 21,22 .
These early calculations were so-called "one-shot". That is following a DFT calculation, the relevant correlated orbitals and the corresponding single-particle Hamiltonian were identified. This DFT Hamiltonian was supplemented by local Coulomb interactions for the dor f -orbitals and solved with DMFT. Physical proper-ties such as the spectral function, susceptibility or magnetization were calculated from this "one-shot" DMFT solution.
Since the DMFT correlations change the site and orbital occupation and consequently the charge density, a natural next step is to do a "charge self-consistent" (CSC) DFT+DMFT calculation. [23][24][25][26][27][28][29] That is, from the DMFT Green's function, a new charge distribution is calculated which in turn serves as input for the DFT potential. This leads to a new DFT Kohn-Sham Hamiltonian, subsequently a new DMFT Green's function etc. This cycle is repeated until convergence. While it has been pointed out in the literature [23][24][25][26][27][28][29] how the DMFT spectral function, the double counting and the d (f ) energy level changes due to CSC for specific materials where charge is moved from one site to another, little attention has been paid to the redistributed charge itself, its spatial arrangement, and more.
The aforementioned change of double counting and d level shift can be understood as follows: In a typical situation, say for a transition metal oxide, the dominant d states crossing the Fermi energy have some oxygen p admixture; conversely, the oxygen states below the Fermi level have some d contribution. Including electronic correlations in a so called d+p DMFT calculation will reduce the d occupation somewhat, and increase the p occupation on the oxygen sites. In the next DFT step, the larger p occupation will increase the p (Hartree) energy and decrease the d (Hartree) energy. This counteracts the first shot DMFT to have less d and more p electrons, dampening the charge redistribution of the "first-shot" DFT+DMFT.
In this paper we study effects of CSC beyond this gross effect of a p-d orbital and site reoccupation. In Section II, we recapitulate the CSC DFT+DMFT approach and outline our implementation thereof. In Sections III A and III B, we show that even in a single-orbital, d-only DMFT calculation, there is a charge redistribution akin the d-p reoccupation effect mentioned above. This runs counter to the naive expectation that there can be no charge redistribution in this situation since the number of elec-trons in the single, predominately d-like orbital centered around the transition metal site is fixed. The two materials studied, where a restriction to a single d band is justified, are Sr 2 CuTeO 6 and HgBa 2 CuO 4 in Section III A and III B, respectively. In Section III C, we study the effect of correlation-induced orbital order on the charge redistribution and self-consistent DFT+DMFT results. Specifically, we consider an ultra-thin layer of the cubic perovskite material SrVO 3 , where breaking of the cubic symmetry stabilizes the in-plane xy orbital against the xz and yz orbitals. This orbital ordering is strongly enhanced in DMFT because of electronic correlations. Finally, Section IV summarizes our main findings.

II. METHODOLOGY
We now present the formalism and our implementation of CSC DFT+DMFT which is in a basis of maximally localized Wannier functions (MLWF). For these, the measure of localization introduced by Marzari and Vanderbilt 7 is the spread in real space. This provides for a very flexible approach portable to any bandstructure method. Moreover, the methodology allows bondcentered or molecular Wannier functions. Our starting point is the wien2wannier 10 interface between Wien2k 9 and Wannier90 6 , and the w2dynamics 30 continuous-time quantum Monte Carlo 31 DMFT implementation.
We combine and extend these methods to include CSC. Let us, for the sake of completeness and given the sparse presentation in the literature, recapitulate here the CSC DFT+DMFT approach, and discuss the peculiarities of our implementation. Readers only interested in the physical applications and effects of CSC can safely skip the rest of this Section.
The CSC DFT+DMFT method relies on the simultaneous convergence of two local observables: the electronic density as the central quantity of DFT and the local Green's function as the central quantity of DMFT. Both mutually affect each other in the CSC cycle. The charge density at position r is given by while the local DMFT Green's function is in the basis of localized Wannier orbitals χ m . Here, m, m enumerate orbitals on a site, β is the inverse temperature, and the factor e iωn0 + ensures the convergence of the sum over Matsubara frequencies ω n = (2n + 1)π/β. In both expressions there appears the full Green's function of the solid, which can be written as with − ∇ 2 2 , V KS , and µ being the kinetic energy operator, Kohn-Sham (KS) effective potential, and chemical potential, respectively. The effective local self-energy ∆Σ =Σ −Σ dc is determined from the DMFT self-energŷ Σ by subtracting a double counting correction termΣ dc , which, as far as possible, accounts for electronic correlations already included in DFT. The KS potential depends on r and consists of the external potential V ext due to the nuclei, a Hartree potential V H describing part of the electron-electron Coulomb repulsion and an exchangecorrelation potential V xc . The latter is obtained here within the generalized gradient approximation (GGA) 11 , but other functionals are possible as well, e.g. hybrid functionals to improve on the exchange contribution.
In DFT, the effective potential is obtained from a charge self-consistent procedure, shown in the upper left part of Fig 1. This DFT cycle starts with an initial choice for the electron density, from which the effective potential V KS is constructed. Incorporating V KS , the Kohn Sham equation is solved to obtain a new density and so forth until convergence. The DFT cycle closes with a converged charge density and provides a reasonable electronic structure as a starting point for DMFT calculations.
There is however an important step between DFT and DMFT, identifying a localized basis (upper right part in Fig 1) since DMFT treats only local correlations. To this end, we employ Wannier functions that are constructed by Fourier transform of the DFT Bloch waves |ψ νk : Here,Û (k) is a unitary matrix, Ω denotes the volume of the unit cell, ν and α are the band indices of the Bloch waves and Wannier functions, respectively. We assume here that we can restrict ourselves to a band window with only C Bloch waves. In the scheme of maximally localized Wannier functions 7 ,Û (k) is obtained by minimizing the spread of the Wannier functions.
Eq. (4) works for isolated bands. However, in most cases, the target bands are "entangled" with further bands at least at some k-points. These additional bands might be less important for the physics but need to be projected out by a so-called "disentanglement" procedure. At each k-point, there is a set of C o (k) Bloch functions which is larger than or equal to the number of target bands, i.e., C o (k) ≥ C. The disentanglement transformation takes the form Here, the band index ν belongs to the "outer window" with C o (k) Bloch wave functions, while ν , α label the C target bands. Hence,V (k) is a rectangular C o (k) × C matrix. A Fourier transformation of |w αR leads to the FIG. 1: Schematic representation of the DFT+DMFT approach. In a non-CSC or "one shot" DFT+DMFT calculation, the DFT Hamiltonian is not updated and both the DFT and DMFT cycle close individually, i.e., we have the orange and green arrows in the schematic, but not the blue ones. In a CSC DFT+DMFT, neither DFT nor DMFT is iterated individually. Instead, both of them are closed together, i.e., we have the green and blue arrows in the schematic, but not the orange ones.
Wannier orbitals in k-space whose occupation will be at the focus of the physics discussed below: The Hamiltonian in Wannier space W is defined in terms of the |w αk and obtained by an unitary transformation for isolated bands and with an additional projection ("downfold") in case of entangled bands, i.e., respectively. In DMFT, this Hamiltonian is now supplemented with the local Coulomb interactions, and the lattice problem defined this way is mapped onto an auxiliary impurity problem which is solved self-consistently. 13,14 The noninteracting Green's functionĜ(iω n ) of the impurity problem can be considered as a dynamical mean field. The DMFT algorithm (see the lower right part of Fig. 1) consists of: (i) Applying the lattice Dyson equation for the local interacting Green's functionĜ(iω n ) In order to enhance convergence, one normally starts withΣ =Σ dc , i.e., using the Hartree-energy as a first guess for the self-energy. A total number of k-points, n k , is considered in the reducible Brillouin Zone.
(ii) Applying the impurity Dyson equation, which relates the non-interacting impurity Green's function to the (lattice and impurity) self-energy and interacting Green's functionĜ (iii) Solving the Anderson impurity problem (AIM) defined by the non-interacting Green's function and the local Coulomb interaction U , i.e., calculating its interacting Green's functionĜ imp (iω n ) This is numerically the most involved step; we employ the continuous-time quantum Monte-Carlo method 31 in the w2dynamics implementation. 30 (iv) Applying the impurity Dyson equation once again, this time to calculate the self-energy as the difference between the inverse non-interacting impurity Green's func-tionĜ(iω n ) and the interacting (lattice and impurity) Green's functionĜ(iω n ) This self-energy is now used again in step (i) to calculate a new local Green's function. This procedure is referred to as "DMFT cycle" in Fig. 1. In a oneshot DFT+DMFT calculation we would stop after this DMFT calculation, extracting the interacting Green's function and further physical quantities from the converged DMFT solution. By contrast, in CSC DFT+DMFT calculations, we need to determine the DMFT-modified electron density (lower left segment in Fig. 1); recalculate from this the Kohn-Sham potential and the Bloch waves without DFT self-consistency; redo for these a Wannier function projection; which is the starting point for another DMFT step (see green and blue arrows in Fig. 1).
We still need to discuss how we calculate the DMFTmodified electron density ρ(r) = ρ DF T (r) + ∆ρ(r), which we defined in terms of the Kohn-Sham or DFT ρ DF T (r) and the correlation-induced difference ∆ρ(r). The latter can be calculated as: 26 where µ DF T and µ are the DFT and DMFT chemical potentials, respectively andĜ DF T (iω n ) = k [iω n + µ DF T −Ĥ W KS (k)] −1 is the DFT Green's function. It is computationally convenient to express ∆ρ(r) in momentum space, which can be deduced from Eq. (13) as It is to be noted that no convergence factor in the frequency summation needs to be used for ∆N W (k) because both Green's functions asymptotically decay as 1/ω n . Note that the change of occupation in Wannier space ∆N W α,α has an explicit k-dependence, which will have significant consequences in the following section.
In order to update the DFT charge density we need to transform ∆N W (k) from the Wannier to the Bloch basis using the unitary and disentanglement matrices,Û (k) andV (k), that define this transformation: Knowing the correlation-induced change of occupation in the Bloch or Kohn-Sham basis we can finally calculate the modified density since we know the spatial density D k ν ν (r) = ψ kν (r)ψ * kν (r) of each Bloch wave: The full CSC DFT+DMFT hence consists of the following workflow, schematically depicted in Fig. 1: • A converged charge density is obtained within DFT to have a reasonable electronic structure to start with (upper left part of Fig. 1). The target bands are identified as a prelude for the Wannier projection. In the following CSC DMFT cycle (green and blue arrows in Fig. 1), a single DFT iteration is performed to update the DFT Kohn-Sham Hamiltonian (i.e., without the orange arrow in the upper left part). We employ the Wien2k program package here.
• Maximally localized Wannier functions are computed within the target subspace as explained in Eqs. (4)-(6) (upper right section of Fig. 1). The DFT Kohn-Sham Hamiltonian is transformed into the Wannier basis following Eq. (7). We employ wien2wannier 10 and Wannier90 6 to this end.
• A single DMFT cycle is performed using w2dynamics 30 (lower right part of Fig. 1). This provides the self-energyΣ, local Green's function G, and the DMFT chemical potential µ, which is fixed to the particle number. Let us note that, for practical purposes, it is beneficial to start with a converged "one-shot" DFT + DMFT calculation. Moreover, a mixing (under-relaxation) between old and new DMFT self-energy is employed.
• For the correlated charge distribution (lower left part of Fig. 1), first ∆N W (k) is calculated taking the difference between DMFT and DFT Green's functions,Ĝ andĜ DF T , as in Eq. (14). As described in Eqs. (16)- (18), ∆N W (k) is transformed back to the DFT eigenbasis and used to obtain the correlation-induced change of density ∆ρ(r) and the total density ρ(r) of the correlated solution.
• The DFT+DMFT charge density, ρ(r), is finally compared with the old density. If the difference does not satisfy the convergence criteria, the new density is mixed with the old density and the result serves as the new density for a new V DF T and a new solution of the Kohn-Sham equation etc. until convergence. At the same time, a convergence of G(τ ) is also checked.

III. APPLICATIONS
In the following, fully CSC DFT+DMFT calculations are employed to shed light on correlation-induced charge redistribution beyond the gross effect of moving electrons from a WF centered at one atom to a WF centered at another atom. Two cuprates, Sr 2 CuTeO 6 and HgBa 2 CuO 4 , whose physics is dominated by a single band, are studied. The systems are different in several aspects. First Sr 2 CuTeO 6 exhibits a single isolated band around Fermienergy, while in HgBa 2 CuO 4 the single d band is entangled with other bands crossing it. On the technical side this requires disentanglement to project onto a single Wannier d orbital for HgBa 2 CuO 4 as discussed in the previous section.
Next, a multi-orbital situation is considered with a single, free-standing layer of SrVO 3 and t 2g orbitals at the Fermi energy that are well isolated from the other orbitals. Here, the interplay between structural confinement, orbital ordering, electronic correlations and CSC is discussed in detail.

A. Sr2CuTeO6
To describe the physics of cuprates, an effective single band model can be derived where the contributing orbital is predominantly of Cu-d x 2 −y 2 character, with some admixture of O-p x/y . The compound, Sr 2 CuTeO 6 , exhibits square lattice Heisenberg antiferromagnetism 33 in a quasi two-dimensional plane, consisting of Cu and O atoms, see Fig. 2 (left). It is quite unique in the cuprate group in having a completely isolated and weakly dispersing band around the Fermi-energy, see the white band in Fig. 2 (right). Thus, no disentanglement is needed in this material.
We take for our calculations the I4/m symmetry of the lattice with the experimental lattice parameters 32 , i.e., in-plane lattice constant a= 5.4308Å, out of plane lattice constant c = 8.4664Å. A slight complication of the lattice structure is that the CuO 6 octahedra in Sr 2 CuTeO 6 are rotated around the z-direction. In contrast to the CuO 2 planes of other cuprates, cf. Section III B below, Sr 2 CuTeO 6 has planes with four O per Cu; no oxygen is shared, which explains the low itinerancy.
In DFT, Sr 2 CuTeO 6 is a metal with a single half-filled band crossing the Fermi energy, predominantly of Cud x 2 −y 2 character. Electronic correlations result, however, in an insulating phase with two Hubbard bands separated by U . This is captured in our DMFT calculations, performed with U = 6.5 eV at inverse temperature β = 40 eV. During the charge self-consistency cycle, the self-energy and density are under-relaxed; 1000 k-points are considered in the reducible Brillouin zone for all the calculations.
The half-filled DFT band remains half-filled in DMFT. Naively one might expect that for an unchanged delectron occupation (half filling) there can be no CSC effect. However, the occupation in k-space is altered. In DFT (white band in Fig. 2) some k-points (in-between X-N -Γ) are below the Fermi level, and hence filled with one electron, whereas for all other k-points the occupation is zero.
In DMFT this half-filled band is split into two Hubbard bands that are broadened because of the imaginary part of the self-energy, the lifetime. This splitting means that now every k-point is occupied with half an electron (lower Hubbard band), whereas the remaining half electronic state (upper Hubbard band) remains unoccupied. That is, we have a major change of the occupation ∆N (k) in k-space, as calculated from the differences between G and G DF T at each k-point in Eq. (14). For the orbital occupation the sum over the entire Brillouin zone is taken, preserving the number of electrons in the d x 2 −y 2 orbital.
For the change of charge ∆ρ(r) in real space, however, each ∆N (k) in Eq. (18) is weighted with the spatial distribution of the corresponding Wannier functions. Hence, the splitting into Hubbard bands results in a charge redistribution: the Wannier functions have a different spatial dependence at each k-point. This correlation-induced correction to the charge density within the Cu-O plane is shown in Fig. 3. Here, the yellow (cyan) color corresponds to a gain (loss) of electron density in real space. As the single band of our consideration is predominantly of d x 2 −y 2 character, the contribution of each sign has the same orbital symmetry; the total change in density within the unit cell (shown as the black dashed box) is zero. The charge redistribution around each Cu ion can be understood easily for cubic crystal symmetry; see for HgBa 2 CuO 4 in section III B. Here, with lower symmetry, the charge redistribution around each Cu ion shows eight lobes with positive and negative contributions. Each Cu is surrounded by 4 oxygen atoms at the edges of the dotted box. As one can clearly see the positive contribution at these O sites is larger than the negative one. That means that even in our d-only model calculation there is some charge redistribution from Cu d to oxygen p. This is akin to the situation in d-p models where charge is moved from d to p orbitals as well. However in our calculation this effect occurrs even though we have only a single orbital in the DMFT calculation. This Wannier orbital is centered around the Cu sites and predominantly of d x 2 −y 2 character. But it has some admixture of oxygen p, i.e., it has some charge density at the neighbouring oxygen sites as well. This admixture requires some k-dependence of the Wannier functions [Eq. (6)]; and the occupation is reduced, e. g., around the N point, while increased in the remainder, eventually leading to the charge distribution pattern of Fig. 3.

B. HgBa2CuO4
Let us now turn to HgBa 2 CuO 4 , a prototype of the high temperature cuprate superconductors. 34 The ar-rangement of the CuO 6 octahedra is distinctly different in undoped HgBa 2 CuO 4 compared to that of Sr 2 CuTeO 6 , see the crystal structure in Fig. 4(a). The system belongs to the space group P4/mmm, with Hg, BaO, CuO 2 , and BaO layers stacked vertically along the c-axis of a tetragonal unit cell. Each oxygen atom in the CuO 2 plane is shared by two Cu atoms, resulting in a more direct hopping and a larger bandwidth of the Cu-d x 2 −y 2 band compared to that of Sr 2 CuTeO 6 . But the d x 2 −y 2 bands of HgBa 2 CuO 4 are no longer isolated. This requires disentanglement for constructing an effective single band model. Like Sr 2 CuTeO 6 , HgBa 2 CuO 4 is metallic in DFT, but is insulating if electron correlations are included as well as in experiment. For all the calculations we use 845 kpoints in the full Brillouin zone; and for the DMFT at β = 40 we employ U = 6.5 eV. This splits the DFT band into two Hubbard bands (see Fig. 5) and redistributes the k-space occupation of the Wannier orbitals as in the case of Sr 2 CuTeO 6 [ Fig. 2 (a)]. The most remarkable difference to that material is the much larger bandwidth in both DFT (white line) and DMFT (color).
Hence, for the very same reason as in the previous Section, ∆N W has a strong k-dependence. Since we have an effectively two-dimensional model, we plot in Fig. 4(b) ∆N W (k) in the plane k z = 0 of the Brillouin zone. The yellow section of the plane represents the set of k-points that have positive ∆N W . These states were unoccupied in DFT but get half-occupied in DMFT due to the lower Hubbard band dispersing throughout the Brillouin zone. The negative counterpart (blue) is around the Γ point where all states where occupied in DFT. The boundary between these two regions is exactly the DFT Fermi-surface marked with a cyan line. In spite of the fact that both cuprates can be modeled using only a single band, there is a significant correction to the charge density by DMFT. In Fig. 5(b), the correlation induced charge redistribution ∆ρ(r) is depicted as an iso-surface plot. A negative sign of ∆ρ(r) (cyan) suggests electron loss from the single band comprised of Cu-d x 2 −y 2 and O-p orbitals, a positive sign (yellow) electron gain. There are gains as well as losses around both Cu and O sites; altogether some net charge is transferred from Cu to O.
Let us hence focus on the O charge gain and Cu loss in the following. Around the O sites the gain has the form of a p x -and p y -orbital density pointing to the Cu sites. It stems from the admixture of these orbitals to our single band. On the Cu site in turn there are blue d x 2 −y 2 -like lobes of removed charge pointing towards the neighboring oxygen sites. This indicates that the redistribution ∆N W of Wannier orbitals in k-space, effectively reduces the level of admixture between these orbitals when moving from the DFT metal to the DMFT insulator. Even though one might naively assume that in a single Wannier band with fixed occupation CSC effects are minor, the density correction is significant for both cuprates.

C. SrVO3
SrVO 3 crystallizes in a cubic perovskite lattice structure and has been the testbed material for DFT+DMFT [36][37][38][39][40] and GW+DMFT 41-44 method development. A strong interplay between the octahedral crystal field in VO 6 and electron correlation determines the properties of this material. Bulk SrVO 3 exhibits robust metallicity also upon chemical substitution such as in Ca 1−x Sr x VO 3 35 . However, the material undergoes a metal-insulator transition if manipulating its dimensionality 45 . Ultra thin layers (up to 3 monolayers) of SrVO 3 are insulating, which opens the possibility to control the metal-insulator transition by applied electric field or strain, paving the way for a Mott transistor 46 .
In ultra-thin layers the bulk t 2g symmetry is broken: the out-of-plane d xz/yz orbitals have a reduced bandwidth, while the in-plane d xy bandwidth remains almost unchanged. Given the 3d 1 electronic configuration of vanadium, ultra-thin layers hence favor the electrons to be placed in the d xy orbital. An orbital polarization develops. The orbital reoccupation is quite dramatic in DFT+DMFT: from 1/3 for all t 2g Wannier orbitals in the metallic bulk to almost a occupation of 1 electron in the d xy orbital for a on-layer film. Let us note that DFT underestimates the orbital polarization which is strongly enhanced by electronic correlations: in DFT d xy and d xz/yz orbitals have 0.6 and 0.2 electrons, respectively, for the ultra-thin film; and it is metallic. Hence a freestanding monolayer of SrVO 3 is ideally suited to study the CSC DFT+DMFT charge redistribution caused by an orbital polarization. In our calculation we set the lattice constant (a=3.92Å) to that of a single SrVO 3 layer on SrTiO 3 as this is the experimental substrate 45 . Fig. 6 (lower panel) shows the DFT density of states (DOS); and the dashed lines represent the DOS of the V-t 2g orbitals in DFT. Let us note that the position and width of the d xy band (red) is significantly larger compared to that of the d xz/yz bands (blue). This leads to a first tendency towards an orbital polarization already in DFT, giving the aforementioned occupations.
For the DMFT calculations at inverse temperature β = 40, we employ the Kanamori interaction parameters U = 4.0 eV, J = 0.75 eV from the literature. 46 The effect of electron correlations in Fig. 6 (upper panel) is twofold: (i) the system becomes a Mott insulator, and (ii) in the insulating phase, d xy is half-filled while the other orbitals are essentially empty. This kind of physics has been observed before 46 but let us now turn to the effect of CSC.
The dashed curves represent the spectral function corresponding to the one-shot DFT+DMFT calculation without CSC, while the solid curves represent the full CSC results. Let us note that the insulating energy gap is slightly reduced by the CSC, which can be attributed to the charge redistribution within the t 2g manifold shown as an inset in Fig. 6. As is to be expected, the charge redistribution has a positive d xy -like shape as these orbitals become more occupied. Perpendicular to the plane, shown in the inset, there is a reduced charge in a d xz -and d yz -like shape.
This changed orbital occupation influences, in turn, the DFT electronic structure, see solid curves in Fig. 6 (lower panel): The d xy orbital that is more occupied in DMFT, is shifted upwards to higher energies in DFT and vice versa for d xz/yz , as is to be expected already from the Hartree term. That is, DFT partially compensates the correlation effect of DMFT, but a large net effect remains. This net effect is shown in Fig. 6: the charge redistribution in the inset, the CSC DFT and DMFT results as solid lines in the main panel. One shot DFT+DMFT has almost filled d xy orbital and almost empty d xz -and d yz -orbitals, while full CSC results in a slight reduction of d xy -orbital occupancy, compared to that in one-shot calculation and vice versa for d xz -and d yz -orbitals. The effect of full CSC is not negligible, also regarding the reduction of the band-gap.

IV. SUMMARY AND CONCLUSIONS
We have implemented a fully charge self-consistent DFT+DMFT method, using maximally localized Wannier functions constructed with Wannier90 6 , the Wien2k program package, wien2wannier as an interface, and w2dynamcis as an impurity solver. We applied the method to strongly correlated electron systems and discussed different physical and technical aspects.
The cuprates, Sr 2 CuTeO 6 and HgBa 2 CuO 4 , can be modeled by a single Wannier orbital. In this situation, one might assume that the charge self-consistency has no effect since this single orbital must remain halffilled; there is no charge redistribution to other orbitals. Nonetheless, the real space charge density is changed with full CSC DFT+DMFT.
In both cuprates, charge is removed from around the Cu site and added around the O sites. Note that oxygen p-states are mixed into the single, predominantly d x 2 −y 2 orbital. The reason for this change is a change of occupation of the Wannier orbitals in k-space. While for the metallic DFT solution, Wannier functions in some part of the the Brillouin Zone are singly occupied, in DMFT the band splits into two Hubbard bands and all k-points are occupied equally with half an electron.
Besides this common ground, there are some differences between Sr 2 CuTeO 6 and HgBa 2 CuO 4 . The former has much weaker p-d hybridization and itinerancy. Hence we do not need disentanglement to Wannier project onto the single orbital. As for the CSC, the changes at the oxygen are much less pronounced because the O states admix to a much lesser extent in Sr 2 CuTeO 6 . It also has a lower symmetry, which results in a more complicated charge redistribution pattern.
A significant correlation-induced occupation redistribution within the V-t 2g manifold is observed in a single layer of SrVO 3 . Here, the interplay between crystal field and electron correlation results in a pronounced orbital polarization. The orbital polarization can be clearly identified in the charge redistribution. The CSC has the tendency to counteract the DMFT orbital polarization, which is however hardly reduced at self-consistency with respect to one-shot DFT+DMFT.
In all the cases, using single or multi-band models, ∆N W (k) has a significant k-dependence which translates to an r-dependence of ∆ρ(r). This shows that there are more profound effects of CSC in DFT+DMFT than the gross effect of charge redistribution from one site to another found in previous studies.