Calculation of screened Coulomb interaction parameters for the charge-disproportionated insulator CaFeO 3

,


I. INTRODUCTION
Strongly correlated materials exhibit many interesting physical phenomena, such as high-temperature superconductivity, colossal magnetoresistance, and metalinsulator transitions [1][2][3], which makes them also very attractive candidates for a variety of technological applications, such as, e.g., Mott-transistors [4][5][6][7].Here, the idea is to exploit metal-insulator transitions for achieving higher carrier densities, larger on-off ratios, lower switching voltage, and faster switching times.The emergence of a metal-insulator transition in strongly correlated materials is typically associated with a large on-site Coulomb repulsion, represented by the Hubbard parameter U , which forces the electrons to localize and thereby produces an insulating state.A metal-insulator transition can then be triggered by variations in temperature, strain, doping, applied voltage, etc.
However, for certain transition-metal oxides exhibiting metal-insulator transitions, it was suggested that their basic physics can be explained within a minimal model where the Hubbard U is in fact strongly screened and therefore rather small, and instead the metal-insulator transition is controlled by the strength of the Hund's interaction J [8][9][10][11][12].In these cases, the metal-insulator transition is accompanied by a charge disproportionation of the transition-metal cations, resulting in crystallographically inequivalent sites and a breathing distortion of the surrounding oxygen polyhedra (see Fig. 1).
Such strong screening of the Coulomb repulsion between the transition-metal d electrons can be caused by the O-p electrons and is the more efficient the closer in Figure 1.Experimental crystal structure of CaFeO3 below the metal-insulator transition [22].The breathing distortion creates long-bond (light orange) and short-bond (dark orange) FeO6 octahedra, which are rotated and tilted by other distortions.
energy the O-p states are to the correlated d states.This makes materials with a small or even negative chargetransfer energy good candidates for this type of behavior [8,13], such as oxides containing late transition-metal cations with high oxidation states.Indeed, rare-earth nickelates [14][15][16] and certain ferrites [17][18][19] have been identified as charge-disproportionated insulators.For the case of the well-studied rare-earth nickelates, RNiO 3 , the strong screening of the local Coulomb repulsion was confirmed by ab initio calculations of the screened interaction parameters [20,21].
Long before the relevance of charge disproportionation for the metal-insulator transition in the rare-earth nickelates has been recognized, a metal-insulator transition around room temperature involving charge disproportionation was reported for paramagnetic CaFeO 3 arXiv:2310.16541v2[cond-mat.str-el]12 Feb 2024 [17,23].Here, the Fe cation has a nominal average oxidation state of 4+ and disproportionates according to 2d 4 → d 3 + d 5 , with a high-spin d 5 configuration at ambient pressure [24].The charge disproportionation in CaFeO 3 also couples to a breathing distortion of the FeO 6 octahedra [22,25,26], resulting in a symmetry lowering from P bnm to P 2 1 /n and a transition from metal to insulator [23,27].The low-temperature P 2 1 /n crystal structure of CaFeO 3 is depicted in Fig. 1.
Several computational studies based on densityfunctional theory (DFT) have addressed the charge disproportionation and metal-insulator transition in CaFeO 3 in order to clarify the underlying physics [28][29][30][31][32][33][34][35][36].Thereby, to obtain a realistic description of the insulating state, the local electron-electron interaction is typically explicitly treated within DFT+U [37] or DFT+dynamical mean-field theory (DMFT) [38,39].Different empirically chosen interaction parameters have been used in these calculations, ranging from U = 3 eV to more than 7 eV for the Hubbard parameter, and from J = 0 eV to as much as 2 eV for the Hund's interaction.A first-principles-based calculation of screened interaction parameters using cRPA can also give insights on whether these choices are physically reasonable.In our recent tight-binding+DMFT study [12], we have found a very rich phase diagram for paramagnetic CaFeO 3 as a function of U and J, which shows that using suitable interaction parameters is crucial for obtaining the correct phase.However, we also note that different studies often use different definitions of the correlated local Fe-d orbitals, which hampers a systematic comparison of the corresponding interaction parameters.
Here, we calculate screened interaction parameters for CaFeO 3 using the constrained random-phase approximation (cRPA) [40][41][42].Our results confirm a strong screening of the Hubbard parameter U , while the Hund's interaction J is less affected by the screening.We also compare different choices of the correlated d orbitals, corresponding either to an atomic-orbital-like basis, as commonly used in DFT+U , or to more molecular-orbitallike hybridized Fe-d/O-p frontier orbitals, which closely resemble the nominal oxidation states and are more commonly used in model Hamiltonian or DFT+DMFT studies.Furthermore, we assess the validity of the typicallyused spherical approximation for the local electronelectron interaction and discuss the applicability of the calculated screened interaction parameters in DFT+U and DFT+DMFT calculations of CaFeO 3 .To that end, we also revisit the U (and J) dependence of the relaxed high-and low-temperature structure of CaFeO 3 within DFT+U and compare the interaction parameters used for these structural relaxations with the ones obtained within cRPA.

