Effect of charge self-consistency in DFT + DMFT calculations for complex transition metal oxides

We investigate the effect of charge self-consistency (CSC) in density-functional theory plus dynamical mean-ﬁeld theory calculations compared to simpler “one-shot” calculations for materials where interaction effects lead to a strong redistribution of electronic charges between different orbitals or between different sites. We focus on two systems close to a metal-insulator transition (MIT), for which the importance of CSC is currently not well understood. Speciﬁcally, we analyze the strain-related orbital polarization in the correlated metal CaVO 3 and the spontaneous electronic charge disproportionation in the rare-earth nickelate LuNiO 3 . In both cases, we ﬁnd that the CSC treatment reduces the charge redistribution compared to cheaper one-shot calculations. However, while the MIT in CaVO 3 is only slightly shifted due to the reduced orbital polarization, the effect of the site polarization on the MIT in LuNiO 3 is more subtle. Furthermore, we highlight the role of the double-counting correction in CSC calculations containing different inequivalent sites.


I. INTRODUCTION
During recent years, the combination of density-functional theory (DFT) and dynamical mean-field theory (DMFT) has become a widespread tool to calculate properties of "correlated materials," i.e., materials where the strong Coulomb repulsion between electrons in partially filled d or f shells leads to effects that cannot easily be treated within effective noninteracting electron theories [1].The basic idea in combining DFT and DMFT is the assumption that for the relevant materials the electronic degrees of freedom can be separated into a "weakly interacting" part, for which a standard DFT treatment is adequate, and a "correlated subspace," which requires a more elaborate treatment of the electron-electron interaction.The latter leads, in general, to a redistribution of electrons within the correlated subspace compared to the DFT result.This change should then enter, in a self-consistent way, the effective potential felt by the weakly interacting electrons, which is achieved by iterating between DFT and DMFT steps.However, such a charge self-consistent (CSC) DFT+DMFT calculation leads to a higher computational cost compared to simpler "one-shot" (OS) calculations, where the charge rearrangement within the correlated subspace is neglected in the DFT calculation.
Thus, as DFT+DMFT develops further toward a standard ab initio-based computational method for materials science [2,3], it becomes essential to know in which cases it is possible to reduce the required computational effort by using more approximate variants of the method, e.g., by neglecting charge self-consistency.While CSC DFT+DMFT calculations have become more common recently, the DFT+DMFT method has also been applied to larger and more complex systems, such as, e.g., oxide heterostructures [4][5][6][7], defective systems [8][9][10], or large molecules [11,12].Therefore, a detailed understanding of the effect of charge self-consistency is desirable in order to better assess in which cases a CSC calculation is crucial or, more importantly, in which circumstances a one-shot calculation is sufficient.Unfortunately, there are currently very few studies available that provide a systematic quantitative comparison between CSC and one-shot calculations.It is the purpose of the present work to start filling this gap.
It can be assumed that charge self-consistency is particularly relevant for systems where correlation effects lead to a redistribution of electrons, e.g., for systems with charge and/or orbital ordering.For example, existing studies of epitaxially strained SrVO 3 monolayers demonstrate a reduced orbital polarization in CSC calculations compared to OS [13,14].Here we therefore analyze the effect of charge selfconsistency for two specific but representative cases.First, strained CaVO 3 , where orbital polarization is considered relevant for a reported strain-induced metal-insulator transition [15,16].And, second, LuNiO 3 , which is representative for a whole series of rare-earth nickelates that exhibit a transition to a charge-ordered (or charge-disproportionated) insulating state, which is also strongly coupled to a structural distortion [17].
Most previous work addressing the influence of charge self-consistency in DFT+DMFT calculations in transition metal (TM) oxides employed a so-called p-d model to define the correlated subspace [18][19][20][21][22][23], i.e., using a basis of rather localized, atomic-like orbitals constructed from a broad energy window that includes the TM d as well as all oxygen p bands.This appears conceptually appealing, since a wider energy window corresponds to a larger, and thus more complete, basis set, and since the use of rather localized orbitals provides better justification for the DMFT assumption of a purely local self-energy and Coulomb interaction [22].On the other hand, this also increases the computational load compared to using a minimal correlated subspace of "frontier" orbitals, corresponding to only a narrow energy window around the Fermi level.In TM oxides, the latter typically includes either t 2g or e g bands.
In the present work, we focus on DFT+DMFT calculations that employ such a minimal correlated subspace, corresponding to only a small number of near-Fermi-surface bands.These are expressed in a localized basis through a suitable transformation in terms of Wannier functions [24].By including only the minimal number of orbitals needed to describe the dominant low-energy physics within DMFT, this scheme requires a comparatively small computational cost.Furthermore, it often allows for an intuitive interpretation of Wannier occupations in terms of formal charge states, since the corresponding Wannier functions include the hybridization with the oxygen p states as "tails" located on the oxygen sites.
Another critical point arising in DFT+DMFT calculations using a p-d-type orbital subspace, is that the physically very important charge transfer energy, p-d , which describes the energy difference between oxygen p and transition metal d states, effectively becomes controlled by the so-called double-counting correction [18,23,25].The latter is required to account for the electron-electron interaction within the correlated subspace that is already included on the DFT level and is notoriously ill defined [26].Different expressions to account for the double counting (DC) have been suggested [23,27], but in some cases the double-counting needs to be adjusted manually, in order to obtain satisfactory results [18,19,28].It was shown that CSC calculations for such p-d-type calculations produce essentially the same spectral properties as OS calculations, if one tunes the DC correction to yield the same d-state occupancy [18].It is, however, not clear a priori, that more complex observables, e.g., the total energy, need to agree within both approaches.We note that the use of a minimal correlated subspace avoids the problem that the DC correction critically affects the important charge transfer energy, because charge neutrality between DFT and DMFT is ensured, and thus the DC potential shift can be absorbed in the chemical potential in DMFT [14].However, as we show in the following, the DC correction can still have a strong effect for systems with multiple inequivalent correlated sites.
In the next section (Sec.II), we provide a detailed description of the DFT+DMFT method as used in this work, specifying also all important computational parameters.We then discuss the two specific cases of CaVO 3 and LuNiO 3 in Sec.III, where we also provide a brief introduction in the relevant physical background for each of these materials.Finally, our main conclusions are summarized in Sec.IV.

