Dynamical Mean-Field Theory within the Full-Potential Methods: Electronic structure of Ce-115 materials

# Dynamical Mean-Field Theory within the Full-Potential Methods: Electronic structure of Ce-115 materials

Kristjan Haule, Chuck-Hou Yee, Kyoo Kim Department of Physics, Rutgers University, Piscataway, NJ 08854, USA
July 3, 2019July 3, 2019
July 3, 2019July 3, 2019
###### Abstract

We implemented the charge self-consistent combination of Density Functional Theory and Dynamical Mean Field Theory (DMFT) in two full-potential methods, the Augmented Plane Wave and the Linear Muffin-Tin Orbital methods. We categorize the commonly used projection methods in terms of the causality of the resulting DMFT equations and the amount of partial spectral weight retained. The detailed flow of the Dynamical Mean Field algorithm is described, including the computation of response functions such as transport coefficients. We discuss the implementation of the impurity solvers based on hybridization expansion and an analytic continuation method for self-energy. We also derive the formalism for the bold continuous time quantum Monte Carlo method. We test our method on a classic problem in strongly correlated physics, the isostructural transition in Ce metal. We apply our method to the class of heavy fermion materials CeIrIn, CeCoIn and CeRhIn and show that the Ce electrons are more localized in CeRhIn than in the other two, a result corroborated by experiment. We show that CeIrIn is the most itinerant and has a very anisotropic hybridization, pointing mostly towards the out-of-plane In atoms. In CeRhIn we stabilized the antiferromagnetic DMFT solution below K, in close agreement with the experimental Néel temperature.

###### pacs:
71.27.+a,71.30.+h

## I Introduction

One of the most active areas of condensed matter theory is the development of new algorithms to simulate and predict the behavior of materials exhibiting strong correlations. Recent developments in the dynamical mean-field theory (DMFT)Antoine (), a powerful many-body approach, hold great promise for more accurate and realistic descriptions of physical properties of this challenging class of materials.

The crucial step towards realistic description of strongly correlated materials was the formulation of DFT+DMFTAnisimov (); Lichtenstein (); PhysicsToday (), a method formed by the combination of density functional theory (DFT) and DMFT (for a review see Ref. our-rmp, ). To date, this method already has substantially advanced our understanding of the physics of the Mott transition in real materials and demonstrated its ability to explain phenomena including the structural phase diagrams of actinides V6 (); V7 (); PuAm (), phonon response V8 (), optical conductivity V9 (); FeAs (), valence and x-ray absorption V10 (); PuCris (); Xray () and transport V11 () of archetypal strongly correlated materials.

At present, much effort is devoted to the development of a robust and precise implementation of DFT+DMFT using state of the art DFT electronic structure codesSavrasov04 (); Lichtenstein_new (); KKR (); Wannier (); Wannier2 () and advanced impurity solversWerner (); Rubtsov (); CTQMC (); Werner2 (). This article describes in detail the implementation of this method within full-potential codes. There are three major issues that arise in DFT+DMFT implementations: i) quality of the basis set, ii) quality of the impurity solvers, and iii) choice of correlated orbitals onto which the full Green’s function is projected. Modern DFT implementations largely resolve the first issue, recent development of new impurity solvers Rubtsov (); Werner (); CTQMC (); SUNCA (); XDai (); SergejIPT (); LichSolver () have focused attention on the second, while the third is rarely discussed in the literature. Many DFT+DMFT proposals in the literature are based on downfolding to low energy model Hamiltonians Anisimov (); Anisimov2 (); Wannier (); Wannier2 (), which requires an atomic set of orbitals and treats the kinetic operator on the level of an effective tight binding model. In contrast, we avoid the ambiguities of downfolding and instead keep the kinetic part of the Hamiltonian and electronic charge expressed in a highly accurate full potential basis set. The advantage of our method is its ability to perform fully self-consistent electronic charge calculations. We concentrate here on the Linear Augmented Plane Wave basis (LAPW) LAPW-book () as implemented in the Wien2K code Wien2K () and the LMTO basis as implemented in LmtArt Savrasov96 (), in combination with the impurity solvers based on the hybridization expansion SUNCA (); Werner (); CTQMC (); Werner2 ().

