Dynamical MeanField Theory within the FullPotential Methods: Electronic structure of Ce115 materials
Abstract
We implemented the charge selfconsistent combination of Density Functional Theory and Dynamical Mean Field Theory (DMFT) in two fullpotential methods, the Augmented Plane Wave and the Linear MuffinTin 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 selfenergy. 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 outofplane 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.+hI 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 meanfield theory (DMFT)Antoine (), a powerful manybody 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. ourrmp, ). 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 xray 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 fullpotential 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 selfconsistent electronic charge calculations. We concentrate here on the Linear Augmented Plane Wave basis (LAPW) LAPWbook () 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 noncausal DFT+DMFT equations, while the projection on to the solution of the Schrödinger equation within the MuffinTin (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 BaymKadanofflike 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 selfenergy 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 Ce155 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 Ce115 compounds from the electronic structure point of view.
Ii Projection on to correlated orbitals within fullpotential methods
DFT+DMFT contains some aspects of band theory, adding a “frequencydependent local potential” to the KohnSham 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 “atomiclike” 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 timedependent potential , derived from the solution of the atomic problem in the presence of a meanfield 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 projectionembedding step. These ingredients are sketched schematically in Fig. 1. The projectionembedding 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+DMFTourrmp (); 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 tightbinding LMTO’s Anisimov (); Lichtenstein (), nonorthogonal LMTO’s Savrasov04 (), Nthorder MuffinTin orbitals Pavarini (), numericallyorthogonalized LMTO’s Purovskii (), and maximallylocalized Wannier orbitals Wannier2 (); Korotin (). The basis functions must fully respect the symmetries of the problem and be atomcentered, rather than bondcentered. Hence maximallylocalized Wannier functions MLWO () are not a good starting point for DMFT.
Localized basis sets are a better starting point for our purposes, but the nonorthogonality 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 fullpotential 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 atomiclike 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 noncausal DMFT equations, which result in an unphysical auxiliary impurity problem. The second oftenemployed 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 selfenergy of the orbitals within the correlated subset. We specify the projection scheme by the projection operator , which defines the mapping between realspace 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 (); ourrmp ()
(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 offdiagonal elements of the correlated Green’s function in order to reduce the minussign problem in MonteCarlo 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 KohnSham basis, where and are band indices.
The inverse process of embedding , i.e. the mapping between the correlated orbitals and realspace , is defined by the same fourindex tensor. However, instead of integrals over realspace, its application is through a discrete sum over the local degrees of freedom,
(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.,
(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
(4)  
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
(5)  
(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 noncausal 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,
(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.

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 atomiclike basis set, hence condition (4) is violated. Therefore, we will focus our discussion on DFT+DMFT implemented within fullpotential 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 selfconsistently, 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 KohnSham bands.
To be more concrete, we will give the proofs of the “weight loss problem” and “causality problem” within the fullpotential LAPW basis. The equivalent derivation is possible for the fullpotential LMTO basis. Inside the MT spheres, the fullpotential LAPW basis functions can be written LAPWbook ()
(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 KohnSham states are superpositions of the basis functions
(10) 
and take the following form inside the MT spheres:
(11) 
where , or equivalently, .
The projectors (5) and (6) can be expressed in the KohnSham basis:
(12) 
Hence, projector takes the form
(13) 
Using projector , we get the following expression for the partial density of states
(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 selfenergy, , 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 selfconsistency condition (8) takes the form
(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,
(16) 
This equation must be satisfied for any matrix form of the selfenergy . 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 nonseparable 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 selfenergy 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 HubbardI (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 selfenergy 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 selfenergy. 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 KohnSham basis:
(17) 
The partial density of states computed from the correlated Green’s function using is
(18)  
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
(19) 
The condition Eq. (16) can then be expressed as
(20) 
which is clearly satisfied when is invertible matrix because . This is satisfied when the KohnSham 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
(21) 
Here index runs over the local basis in which the green’s function is minimally offdiagonal (cubic harmonics or relativistic harmonics).
The projector is separable, as postulated in Eq. (19), and the transformation is
(22) 
with
(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 nonlocal 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 nonorthonormal 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 KohnSham basis, restricted to some subset of low energy bands. The local orbitals used for the projection were either allelectron atomic partial waves in the PAW framework, or pseudoatomic wave functions in mixedbasis 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 pseudoatomic 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 selfconsistent 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 selfenergy. This simplifies the selfconsistent DMFT problem, but makes it impossible to implement the charge selfconsistency. 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
(24) 
where runs over all space (orbitals,momenta) and time (frequency). The quantities apprearing in the above functional are
(25)  
(26)  
(27)  
where is trace over time only (not space), is the potentials due to ions, are the Hartree, and exchangecorrelation 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.,
(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
With the above functional dependence in mind, minimization with respect to gives
and minimization with respect to leads to
Hence the Hartree and exchangecorrelation potential are computed in the same way as in DFT method (note however is electron density in the presence of DMFT selfenergy), while the DMFT selfenergy 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 selfenergy . The impurity Green’s function is , hence the DMFT selfconsistency condition reads
(29) 
where , and is the interaction included in DFT (double counting). The selfconsistency condition takes the explicit form
(30) 
where and is the muffintin radius.
For efficient evaluation of the DMFT selfconsistency condition Eq. (30), we choose to work in the KohnSham (KS) basis. At each DFT+DMFT iteration, we first solve the KSeigenvalue problem
(31) 
Then we express the projection in KS basis, , where run over all bands. We then perform the embedding of the selfenergy, i.e., transforming it from DMFT base to the KS base
(32) 
In KSbase, we can invert the Green’s function Eq. (30), to obtain the practical form of the selfconsistency condition
(33)  
(34) 
This is of course equivalent to Eq. (30). Finally we solve this selfconsistency equation for a given selfenergy to obtain the hybridization function and the impurity levels .
We note in passing that the selfenergy is a complex function, and its imaginary part is related to the electronelectron 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 clusterDMFT is very straightforward. One needs to increase the unit cell to include more sites of the same atom type. The selfenergy 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 clusterDMFT case, the selfenergy in KSbasis Eq. (32) has to be summed over both and , and selfconsistency condition Eq. (34) becomes a matrix equation in . The challenging part of the clusterDMFT formalism is in solving the clusterimpurity 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 powerlaw, however, these techniques usually can not reach the interesting regime of strong correlations and low temperatures.
The major bottleneck in evaluating the DMFT selfconsistency 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. Nonnegligible 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 selfconsistency condition Eq. (34) can be inverted.
To optimize the sum in Eqs. (32) and (33), one can notice that local quantities like selfenergy and local green’s function possess a large degree of symmetry when written in proper basis (real harmonics, relativistic harmonics): many offdiagonal matrix elements vanish, and many matrix elements are equivalent. For example, in a system with cubic symmetry, one has only two types of selfenergy and . Hence, instead of summing over matrix elements in Eq. (32), one can rewrite the sum over two matrix elements , i.e.,
(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 nonspin polarized solution as starting point. In the ordered state, the DMFT selfenergy is allowed to break the symmetry, while typically the exchangecorrelation 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 KSeigenvalue problem
to obtaine KS eigenvectors, core, and semicore charge, and linearization energies .

We start with a guess for the lattice selfenergy correction (here is the dynamic part of the selfenergy 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 selfenergy (shifted by double counting) to KohnSham base by the transformation Eq. (32) to obtain .

: Using the current DMFT selfenergy , and the current DFT KSpotential , 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
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
and on real axis we solve
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
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 fourdimensional integral is evaluated analytically (see chapter V)
When DMFT is done on imaginary axis (using imaginary time impurity solvers), we evaluate
(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 selfconsistent solution. We devised the following method to remove this instability:

The diagonal components of the selfenergy were fitted by a polelike 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.,
(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 KohnSham 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.

: Impurity solver uses , , and Coulomb repulsion (which is represented by Slater integrals , , and ) as the input and gives the new selfenergy as the output.
Currently we integrated the following impurity solvers: OCA (see chapter VII.2), Noncrossing 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