II. METHODS
In this section, we introduce the interaction Hamiltonian acting on the correlated subspace of Fe-d electrons and discuss the determination of the corresponding screened interaction parameters using cRPA.Then, we describe the details and all parameters used in our DFT(+U ) calculations and in the construction of the Wannier basis representing the correlated orbitals.

A. Interaction Hamiltonian and constrained random-phase approximation
The general local two-particle interaction Hamiltonian can be written as follows (see, e.g., [43]): where c mσ is the annihilation operator of an electron with spin σ in orbital m and the matrix elements U mm ′ m ′′ m ′′′ are defined through Here, ϕ m (r) are the orbitals and v(r, r ′ ) is the interaction between these orbitals, which in the following can be either the bare or the screened Coulomb interaction.For practical calculations, the full four-index form of U mm ′ m ′′ m ′′′ is sometimes reduced to a two-index form according to U mm ′ = U mm ′ mm ′ and J mm ′ = U mm ′ m ′ m .These two-index terms are assumed to dominate the interaction.The third possible two-index matrix I mm ′ = U mmm ′ m ′ equals J mm ′ if the orbitals ϕ m (r) are realvalued, which can be directly seen from Eq. (2).In the following, we always use real-valued, maximally localized Wannier functions [44], and thus we only discuss U mm ′ and J mm ′ .
The interpretation of these two-index terms becomes clear when considering only the density-density terms of the Hamiltonian with the density operators n mσ = c † mσ c mσ and the convention σ = −σ.Here, U mm ′ = U σ σ mm ′ represents the interaction between electrons with opposite spins and U mm ′ − J mm ′ = U σσ mm ′ the interaction between electrons with the same spin.
To obtain screened interaction parameters within cRPA [40][41][42]45], the Hilbert space is divided into the correlated subspace, represented by the orbitals ϕ m (r) for which H int is defined, and the rest.Then, the bare electron-electron interaction v(r, r ′ ) = |r − r ′ | −1 is screened by the polarization function P rest including all electronic transitions that involve the rest of the Hilbert space, but not those purely within the correlated subspace itself.The screened Coulomb interaction U is then given by U Typically, P rest is calculated as P total − P corr , where P total is the total polarization, which can easily be calculated from the Kohn-Sham wave-functions, and P corr is the polarization only within the correlated subspace.If the correlated bands are not entangled with the rest, P corr is well-defined and can also be calculated from the Kohn-Sham wave-functions of the corresponding bands.If there is entanglement between the correlated and uncorrelated subspace, different approaches have been introduced to calculate the polarization function.
One option is to neglect the hybridization between the correlated subspace and the rest subspace in the Hamiltonian, and diagonalize these two subspaces individually.This leads to two independent sets of wave functions and eigenenergies, which for small hybridization differ only slightly from the original Kohn-Sham bands.The polarizations P total and P corr can then be calculated from these new wave functions and eigenenergies [41].This method is called the "disentanglement method" in the following.
Alternatively, Ref. 42 suggested to obtain P corr by assigning weights to each transition in the polarization function that measure the probability of an electron residing in the correlated suspace both before and after the transition.These weights are based on the projections of the Bloch functions onto the correlated orbitals.We refer to this method as the "weighted method".
Finally, Ref. 45 proposed a slightly different expression for P corr , derived from the general Kubo-Nakano formula for the response function corresponding to density fluctuations within the correlated subspace.We call this the "projector method" in the remainder of this work.It also is the default option implemented in the Vienna Ab-initio Simulation Package (VASP) [46,47] and is used throughout this work, except in Sec.III C 3, where we present an explicit comparison between all three methods for obtaining the polarization function.

B. Spherical parametrization and symmetrization of the interaction
For DFT+U and DFT+DMFT calculations, the interaction matrix U mm ′ m ′′ m ′′′ is typically constructed assuming spherical symmetry, i.e, v(r, r ′ ) = v(|r − r ′ |) and atomic-like orbitals ϕ m (r), where m is the magnetic quantum number corresponding to an electronic shell with well-defined orbital momentum l.This allows to expand the interaction v in terms of complex spherical harmonics Y kq and to arrive at a unique parametrization of the interaction matrix in terms of l + 1 independent parameters, the Slater integrals F k , where k is an even integer with 0 ≤ k ≤ 2l (see, e.g., [43,48,49]).The resulting interaction matrix in the basis of complex spherical harmonics is where the "2 × 3 matrices" are Wigner 3-j symbols.
The Slater integrals are in principle defined as the radial integrals appearing in this expansion.Conversely, they can be extracted from the interaction matrix by inverting Eq. ( 4) and using the orthogonality relation of the α tensor as specified in Eq. (A.1): Eq. ( 6) allows us to reduce any interaction matrix, e.g., one obtained from cRPA to a set of Slater parameters, which can then be used in Eq. ( 4) to construct a spherically symmetrized interaction matrix 1 .Note that while these formulas are specified in the basis of complex spherical harmonics, in Sec.III we represent the interaction matrices in the basis of real cubic harmonics.These two types of functions are related by a simple linear transformation.
Often, the parametrization of the interaction is further simplified to only two parameters U and J, which are defined in terms of the two-index matrices: From Eq. ( 6) and Eq. ( 7), it follows directly that F 0 ≡ U and, for l = 2, J ≡ (F 2 + F 4 )/14.(The general relation between Slater integrals and J for all l is given in Eq. (A.2)).However, the inverse operation of obtaining the interaction matrix from only U and J is underdefined for l ≥ 2 since this requires knowledge of all l + 1 Slater integrals.Therefore, a fixed ratio of F 4 /F 2 ≈ 0.63 is usually assumed for d electrons (see, e.g., [43]), which is based on calculations of the unscreened interaction using atomic-orbital-like orbitals [50,51].
Within cRPA, it is straightforward to not only calculate the local elements of the interaction matrix between orbitals on the same site but also between orbitals on different sites in the same unit cell.For practical purposes, it is often assumed that the intersite terms are not strongly orbitally dependent [52].Therefore, we only discuss the orbitally averaged nearest-neighbor (NN) interaction where now m and m ′ are on NN sites.