The first half of the article introduces the basic steps of implementing the DFT+DMFT algorithm and provides a pedagogical introduction to the method. Section II is devoted to a crucial element of the DFT+DMFT formalism, namely the projection of the full electronic Green’s function to the correlated subset. We show that the projection used in the LDA+U method leads to non-causal DFT+DMFT equations, while the projection on to the solution of the Schrödinger equation within the Muffin-Tin (MT) spheres misses electronic spectral weight. We propose a new projection that leads to causal DMFT equations and captures all electronic spectral weight. Section III derives the DFT+DMFT equations from a Baym-Kadanoff-like functional formalism. Section IV provides a detailed flowchart of all the steps of the algorithm. In section V we discuss the necessary changes to the tetrahedron method when used in the context of DMFT. Section VI described the algorithm to compute transport properties within DFT+DMFT. Section VII describes the impurity solvers based on the hybridization expansion, the One Crossing Approximation (OCA), and the Bold continuous time quantum Monte Carlo algorithm (-CTQMC). Finally, section VIII discusses a new algorithm for analytic continuation of the self-energy from the imaginary to real axis.

In the second half of the article, we describe the results obtained by applying our new implementation of DFT+DMFT to several correlated materials. As a first test of the algorithm, in section IX we present its application to elemental cerium. Section X is devoted to a class of heavy fermion materials, CeRhIn, CeCoIn and CeIrIn, dubbed Ce-155 materials. We show the difference in the electronic structure among these three materials and demonstrate that the Ce electrons are most localized in CeRhIn and order antiferromagnetically below K, in agreement with experiment, while the Ce electrons are most itinerant in CeIrIn5. We explain the origin of the subtle difference between the three Ce-115 compounds from the electronic structure point of view.

## Ii Projection on to correlated orbitals within full-potential methods

DFT+DMFT contains some aspects of band theory, adding a “frequency-dependent local potential” to the Kohn-Sham Hamiltonian. It also contains some aspects of quantum chemistry, carrying out an exact local configuration interaction procedure by summing all local diagrams, which requires the definition of an “atomic-like” or “local” Green’s function. The operation of extracting the local Green’s function from the full Green’s function is called projection (or truncation). The reverse operation of expressing the local time-dependent potential , derived from the solution of the atomic problem in the presence of a mean-field environment, is called embedding. The various DFT+DMFT implementations differ not only in the choice of basis set, but also in the choice of the projection-embedding step. These ingredients are sketched schematically in Fig. 1. The projection-embedding step connects the atomic and solid state physics, and its proper definition is a conceptual issue of DFT+DMFT method.

In the current formulation of DFT+DMFTour-rmp (); Held (); Lich (), one must define the correlated orbitals to which the Coulomb correlation is applied, i.e., , where is a localized orbital. Usually, this is achieved by transforming the DFT Hamiltonian to a set of localized Wannier orbitals. These Wannier orbitals are then identified as the local correlated orbitals of DMFT. Various choices of these orbitals were proposed in the literature, including tight-binding LMTO’s Anisimov (); Lichtenstein (), non-orthogonal LMTO’s Savrasov04 (), Nth-order Muffin-Tin orbitals Pavarini (), numerically-orthogonalized LMTO’s Purovskii (), and maximally-localized Wannier orbitals Wannier2 (); Korotin (). The basis functions must fully respect the symmetries of the problem and be atom-centered, rather than bond-centered. Hence maximally-localized Wannier functions MLWO () are not a good starting point for DMFT.

Localized basis sets are a better starting point for our purposes, but the non-orthogonality of these sets pose a serious challenge. Straighforward orthogonalization mixes the character of the orbitals, resulting in mixed the partial occupancies and partial density of states, leading to incorrect partial electron counts. For example, within modern DFT implementations, cerium metal has approximately one electron. Naïve orthogonalization results in a considerably higher electron count, leading to an unphysical DMFT solution.

Even more challenging is the formulation of the good localized orbitals in full-potential basis sets. Here, multiple basis functions are used to obtain more variational freedom. To implement DMFT in such basis sets, the group of orbitals representing the correlated electrons in the solid must be contracted to form a single set of atomic-like heavy orbitals, i.e., one orbital per Ce atom, one orbital per Fe atom, etc.

A straighforward projection on to the orbital angular momentum eigenfunctions leads to non-causal DMFT equations, which result in an unphysical auxiliary impurity problem. The second often-employed choice is the projection on to the solution of the Schrödinger (Dirac) equation inside the MT sphere . While this choice is certainly superior to the straighforward projection, it does not take into account the contributions due to the energy derivative of the radial wave function and the localized orbitals (LO) at other energies , and hence misses some electronic spectral weight of the correlated orbital. Alternative choices are possible which simultaneously capture all spectral weight and obey causality. We implemented one of them and we believe it is superior to other choices in the literature.