A. DFT calculations
Structural relaxations for CaVO 3 within the 20 atom unit cell in Pbnm space group symmetry are performed using the QUANTUM ESPRESSO package [29].We employ scalarrelativistic ultrasoft pseudopotentials, with the 3s and 3p semicore states included in the valence for both V and Ca, together with the exchange-correlation functional according to Perdew, Burke, and Ernzerhof [30].Cell parameters and internal coordinates are relaxed until all force components are smaller than 0.1 mRy/a 0 (a 0 : Bohr radius) and all components of the stress tensor are smaller than 0.1 kbar.The plane-wave energy cutoff is set to 70 Ry for the wave functions and 840 Ry for the charge density.A 6 × 6 × 4 Monkhorst-Pack k-point grid is used to sample the Brillouin zone, and the Methfessel-Paxton scheme with a smearing parameter of 0.02 Ry is used to broaden electron occupations.For the calculation of epitaxially strained CaVO 3 , the in-plane lattice parameters are increased by 4% and kept fixed, while the c component of the cell and all atomic positions are relaxed.
All DFT calculations for LuNiO 3 as well as the DFT parts of all our CSC DFT+DMFT calculations are performed using the projector augmented wave (PAW) method [31], implemented in the Vienna Ab initio Simulation Package (VASP) [32][33][34], and also using the exchange correlation functional according to Perdew, Burke, and Ernzerhof [30].For Ni, we use the PAW potential where the 3p semicore states are included as valence electrons, while for Lu, we use the PAW potential corresponding to a 3+ valence state with f electrons frozen into the core.For the CaVO 3 calculations with VASP, we use the PAW potentials where the s and p semicore states are included as valence electrons for both Ca and V. Furthermore, a k-point mesh with 9 × 9 × 7 grid points along the three reciprocal lattice directions is used and a plane-wave energy cut-off of 550 eV is chosen for LuNiO 3 and 600 eV for CaVO 3 .The structure of LuNiO 3 is fully relaxed within Pbnm symmetry, both internal and lattice parameters, until the forces acting on all atoms are smaller than 10 −4 eV/Å.

