On the effect of charge self-consistency in DFT+DMFT calculations for complex transition metal oxides
Alexander Hampel, Sophie Beck, Claude Ederer

TL;DR
This paper examines how charge self-consistency in DFT+DMFT calculations influences electronic charge redistribution and metal-insulator transitions in complex transition metal oxides, revealing nuanced effects on material properties.
Contribution
It provides a detailed comparison of charge self-consistent and one-shot DFT+DMFT calculations for materials near a metal-insulator transition, highlighting the impact on charge redistribution and transition behavior.
Findings
Charge self-consistency reduces charge redistribution in studied materials.
The metal-insulator transition in CaVO3 is only slightly affected by CSC.
Charge disproportionation in LuNiO3 is subtly influenced by CSC.
Abstract
We investigate the effect of charge self-consistency (CSC) in density functional theory plus dynamical mean-field theory (DFT+DMFT) 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, for which the importance of CSC is currently not well understood. Specifically, we analyze the strain-related orbital polarization in the correlated metal CaVO and the spontaneous electronic charge disproportionation in the rare-earth nickelate LuNiO. In both cases, we find that the CSC treatment reduces the charge redistribution compared to cheaper one-shot calculations. However, while the MIT in CaVO is only slightly shifted due to the reduced orbital polarization, the effect of…
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
On the effect of charge self-consistency in DFT+DMFT calculations for complex transition metal oxides
Alexander Hampel
Materials Theory, ETH Zürich, Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland
Center for Computational Quantum Physics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA
Sophie Beck
Materials Theory, ETH Zürich, Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland
Claude Ederer
Materials Theory, ETH Zürich, Wolfgang-Pauli-Strasse 27, 8093 Zürich, Switzerland
Abstract
We investigate the effect of charge self-consistency (CSC) in density functional theory plus dynamical mean-field theory (DFT+DMFT) 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, for which the importance of CSC is currently not well understood. Specifically, we analyze the strain-related orbital polarization in the correlated metal CaVO3 and the spontaneous electronic charge disproportionation in the rare-earth nickelate LuNiO3. In both cases, we find that the CSC treatment reduces the charge redistribution compared to cheaper one-shot calculations. However, while the MIT in CaVO3 is only slightly shifted due to the reduced orbital polarization, the effect of the site polarization on the MIT in LuNiO3 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 so-called “correlated materials”, i.e., materials where the strong Coulomb repulsion between electrons in partially filled or shells leads to effects that cannot easily be treated within effective non-interacting electron theories Held (2007). 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 towards a standard ab initio-based computational method for materials science Grieger et al. (2012); Adler et al. (2019), 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, Ishida and Liebsch (2009); Zhong et al. (2015); Lechermann (2018); Beck and Ederer (2019) defective systems, Backes et al. (2016); Sing et al. (2017); Souto-Casares et al. (2019) or large molecules Jacob et al. (2010); Turkowski et al. (2012). 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 SrVO3 monolayers demonstrate a reduced orbital polarization in CSC calculations compared to OS Bhandary et al. (2016); Schüler et al. (2018). Here, we therefore analyze the effect of charge self-consistency for two specific but representative cases. First, strained CaVO3, where orbital polarization is considered relevant for a reported strain-induced metal-insulator transition Gu et al. (2013); Beck et al. (2018). And second, LuNiO3, 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 Catalano et al. (2018).
Most previous work addressing the influence of charge self-consistency in DFT+DMFT calculations in transition metal (TM) oxides employed a so-called “”-model to define the correlated subspace, Wang et al. (2012); Dang et al. (2014); Park et al. (2014a); Aichhorn et al. (2009, 2011); Karolak et al. (2010) i.e., using a basis of rather localized, atomic-like orbitals constructed from a broad energy window that includes the TM as well as all oxygen 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 Aichhorn et al. (2011). 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 or 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 Lechermann et al. (2006). 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 states as “tails” located on the oxygen sites.
Another critical point arising in DFT+DMFT calculations using a -type orbital subspace, is that the physically very important charge transfer energy, , which describes the energy difference between oxygen and transition metal states, effectively becomes controlled by the so-called double-counting correction Wang et al. (2012); Karolak et al. (2010); Solovyev and Terakura (1998). 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 Kotliar et al. (2006). Different expressions to account for the double counting (DC) have been suggested Karolak et al. (2010); Haule (2015), but in some cases the double-counting needs to be adjusted manually, in order to obtain satisfactory results Dang et al. (2014); Wang et al. (2012). It was shown that CSC calculations for such -type calculations produce essentially the same spectral properties as OS calculations, if one tunes the DC correction to yield the same -state occupancy Wang et al. (2012). 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 Schüler et al. (2018). 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 CaVO3 and LuNiO3 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.
II Theoretical framework
II.1 DFT calculations
Structural relaxations for CaVO3 within the 20 atom unit cell in space group symmetry are performed using the Quantum ESPRESSO package Giannozzi et al. (2009). We employ scalar-relativistic ultrasoft pseudopotentials, with the and semicore states included in the valence for both V and Ca, together with the exchange-correlation functional according to Perdew, Burke, and Ernzerhof Perdew et al. (1996). Cell parameters and internal coordinates are relaxed until all force components are smaller than 0.1 mRy/ (: 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 wavefunctions and 840 Ry for the charge density. A Monkhorst-Pack -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 CaVO3, the in-plane lattice parameters are increased by 4% and kept fixed, while the -component of the cell and all atomic positions are relaxed.
All DFT calculations for LuNiO3 as well as the DFT parts of all our CSC DFT+DMFT calculations are performed using the projector augmented wave (PAW) method Blöchl (1994), implemented in the “Vienna Ab initio Simulation Package”(VASP) Kresse and Hafner (1993); Kresse and Furthmüller (1996); Kresse and Joubert (1999), and also using the exchange correlation functional according to Perdew, Burke, and Ernzerhof Perdew et al. (1996). For Ni, we use the PAW potential where the 3 semi-core states are included as valence electrons, while for Lu, we use the PAW potential corresponding to a valence state with -electrons frozen into the core. For the CaVO3 calculations with VASP, we use the PAW potentials where the and semi-core states are included as valence electrons for both Ca and V. Furthermore, a -point mesh with grid points along the three reciprocal lattice directions is used and a plane wave energy cut-off of 550 eV is chosen for LuNiO3 and 600 eV for CaVO3. The structure of LuNiO3 is fully relaxed within symmetry, both internal and lattice parameters, until the forces acting on all atoms are smaller than eV/Å.
II.2 DFT+DMFT calculations
II.2.1 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 , then a local Coulomb interaction is added, and the resulting Hubbard-like lattice Hamiltonian is solved via the DMFT approximation Georges et al. (1996); Held (2007). Without feedback to the DFT part, this corresponds to a OS calculation. To perform CSC calculations, one computes a correction to the charge density, , which is then passed back to the DFT code (here VASP) to calculate new KS wave-functions and hence, update the correlated subspace Amadon et al. (2008); Lechermann (2018). 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:
[TABLE]
where is the chemical potential and is the Kohn-Sham (KS) Hamiltonian. The lattice self-energy is obtained by solving the effective DMFT impurity problem (see next sub-section).
The lattice Green’s function in Eq. (1) is expressed in the KS (Bloch) basis. To achieve the up/down-folding between the quantities defined within the correlated subspace and the Green’s function in the KS basis,
[TABLE]
projector functions are introduced. The projector functions are defined as projections of the KS eigenstates onto localized orbitals
, . Here, serves as compound index for all local quantum numbers (site, orbital, and spin-character).
In our VASP-based OS and CSC calculations, the local basis functions are constructed from projection to localized orbitals (PLO) Amadon et al. (2008); Schüler et al. (2018). To construct optimal projector functions, we apply the scheme introduced in Ref. Schüler et al., 2018, 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 . We use initial projections on - or -like orbitals within the energy window of the correlated subspace . The resulting projectors are in general not orthogonal to each other, and need to be orthonormalized:
[TABLE]
[TABLE]
to define an orthonormal set of PLO-based Wannier functions describing the correlated subspace . 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 Aichhorn et al. (2016); Parcollet et al. (2015).
The projectors 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 LuNiO3 is smaller than %.
We note that the strong octahedral rotations present within symmetry lead to large off-diagonal crystal-field terms in the KS Hamiltonian, and the non-interacting 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 Gull et al. (2011) used to solve the effective impurity problem (see next sub-section), we perform a local unitary transformation of each impurity Green’s function after the down- respectively before the up-folding, which diagonalizes the initial non-interacting 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 CaVO3 we also perform OS DFT+DMFT calculations based on the electronic structure obtained with Quantum 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) Marzari et al. (2012) using the Wannier90 code Mostofi et al. (2014). 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 Mostofi et al. (2014).
The code used for all DFT+DMFT calculations in this paper is publicly available on github Hampel et al. (2019a).
II.2.2 Solving the impurity problem
For both CaVO3 ( subspace) and LuNiO3 ( subspace) the effective impurity problem within the DMFT cycle is solved with a continuous-time QMC hybridization-expansion solver Gull et al. (2011) implemented in TRIQS/cthyb Seth et al. (2016), 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 Vaugier et al. (2012),
[TABLE]
including all spin-flip and pair-hopping terms. Here, the operator creates an electron in the atom-centered Wannier orbitals of type and spin . The interaction parameters are given by the local intra-orbital Coulomb repulsion , and the Hund’s coupling . To reduce the QMC noise in the high-frequency regime of the impurity self-energy and , we represent both quantities in the Legendre basis Boehnke et al. (2011) and sample the Legendre coefficients directly within the TRIQS/cthyb solver.
II.2.3 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 Solovyev et al. (1994); Anisimov et al. (1997). Specifically, we use the parameterization given in Ref. Held, 2007 for the DC potential,
[TABLE]
where is the occupation of impurity site , and the average Coulomb interaction between orbitals, , is defined as Held (2007)
[TABLE]
The potential shift of Eq. (6) is added to the impurity self-energy; its form is directly tailored to the Hubbard-Kanamori interaction Hamiltonian in Eq. (5) for a - or -model resulting from an octahedral crystal-field environment of interacting orbitals ( for CaVO3 and for LuNiO3).
In this work, we draw particular attention on how the occupations 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 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 can have a strong influence on the results, and that one should be aware of this issue when evaluating the DC correction.
II.2.4 Calculation of observables
From the imaginary-time Green’s function, we calculate the spectral weight around the Fermi level, , which indicates whether the system is metallic () or insulating () Fuchs et al. (2011). For (), is identical to the spectral function at . For finite temperatures, it represents a weighted average around with a width of Fuchs et al. (2011).The full real-frequency spectral function is obtained via analytic continuation using the maximum entropy method Jarrell and Gubernatis (1996). 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 Lechermann et al. (2006):
[TABLE]
where are the KS eigenvalues with corresponding weights within the correlated subspace , and denotes quantities evaluated from the DMFT solution. The interaction energy is calculated using the Galitskii-Migdal formula Abrikosov et al. (2012); Galitskii and Migdal (1958), 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 . This way, we achieve an accuracy in the total energy of meV. All DMFT calculations are performed for eV*-1*, which corresponds to a temperature of 290 K.
III Materials & 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 CaVO3. While in the former case this material is a correlated metal Nekrasov et al. (2005); Pavarini et al. (2004), it has recently been demonstrated that tensile epitaxial strain leads to a transition towards the Mott insulating state within OS DFT+DMFT calculations Beck et al. (2018). An important aspect in this transition is the strain-induced crystal-field splitting between the partially filled 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 CaVO3, all correlated sites are symmetry-equivalent and thus the DC correction is irrelevant when using a minimal “-only” correlated subspace.
Second, we consider the rare earth nickelate LuNiO3, which exhibits a complex interplay between a specific structural distortion and an associated charge ordering, resulting in a transition from a paramagnetic metallic towards an also paramagnetic but insulating phase Catalano et al. (2018). In previous works, it was shown that DFT+DMFT is well suited to describe this correlated insulating state Park et al. (2012, 2014a); Subedi et al. (2015). Since two symmetry-inequivalent types of Ni sites appear in the insulating state, this case allows to analyze the effect of a site-dependent DC correction within CSC DFT+DMFT calculations.
Both materials, CaVO3 and LuNiO3, exhibit a distorted perovskite structure with space group (in the case of LuNiO3 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 GdFeO3-type distortion ( tilt system in Glazer notation), as depicted in Fig. 1. The -levels of the TM ions are split into and 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).
III.1 CaVO3 - orbital polarization
As stated above, bulk CaVO3 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 ultra-thin films Beck et al. (2018); Gu et al. (2013); McNally et al. (2019). As has been pointed out in Ref. Pavarini et al., 2005, the orbital polarization resulting from the orthorhombic distortion of the perovskite structure (related to the tilts and rotations of the octahedral network) is an 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 Gu et al. (2013); Sclauzero et al. (2016); Beck et al. (2018). For example, as depicted in Fig. 1, tensile epitaxial strain will lift the degeneracy of the -states, lowering the energy of one orbital compared to the other two. Since the orbital polarization in CaVO3 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 CaVO3 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 . In all cases, the spectral weight is finite for small values of , where the system is metallic, and then becomes zero in the insulating phase for large , with a rather sharp transition at . For the unstrained bulk system, all three approaches give identical results for the spectral weight as function of , with a critical value of =5.5 eV. Thus, at eV, which is typically considered as realistic value for transition metal oxides Pavarini et al. (2004), we find a finite weight corresponding to metallic behaviour, 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 (top left), it can be seen that the orbital polarization is weak in the metallic regime, but is significantly enhanced above , where the occcupation 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 Beck et al. (2018). 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 .
Under 4% tensile strain (right panels in Fig. 2), the MIT is shifted to lower values, below the realistic value of 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 eV. The CSC calculation, however, places the MIT at a slightly higher value of 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 Sclauzero et al. (2016); Beck et al. (2018) (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 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 Pavarini et al. (2004), 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 compared to the OS case.
To illustrate the difference between OS and CSC calculations in the strained case, we plot the spectral function at eV for both cases in Fig. 3. Here, the three different line-styles correspond to the three different -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 SrVO3 Bhandary et al. (2016).
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 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 , 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 method, 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 models due to larger hybridization effects Dang et al. (2014).
III.2 LuNiO3 — charge-ordering and structural energetics
The second case that we analyze is LuNiO3. This material belongs to the family of rare-earth nickelates, NiO3, where 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 to . The corresponding structural distortion results in a three dimensional checkerboard-like arrangement of long bond (LB) and short bond (SB) NiO6 octahedra, referred to as breathing mode distortion Medarde et al. (2008), and schematically shown on the right side of Fig. 1. Recent theoretical work indicates that this transition is related to an electronic instability towards 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) Peil et al. (2019); Hampel et al. (2019b). Furthermore, the choice of the site cation determines the degree of octahedral rotations in the corresponding high symmetry structure, and thus the bandwidth. The latter then controls how close the system is to the electronic instability, driving trends across the series Peil et al. (2019); Hampel et al. (2019b); Mercy et al. (2017); Varignon et al. (2017); Hampel and Ederer (2017).
Here, we use the case of LuNiO3 to analyze if, and how, the charge disproportionation, as a specific example for charge-ordering phenomena in general, is affected by the inclusion of charge self-consistency in DFT+DMFT. Earlier studies by Park et al. (2014b) also investigated the effect of CSC and DC for LuNiO3 using a -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 -like frontier orbitals per Ni site for our DFT+DMFT calculations. As shown in Ref. Subedi et al., 2015, the electronic instability towards charge disproportionation and the resulting site-selective Mott transition Park et al. (2012) occurring in the paramagnetic state is well described within DFT+DMFT using such a minimal subspace.
Subedi et al. (2015) found that the CDI state emerges in the frontier model for nickelates for rather large values of the Hund’s coupling , and is very sensitive to its variation. The fact that the Hund’s coupling is the critical ingredient in the occurence of the CDI state was first proposed by Mazin et al. (2007). They showed in an atomic picture that when becomes small and is overcome by the potential energy difference between the Ni sites, , 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 bands, whereas the Hund’s coupling is less sensitive to screening Mazin et al. (2007). A strong screening of in nickelates has been confirmed by recent studies using the constrained random phase approximation (cRPA) Seth et al. (2017a); Hampel et al. (2019b). Moreover, in Ref. Isidori et al., 2019 it is shown, that such a CDI regime for small or negative 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 LuNiO3 by employing a symmetry-based mode decomposition Perez-Mato et al. (2010), as outlined in Refs. Balachandran and Rondinelli, 2013; Hampel and Ederer, 2017; Hampel et al., 2019b. This allows to add the breathing mode distortion, with symmetry label , on top of the relaxed structure, and systematically vary its amplitude without changing any other parameter of the unit cell. We use the software ISODISTORT Campbell et al. (2006) to perform the mode decomposition.
III.2.1 Results for fixed structure
First, we calculate the properties of LuNiO3 for a fixed structure, using the experimentally observed breathing mode amplitude, Å Alonso et al. (2001), thereby varying the strength of the Hund’s coupling . As discussed above and shown in Ref. Subedi et al., 2015, the charge disproportionation and the resulting MIT depend sensitively on , which thus allows us to critically examine the influence of CSC on the most crucial system properties. We use a fixed value of 1.85 eV, which corresponds to the value calculated for LuNiO3 using the cRPA Hampel et al. (2019b); Seth et al. (2017b). The results are depicted in Fig. 4, where in the top panel the charge disproportionation, , i.e. the difference of the occupation between the LB and SB Ni sites, is shown as function of . The bottom panel shows the corresponding value for , indicating whether the system is metallic or insulating. The dashed vertical line corresponds to the value obtained within cRPA Hampel et al. (2019b); Seth et al. (2017b). 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 ” (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.2.3, this is the correct way to perform such CSC DFT+DMFT calculations, since the converged 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 eV, indicated by clear jumps in and . The jump in is related to a drastic change in the DC potential difference between the Ni sites, since, for not too large (see also below), the DC correction tends to increase the charge disproportionation by further lowering the states on the more occupied LB site compared to the less occupied SB site. This is discussed and analyzed in more detail in Appendix A.
For further increasing , stays almost constant until eV, where decreases again. Finally, at around eV, the system becomes metallic again. This can be explained by the fact that for increasing , the DC potential, proportional to (see Eq. (7) for ), decreases, and eventually changes sign for eV where . Thus, above eV the DC correction opposes the charge disproportionation by lowering the 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.2.3, it is ambiguous whether to use the DMFT impurity occupations or the occupations of the Wannier functions obtained within DFT, , to evaluate the DC correction. We first compare to the OS calculations where 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 eV. In addition, a small shift to larger can be observed compared to the CSC case. Thus, the tendency towards the CDI state is slightly stronger than in the CSC calculations.
In contrast, the OS calculations using (shown in orange) leads to a significantly reduced , which increases slowly with increasing . Moreover, for small eV, clear metallic behavior is observed, while from to 1.0 eV, the system undergoes the MIT, where eventually at eV the system is completely in the CDI state with . The occupations obtained in the initial DFT step are and .
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 , albeit with a small further reduction of .
The fixed structure calculations for LuNiO3, 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 towards 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 LuNiO3, reducing by approx. 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 towards 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 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 CaVO3 discussed in Sec. III.1. 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 towards the CDI state.
III.2.2 Influence on energetics
Another important aspect is the influence of charge self-consistency in total energy calculations for different amplitudes, i.e., for different amplitudes of the structural breathing mode distortion. As the 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 Kotliar et al. (2006), but nevertheless we discuss the difference here to better understand the influence of performing full CSC calculation Park et al. (2014b); Pourovskii et al. (2007). To analyze the resulting effects, we again use eV and two different values for , 0.42 eV (the cRPA value) and 1.1 eV (where 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, where the top panels show the total energy as function of the amplitude, and the bottom panels show the corresponding .
For the smaller value, eV, both the OS (green) and CSC (red) result in an energy minimum at a finite amplitude close to the experimental value (indicated by the vertical line). However, the OS calculation exhibits a much stronger response on the amplitude, and hence shows a significantly deeper energy minimum. In contrast, the “artificial” CSC calculation using for the DC correction (purple), exhibits no energy minimum for . Furthermore, the “correct” CSC calculation using undergoes a MIT to the CDI between and Å, while the corresponding OS calculation is already insulating without structural distortion and the CSC calculation with remains metallic for any calculated amplitude.
For eV, both CSC calculations, done either with DFT (purple) or DMFT occupations (red), agree very well (due to in the DC) and do not result in a stable finite breathing mode amplitude, even though both undergo a MIT at around Å 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 Å. Note that here we used for the DC correction, but the same result would be obtained using , due to .
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 .
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 CaVO3, 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 SrVO3 in Refs. Bhandary et al., 2016; Schüler et al., 2018. This has a small but noticeable effect on , the critical 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 CaVO3 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, LuNiO3, 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 %. Thus, similar to reducing the orbital polarization for strained CaVO3, including charge self-consistency 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 LuNiO3, OS calculations overestimate the response on the 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.
Acknowledgements.
This work was supported by ETH Zurich and the Swiss National Science Foundation through NCCR-MARVEL. Calculations have been performed on the cluster “Piz Daint” hosted by the Swiss National Supercomputing Centre.
Appendix A Influence of the DC on the effective inter-site splitting
In this appendix we explicitly show, how the DC corrections affects the level splitting between the two inequivalent Ni sites in the charge disproportionated state, which in turn controls the tendency to form a CDI state in the rare earth nickelates.
As outlined in the main text, Subedi et al. (2015) found that the CDI state emerges in the frontier model for nickelates when the following inequality is satisfied (derived from the the atomic limit):
[TABLE]
Here, is the “bare” site splitting, i.e., the difference in the average orbital energy between the SB and LB Ni sites, and is given as:
[TABLE]
where the first term, , denotes the corresponding splitting obtained within DFT from the on-site energies of the Wannier functions, and is found to be eV for Å in LuNiO3. The second term, , arises from the difference in the DC potential between the SB and LB sites:
[TABLE]
To further elucidate the interplay of the site-dependent DC potential in our DFT+DMFT calculations for LuNiO3, we analyzed the behavior of the critical quantity for obtaining a CDI state as function of for the different DC schemes applied in our study. The behavior of is depicted in Fig. 6 for the different flavors of DC, and for OS and CSC calculations at fixed eV and eV. The parameter regime which corresponds to the CDI in the atomic limit is highlighted in magenta. To directly compare our calculations with Ref. Subedi et al., 2015 we also performed OS calculations without applying any DC correction.
For the OS calculation without DC, (blue circles), becomes negative for eV. Then, for the calculations with DC evaluated with DFT occupations (OS: orange squares, CSC: purple crosses) this happens at a slightly smaller value, leading to an increased response in compared to the calculations without DC correction. This shows, that the DC correction enhances the CDI state for positive values of ( eV), since then is negative thus .
The strong tendency to form the CDI state in the calculations with 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 eV, and in the OS calculations (green stars) for even smaller . 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 pciture of the underlying physics.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1Held (2007) K. Held, Advances in Physics 56 , 829 (2007) . · doi ↗
- 2Grieger et al. (2012) D. Grieger, C. Piefke, O. E. Peil, and F. Lechermann, Physical Review B 86 , 155121 (2012) . · doi ↗
- 3Adler et al. (2019) R. Adler, C.-J. Kang, C.-H. Yee, and G. Kotliar, Reports on Progress in Physics 82 , 012504 (2019) . · doi ↗
- 4Ishida and Liebsch (2009) H. Ishida and A. Liebsch, Physical Review B 79 , 045130 (2009) . · doi ↗
- 5Zhong et al. (2015) Z. Zhong, M. Wallerberger, J. M. Tomczak, C. Taranto, N. Parragh, A. Toschi, G. Sangiovanni, and K. Held, Physical Review Letters 114 , 246401 (2015) . · doi ↗
- 6Lechermann (2018) F. Lechermann, in Handbook of Materials Modeling (Springer International Publishing, 2018) pp. 1–20. · doi ↗
- 7Beck and Ederer (2019) S. Beck and C. Ederer, Phys. Rev. Materials 3 , 095001 (2019) . · doi ↗
- 8Backes et al. (2016) S. Backes, T. C. Rödel, F. Fortuna, E. Frantzeskakis, P. Le Fèvre, F. Bertran, M. Kobayashi, R. Yukawa, T. Mitsuhashi, M. Kitamura, K. Horiba, H. Kumigashira, R. Saint-Martin, A. Fouchet, B. Berini, Y. Dumont, A. J. Kim, F. Lechermann, H. O. Jeschke, M. J. Rozenberg, R. Valentí, and A. F. Santander-Syro, Physical Review B 94 , 241110 (2016) . · doi ↗