The central objects of DMFT are the local Green’s function and the local self-energy of the orbitals within the correlated subset. We specify the projection scheme by the projection operator , which defines the mapping between real-space objects and their orbital counterparts (see Fig. 1). The operator acts on the full Green’s function and gives the correlated Green’s function  Savrasov04 (); our-rmp ()

 GτLL′=∫drdr′P(rr′,τLL′)G(rr′). (1)

The integrals over and are performed inside the sphere of size around the correlated atom at position . The subscript can index spherical harmonics , cubic harmonics, or relativistic harmonics , depending on the system symmetry. We always choose the basis which minimizes the off-diagonal elements of the correlated Green’s function in order to reduce the minus-sign problem in Monte-Carlo impurity solvers. In general, is a multidimensional tensor with one pair of indices in the space of local correlated orbitals and the other pair in the space of the full basis set, which can be expressed in a real space or Kohn-Sham basis, where and are band indices.

The inverse process of embedding , i.e. the mapping between the correlated orbitals and real-space , is defined by the same four-index tensor. However, instead of integrals over real-space, its application is through a discrete sum over the local degrees of freedom,

 Σ(r,r′)=∑τLL′∈HP(r′r,τL′L)ΣτLL′ (2)

Here means to only sum over correlated orbitals. In actinides, the sum would run over orbitals, in lanthanides over and in transition metals over orbitals. runs over all atoms in the solid and over the full space. Note that within the correlated Hilbert subspace, the embedding and projection should give unity , i.e.,

 ∫drdr′P(rr′,τL1L2)P(r′r,τ′L3L4)=δL1L4δL2L3δττ′, (3)

while the projection from the full Hilbert space to the correlated set, followed by embedding, gives the correlated local Green’s function in real space

 G(r,r′)= (4) ∑τLL′∈HP(r′r,τL′L)∫dr1dr2P(r1r2,τLL′)G(r1r2)

which is the central object of the functional definition of the DMFT described below. In general, the two operators and could be different, but they must satisfy the condition Eq. (3).

The two simplest projections, namely, the projection on to the orbital angular momentum functions , and the projection on to the solution of the Schrödinger equation, can be explicitly written as

 P0(rr′,τLL′) = YL(^rτ)δ(r−r′)Y∗L′(^r′τ) (5) P1(rr′,τLL′) = YL(^rτ)u0l(rτ)u0l′(r′τ)Y∗L′(^r′τ) (6)

where is the vector defined with the origin placed at the atomic position , and is the solution of the radial Schrödinger equation for angular momentum at a fixed energy .

In the following, we will show that the projection , used in some implementations of DMFT Lichtenstein_new (), captures the full spectral weight of the correlated character , but leads to non-causal DMFT equations. On the other hand gives causal DMFT equations, but misses some spectral weight.

In our view, a good DFT+DMFT implementation should satisfy the following conditions

• Correct correlated spectral weight: The projected density of states, computed from the projected Green’s function,

 ρL(ω)=12πi[G†LL(ω)−GLL(ω)], (7)

should capture the partial electronic weight inside a given MT sphere at all frequencies, i.e., . In particular, must include the electronic weight contained in and local orbitals. Projection should not include any weigh of other character, nor miss correlated weight.

• DMFT equations are causal: For any causal self-energy , the DMFT self-consistency condition

 1ω−Eimp−Σ−Δ=∑kPk[(ω+μ−HDFTk−EkΣ)−1] (8)

should give a causal hybridization function . Here we used projections in momentum space as opposed to their real-space definitions in Eqs. (1) and (5), (6).

• Sufficient accuracy of the hybridization function: The hybridization function is usually very sensitive to the choice of the projector. Therefore, we require that in the relevant low energy region, the hybridization function is similar to its DFT counterpart. Explicitly, must be sufficiently close to its DFT estimate, . Here stands for the full Green’s function when . The choice of is dictated by the fact that the hybridization , computed by is not well behaved for , as we will show below. The motivation for using P0 in the above equation is that we want to project the full Hilbert space to a correlated subset with pure angular momentum, either f or d, but not to a mixure of characters.

• Good representation of kinetic energy and electronic density: Finally, it is crucial to faithfully represent the kinetic energy operator and electronic density in real space, a feat most modern DFT implementations achieve. The DFT+DMFT implementation should not reduce the precision already achieved in DFT underlying code.