C. DFT calculations and Wannier orbitals
For all calculations, we use the VASP code (version 6.3.0)[46,47] in combination with the PBEsol exchangecorrelation functional [53].We use standard projector augmented wave (PAW) potentials [54,55] with valence states 3s 2 3p 6 4s 2 for Ca, 3s 2 3p 6 3d 7 4s 1 for Fe, and 2s 2 2p 4 for O.We treat CaFeO 3 in a 20-atom unit cell, which allows us to describe both the P bnm and P 2 1 /n structures, using a 7 × 7 × 5 k-point mesh, a plane-wave energy cutoff of 600 eV, and an energy tolerance of 10 −8 eV between electronic steps, which leads to good convergence.
Structural relaxations are performed up to an energy tolerance of 10 −7 eV between ionic steps.Except where explicitly noted, relaxations of the unit cell and atomic positions are done using DFT+U [56] with U = 4 eV and J = 1 eV, similar to previous studies [29,32,34], but we also present results where we systematically vary U over a wider range.For the structural relaxation, we emulate the complex low-temperature helical magnetic order in CaFeO 3 [22,57] with the much simpler A-type antiferromagnetic ordering, which has the lowest energy out of A-, C-and G-type orderings in our calculations.For the high-temperature P bnm structure, we also perform calculations where we relax only the atomic positions using spin-degenerate DFT without +U correction, which we call "bare DFT" here.
We compare the relaxed to the experimental structures using the mode-decomposition tool ISODISTORT [58,59].All amplitudes are normalized to the 20-atom unit cell, and the symmetry-labeling convention uses a cubic perovskite parent cell with Fe at the origin.The most important modes are the R + 1 breathing mode in the P 2 1 /n structure, which splits the Fe into two inequivalent sites, and the M + 2 Jahn-Teller mode in the P bnm structure, which locally splits the Fe-e g orbitals but experimentally is not relevant in CaFeO 3 .Also present in P bnm are the modes related to octahedral rotations and tilting, M + 3 , R + 4 , and X + 5 , which mainly decrease the hopping and therefore the band width [29].Additionally, there are three more symmetry-allowed modes with very small amplitudes: R + 5 (P bnm), R + 3 , and M + 5 (both We use Wannier90 (version 3.1.0)[60] to construct the Wannier orbitals defining the correlated subspace for our cRPA calculations.We start with initial projections corresponding to d-shell cubic harmonics defined in a local coordinate system aligned along the Fe-Fe distances.Thereby, we employ a disentanglement (outer) window from −10 to 5 eV relative to the Fermi energy, which is chosen to be large enough to encompass all 56 O-p and Fe-d derived bands, and, optionally, a frozen (inner) window from −1.1 to 5 eV.We always first perform the disentanglement, where convergence is defined as a relative change of less than 10 −9 of the gauge-invariant part of the spread, and then the localization (Wannierization) down to an absolute change of less than 10 −9 Å 2 of the total spread.Note that the resulting Wannier functions still resemble the cubic harmonics, which is required to map the cRPA results onto Eq. ( 6).
For the cRPA calculations, we first converge a nonmagnetic DFT calculation with 512 bands and otherwise the same parameters as listed above.Then we perform one exact-diagonalization step in DFT to obtain the Kohn-Sham band energies.We perform cRPA calculations only for the zero frequency component of the screened interaction and with the correlated subspace defined by the Wannier orbitals.

A. Structural relaxations
We first perform full structural relaxations of CaFeO 3 within both P bnm and P 2 1 /n symmetries using DFT+U and A-type antiferromagnetic order, as described in Sec.II C. The corresponding results are listed in Table I, together with experimental data from Ref. 22.
It can be seen that the lattice parameters obtained for the P bnm and P 2 1 /n structures differ only very little, Table I.Lattice parameters and distortion mode amplitudes (with respect to the cubic P m 3m high symmetry reference structure) obtained from DFT+U structural relaxations within P bnm and P 21/n symmetries.The third line corresponds to a calculation where lattice parameters are fixed to those obtained from DFT+U and only the atomic positions are relaxed using spin-degenerate DFT without a "+U " correction.The last two lines correspond to experimental data from Ref. 22 obtained at 15 K (P 21/n) and 300 K (P bnm).The M + 5 , R + 1 and R + 3 modes are not allowed in P bnm symmetry.

Unit cell
Mode amplitudes (Å) with only a small deviation from β = 90 • in the monoclinic case.Furthermore, the agreement with the experimental lattice parameters is very good (there is only a slight underestimation of the c lattice parameter by 2.1 % and the monoclinic distortion is slightly lower than 90 • rather than larger).
The calculated mode amplitudes also agree very well with the experimental data, except for the M + 2 Jahn-Teller mode in the P bnm case, which is essentially zero in the room-temperature structure of Ref. 22, whereas a noticeable Jahn-Teller distortion develops in the zero temperature DFT+U A-type antiferromagnetic case.This is due to the nominally Jahn-Teller-active d 4 electron configuration of the Fe 4+ cation, and has also been found with a very similar amplitude of 0.26 Å in a previous DFT+U study of CaFeO 3 if the symmetry is constrained to P bnm [35].
In the P 2 1 /n structure, which is about 11 meV per formula unit lower in energy than the P bnm structure, the M + 2 mode nearly vanishes and instead the R + 1 breathing mode appears with a magnitude of 0.16 Å, in very good agreement with the experimental data.Furthermore, it can be seen that the modes related to octahedral tilts and rotations, R + 4 and M + 3 (and also X + 5 ), do not change significantly between P bnm and P 2 1 /n, and also agree well with the experimental data.
Fig. 2 shows the evolution of the relaxed R + 1 breathing mode amplitude and the local magnetic moments of the two inequivalent Fe sites as function of the Hubbard U parameter for fixed J = 1 eV in the P 2 1 /n structure.We note that within the DFT+U method, U and J refer to rather localized atomic-like orbitals (similar to what we refer to as the localized basis in Sec.III B).This basis does not correspond to the minimal models for charge disproportionation discussed, e.g., in Refs.8-12 where the Hund's interaction drives the transition to the chargedisproportionated insulating state.Instead, as seen in Fig. 2, the Hubbard U remains the main control parameter in these DFT+U calculations.The R + 1 amplitude becomes nonzero for U ≥ 2 eV and is accompanied by a sizeable difference in the local magnetic moments, indicating the disproportionation of the Fe sites into nominal d 3 and d 5 valences.For 4 eV ≤ U ≤ 7 eV the system is in-  sulating and the R + 1 amplitude stays nearly constant.For larger U > 7 eV, the breathing mode amplitude decreases and the system becomes metallic again.This is similar to what has been observed for the rare-earth nickelates in Ref. 61, where a large U shifts the transition-metal d levels below the O-p states in energy, which then disfavors the charge disproportionation.Thus, even though a small or maybe even slightly negative charge-transfer energy seems to be necessary for this type of charge disproportionation, it appears unfavorable if the charge-transfer energy becomes too negative.
Overall, our results show that structural parameters in good agreement with available experimental data can be obtained within DFT+U using interaction parameters of U = 4 eV and J = 1 eV and an appropriate magnetic order, consistent with previous work [32,34].
Finally, to asses the effect of the imposed A-type magnetic order on the high-temperature P bnm structure, we perform an additional relaxation for the spin-degenerate case using U = 0 eV, where we fix the lattice parameters to the ones obtained from the P bnm DFT+U relaxation and then relax only the internal atomic positions.We note that fixing the lattice parameters to the DFT+U values is necessary, since the spin-degenerate DFT cal- culation leads to a low-spin state of the Fe cation which would in turn lead to a significantly lower unit cell volume.Table I shows that, in this case, the M + 2 Jahn-Teller mode disappears almost completely, while all other mode amplitudes are rather similar to the magnetic DFT+U calculation.This shows that the Jahn-Teller mode is indeed stabilized by the +U correction in combination with the magnetic order, while the octahedral rotation modes are unaffected by this.Throughout the remainder of this work, we use the so-obtained structure as the hightemperature P bnm structure unless otherwise noted.

B. Band structure and the correlated subspace
Next, we discuss the electronic structure obtained from spin-degenerate DFT calculations without a +U correction in the P bnm structure.Fig. 3 (a) and (b) show the total density of states (DOS), its projection on different atoms and orbitals, as well as the band structure around the Fermi level.The DOS shows three (partially overlap-ping) groups of bands with mixed atomic/orbital character due to hybridization.Bands with predominant Fee g character above 0.5 eV are separated from lower-lying bands with predominant Fe-t 2g character between about −1 to 0.5 eV due to a strong octahedral crystal field splitting.The Fe-t 2g dominated bands overlap slightly with O-p dominated bands below −1 eV.The fact that the Fe-d dominated bands are higher in energy than the O-p dominated ones but also overlap slightly indeed indicates a small positive but nearly vanishing charge-transfer energy within the classification scheme according to Zaanen, Sawatzky, and Allen [62].
We now construct two different sets of Wannier orbitals, starting from initial projections on Fe-centered d orbitals.The first set represents a so-called frontier basis that describes only the Fe-d dominated bands in the energy range above approximately −1 eV and directly maps on minimal models describing the basic physics of charge disproportionation [9,11,12,63].The second set corresponds to a localized basis, resembling an atomic-orbitallike basis similar to what is typically used to define the correlated orbitals within DFT+U calculations.
In both cases, we use an outer energy window comprising all bands shown in Fig. 3 (a) and (b).For the frontier basis, we use an additional inner (frozen) energy window for the disentanglement procedure which comprises all of the Fe-e g dominated bands as well as the upper part of the Fe-t 2g dominated bands down to −1.1 eV, as indicated in Fig. 3 (b) by the green shading.We note that since this inner energy window encompasses most of the Fe-d dominated bands, except for the narrow overlapping region immediately below −1 eV, the specific choice of the outer energy window does not have a strong influence on the resulting Wannier functions.The corresponding Wannier bands [green lines Fig. 3 (b)] remain restricted to the energy range of the Fe-d dominated bands, and perfectly reproduce the DFT bands [black lines in Fig. 3  (b)] within the inner energy window.Since this inner energy window contains the Fermi level and there is no entanglement with higher lying bands, the Wannier orbitals recover an occupation of exactly four electrons per Fe site, corresponding to the nominal valence of the Fe 4+ cation.
The resulting Wannier orbitals are centered on the Fe sites but also exhibit strong "tails" on the surrounding oxygen ligands, representing the strongly hybridized character of the corresponding bands.Fig. 3(c) depicts the (x 2 − y 2 )-type Wannier function as an example.One can recognize the antibonding character of the underlying hybridization between Fe-d and O-p orbitals.Such frontier orbitals are often used in DFT plus dynamical mean-field theory calculations or other more model-based approaches that consider only the effective d bands (see, e.g., [12]).
For the construction of the localized basis, we do not use a frozen inner window during the disentanglement.This results in more localized orbitals that closely resemble atomic orbitals and, as already mentioned, are similar to the orbitals used for the projected DOS shown in Fig. 3(a) and for the +U correction applied in the VASP code [54].Fig. 3(d) shows the resulting (x 2 − y 2 )type Wannier function, which does not exhibit the pronounced p-like tails on the surrounding oxygen ligands (small s-like tails remain, due to a corresponding weak hybridization and the orthogonalization between Wannier functions centered on different sites).
We note that the additional structural distortions in the low-temperature P 2 1 /n structure do not lead to substantial changes of the band structure around the lower bound of the frozen energy window, so that we can use the same energy window relative to the Fermi energy for both crystal structures.
C. Screened-interaction parameters from cRPA

Spherical symmetrization and effective interaction parameters
We first discuss the effect of the spherical symmetrization of the interaction parameters by comparing the matrix elements obtained directly from cRPA with their spherically averaged form, i.e., by using Eq. ( 6) and Eq. ( 4).Thereby, we focus on the two-index interaction matrices U σσ mm ′ and U σσ mm ′ between electrons with opposite and same spins, respectively.We note that the largest interaction that is not contained in the two-index matrices is U yz,z 2 ,yz,x 2 −y 2 = 0.19 eV (and equivalent matrix elements) in the frontier basis, and U yz,z 2 ,yz,x 2 −y 2 = 0.29 eV in the localized basis, which is clearly smaller than the dominant components of the two-index matrices (see below).
In the frontier basis, the two-index matrices obtained directly from cRPA are The order of the orbitals in these matrices is xy, yz, z 2 , xz, and It can be seen that there is a noticeable difference of around 0.21 to 0.26 eV between the intra-orbital terms (on the diagonal of U σσ cRPA ) corresponding to the t 2g -like orbitals and that of the e g -like orbitals.However, this orbital dependence of the intra-orbital interaction is weaker than what has been found for other materials [48,64].The deviations of the various (off-diagonal) inter-orbital terms from the symmetrized form are all smaller than 0.1 eV and even slightly smaller for the same-spin case.
In the localized basis, the corresponding interaction matrices are Here, there is almost no orbital dependence of the intraorbital entries in the cRPA matrices, and also the deviations of the inter-orbital matrix elements from their symmetrized form is rather small.Notably, the interaction for parallel spins has some relatively small, but negative inter-orbital entries, both before and after symmetrization.This could be a sign of overscreening, which is discussed further below.
Table II.Deviations in eV between the cRPA matrices and the spherically symmetrized two-index interaction matrices for same-spin (σσ) and opposite spin (σσ), both for the frontier and the localized basis.The symmetrization is performed either without or with the additional assumption F4/F2 = 0. To better quantify the deviations between the unsymmetrized and symmetrized two-index matrices, we determine the maximum absolute deviation of a single matrix element, ∆ max , and we also compute the Frobenius norm These deviations are summarized in Table II.Additionally, we also quantify the deviation between the cRPA values and an alternative spherically symmetrized interaction matrix, which is obtained by first calculating parameters U and J according to Eq. ( 7) and then assuming F 4 /F 2 = 0.63 to reconstruct the interaction matrix, as done in many DFT+U and DFT+DMFT implementations.
For the frontier basis, both ∆ max and ∆ Frob are noticeably larger for the opposite-spin interaction than for the same-spin case, due to the larger variations of the intra-orbital entries, as already discussed above.In the localized basis, the deviations from spherical symmetry are less pronounced for the opposite-spin interaction, and even slightly smaller than for the same-spin case.On the other hand, assuming F 4 /F 2 = 0.63 in the symmetrization has a stronger effect in the localized basis, while in the frontier basis this barely increases the overall deviation from spherical symmetry.The effect of the F 4 /F 2 ratio discussed further in Sec.III C 2.
Compared to the interaction matrices reported in Ref. 64 for LaFeAsO, the deviations from spherical symmetry are less pronounced in CaFeO 3 .This could be due to the different structure or coordination of the Fe cation in both systems (octahedral in CaFeO 3 versus tetrahedral in LaFeAsO), or also due to different degrees of covalency in the two systems.All in all, a spherical approximation of the local interaction appears to be relatively unproblematic for CaFeO 3 .
Next, we discuss the effective parameters obtained for the screened and unscreened interaction, which are listed in Table III.We focus on the Hubbard U and Hund J obtained directly by averaging according to Eq. ( 7) and start the discussion with the values for the frontier basis.The comparison of screened and unscreened values shows that U is strongly screened to U = 1.75 eV, which corresponds to only 12 % of its unscreened value, while J, as expected, is less screened to J = 0.53 eV, i.e., 78 % of its unscreened value.The rather strong screening of U in CaFeO 3 is indeed comparable to what has been found in cRPA calculations for rare-earth nickelates [20,21], and can be explained by the close vicinity between the correlated bands and the screening bands, which even overlap slightly [see Fig. 3(b)].
A systematic comparison with literature data is hampered by the use of different types of correlated subspaces (full d-shell vs. t 2g -or e g -only models, frontier vs. localized orbitals) or different parametrizations of the interaction Hamiltonian.The values obtained in Ref. 21 for LuNiO 3 , if transformed from the e g -only Kanamori parametrization used in [21] to the Slater parametrization used here, correspond to U = 1.3 eV and J = 0.5 eV, i.e., an even smaller U than in CaFeO 3 and a fairly similar J.Note that, due to the use of an e g -only correlated subspace in [21], the t 2g bands will also contribute to the screening.In cRPA calculations for other transitionmetal oxides, U is typically reduced to around 15-25 % of its bare value [48,65,66].The overall magnitude of the screened U and J in CaFeO 3 is also comparable to values obtained for Fe in LaFeAsO, with J = 0.59 eV, almost equal to CaFeO 3 , and a slightly larger U = 2.14 eV [64].
In the localized basis, the unscreened values of U and J increase significantly compared to the frontier basis, which reflects the higher Coulomb repulsion due to the stronger localization of the corresponding orbitals.However, and perhaps surprisingly, the screened value of U = 1.25 eV is noticeably smaller than for the frontier basis and corresponds to only 5 % of the unscreened value.In contrast, the screening of J in the localized basis is similar (in percent) to that in the frontier basis, resulting in a screened J = 0.92 eV, i.e., larger than in the frontier basis.
The extremely strong screening in the localized basis in CaFeO 3 can be understood from the strong band entanglement and hybridization.As a result, virtually every band inside the energy window has a mixed character, with contributions from both the correlated subspace and the screening subspace.This results in many potential screening channels close to the Fermi energy.To avoid this, some previous studies have used a screening subspace similar to that for the frontier orbitals also for the localized basis [48,64,67].While this reduces the screening on U , it also introduces an inconsistency between correlated and screening subspace, and thus we have not used this approach in our work.A small decrease in the screened U (and an increase in the screened J) in a localized basis compared to a frontier basis has also been observed for SrMoO 3 [66].However, in that case, the effect is much less extreme than here.Generally, cRPA is known to have a tendency to overscreen, which becomes more pronounced the less the correlated and screening subspaces are separated [68,69].Consequently, one can expect that the overscreening affects both basis sets, but is particularly strong in the localized basis.Finally, we note that our cRPA calculations are based on the low- spin state obtained within bare DFT, which might also influence the degree of screening.
Since the localized basis is similar to the correlated subspace used in our DFT+U calculations, the corresponding screened interaction parameters can also be compared to the values of U = 4 eV and J = 1 eV used for our structural relaxations in Sec.III A. Using the cRPA value of U = 1.25 eV in these calculations would lead to a rather poor description of the P 2 1 /n phase with a low-spin state of the Fe cations and essentially zero R + 1 mode, as shown in Fig. 2. On the other hand, the cRPA value of J = 0.92 eV is very close to our initial choice of J = 1 eV and (in combination with an appropriate U value) appears quite reasonable.We further discuss this partial mismatch between the interaction parameters obtained from cRPA and the values required in DFT+U for a good description of the charge-disproportionated state in CaFeO 3 in Sec.IV.
Additionally, we compute the orbitally averaged nearest-neighbor intersite interaction V NN and J NN , also shown in Table III.V NN assumes a moderate value of 0.59 eV in the frontier basis, screened to 16 % of its bare value, and comparable to LuNiO 3 [20] and LaFeAsO [64].In the localized basis, the unscreened V NN is comparable to that of the frontier basis, but is then completely screened, to only 1 % of the bare value, and thus becomes negligibly small.However, it is unclear whether this indeed indicates a highly efficient inter-site screening, and thus a good quality of the locality assumption for the screened interaction, or whether this values is affected by a potential overscreening.J NN is negligible in all cases (both screened and unscreened).
Finally, we also compute screened interaction parameters in the frontier basis for the relaxed structure with P 2 1 /n symmetry, and obtain only a negligible difference of 3 meV for U and less than 1 meV for J compared to the P bnm structure.

F4/F2 ratio in cRPA
Next, we analyze the ratio F 4 /F 2 obtained from the spherically averaged cRPA interaction parameters.Previous work has indicated that this ratio can deviate significantly from the value of 0.63 typically used to compute the three Slater parameters from U and J, and can reach up to 0.86 [48].As shown in Table III, this also applies to the case of CaFeO 3 .The unscreened F 4 /F 2 ratios are slightly smaller than the "atomic" value of 0.63, in particular for the localized basis, while in the screened case this ratio is increased to 0.80 and 0.76 for the localized and frontier basis, respectively.Thus, while a constant ratio of F 4 /F 2 = 0.63 still appears as a reasonable approximation for the unscreened interaction, this is no longer guaranteed in the presence of strong screening.
To assess whether a variation in F 4 /F 2 can also affect physical properties, we estimate the energy difference between the low-spin (LS) and high-spin (HS) states of the Fe 4+ or other d 4 cations in the zero band-width limit by evaluating the difference in the local interaction and crystal field energies for a t 4 2g e 0 g and a t 3 2g e 1 g configuration (see also the Appendix of Ref. 12) where ∆ CF is the local crystal-field splitting between e g and t 2g states.Thus, for the standard value of F 4 /F 2 = 0.63, a slightly larger J is necessary to obtain a high spin state.For a realistic crystal-field splitting ∆ CF = 2.5 eV, this corresponds to a difference of about 0.05 eV in the critical J.While this difference does not appear very large, it could be important for systems that are very close to a transition between the high-spin and low-spin states.

Different methods for disentangling the correlated and screening subspaces
Finally, as discussed in Sec.II A, for band-structures where the correlated bands are entangled with the rest, different methods for calculating the corresponding cRPA polarization functions have been suggested.Here, we compare the interaction parameters obtained from the "projector method" [45], which is the default method in this paper, the "weighted method" [42], and the "disentanglement method" [41], as implemented in VASP.The comparison of different interaction parameters is shown in Table IV.In the frontier basis, the projection and weighted methods give almost identical results for all interaction parameters, while the disentanglement method predicts larger values, and thus less screening, especially for the on-and intersite Coulomb repulsion U and V NN .The former two methods likely agree well because the Fe-dominated bands captured by the frontier basis only slightly overlap with the O-dominated bands.This in turn means that almost all bands clearly belong to either the correlated or uncorrelated subspace and only a few bands right below the frozen window at −1.1 eV partially belong to both subspaces.Furthermore, both methods work on the unaltered DFT band-structure.In contrast, neglecting the hybridization between correlated and uncorrelated subspaces, will slightly affect the band structure, in particular in the overlapping region below −1 eV, which can explain the more significant deviation of the disentanglement method compared to the other two method.
In the localized basis, the differences in interaction parameters, in particular U , are much more pronounced between the different methods.Again, the disentanglement method results in a significantly higher value of U and thus less screening.However, also the projection and weighted methods now differ quite drastically, with a very small U = 0.46 eV obtained from the latter.The large spread of the calculated interaction parameters for the localized basis is most likely a result of the strong band entanglement, such that small differences in the way the two subspaces are separated can lead to rather large differences in the corresponding polarization functions.
On first view, the interaction parameters obtained from the disentanglement method appear quite reasonable and compatible with the values used in our DFT+U calculations in Sec.III A. However, ignoring the hybridization between correlated and uncorrelated subspaces in the localized basis reduces the band width of the correlated bands to only 0.6 eV for the Fe-e g and 0.5 eV for the Fe-t 2g bands, and thus leads to severe changes compared to the DFT bands.This is very different from the situation for which this method was introduced [41].Therefore, this match of the interaction values seems likely to be a coincidence given the strong modification of the disentangled band structure.
The intersite terms are always negligibly small for the localized basis.

IV. SUMMARY AND DISCUSSION
In this work, we have presented cRPA calculations of the screened interaction parameters for the chargedisproportionated insulator CaFeO 3 .Thereby, we have compared results for two different types of Wannier functions used to describe the Fe-d states.First, a frontier orbital basis, incorporating the hybridization with the oxygen ligands, which describes only the bands immediately around the Fermi level and maps well on simple models involving only Fe cations with formal charge states.Second, a more localized basis resembling an atomic orbitallike basis, similar to the basis typically used in DFT+U calculations.
We have shown that cRPA predicts CaFeO 3 to have a strongly screened Hubbard U and a relatively large Hund's interaction J, consistent with other chargedisproportionated insulators such as the rare-earth nickelates [20,21].This supports the picture that the basic physics of these systems can be described within a minimal model where the insulating state is driven by the Hund's interaction rather than a strong Hubbard U [9,11,12,63].Additionally, the DFT band structure indicates that CaFeO 3 has a relatively small but positive charge-transfer energy.
We have also tested the quality of certain approximations typically applied to parametrize the interaction matrices in first-principles calculations, such as the assumption of spherical symmetry and a fixed ratio between the Slater integrals F 2 and F 4 .Indeed, the spherical approximation appears relatively uncritical, both for the more localized and the frontier basis, with the largest deviations from spherical symmetry observed for the intra-orbital interactions in the frontier basis, which vary by around 0.2 to 0.25 eV.The calculated F 4 /F 2 ratio was found to deviate noticeably from the commonly used fixed ratio of 0.63, in agreement with previous cRPA calculations for other transition-metal oxides [48].This can potentially become important when simulating a material close to a transition, e.g., between the high-and low-spin states.
In addition, we have discussed the choice of interaction parameters in DFT+U calculations, to obtain structural parameters for both the high-and low-temperature structures that are in good agreement with experimental data.We have shown that a simple A-type antiferromagnetic order and moderate electronic interactions suffice to reproduce the experimentally observed low-temperature structure.By subsequently neglecting the local electronelectron interactions and magnetic order, a good description of the experimental high-temperature structure can also be obtained.
However, our study also highlights some difficulties of directly using interaction parameters obtained from cRPA in Hubbard-corrected DFT-based methods such as DFT+U or DFT+DMFT, in particular for systems with strong entanglement between the correlated and the uncorrelated subspaces.The cRPA values obtained for U in the localized basis are clearly too small to obtain a good description of the charge-disproportionated state of CaFeO 3 in DFT+U calculations.While this apparent mismatch can, at least partially, be attributed to the well-established tendency for overscreening within the random-phase approximation [68,69], in the present case it is also, to a large extent, a result of the difficulty to disentangle the electronic bands between the correlated and uncorrelated subspaces.From our comparison between different methods in Sec.III C 3 it becomes clear that, in situations with strong entanglement, the results can vary widely and none of the methods is likely to produce interaction parameters that can reliably be used in firstprinciples calculations.Furthermore, screening is fundamentally a frequency-dependent phenomenon, while DFT+U (and also standard DFT+DMFT) uses constant, frequency-independent interaction parameters.It thus remains unclear whether using the zero-frequency value of the screened interaction is necessarily the best choice to be used in such calculations (see also Ref. 70).The further simplifications within the DFT+U method, i.e., the Hartree-Fock-like approximation for the local interaction and the treatment of intersite effects only within the standard semi-local DFT approximations, can also affect the optimal choice for the effective interaction parameters.Finally, it would be desirable to perform cRPA calculations for CaFeO 3 also in the high-spin state, either by constraining the orbital occupations to a t 3 2g e 1 g configuration or by performing cRPA calculations for a magnetically ordered state.While a solution of these issues goes beyond the scope of the present work, our results can serve as starting point for future attempts aimed at improving the currently available methods for calculating effective interaction parameters.

DATA AVAILABILITY
The supporting data for this article are openly available from the Materials Cloud Archive [71].

Figure 2 .
Figure 2. (a) R + 1 breathing mode amplitude and (b) local magnetic moments of the two inequivalent Fe sites obtained from full structural relaxation within P 21/n symmetry as a function of U at J = 1 eV.

Figure 3 .
Figure 3. (a) Total and projected density of states (DOS) and (b) band structure along high symmetry lines of the orthorhombic Brillouin zone of CaFeO3, obtained from the spindegenerate DFT calculation for the P bnm structure.The green shaded region in (b) indicates the frozen energy window used to construct the frontier-type Wannier orbitals, which result in the dispersion shown by the green lines.(c, d) Isosurfaces of the (x 2 − y 2 )-type Wannier orbital (c) in the frontier basis and (d) in the localized basis.Brown (red) spheres indicate Fe (O) atoms.

Table IV .
Screened interaction parameters obtained with different band disentanglement methods for the two different basis sets.We compare the parameters U and J, the ratio F4/F2, and the averaged intersite interaction VNN.