Construction of the correlated subspace
In the DFT+DMFT method, the Kohn-Sham (KS) Hamiltonian within the chosen energy window is mapped onto a basis of localized states, spanning the correlated subspace C, then a local Coulomb interaction is added, and the resulting Hubbard-like lattice Hamiltonian is solved via the DMFT approximation [1,35].Without feedback to the DFT part, this corresponds to a OS calculation.To perform CSC calculations, one computes a correction to the charge density, ρ = ρ DMFT − ρ DFT , which is then passed back to the DFT code (here VASP) to calculate new KS wave functions and hence, update the correlated subspace [6,36].In a fully CSC calculation, this is repeated until ρ does not change compared to the previous iteration.
For the DMFT calculation, the electronic degrees of freedom within the chosen energy window are described via the interacting lattice Green's function: where μ is the chemical potential and ĤKS (k) is the KS Hamiltonian.The lattice self-energy ˆ (k, iω n ) is obtained by solving the effective DMFT impurity problem (see next subsection).
The lattice Green's function in Eq. ( 1) is expressed in the KS (Bloch) basis.To achieve the up-or down-folding between the quantities defined within the correlated subspace and the Green's function in the KS basis,  [14,36].To construct optimal projector functions, we apply the scheme introduced in Ref. [14], choosing a linear combination of the PAW partial wave augmentation channels that maximizes the overlap with the KS states inside a chosen energy window, which matches that of the correlated subspace C. We use initial projections on t 2g -or e g -like orbitals within the energy window of the correlated subspace C. The resulting projectors PL,ν (k) are in general not orthogonal to each other, and need to be orthonormalized: to define an orthonormal set of PLO-based Wannier functions describing the correlated subspace C. The orthonormalization of these PLO-based Wannier functions, as well as the whole DFT+DMFT self-consistency cycle has been implemented using the TRIQS/DFTTOOLS software package [37,38].The projectors P Lν (k) are updated after each DMFT cycle from the new KS states.Thereby, the energy window defining the correlated subspace is kept fixed, while monitoring that the change in the KS energies due to the charge density correction does not move the relevant bands outside of this energy window.The observed change of the KS eigenvalues is relatively small for all cases considered in this work, e.g., the maximum bandwidth reduction in LuNiO 3 is smaller than ∼5%.
We note that the strong octahedral rotations present within Pbnm symmetry lead to large off-diagonal crystal-field terms in the KS Hamiltonian, and the noninteracting Green's function for the effective impurity problem is no longer diagonal.Since this can induce a severe sign problem in the quantum Monte Carlo (QMC) method [39] used to solve the effective impurity problem (see next subsection), we perform a local unitary transformation of each impurity Green's function after the down-folding (respectively, before the up-folding), which diagonalizes the initial noninteracting local Hamiltonian on each site transforming the system into the crystal field basis.We note that the transformation is optimized in the first CSC cycle and is kept fixed in consecutive iterations to maintain a consistent orbital basis.
For CaVO 3 we also perform OS DFT+DMFT calculations based on the electronic structure obtained with QUAN-TUM ESPRESSO.In this case, the low-energy tight-binding Hamiltonian, used as input for the OS DMFT calculation, is formulated in the basis of maximally localized Wannier functions (MLWFs) [40] using the WANNIER90 code [41].Note that the PLO basis functions used in our VASP-based DFT+DMFT calculations essentially correspond to the initial Wannier functions constructed by WANNIER90 before the spread minimization, which are based on orthogonalized projections of (pseudo-) atomic orbitals on the Bloch bands [41].
The code used for all DFT+DMFT calculations in this paper is publicly available on github [42].

Solving the impurity problem
For both CaVO 3 (t 2g subspace) and LuNiO 3 (e g subspace) the effective impurity problem within the DMFT cycle is solved with a continuous-time QMC hybridization-expansion solver [39] implemented in TRIQS/CTHYB [43], taking into account all off-diagonal elements of the local Green's function in the crystal-field basis.For each impurity we add a local Coulomb interaction in the form of the Hubbard-Kanamori Hamiltonian [44], including all spin-flip and pair-hopping terms.Here the operator ĉ † mσ creates an electron in the atom-centered Wannier orbitals of type m and spin σ .The interaction parameters are given by the local intraorbital Coulomb repulsion U and the Hund's coupling J.To reduce the QMC noise in the high-frequency regime of the impurity self-energy imp and G imp , we represent both quantities in the Legendre basis [45] and sample the Legendre coefficients G l directly within the TRIQS/CTHYB solver.

Double counting correction
To correct the electron-electron interaction within the correlated subspace already accounted for within VASP, we use the fully localized limit DC correction scheme [46,47].Specifically, we use the parametrization given in Ref. [1] for the DC potential, where n α is the occupation of impurity site α and the average Coulomb interaction between M orbitals, Ū , is defined as The potential shift of Eq. ( 6) is added to the impurity selfenergy; its form is directly tailored to the Hubbard-Kanamori interaction Hamiltonian in Eq. ( 5) for a t 2g or e g model resulting from an octahedral crystal-field environment of M interacting orbitals (M = 3 for CaVO 3 and M = 2 for LuNiO 3 ).
In this work, we draw particular attention on how the occupations n α used for the DC correction are evaluated, i.e., whether they correspond to (a) the occupations of the Wannier functions as obtained from DFT or (b) the occupations corresponding to the impurity Green's function G imp calculated by the QMC solver within the DMFT step.It can be misleading to assume that these quantities are the same, even within a CSC calculation.Indeed, when the system is in a charge-ordered phase, such as, e.g., in heterostructures or nickelates, or in any other case with several inequivalent impurity problems, different impurities can exchange charge within the DMFT loop, potentially leading to drastic changes of the local occupations compared to the ones calculated within the DFT step.In principle, only the occupations evaluated for the impurity problem within DMFT that are used to define the charge density correction, have physical meaning within a CSC DFT+DMFT calculation.By contrast, the occupations obtained in the DFT part, by filling up the lowest energy KS states, do not correspond to the charge density that is used to evaluate the Kohn-Sham potential in a CSC calculation.However, in the case of a OS DFT+DMFT calculation, the question of whether to use DFT or DMFT occupations for the DC correction is ambiguous.An informal (and perhaps unrepresentative) community survey conducted by us, has shown that both variants are currently used in different studies.Here, we show that in certain systems the question of how to extract n α can have a strong influence on the results, and that one should be aware of this issue when evaluating the DC correction.