Downfolding to only a few low energy bands clearly violates the condition number (3), since the hybridization outside the downfolded window vanishes. A more severe problem is that downfolding approximates the kinetic energy operator by expressing it in a small atomic-like basis set, hence condition (4) is violated. Therefore, we will focus our discussion on DFT+DMFT implemented within full-potential basis sets where all bands are kept at each stage of the calculation. Downfolding to a sufficiently large energy window may sometimes be helpful due to its conceptual simplicity, but this approach can not compute the electronic charge self-consistently, as is possible in our implementation. Moreover, the localized orbitals chosen in the downfolding procedure combined with the limited number of hoppings retained often cannot faithfully represent the original Kohn-Sham bands.

To be more concrete, we will give the proofs of the “weight loss problem” and “causality problem” within the full-potential LAPW basis. The equivalent derivation is possible for the full-potential LMTO basis. Inside the MT spheres, the full-potential LAPW basis functions can be written LAPW-book ()

 χk+K(r)=∑LτκAτκk+K,Luτκl(rτ)YL(^rτ) (9)

where corresponds to the solution of the Schrödinger equation at a fixed energy , to the energy derivative of the same solution , and to a localized orbitals at additional linearization energies . Here runs over the atoms in the unit cell.

The Kohn-Sham states are superpositions of the basis functions

 ψik(r)=∑KCkiKχk+K(r) (10)

and take the following form inside the MT spheres:

 ψik(r)=∑τLκAτκiL(k)uτκl(rτ)YL(^% rτ) (11)

where , or equivalently, .

The projectors (5) and (6) can be expressed in the Kohn-Sham basis:

 Pk(ij,τLL′)=∫drdr′ψ∗ik(r)P(rr′,τLL′)ψjk(r′). (12)

Hence, projector takes the form

 P0k(ij,τLL′) =∫drdr′ψ∗ik(r)YL(^rτ)δ(r−r′)Y∗L′(^r′τ)ψjk(r′) (13)

Using projector , we get the following expression for the partial density of states

 DτL(ω)=∑κκ′kiAτκ∗iL(k)Aτκ′iL(k)⟨uτκl|uτκ′l⟩δ(ω+μ−εki) (14)

which exactly coincides with the DFT partial DOS. Hence satisfies the condition number (1). However, it does not lead to causal DMFT equations.

To show that, consider the limit of a diverging self-energy, , as is relevant for the Mott insulators. Despite the diverging , the projection must still produce a finite hybridization. In the case when all the bands at the energy of the pole are correlated, the hybridization should vanish. In this limit, the DMFT self-consistency condition (8) takes the form

 (Στ+Δ)−1LL′ =∑kijPk(ji,τLL′)⎡⎣∑L2L3τ′Pk(::,τ′L2L3)Στ′L3L2⎤⎦−1ij (15)

where stands for the two band indices constituting a matrix in to be inverted. Since is finite while diverges, we neglect to obtain the condition for causal projection,

 δLL′′=∑kij,τL′Pk(ji,τLL′)ΣτL′L′′× ×⎡⎣∑L2L3τ′Pk(::,τ′L2L3)Στ′L3L2⎤⎦−1ij. (16)

This equation must be satisfied for any matrix form of the self-energy . Moreover, it has to be satisfied for each and . We will show below that Eq. (16) is satisfied for a separable projection (see Eq. 19 for a definition), while for a non-separable projection, it likely is not. One can check explicitely that violates the condition Eq. (16). Only after applying an additional trace over will the two matrices cancel. However, for any given choice of , does not satisfy the causality condition. Instead a pole in the self-energy results in a diverging , with the imaginary part having the wrong sign. The projection is implemented in the qtl package NovakQTL () of Wien2KWien2K (). The LDA+U implementation within Wien2K Shick () also uses , but this does not cause any causality issues since the problem is unique to DFT+DMFT. Additionally, simple impurity solvers such as Hubbard-I (Ref. Lichtenstein_new, ) do not incorporate a true hybridization so they also avoid issues with causality.

Finally, let us mention an attractive feature of . Within this scheme, the self-energy is independent of the radial distance from the atom , having only angular dependence in the form . This matches the conceptual fact that the impurity solver within the DMFT framework can not determine the radial dependence of the self-energy. The impurity solver can only be used to obtain the angular dependence of by determining the expansion coefficeints . In the absence of any knowledge of the radial dependence of , the natural choice is a constant function, independent of radius . Since is a function of two vectors, a radial delta function would be an obvious choice. However, issues with causality preclude the use of this projection.

The second projection of Eq. (6) takes the following form in the Kohn-Sham basis:

 P1k(ij,τLL′)=∑κκ′AτκiL(k)Aτκ′∗jL′(k)⟨uτκl|u0l⟩⟨u0l′|uτκ′l′⟩. (17)

The partial density of states computed from the correlated Green’s function using is

 DτL(ω)= (18) ∑κκ′kiAτκiL(k)Aτκ′∗iL(k)⟨uτκl|u0l⟩⟨u0l|uτκ′l⟩δ(ω+μ−εki)

Comparing Eq. (18) with (14), we notice that is replaced by , which leads to incorrect spectral weight. In particular, for , the original overlap in Eq. (14) is , while the overlap obtained by , vanishes.

Causality is not violated for any projection , which is separable, i.e., can be cast into the form

 Pk(ij,τLL′)=UkτiLUkτ∗jL′. (19)

The condition Eq. (16) can then be expressed as

 1=∑kUkτ†(UkτΣτUkτ†)−1UkτΣτ (20)

which is clearly satisfied when is invertible matrix because . This is satisfied when the Kohn-Sham Hilbert space is of larger dimension than the correlated Hilbert space. The projection leads to causal DMFT equations, and therefore is a better choice than . However, some spectral weight is lost at energies away from the linearization energy . To this end, we also implemented an alternative projection within Wien2K package Wien2K (), which preserves both causality and spectral weight. This projector is given by

 P2(rr′,τLL′)=∑ijkκκ′ψik(r)AτκiL(k)⟨uτκl|uτ0l⟩⟨uτ0l′|uτκ′l′⟩Aτκ′∗jL′(k)ψ∗jk(r′)×   ⎷(∑κ1κ2Aτκ1iLAτκ2∗iL⟨uτκ1l|uτκ2l⟩∑κ1κ2Aτκ1iLAτκ2∗iL⟨uτκ1l|uτ0l⟩⟨uτ0l|uτκ2l⟩)⎛⎝∑κ1κ2Aτκ1∗jL′Aτκ2jL′⟨uτκ1l′|uτκ2l′⟩∑κ1κ2Aτκ1∗jL′Aτκ2jL′⟨uτκ1l′|uτ0l′⟩⟨uτ0l′|uτκ2l′⟩⎞⎠. (21)

Here index runs over the local basis in which the green’s function is minimally off-diagonal (cubic harmonics or relativistic harmonics).

The projector is separable, as postulated in Eq. (19), and the transformation is

 UkτiL=∑κAτκiL(k)⟨uτκl|uτ0l⟩SτiL (22)

with

 SτiL= ⎷∑κ1κ2Aτκ1iLAτκ2∗iL⟨uτκ1l|uτκ2l⟩∑κ1κ2Aτκ1iLAτκ2∗iL⟨uτκ1l|uτ0l⟩⟨uτ0l|uτκ2l⟩ (23)

Hence the DMFT equations are causal. Moreover, is identical to and hence the partial density of states , obtained by , is identical to Eq. (14). Hence the projection correctly captures the partial spectral weight. Knowledgeable reader would notice that the projection is slightly non-local because is weakly momentum dependent. At energies where or local orbital substantially contribute to the spectral weight (away from the Fermi level), we give up locality in expense of correctly capturing the spectral weight.

All projection schemes lead to slightly non-orthonormal correlated Green’s function. This is because the interstitial weight is not taken into account and because the full potential basis is overcomplete. To have an orthonormal impurity problem, we compute the overlap and renormalize .

Finally, we remark that the segment of our code which builds projections , and within Wien2K Wien2K () is based on the qtl package of Pavel Novak NovakQTL ().

Similar projections within LDA+DMFT method were proposed before. In particular the method by B. Amadon et.al. Wannier2 () proposed to construct the Wannier functions for the correlated subset only, while the DMFT equations were solved in the Kohn-Sham basis, restricted to some subset of low energy bands. The local orbitals used for the projection were either all-electron atomic partial waves in the PAW framework, or pseudo-atomic wave functions in mixed-basis pseudopotential code. Hence, in the language of projectors, the method was similar to choosing the projector to be , where is the the partial waves or pseudo-atomic wave function. While this method is clearly causal, it looses spectral weight of the correlated angular momentum character. Moreover, the implementation of the method did not allow the self-consistent evaluation of the electronic charge. The method of Anisimov et.al. Anisimov2 () also proposed a construction of the Wannier functions using an arbitrary set of localized orbitals. In their work, the LDA Hamiltonian was truncated to Wannier representation for the purpose of obtaining the DMFT self-energy. This simplifies the self-consistent DMFT problem, but makes it impossible to implement the charge self-consistency. Finally, Savrasov et.al. Savrasov04 () proposed a projector particular to LMTO basis set, for which causality was not proven.