Calculation of observables
From the imaginary-time Green's function, we calculate the spectral weight around the Fermi level, Ā(ω = 0) = − β π Tr G C LL (β/2), which indicates whether the system is metallic [ Ā(0) > 0] or insulating [ Ā(0) ≈ 0] [48].For T = 0 (β → ∞), Ā is identical to the spectral function at ω = 0.For finite temperatures, it represents a weighted average around ω = 0 with a width of ∼k B T [48].The full real-frequency spectral function A(ω) is obtained via analytic continuation using the maximum entropy method [49].The on-site density matrix can be obtained directly from the local Matsubara Green's function as To extract the total energy of the system we use the following formula [24]: where KS ν,k are the KS eigenvalues with corresponding weights f νk within the correlated subspace C and • DMFT denotes quantities evaluated from the DMFT solution.The interaction energy Ĥint DMFT is calculated using the Galitskii-Migdal formula [50,51], and the last term in Eq. ( 8) subtracts the DC energy.To ensure high accuracy, we sample the total energy over a minimum of additional 60 converged DMFT iterations after the CSC DFT+DMFT loop is already converged.Convergence is reached when the standard error of the site occupation during the last 10 DFT+DMFT loops is smaller than 1.5 × 10 −3 .This way, we achieve an accuracy in the total energy of <5 meV.All DMFT calculations are performed for β = 40 eV −1 , which corresponds to a temperature of 290 K.

III. MATERIALS AND RESULTS
To analyze the effect of CSC within DFT+DMFT, we study two representative examples of TM oxides with different levels of complexity.First, we consider the case of unstrained and strained CaVO 3 .While in the former case this material is a correlated metal [52,53], it has recently been demonstrated that tensile epitaxial strain leads to a transition toward the Mott insulating state within OS DFT+DMFT calculations [16].An important aspect in this transition is the strain-induced crystal-field splitting between the partially filled t 2g orbitals, leading to a strong orbital polarization, and thus a local charge redistribution, which can potentially affect the result of a CSC compared to a OS DFT+DMFT calculation.However, in CaVO 3 , all correlated sites are symmetry equivalent and thus the DC correction is irrelevant when using a minimal "t 2g -only" correlated subspace.
Second, we consider the rare-earth nickelate LuNiO 3 , which exhibits a complex interplay between a specific structural distortion and an associated charge ordering, resulting in a transition from a paramagnetic metallic toward an also paramagnetic but insulating phase [17].In previous works, it was shown that DFT+DMFT is well suited to describe this correlated insulating state [20,54,55].Since two symmetryinequivalent types of Ni sites appear in the insulating state, this case allows us to analyze the effect of a site-dependent DC correction within CSC DFT+DMFT calculations.
Both materials, CaVO 3 and LuNiO 3 , exhibit a distorted perovskite structure with Pbnm space group (in the case of LuNiO 3 this corresponds to the high symmetry metallic phase).The corresponding unit cell contains four TM atoms surrounded by edge-connected oxygen octahedra, that are tilted and rotated around the Cartesian axes, corresponding to the so-called GdFeO 3 -type distortion (a − a − c + tilt system in Glazer notation), as depicted in Fig. 1.The d levels of the TM ions are split into e g and t 2g manifolds by the octahedral crystal field, and the remaining degeneracies can be further lifted by additional distortions of the oxygen octahedra (also shown schematically in Fig. 1).