## Iii DFT+DMFT Formalism

To derive the DFT+DMFT equations, we define a functional of the correlated Green’s function and extremise it. The correlated Green’s function is defined by Eq. (4), and the functional to be extremise is

 Γ[G,ρ]=−Trln(G−1)−Tr[ΣtotG]+Φ[G,ρ], (24)

where runs over all space (orbitals,momenta) and time (frequency). The quantities apprearing in the above functional are

 (25) Σtotω(r,r′)=[VH(r)+Vxc(r)]δ(r−r′)+[Σω(r,r′)−EDCδ(r−r′)]Θ(r

where is trace over time only (not space), is the potentials due to ions, are the Hartree, and exchange-correlation potential, respectively. is the sum of all local two particle irreducible skeleton diagrams constructed from , and the Coulomb repulsion (screened by orbitals not contained in ), and is the double counting functional.

We assume that the Coulomb interaction has the same form as in the atom, i.e.,

 ^U=∑La,..Ld,,m,σσ′2l∑k=04πFk{l}2k+1⟨YLa|Ykm|YLc⟩⟨YLb|Y∗km|YLd⟩f†Laσf†Lbσ′fLdσ′fLcσ (28)

however, the Slater integrals are reduced due to screening effects. Typically, we renormalize by 30%, from their atomic values, while , being renormalized more, can be estimate by constraint LDA or constraint RPA Ferdi ().

To extremize the functional Eq. (24), we take and as independent variables, and use the following functional dependence: , , , are functionals of . Consequently, is also a functional of , i.e., . On the other hand, , , , are functionals of the total electron density, hence is also a functional of since . Finally it is easy to check that

 Tr[ΣtotG]=Tr[(VH+Vxc)ρ]+Tr[(Σ−EDC)G].

With the above functional dependence in mind, minimization with respect to gives

 Σ−EDC=δΦDMFT[G]δG−δΦDC[G]δG,

and minimization with respect to leads to

 VH+Vxc=δΦH[ρ]δρ+δΦxc[ρ]δρ.

Hence the Hartree and exchange-correlation potential are computed in the same way as in DFT method (note however is electron density in the presence of DMFT self-energy), while the DMFT self-energy is the sum of all local Feynman diagrams, constructed from and Coulomb interaction .

To sum up all local diagrams, constructed from and screened Coulomb interaction , we solve an auxiliary quantum impurity problem, which has as the impurity green’s function, and as the impurity self-energy . The impurity Green’s function is , hence the DMFT self-consistency condition reads

 ^P(iω+μ−HDFT−^E¯¯¯¯Σ))−1=(iω−Eimp−Σimp−Δ)−1. (29)

where , and is the interaction included in DFT (double counting). The self-consistency condition takes the explicit form

 ∫(r,r′)

where and is the muffin-tin radius.

For efficient evaluation of the DMFT self-consistency condition Eq. (30), we choose to work in the Kohn-Sham (KS) basis. At each DFT+DMFT iteration, we first solve the KS-eigenvalue problem

 [−∇2+VKS(r)]ψki(r)=ϵkiψki. (31)

Then we express the projection in KS basis, , where run over all bands. We then perform the embedding of the self-energy, i.e., transforming it from DMFT base to the KS base

 ¯¯¯¯Σk,ij(ω)=∑τ,L1L2Pkτ(ji,τL2L1)¯¯¯¯ΣτL1L2(ω) (32)

In KS-base, we can invert the Green’s function Eq. (30), to obtain the practical form of the self-consistency condition

 GτLL′ = ∑kijPkτ(ij,LL′)[(iω+μ−ϵk−¯¯¯¯Σk(ω))−1]ji (33) GτLL′ = [1iω−Eτimp−Στ(ω)−Δτ(ω)]LL′ (34)

This is of course equivalent to Eq. (30). Finally we solve this self-consistency equation for a given self-energy to obtain the hybridization function and the impurity levels .