A. CaVO 3 : Orbital polarization
As stated above, bulk CaVO 3 is a moderately correlated metal with weak orbital polarization that can undergo a transition to the Mott-insulating state under tensile epitaxial strain or in ultrathin films [15,16,56].As has been pointed out in Ref. [57], the orbital polarization resulting from the orthorhombic distortion of the perovskite structure (related to the tilts and rotations of the octahedral network) is an FIG. 1. Exemplary Pbnm crystal structure, as well as the main distortion modes relevant for the discussion of CaVO 3 and LuNiO 3 , respectively, i.e., tetragonal strain (left) and octahedral breathing mode distortion (right).The d-orbital energy levels that result from these distortions and the nominal occupations for each compound are also depicted.Note that for simplicity we omitted octahedral rotations in the simplified lattice structures and the corresponding d-level crystal-field splittings and that the local spin direction fluctuates in the paramagnetic case.
important factor in the metal-insulator transition (MIT).Several examples suggest that by an appropriate tuning of the bandwidth and the crystal-field splitting via, for example, strain or dimensional confinement, the resulting charge redistribution enhances the orbital polarization, eventually leading to a MIT [15,16,58].For example, as depicted in Fig. 1, tensile epitaxial strain will lift the degeneracy of the t 2g states, lowering the energy of one orbital compared to the other two.Since the orbital polarization in CaVO 3 can be seen as a measure for the likelihood of the Mott-insulating state, it is clear that describing this quantity accurately is essential for the success of the chosen method.
As described in Sec.II, we perform DFT+DMFT calculations for the bulk structure of CaVO 3 using three different schemes, i.e., OS calculations using either MLWFs (magenta line in Fig. 2) or PLOs (blue lines in Fig. 2) to represent the correlated subspace, as well as CSC calculations using PLOs (green lines in Fig. 2).From this we obtain the orbital occupations and spectral weight at the Fermi level, shown in Fig. 2, as a function of the Coulomb interaction parameter U .In all cases, the spectral weight is finite for small values of U , where the system is metallic and then becomes zero in the insulating phase for large U , with a rather sharp transition at U MIT .For the unstrained bulk system, all three approaches give identical results for the spectral weight as function of U , with a critical value of U MIT = 5.5 eV.Thus, at U ≈ 5 eV, which is typically considered as realistic value for 3d 1 transition metal oxides [53], we find a finite weight corresponding to metallic behavior, in agreement with experimental observations.This shows that the obtained results do not depend on details of the implementation, such as small differences in the basis used to represent the correlated subspace.
From the occupations shown in Fig. 2(a), it can be seen that the orbital polarization is weak in the metallic regime but is significantly enhanced above U MIT , where the occupation of one orbital is decreased compared to the other two orbitals.This is in line with the crystal-field splitting of the bulk on-site Wannier energies, where one orbital is energetically higher than the other two, with only a small difference between the latter [16].Here the two different OS calculations agree extremely well, while the orbital polarization is slightly reduced in the CSC calculation, however, with no apparent effect on the predicted U MIT .
Under 4% tensile strain [Figs.2(c) and 2(d)], the MIT is shifted to lower U values, below the realistic value of U ≈ 5 eV.Here both the MLWF-and PLO-type OS calculations agree within the accuracy of the method and give exactly the same value for the critical interaction parameter of U MIT = 4.7 eV.The CSC calculation, however, places the MIT at a slightly higher value of U MIT = 4.9 eV.
An even stronger difference between OS and CSC calculations can be seen in the orbital polarization, which is generally strongly enhanced compared to the unstrained case, due to a large strain-induced crystal-field splitting [16,58] (see Fig. 1).
Within the OS calculations, both PLO and MLWF based, we find that in the insulating regime two orbitals become completely empty, while the third one is essentially fully occupied by a single electron, i.e., the system exhibits full orbital polarization.In the CSC calculation this orbital polarization is significantly reduced, with a maximal occupation of ∼0.7 in the preferential orbital.The crystal-field-induced orbital polarization, enhanced by electronic interaction effects, has previously been suggested to be an important factor supporting the insulating phase [53], since the resulting effective half-filling of only one orbital promotes the MIT as opposed to fractional occupation of three degenerate levels.This is consistent with our results, since the lower orbital polarization in the CSC calculation correlates with a higher U MIT compared to the OS case.To illustrate the difference between OS and CSC calculations in the strained case, we plot the spectral function A(ω) at U = 5.0 eV for both cases in Fig. 3.Here the three different line styles correspond to the three different t 2g -like orbitals.As discussed previously, in the OS calculation one of the orbitals is essentially completely occupied, while the remaining two are empty.In contrast to this, the CSC calculation shows a correlation-induced charge redistribution from the occupied orbital to the previously empty orbitals.Furthermore, comparing the gap sizes of both cases, it is clearly visible that in the CSC case the gap is reduced compared to OS, similar to what has been reported in earlier studies on SrVO 3 [13].
Overall, we conclude that charge self-consistency plays only a minor role for systems with weak or vanishing orbital polarization, where the corresponding charge redistribution obtained within DMFT compared to the DFT calculation is small.In contrast, for systems with strong differences in orbital occupations, the OS calculation can lead to an overestimation of the orbital polarization, which in turn can affect the tendency of the system to undergo a MIT.While the effect on U MIT is not too strong in the present case, the corresponding differences in spectral properties can be more pronounced.Nevertheless, it appears that for the present case, OS calculations can at least give reliable qualitative information about the overall system behavior, such as, e.g., the effect of tensile epitaxial strain on U MIT , favoring the insulating state.
Furthermore, we note that in our calculations using frontier orbitals, we find very good agreement between the PLO-and MLWF-based methods both in the spectral properties and for the orbital occupations.This is in contrast to previous studies reporting that projector-based methods require a larger U in p-d models due to larger hybridization effects [19].