We note in passing that the self-energy is a complex function, and its imaginary part is related to the electron-electron scattering rate, which is very large in correlated materials. In Mott insulators, it is even diverging. Hence the DMFT ”effective Hamiltonian” can not be diagonalized by standard methods to obatin a set of eigenvalues, i.e., bands. The eigenvalues are complex and hence only the spectral weight is a well defined quantity. The absence of well defined bands in correlated materials makes computational techniques more challenging. For example, the calculation of the chemical potential is far more demanding because one can not assign a unity of charge to each fully occupied band. Rather all complex eigenvalues, even those which are far from the Fermi level, need to be carefully considered. This point will be addressed below in section IV, item 5. Further, the tetrahedron method Tetra (), a very useful technique to reduce the number of necessary momentum points in practical calculation, is not applicable since it needs real eigenvalues. We address the necessary generalization of this method is chapter V.

Note that generalization of the projector and the LDA+DMFT formalism to cluster-DMFT is very straightforward. One needs to increase the unit cell to include more sites of the same atom type. The self-energy and the Green’s function become matrices in index , i.e., , . The transformation is also straightforwardly generalized to matrix form . The only difference in the definition of the projector Eq. (21) is that is replaced by ( remains unchanged), which amounts to the integral over two different spheres around two atoms of the same type. Finally, in cluster-DMFT case, the self-energy in KS-basis Eq. (32) has to be summed over both and , and self-consistency condition Eq. (34) becomes a matrix equation in . The challenging part of the cluster-DMFT formalism is in solving the cluster-impurity problem. In combination with impurity solvers based on the hybridization expansion (discussed below) the computational effort grows exponentially with the number of correlated sites. In the weak coupling impurity solvers, the computational effort grows as a power-law, however, these techniques usually can not reach the interesting regime of strong correlations and low temperatures.

The major bottleneck in evaluating the DMFT self-consistency condition in our method is the multiplication of the projector with in Eq. (32) and multiplication of projection with Green’s function in Eq. (33). Since projection is separable, one can write the operation in terms of matrix products. Still, these sums run over all -points (typically few thousands) and all frequency points (typically few hundreds).

For the efficient implementation of the set of Eqs. (32) and (33), we first notice that the transformation (or its separable part ) is very large and is not desirable to be written to the computer hard disc. Hence we generate it only for one -point at a time, and evaluate both products at this particular -point. Non-negligible amount of time is necessary to generate the transformation Eq. (21), and because this transformation does not depend on frequency, it needs to be used for all frequencies in Eqs. (32) and (33). Hence paralization over frequency is not implemented, while paralization over -points is.

Note that because of the sum over atoms () in Eq. (32), the transformation for all atoms needs to be computed first, and only then the sum in Eq. (32) can be evaluated and the self-consistency condition Eq. (34) can be inverted.

To optimize the sum in Eqs. (32) and (33), one can notice that local quantities like self-energy and local green’s function possess a large degree of symmetry when written in proper basis (real harmonics, relativistic harmonics): many off-diagonal matrix elements vanish, and many matrix elements are equivalent. For example, in a system with cubic symmetry, one has only two types of self-energy and . Hence, instead of summing over matrix elements in Eq. (32), one can rewrite the sum over two matrix elements , i.e.,

 Σk,ij(ω)=∑τ,tΣ(τ)t(ω)Pkτ(ji,t) (35)

where and the indices here stand for the real harmonics rather than spheric harmonics. The later transformation is independent of frequency, while the sum Eq. (35) needs to be performed for all frequencies, hence the compact form of the transformation saves a lot of computer time.

## Iv The algorithm

The implementation of the DFT+DMFT algorithm is done in the following few steps:

• : We converge the LDA/GGA equations to get the starting electronic charge . We use the non-spin polarized solution as starting point. In the ordered state, the DMFT self-energy is allowed to break the symmetry, while typically the exchange-correlation potential is not allowed to break the symmetry (LDA rather than LSDA).

In this preparation step we also obtain good estimates for the Coulomb repulsion (which is represented by Slater integrals , , and ). Slater integrals are computed by the atomic physics program of Ref. Cowan, , and they are scaled down by 30% to account for the screening in the solid. The terms is very different from the atomic and is obtained by constraint LDA calculation, or constraint RPA calculation Ferdi ().

• : We solve the DFT KS-eigenvalue problem

 (−∇2+VKS(r))ψki(r)=ψki(r)εDFTik

to obtaine KS eigenvectors, core, and semicore charge, and linearization energies .