B. LuNiO 3 : Charge-ordering and structural energetics
The second case that we analyze is LuNiO 3 .This material belongs to the family of rare-earth nickelates, RNiO 3 , where R can be any rare-earth ion ranging from Lu to Pr, including Y.All members of the series exhibit a MIT, which is accompanied by a structural transition, lowering the space group symmetry from Pbnm to P2 1 /n.The corresponding structural distortion results in a three dimensional checkerboard-like arrangement of long bond (LB) and short bond (SB) NiO 6 octahedra, referred to as breathing mode distortion [59], and schematically shown on the right side of Fig. 1.Recent theoretical work indicates that this transition is related to an electronic instability toward spontaneous charge disproportionation on the Ni sites, which couples to the breathing mode, leading to a first-order coupled structural-electronic transition into a paramagnetic charge-disproportionated insulator (CDI) [60,61].Furthermore, the choice of the R site cation determines the degree of octahedral rotations in the corresponding high symmetry Pbnm structure, and thus the bandwidth.The latter then controls how close the system is to the electronic instability, driving trends across the series [60][61][62][63][64].
Here we use the case of LuNiO 3 to analyze if, and how, the charge disproportionation, as a specific example for chargeordering phenomena in general, is affected by the inclusion of charge self-consistency in DFT+DMFT.Earlier studies by Park et al. [65] also investigated the effect of CSC and DC for LuNiO 3 using a p-d-type subspace.They found only a small effect due to CSC on total energy calculations but had to adjust the DC correction to obtain a stable finite equilibrium breathing mode distortion.Here we use only the two e g -like frontier orbitals per Ni site for our DFT+DMFT calculations.As shown in Ref. [55], the electronic instability toward charge disproportionation and the resulting site-selective Mott transition [54] occurring in the paramagnetic state is well described within DFT+DMFT using such a minimal subspace.
Subedi et al. [55] found that the CDI state emerges in the frontier e g model for nickelates for rather large values of the Hund's coupling J and is very sensitive to its variation.The fact that the Hund's coupling J is the critical ingredient in the occurence of the CDI state was first proposed by Mazin et al. [66].They showed in an atomic picture that when U − 3J becomes small and is overcome by the potential energy difference between the Ni sites, s , which results from the breathing mode distortion and the charge disproportionation, the CDI state is favored.This regime is accessible in systems with small or negative charge-transfer gap, which results in a strong screening of the Coulomb interaction in the effective d bands, whereas the Hund's coupling is less sensitive to screening [66].A strong screening of U in nickelates has been confirmed by recent studies using the constrained random phase approximation (cRPA) [61,67].Moreover, in Ref. [68] it is shown that such a CDI regime for small or negative U − 3J is also accessible in a general three orbital Hubbard model and is thus not limited to nickelate systems.
To isolate the effect of the structural breathing mode distortion on the electronic charge disproportionation and the total energy of the system, we distinguish the various structural distortions present in LuNiO 3 by employing a symmetry-based mode decomposition [69], as outlined in Refs.[61,64,70].This allows to add the breathing mode distortion, with symmetry label R + 1 , on top of the relaxed Pbnm structure, and systematically vary its amplitude without changing any other parameter of the unit cell.We use the software ISODISTORT [71] to perform the mode decomposition.

Results for fixed structure
First, we calculate the properties of LuNiO 3 for a fixed structure, using the experimentally observed breathing mode amplitude, R + 1 = 0.075 Å [72], thereby varying the strength of the Hund's coupling J.As discussed above and shown in Ref. [55], the charge disproportionation and the resulting MIT depend sensitively on J, which thus allows us to critically examine the influence of CSC on the most crucial system properties.We use a fixed U value of 1.85 eV, which corresponds to the value calculated for LuNiO 3 using the cRPA [61,73].The results are depicted in Fig. 4, where in Fig. 4(a) the charge disproportionation, ν ≡ n LB − n SB , i.e., the difference of the e g occupation between the LB and SB Ni sites, is shown as function of J. Figure 4(b) shows the corresponding value for Ā(0), indicating whether the system is metallic or insulating.The dashed vertical line corresponds to the J value obtained within cRPA [61,73].Different data sets in Fig. 4 correspond to DFT+DMFT calculations with different treatments of the DC correction, both OS and CSC, which we discuss in the following.
We first focus on the data set labeled "CSC n DMFT α " (shown in red), which corresponds to the CSC calculation where the occupations entering the DC correction are calculated from the impurity occupations, and are updated in each DMFT iteration.As discussed in Sec.II B 3, this is the correct way to perform such CSC DFT+DMFT calculations, since the converged n DMFT α give rise to the corrected charge density from which the KS potential is constructed within the DFT step.It can be seen, that the transition to the CDI occurs at J = 0.2 eV, indicated by clear jumps in ν and Ā(0).The jump in ν is related to a drastic change in the DC potential difference between the Ni sites, since, for not too large J (see also below), the DC correction tends to increase the charge disproportionation by further lowering the e g states on the more occupied LB site compared to the less occupied SB site.This is discussed and analyzed in more detail in Appendix.
For further increasing J, ν stays almost constant until J ≈ 0.8 eV, where ν decreases again.Finally, at around J = 1.2 eV, the system becomes metallic again.This can be explained by the fact that for increasing J, the DC potential, proportional to Ū = U − 5  3 J [see Eq. ( 7) for M = 2], decreases, and eventually changes sign for J = 1.11 eV where Ū = 0. Thus, above J = 1.11 eV the DC correction opposes the charge disproportionation by lowering the e g levels of the SB sites relative to the LB sites.
Next we compare the CSC calculations with the simpler OS calculations.As discussed in Sec.II B 3, it is ambiguous whether to use the DMFT impurity occupations or the occupations of the Wannier functions obtained within DFT, n DFT α , to evaluate the DC correction.We first compare to the OS calculations where n DMFT α has been used for the DC correction (shown in green).It can be observed that in these OS calculations the system is already in the CDI state even for J = 0.2 eV.In addition, a small shift to larger ν can be observed compared to the CSC case.Thus, the tendency toward the CDI state is slightly stronger than in the CSC calculations.
In contrast, the OS calculations using n DFT α (shown in orange) leads to a significantly reduced ν, which increases slowly with increasing J.Moreover, for small J < 0.5 eV, clear metallic behavior is observed, while from J = 0.5 to 1.0 eV, the system undergoes the MIT, where eventually at J = 1.0 eV the system is completely in the CDI state with ν > 1.0.The occupations obtained in the initial DFT step are n DFT LB ≈ 1.15 and n DFT SB ≈ 0.85.For comparison, we also perform CSC calculations where the DFT occupations are used for the DC correction (shown in purple).However, one should note, that these calculations are somewhat artificial, since the DFT Wannier orbital occupations loose their physical meaning in a CSC calculation, and are used here just to allow for a more systematic comparison between OS and CSC calculations.One can see that overall the results of these calculations show similar behavior than the corresponding OS calculation using n DFT α , albeit with a small further reduction of ν.
The fixed structure calculations for LuNiO 3 , show that performing CSC calculations leads to a small reduction of the charge disproportionation compared to OS calculations, if in both calculations the DMFT impurity occupations are used to determine the DC potential.Moreover, we find that the DC has a very strong effect, so that a OS calculation with DFT occupations significantly underestimates the tendency toward charge disproportionation compared to the "correct" CSC calculation.
Overall, we conclude that CSC has a small, but certainly not negligible, influence on the DFT+DMFT calculations for LuNiO 3 , reducing ν by approximately 10%.However, this only holds if DMFT occupations are used in the OS calculation to evaluate the DC correction.If DFT occupations are used in the OS calculation, then the tendency toward the CDI state is significantly weakened, indicated by the much smaller ν, which is related to the smaller difference in the DC potential shift.However, compared to a hypothetical CSC calculation also using n DFT α for the DC correction, ν is again slightly enhanced in the OS calculation.Thus, one can clearly distinguish between the effect of the DC correction, and the effect of the charge density correction in the CSC calculation.The latter tends to reduce the charge disproportionation, independently of the chosen DC scheme, and analogous to reducing the orbital polarization in the case of CaVO 3 discussed in Sec.III A. Finally, our results also indicate that the OS calculations using DMFT occupations for the DC correction already provide a good approximation for the CSC calculation, even though they slightly overestimate the SB-LB splitting and thus the tendency toward the CDI state.