• We start with a guess for the lattice self-energy correction (here is the dynamic part of the self-energy with the property ). A reasonable starting point is and . The potential in the first DMFT iteration is thus the DFT potential.

• : Next we embed the DMFT self-energy (shifted by double counting) to Kohn-Sham base by the transformation Eq. (32) to obtain .

• : Using the current DMFT self-energy , and the current DFT KS-potential , we compute the current chemical potential. This is done in the followin steps:

• Complex eigenvalues of the full Green’s function are found in the large enough energy interval (at least [,]) by solving

 ∑j[εDFTkiδij+¯¯¯¯Σkij(ω)]Ckjl(ω)=Ckil(ω)εkl(ω).

Here are DMFT eigenvectors expressed in KS base. The DMFT eigenvalues outside this interval are set to DFT eigenvalues. We need only eigenvalues in this step, but not eigenvectors.

• The chemical potential is determined using precomputed complex and frequency dependent eigenvalues . On imaginary axis we solve

 Nval=T∑kl,ωn1iωn+μ−εkl(iωn)

and on real axis we solve

 Nval=−1πIm∑kl∫f(ω)dωω+μ−εkl(ω)

If enough -points can be afforded, we use special point method, otherwise the “complex tetrahedron method” can be used (see chapter V).

For numerical evaluation of the real axis density, we discretize the integral

 Nval=−1πIm∑if(ωi)∑kl∫biaidωω+μ−εkl(ωi)

with and . When using the special point method, the integral over frequency is evaluated analytically, and the terms of the form are summed up. Alteratively, we sometimes use the complex tetrahedron method, where the four-dimensional integral is evaluated analytically (see chapter V)

When DMFT is done on imaginary axis (using imaginary time impurity solvers), we evaluate

 N=∑klf(ε0kl−μ)+ 2T∑0<ωn<ωN∑kl[1iωn+μ−εkl(iωn)−1iωn+μ−ε0kl] −1πarctan(ε∞kl−μωN)+1πarctan(ε0kl−μωN) (36)

Here is the real part of the eigenvalue at arbitrary frequency. We choose the lowest or the last Mastubara point. Again, the tetrahedron method can be used for momentum sum.

For Mott insulators, the above described method is not very efficient, because even a small numerical error in computing places chemical potential at the edge of the Hubbard band, either upper or lower. This instability usually does not allow one to reach a stable self-consistent solution. We devised the following method to remove this instability:

• The diagonal components of the self-energy were fitted by a pole-like expression .

• Next, we neglected broadening of the pole (), which should be small in the Mott insulating state. We computed a quasiparticle approximation for the Green’s function , i.e.,

 (Gqpk)−1ij=iω−εDFTki−¯¯¯¯Σ∞,ij−UkτiL√WL1iω−PL√WLUkτ∗jL (37)

where is part of the projector defined above.

• The above Green’s function formulae can be cast into a block form

 (38)

Here is the quasiparticle Hamiltonian which can be diagonalized to obtain the quasiparticle bands. We notice that the number of quasiparticle bands of the Mott insulator is larger then the number of Kohn-Sham bands because Mott insulators have at least two Hubbard bands. The quasiparticle bands are not very accurate away from the Fermi level, however they are sufficiently acurate at low energy and allow one to identify gaps at the Fermi level. Once a gap in the spectra of is identified, the charge is computed using the full DMFT density matrix to verify the neutrality of the solid. If the solid is neutral when chemical potential is in the gap, the chemical potential is set to the middle of the gap.

• : Impurity hybridization function and impurity levels are computed in this step.

We use equation (33) to get and we use the high frequency expansion of both equations (33) and (34) to determin impurity levels

 EimpLL′=−EDCδLL′+∑kiPkτ(ii,LL′)εDFTki
• : Impurity solver uses , , and Coulomb repulsion (which is represented by Slater integrals , , and ) as the input and gives the new self-energy as the output.

Currently we integrated the following impurity solvers: OCA (see chapter VII.2), Non-crossing approximation (NCA), Continuous time quantum Monte Carlo (CTQMC) CTQMC (). The latter is implemented on imaginary axis, and the former two on real axis.

Before the impurity solver is run, we exactly diagonalize the atomic problem in the presence of crystal fields, to obtain all atomic energies and the matrix elements of electron creation operator in the atomic basis . Since the impurity levels can change during the iteration, the crystal field of the atomic problem can change as well. In case of -systems, the crystal field splittings are small and one can assume that they do not change substantially from their DFT value. Hence the exact diagonalization can be done only once at the beginning. For the