Influence on energetics
Another important aspect is the influence of charge selfconsistency in total energy calculations for different R + 1 amplitudes, i.e., for different amplitudes of the structural breathing mode distortion.As the R + 1 amplitude, and thus ν, changes, the DC potential and energy correction change accordingly.In addition, within the CSC calculation, the Hartree energy and other DFT energy contributions are evaluated from the corrected, self-consistent charge density.Strictly speaking, only a full CSC calculation gives physical meaningful total energies [26], but nevertheless we discuss the difference here to better understand the influence of performing full CSC calculation [65,74].To analyze the resulting effects, we again use U = 1.85 eV and two different values for J, 0.42 eV (the cRPA value) and 1.1 eV (where Ū ≈ 0 and thus the DC correction vanishes).For both cases, we compare OS with CSC calculations with different treatments of the DC correction, as introduced above.The results are shown in Fig. 5. Figures 5(a For the smaller value, J = 0.42 eV, both the OS (green) and CSC (red) result in an energy minimum at a finite R + 1 amplitude close to the experimental value (indicated by the vertical line).However, the OS calculation exhibits a much stronger response on the R + 1 amplitude, and hence shows a significantly deeper energy minimum.In contrast, the "artificial" CSC calculation using n DFT α for the DC correction (purple), exhibits no energy minimum for R + 1 > 0. Furthermore, the "correct" CSC calculation using n DMFT α undergoes a MIT to the CDI between R + 1 = 0 and R + 1 = 0.03 Å, while the corresponding OS calculation is already insulating without structural distortion and the CSC calculation with n DFT α remains metallic for any calculated R + 1 amplitude.For J = 1.1 eV, both CSC calculations, done either with DFT (purple) or DMFT occupations (red), agree very well (due to Ū ≈ 0 in the DC) and do not result in a stable finite breathing mode amplitude, even though both undergo a MIT at around R + 1 = 0.03 Å and exhibit a large charge disproportionation ν in the insulating state.In contrast, the OS calculation (orange), shows a stronger response, and predicts a breathing mode amplitude of R + 1 = 0.06 Å.Note that here we used n DFT α for the DC correction, but the same result would be obtained using n DMFT α , due to Ū ≈ 0. These results show that, even though the effect of charge self-consistency on ν for fixed crystal structure seems to be relatively minor, the effect on the energetics can be quite drastic, such that one can obtain a finite breathing mode distortion within a OS calculation, while the CSC calculation does not exhibit an energy minimum for R + 1 > 0.

IV. SUMMARY
We have studied the effect of charge self-consistency and the role of the DC correction within CSC DFT+DMFT calculations in two representative examples of transition metal oxides, using only a minimal correlated subspace corresponding to few frontier bands around the Fermi level.Our goal is to better understand in which cases charge self-consistency is really required in order to obtain accurate results, and in which cases a computationally much cheaper OS calculation might be sufficient.
For CaVO 3 , we find that the strong orbital polarization in the insulating phase under tensile strain is significantly overestimated by about 30% in OS compared to CSC calculations, in agreement with similar calculations for SrVO 3 in Refs.[13,14].This has a small but noticeable effect on U MIT , the critical U for the MIT, which is slightly underestimated in the OS calculations.In contrast, for the unstrained system, where the orbital polarization is much smaller, the difference between CSC and OS calculations is nearly negligible, even though also in this case the orbital polarization is slightly overestimated in OS calculations.Furthermore, we also compared OS calculations using PLO-based and MLWF-based schemes for constructing the correlated subspace and found very good agreement between the two methods.
While for CaVO 3 all TM sites are symmetry equivalent, and thus the site-dependent but orbitally independent DC correction does not affect the results, for the second example investigated in this work, LuNiO 3 , the DC correction becomes rather important.Here, we find that if DMFT occupations are used to evaluate the DC correction in the OS calculation, one can obtain results that are in rather good agreement with the CSC calculation, even though the charge disproportionation ν is overestimated by ∼10%.Thus, similar to reducing the orbital polarization for strained CaVO 3 , including charge selfconsistency leads to a somewhat more homogeneous charge distribution compared to a OS calculation.Nevertheless, it appears that in order to obtain qualitative insights or general trends, OS calculations can be a reasonable approximation, even in charge ordered systems, if the DMFT occupations are used for the DC.However, our analysis of the energetics of the breathing mode distortion shows that for certain observables, such as the total energy and resulting structural distortions, charge self-consistency can be crucial.For example in the case of LuNiO 3 , OS calculations overestimate the response on the R + 1 mode, in the most extreme case leading to a stable finite breathing mode amplitude, which is absent in the CSC calculation.In this case it is is inevitable to perform a full CSC calculation to obtain reliable results.
In summary, the effect of charge self-consistency is mainly to reduce a potential site or orbital polarization by favoring a more "homogeneous" distribution of electrons over all sites and/or orbitals.For the cases studied in this work, this results in a weak to moderate charge redistribution, which can be quantitatively relevant, depending on the specific application.In particular for total energy calculations, which depend on a subtle balance between different contributions, charge self-consistency can be crucial to obtain quantitatively and even qualitatively correct results.Nevertheless, it appears that cheaper OS calculations are often sufficient to gain insight into the system properties on a qualitative level, even though the, in principle ambiguous, choice of DFT or DMFT occupations to evaluate the DC correction in the OS calculations can become crucial.In the present examples, the use of DMFT occupations provided better agreement with the full CSC calculation, but in other cases this approach might also severely overestimate the electron transfer between inequivalent sites.
We hope that our detailed analysis of two specifically selected cases, provides useful insights for future DFT+DMFT studies of related material systems, thus allowing the treatment of larger and more complex materials systems by avoiding the higher computational cost of a CSC calculation when possible.
positive values of Ū (J < 1.11 eV), since then DC s is negative thus s > DFT s .The strong tendency to form the CDI state in the calculations with n DMFT α can be explained along the same lines.It can be seen that in the CSC calculations (red crosses), the CDI regime is entered already for J = 0.2 eV, and in the OS calculations (green stars) for even smaller J. Importantly, it can be seen that in these cases the DC potential jumps at the MIT, strongly favoring the CDI state.We note that of course the underlying atomic limit consideration neglects the important effect of the bandwidth, but it nevertheless can give a transparent picture of the underlying physics.

FIG. 2 .
FIG. 2. DFT+DMFT results obtained from OS and CSC calculations with VASP (PLO basis), compared to OS calculations using QUANTUM ESPRESSO (QE, MLWF basis), for bulk [(a) and (b)] and strained [(c) and (d)] CaVO 3 .Panels (a) and (c) show the orbitally resolved occupations as a function of the interaction parameter U .Panels (b) and (d) show the averaged spectral weight at the Fermi level, Ā(0).

FIG. 4 .
FIG. 4. Results of different DFT+DMFT calculations for LuNiO 3 using the experimental R + 1 amplitude, U = 1.85 eV, and varying J. CSC and OS calculations are labeled accordingly.For calculations labeled n DFT α (n DMFT α ) the DFT (DMFT) occupations have been used to evaluate the DC corrcetion.The dashed vertical line marks the cRPA value of J [61].(a) Charge disproportionation ν; (b) corresponding spectral weight at the Fermi level.
) and 5(c) show the total energy as function of the R + 1 amplitude, and Figs.5(b) and 5(d) show the corresponding Ā(0).

FIG. 5 .
FIG. 5. Comparison of energetics from DFT+DMFT for LuNiO 3 as function of the R + 1 amplitude.Calculations without CSC are labeled OS both in combination with DC occupations obtained from DFT (n DFT α ) or with occupations obtained from each DMFT step (n DMFT α ).Panels (a) and (b) show results for the small J = 0.42 eV, and panels (c) and (d) for the large J = 1.1 eV, where the upper panel shows the energy as function of the R + 1 amplitude and the panels at the bottom the corresponding spectral weight at the Fermi level.