Gutzwiller Density Functional Theory: a formal derivation and application to ferromagnetic nickel

We present a detailed derivation of the Gutzwiller Density Functional Theory that covers all conceivable cases of symmetries and Gutzwiller wave functions. The method is used in a study of ferromagnetic nickel where we calculate ground state properties (lattice constant, bulk modulus, spin magnetic moment) and the quasi-particle band structure. Our method resolves most shortcomings of an ordinary Density Functional calculation on nickel. However, the quality of the results strongly depends on the particular choice of the double-counting correction. This constitutes a serious problem for all methods that attempt to merge Density Functional Theory with correlated-electron approaches based on Hubbard-type local interactions.


Introduction
Density Functional Theory (DFT) is the workhorse of electronic structure theory [1]. Based on the Hohenberg-Kohn theorem [1], the ground-state properties of an interacting many-electron system are calculated from those of an effective single-particle problem that can be solved numerically. An essential ingredient in DFT is the so-called exchangecorrelation potential which, however, is unknown and sensible approximations must be devised, e.g., the local (spin) density approximation, L(S)DA. In this way, the electronic properties of metals were calculated systematically [2]. Unfortunately, the L(S)DA leads to unsatisfactory results for transition metals, their compounds, and for heavy-fermion systems. The electrons in the narrow 3d or 4f /5f bands experience correlations that are not covered by current exchange-correlation potentials.
For a more accurate description of electrons in narrow bands, Hubbard-type models [3,4] have been put forward. However, simplistic model Hamiltonians can describe limited aspects of real materials at best, while, at the same time, they reintroduce the full complexity of the many-body problem. Recently, new methods were developed that permit the (numerical) analysis of multi-band Hubbard models, and, moreover, can be combined with DFT, specifically, the LDA+U method [5], the LDA+DMFT (Dynamical Mean-Field Theory) [6,7], and the Gutzwiller variational approach [8,9,10,11]. The LDA+U approach treats atomic interactions on a mean-field level so that it is computationally cheap but it ignores true many-body correlations. The DMFT becomes formally exact for infinite lattice coordination number, Z → ∞, but it requires the self-consistent solution of a dynamical impurity problem that is numerically very demanding. The Gutzwiller DFT is based on a variational treatment of local manybody correlations. Expectation values can be calculated for Z → ∞ without further approximations, and the remaining computational problem remains tractable.
Previously, we used the DFT to obtain the bare band structure from which we calculated the properties of nickel [8,12,13,14] and LaOFeAs [15]. For these studies, we developed a formalism that applies to general Gutzwiller-correlated states for arbitrary multi-band Hubbard Hamiltonians. However, some single-particle properties remained fixed at their DFT values. In contrast, in Refs. [9,10,11] the correlated electron density was fed back into the DFT calculations but the Gutzwiller quasi-particle Hamiltonian was introduced in an ad-hoc manner. In this work, we bridge the gap between the two approaches. We present a formal derivation of the Gutzwiller DFT as a generic extension of the DFT. Our formulae apply for general Gutzwiller-correlated wave functions and reproduce expressions used previously [9,10] as special cases. As an application, we give results for nickel in face-centered cubic structure. The Gutzwiller DFT results for the lattice constant, the magnetic spin-only moment, and the bulk modulus agree very well with experiments. Moreover, the quasi-particle bandstructure from Gutzwiller DFT is in satisfactory agreement with data from Angle-Resolved Photo-Emission Spectroscopy (ARPES). As found earlier [9,10], the Gutzwiller DFT overcomes the limitations of DFT for the description of transition metals.
Our paper is organized as follows. In Sect. 2 we recall the derivation of Density Functional Theory (DFT) as a variational approach to the many-body problem. We extend our concise derivation to Density Functional Theory for many-particle Hamiltonians in Sect. 3. In particular, we formulate the Gutzwiller density functional whose minimization leads to the Gutzwiller-Kohn-Sham Hamiltonian. The theory is worked out in the limit of large coordination number, Z → ∞, where explicit expressions for the Gutzwiller density functional are available. In Sect. 4 we restrict ourselves to lattice systems that are invariant under translation by a lattice vector so that the quasiparticle excitations can be characterized by their Bloch momentum. In Sect. 5 we present results for face-centered cubic (fcc) nickel where Z = 12. A short summary, Sect. 6, closes our presentation. Technical details are deferred to the appendices.

Density Functional Theory
We start our presentation with a concise derivation of Density Functional Theory that can readily be extended to the Gutzwiller Density Functional Theory.

Many-particle Hamiltonian and Ritz variational principle
Our many-particle Hamiltonian for electrons with spin σ =↑, ↓ reads (h ≡ 1) with V (r − r ′ ) = 1 2 The electrons experience the periodic potential of the ions, U(r), and their mutual Coulomb interaction, V (r−r ′ ). The total number of electrons is N = N ↑ +N ↓ . According to the Ritz variational principle, the ground state of a HamiltonianĤ can be obtained from the minimization of the energy functional in the subset of normalized states |Ψ in the Hilbert space with N electrons, Ψ|Ψ = 1.

Levy's constrained search
The minimization of the energy functional (3)  To this end, we consider the subset of normalized states |Ψ (n) with fixed electron densities n σ (r), In the following we accept 'physical' densities only, i.e., those n σ (r) for which states |Ψ (n) can be found. For the subset of states |Ψ (n) we definê H e =Ĥ kin +V xc , Here, we extracted the Hartree terms from the Coulomb interaction H int in eq. (1) so thatV xc contains only the so-called exchange and correlation contributions. In the subset of normalized states |Ψ (n) we consider the functional For fixed densities n σ (r), the HamiltonianĤ e defines an electronic problem where the periodic potential of the ions is formally absent.

Constrained search.
The formal task is to find the minimum of the energy functional F in (6) with respect to |Ψ (n) , Recall that the electron densities n σ (r) are fixed in this step. We denote the resulting optimal many-particle state |Ψ (n) 0 . Thus, we may writē For later use, we define the functionals for the kinetic energy and the exchange-correlation energy so thatF According to the Ritz variational principle, the ground-state energy E 0 is found from the minimization of this functional over the densities n σ (r), The ground-state densities n 0 σ (r) are those where the minimum of D [{n σ (r)}] is obtained.

Single-particle Hamiltonian and Ritz variational principle
We consider the subset of single-particle product states |Φ (n) that are normalized to unity, Φ (n) |Φ (n) = 1. As before, the upper index indicates that they all lead to the same (physical) single-particle densities n sp σ (r), As our single-particle Hamiltonian we consider the kinetic-energy operatorĤ kin , see eq. (5). For fixed single-particle densities n sp σ (r) we define the single-particle kineticenergy functional 2.3.1. Constrained search. As in Sect. 2.2, we carry out a constrained search in the subset of states |Φ (n) . The task is the minimization of the kinetic-energy functional F sp {n sp σ (r)} , |Φ (n) . We denote the optimized single-particle product state |Φ (n) 0 so that we find the density functional for the kinetic energy as 2.3.2. Density Functional, ground-state density and ground-state energy. As the density functional D sp [{n sp σ (r)}] that corresponds to the single-particle problem we define with the kinetic energy term from (17)

2.3.3.
Noninteracting V -representability. In order to link the many-particle and singleparticle approaches we make the assumption of non-interacting V -representability [1]: For any given (physical) densities n σ (r) we can find a subset of normalized single-particle product states |Φ (n) with N electrons such that Moreover, we demand that the density functionals D [{n σ (r)}] (12) for the interacting electrons and D sp [{n σ (r)}] (18) for the single-particle problem agree with each other [17], Then, the single-particle problem leads to the same ground-state density n 0 σ (r) and ground-state energy E 0 as the interacting-particle Hamiltonian because the density variation is done with the same density functional (Hohenberg-Kohn theorem) [1].
The condition (20) is equivalent to because the interaction with the external potential and the Hartree term only depend on the densities. Eq. (21) then leads to an exact expression for the single-particle exchange-correlation energy This is our defining equation for E sp,xc [{n σ (r)}] in eq. (18).

Kohn-Sham Hamiltonian
In the following we address the single-particle energy functional directly, i.e., the Ritz variational problem without a prior constrained search, For the extension to the Gutzwiller Density Functional Theory in Sect. 3, we expand the field operators in a basis, where the index i typically represents a combination of site or crystal momentum index and an orbital index. For a canonical basis we must have completeness and orthogonality, When we insert (24) into (5), we obtain the operator for the kinetic energy in a general single-particle basis, where the elements of the kinetic-energy matrix T are given by with ξ i,σ (r) = r|i, σ .

Energy functional.
We introduce the single-particle density matrixρ. Its elements in the general single-particle basis read Then, the densities are given by Using these definitions, we can write the energy functional in the form The fact that |Φ are normalized single-particle product states is encoded in the matrix relationρ This is readily proven by using a unitary transformation between the operatorsĉ i,σ and the single-particle operatorsb k,σ that generate |Φ , see Appendix A.1.
When we minimize E [{n σ (r)} ,ρ] with respect toρ we must take the condition (31) into account using a matrix Ω of Lagrange multipliers Ω l,m;σ . Moreover, we use the Lagrange multipliers κ σ (r) to ensure eq. (29), i.e., altogether we address the functional

Minimization.
When we minimize G DFT in eq. (32) with respect to n σ (r) we find where V Har (r) is the Hartree interaction and v sp,xc,σ (r) is the single-particle exchangecorrelation potential. The minimization with respect toρ is outlined in Appendix A.2 [18]. It leads to the Kohn-Sham single-particle Hamiltonian where the elements of the Kohn-Sham Hamilton matrix T KS are given by Explicitly, Here, we defined the Kohn-Sham potential V KS σ (r) that, in our derivation, is identical to the Lagrange parameter κ σ (r).
The remaining task is to find the basis in which the Kohn-Sham matrix T KS is diagonal, see Appendix A.3.

Density Functional Theory for many-particle Hamiltonians
The Kohn-Sham potential (35) cannot be calculated exactly because the functionals in eq. (22) are not known. Therefore, assumptions must be made about the form of the single-particle exchange-correlation potential, e.g., the Local Density Approximation [1]. Unfortunately, such approximations are not satisfactory for transition metals and their compounds, and more sophisticated many-electron approaches must be employed.

Hubbard Hamiltonian and Hubbard density functional
3.1.1. Multi-band Hubbard model. A better description of transition metals and their compounds can be achieved by supplementing the single-particle system defined byĤ kin in Sect. 2.3 by a multi-band Hubbard interaction. Then, our new reference system is defined byĤ whereV loc describes local interactions between electrons in Wannier orbitals on the same site R. The local single-particle operatorV dc accounts for the double counting of their interactions in the Hubbard termV loc and in the single-particle exchange-correlation energy E sp,xc . We assume thatV loc andV dc do not depend on the densities n σ (r) explicitly.
For the local interaction we set Note that only electrons in the small subset of correlated orbitals (index c) experience the two-particle interactionV loc : When there are two electrons in the Wannier orbitals φ R,c 3 ,σ 3 (r) and φ R,c 4 ,σ 4 (r) centered around the lattice site R, they are scattered into the orbitals φ R,c 1 ,σ 1 (r) and φ R,c 2 ,σ 2 (r), centered around the same lattice site R; for the definition of basis states, see Appendix A.3. Typically, we consider c = 3d for the transition metals and their compounds. The interaction strengths are parameters of the theory. Later, we shall employ the spherical approximation so that U ··· ··· for d-electrons can be expressed in terms of three Racah parameters A, B, and C. Fixing C/B makes it possible to introduce an effective Hubbard parameter U and an effective Hund's-rule coupling J, see Sect. 5 and Appendix C.

Hubbard density functional.
According to Levy's constrained search, we must find the minimum of the functional in the subset of normalized states with given (physical) density n σ (r), see eq. (4). The H,0 of the HamiltonianĤ H for fixed densities n σ (r). In analogy to Sect. 2.3, we define the Hubbard density functional and E H,xc [{n σ (r)}] is the exchange-correlation energy forĤ H . As in Sect. 2.3, the Hubbard density functional agrees with the exact density functional if we choose Then, the Hubbard approach provides the exact ground-state densities and ground-state energy of our full many-particle Hamiltonian (Hubbard-Hohenberg-Kohn theorem). Of course, our derivation relies on the assumption of Hubbard V -representability of the densities n σ (r).

Hubbard single-particle potential.
When we directly apply the Ritz principle, we have to minimize the energy functional We include the constraints eq. (4) and the normalization condition using the Lagrange parameters κ σ (r) and E 0 in the functional As in Sect. 2.4, see eqs. (33) and (40), the variation of G H with respect to n σ (r) gives the single-particle potential The Hubbard-model approach is based on the idea that typical approximations for the exchange-correlation energy, e.g., the local-density approximation, are suitable for the Hubbard model, Indeed, as seen from eq. (46), in the Hubbard exchange-correlation energy E H,xc the exchange-correlation contributions in the exact E xc are reduced by the Hubbard term , reflecting a more elaborate treatment of local correlations. The minimization of (47) with respect to |Ψ constitutes an unsolvable manyparticle problem. The ground state |Ψ 0 is the solution of the many-particle Schrödinger equation with energy E 0 , with the single-particle Hamiltonian The Schrödinger equation (51) can be used as starting point for further approximations, for example, the Dynamical Mean-Field Theory (DMFT). In the following we will address the functional in eq. (47) directly.

Gutzwiller density functional
In the widely used LDA+U approach [5], the functional in eq. (47) is evaluated and (approximately) minimized by means of single-particle product wave functions. However, this approach treats correlations only on a mean-field level. In the more sophisticated Gutzwiller approach, we consider the functional in eq. (47) in the subset of Gutzwiller-correlated variational many-particle states.

Gutzwiller variational ground state.
In order to formulate the Gutzwiller variational ground state [4,8], we consider the local (atomic) states |Γ R that are built from the correlated orbitals. The local Hamiltonians take the form where |Γ R contains |Γ R | electrons. Here, we introduced and the local many-particle operatorsm R;Γ, The Gutzwiller correlator and the Gutzwiller variational states are defined aŝ Here, |Φ is a single-particle product state, and λ Γ,Γ ′ (R) defines the matrixλ(R) of, in general, complex variational parameters.

Gutzwiller functionals.
We evaluate the energy functional (47) in the restricted subset of Gutzwiller variational states, Note that we work with the orbital Wannier basis, see Appendix A.3, The elements of the Gutzwiller-correlated single-particle density matrix are and the densities become The expectation values for the atomic operators are given by The diagrammatic evaluation of ρ G (R ′ ,b ′ ),(R,b);σ and of m G R;Γ,Γ ′ shows that these quantities are functionals of the non-interacting single-particle density matricesρ, see eq. (28), and of the variational parameters λ Γ,Γ ′ (R). Moreover, it turns out that the local, noninteracting single-particle density matrixC(R) with the elements plays a prominent role in the Gutzwiller energy functional, in particular, for infinite lattice coordination number. Therefore, we may write In the Lagrange functional we shall impose the relation (61) with the help of the Hermitian Lagrange parameter matrixη with entries η (R,b),(R,b ′ );σ . Lastly, for the analytical evaluation of eq. (62) it is helpful to impose a set of (real-valued) local constraints (l = 1, 2, . . . , N con ) which we implement with real Lagrange parameters Λ l (R). In the following, we abbreviate i = (R, b) and j = (R ′ , b ′ ). Consequently, in analogy with Sect. 2.4, we address as our Lagrange functional, cf. eq. (32). Here, we took the condition (59) into account using Lagrange parameters κ σ (r) because the external potential, the Hartree term and the exchange-correlation potential in eq. (56) depend on the densities.
(i) As in the derivation of the exact Schrödinger equation (51), the variation of G G DFT with respect to n σ (r) generates the single-particle potential, see eqs. (40) and (49).
(ii) The minimization with respect toC(R) gives (iii) The minimization with respect to the Gutzwiller correlation parametersλ(R) results in for all λ Γ,Γ ′ (R). Note that, in the case of complex Gutzwiller parameters, we also have to minimize with respect to (λ Γ,Γ ′ (R)) * . Using these equations we may calculate the Lagrange parameters Λ l (R) that are needed in eq. (67).
(iv) The minimization of G G DFT with respect toρ generates the Landau-Gutzwiller quasiparticle Hamiltonian, see Appendix A.2, with the entries where we used eqs. (56) and (69). The single-particle state |Φ is the ground state of the Hamiltonian (70) from which the single-particle density matrixρ follows.
The minimization program outlined in steps (i)-(iv) requires the evaluation of the energy E G in eq. (56). In particular, the correlated single-particle density matrixρ G , eq. (58), must be determined. All equations derived in this section are completely general. They can, at least in principle, be evaluated by means of a diagrammatic expansion method [19,20,21]. The leading order of the expansion corresponds to an approximation-free evaluation of expectation values for Gutzwiller wave functions in the limit of high lattice coordination number. This limit will be studied in the rest of this work.

Gutzwiller density functional for infinite lattice coordination number
For Z → ∞, the Gutzwiller-correlated single-particle density matrix and the Gutzwiller probabilities for the local occupancies can be calculated explicitly without further approximations.

Local constraints.
As shown in Refs. [8,22] it is convenient for the evaluation of Gutzwiller wave functions to impose the following (local) constraints and where we abbreviated Â Φ ≡ Φ|Â|Φ . Note that, for complex constraints, the index l in (63) labels real and imaginary parts separately.

Atomic occupancies.
In the limit of infinite lattice coordination number, the interaction and double-counting energy can be expressed solely in terms of the local variational parametersλ(R) and the local density matrixC(R) of the correlated bands in |Φ , The remaining expectation values m R;Γ 1 ,Γ 4 Φ are evaluated using Wick's theorem. Explicit expressions are given in Refs. [8,23].

3.3.
3. Correlated single-particle density matrix. The local part of the correlated singleparticle density matrix is given by It can be evaluated using Wick's theorem. As can be seen from eq. (75), it is a function of the variational parameters λ Γ,Γ ′ (R) and of the local non-interacting single-particle density matrixC(R). For R = R ′ , we have for the correlated single-particle density matrix with the orbital-dependent renormalization factors q a,σ b,σ (R) for the electron transfer between different sites. Explicit expressions in terms of the variational parameters λ(R) and of the local non-interacting single-particle density matrixC(R) are given in Refs. [8,23].

Implementation for translational invariant systems
For a system that is invariant under translation by a lattice vector, all local quantities become independent of the site index, e.g., λ Γ,Γ ′ (R) ≡ λ Γ,Γ ′ for the Gutzwiller variational parameters. Since k from the first Brillouin zone is a good quantum number, we work with the (orbital) Bloch basis, see Appendix A.3.
As shown in Sect. 3.2.3, the minimization of the Gutzwiller energy functional requires two major steps, namely, the variation with respect to the Gutzwiller parametersλ and the variation with respect to the single-particle density matrixρ that characterizes the single-particle product state |Φ .

Gutzwiller-Kohn-Sham Hamiltonian
The minimization of the energy functional with respect to the single-particle density matrix leads to the Gutzwiller-Kohn-Sham Hamiltonian. In the orbital Bloch basis φ k,b,σ (r), see Appendix A.3, the corresponding quasi-particle Hamiltonian readŝ see eq. (70). In this section, we shall explain how this Hamiltonian can be calculated. Note, however, that the actual numerical implementation within QuantumEspresso is done in first quantization and uses plane-waves. We therefore derive the plane-wave representation of the Gutzwiller-Kohn-Sham equations in Appendix B.

Derivation of matrix elements.
When we apply the general expressions (71), the matrix elements of the quasi-particle Hamiltonian are obtained as where we have from eq. (69) Moreover, are the entries of the single-particle density matrix in the orbital Bloch basis, and ρ G b ′ ,b;σ (k) denotes the corresponding quantities in the Gutzwiller-correlated state. In the limit of infinite lattice coordination number, we may express V G loc/dc in eq. (71) as a function of the Gutzwiller variational parameters λ Γ,Γ ′ and of the local density matrix C. Therefore, V G loc/dc are formally independent of the single-particle density matrixρ so that they do not contribute to h G b,b ′ ;σ (k). Equation (78) shows that, apart from an overall shift of the orbitals through η b,b ′ ;σ , we can still work with the matrix elements h 0 a,a ′ ;σ ′ (p) of the single-particle operatorĤ 0 that enters the many-particle Schrödinger equation (51).
In the orbital Bloch basis, eqs. (75) and (76) take the form When we insert eq. (81) into eq. (78) we thus find for the entries of the Gutzwiller quasi-particle Hamiltonian Recall that the local single-particle densitiesC are treated as independent parameters. Since the q-factors and the correlated local single-particle density matrix C G b,b ′ ;σ = ρ G (R,b),(R,b ′ );σ are solely functions of the variational parametersλ and ofC, they are treated as constants when we take the derivative of ρ G b ′ ,b;σ (k) with respect to ρ b ′ ,b;σ (k) in eq. (81).

Diagonalization of the quasi-particle Hamiltonian. The unitary matrixF
which provides the quasi-particle dispersion ǫ G n,σ (k). We introduce the quasi-particle band operatorŝ in which the Gutzwiller quasi-particle Hamiltonian becomes diagonal, In order to minimize the Gutzwiller density functional G G DFT , we must work with the ground state of the quasi-particle HamiltonianĤ G qp , where the N levels lowest in energy are occupied as indicated by the prime at the product, ǫ G n,σ (k) ≤ E G F,σ . Using eq. (80) we find where the quasi-particle occupancies are unity for occupied quasi-particle levels up to the Fermi energy E G F,σ , and zero otherwise. The particle densities follow from eq. (87), where ρ G b ′ ,b;σ (k) is given by eq. (81). As in DFT, the particle densities must be calculated self-consistently.

Band-shift parameters.
In order to determine the band-shift parameters η b,b ′ ;σ , we must evaluate eq. (67) using the Gutzwiller energy functional in the limit of infinite lattice coordination number. We define the kinetic energy of the Gutzwiller quasiparticles as It is easy to show that Then, eq. (67) gives the effective local hybridizations η b,b ′ ;σ Note that the term in the second line in eq. (92) often vanishes due to symmetry, e.g., in nickel, because C G c,c ′ ;σ = C c,c ′ ;σ .

Minimization with respect to the Gutzwiller parameters
In the ('inner') minimization with respect to the Gutzwiller parameters λ Γ,Γ ′ (which are now independent of R) we assume that the single-particle state |Φ is fixed. Then we have to minimize the function where we introduced Here, the indices c andc denote correlated and non-correlated orbitals, respectively. Note that the sum over b, b ′ in the last line of eq. (93) only contributes in the minimization if (at least) one of these two indices belongs to a correlated orbital. An efficient algorithm for the minimization of (93) has been introduced in Ref. [23]. This minimization also gives us the Lagrange parameters Λ l which enter the outer minimization in eq. (92).

Local Hamiltonian and double-counting corrections
For a Gutzwiller DFT calculation we need to specify the Coulomb parameters in the local Hamiltonian (42) and the form of the double-counting operator in (41).

Cubic symmetry and spherical approximation.
In many theoretical studies one uses a Hamiltonian with only density-density interactions, Here, we introduced↑ =↓ (↓ =↑) and and J(c, c ′ ) are the local Hubbard and Hund's-rule exchange interactions. An additional and quite common approximation is the use of orbital-independent Coulomb parameters, For a system of five correlated 3d orbitals in a cubic environment as in nickel, however, the Hamiltonian (96) is incomplete [24]. The full Hamiltonian readŝ whereV n.dens. loc Here, t = ζ, η, ξ and e = u, v are indices for the three t 2g orbitals with symmetries ζ = xy, η = xz, and ξ = yz, and the two e g orbitals with symmetries u = 3z 2 − r 2 and v = x 2 − y 2 , respectively. The parameters A(t), T (t), S(t, t ′ ; t ′′ , e) in eq. (99) are of the same order of magnitude as the exchange interactions J(c, c ′ ) and, hence, there is no a-priori reason to neglect V n.dens.

loc
. Of all the parameters U(c, c ′ ), J(c, c ′ ), A(t), T (t), S(t, t ′ ; t ′′ , e) only ten are independent in cubic symmetry, see Appendix C.
When we assume that all 3d-orbitals have the same radial wave-function ('spherical approximation'), all parameters are determined by, e.g., the three Racah parameters A, B, C. For comparison with other work, we introduce the average Coulomb interaction between electrons in the same 3d-orbitals, U = c U(c, c)/5 = A + 4B + 3C, the average Coulomb interaction between electrons in different orbitals, U ′ = c<c ′ U(c, c ′ )/10 = A − B + C, and the average Hund's-rule exchange interaction, J = c<c ′ J(c, c ′ )/10 = 5B/2 + C that are related by the symmetry relation U ′ = U − 2J, see Appendix C. Due to this symmetry relation, the three values of U, U ′ , and J do not determine the Racah parameters A, B, C uniquely. Therefore, we make use of the relation C/B = 4 which is a reasonable assumption for metallic nickel [8,24]. In this way, the three Racah parameters and, consequently, all parameters inV full loc are functions of U and J, A = U − 32J/13, B = 2J/13, C = 8J/13. This permits a meaningful comparison of our results for all local Hamiltonians. Later we shall compare our results forV dens loc with orbital-independent values for U, U ′ and J = (U − U ′ )/2, see eq. (97), with those for the full local Hamiltonian,V full loc , for the same values for U and J.

Double counting corrections.
There exists no systematic (let alone rigorous) derivation of the double-counting correction in eq. (41). A widely used form for this operator has first been introduced in the context of the LDA+U method. Its expectation value is given by wheren σ ≡ Nc c=1 C G c,c;σ ,n ≡n ↑ +n ↓ , and N c is the number of correlated orbitals (N c = 10 for nickel). Note that only the two mean values U and J enter this doublecounting operator, i.e., it is the same for all local Hamiltonians introduced above.
The physical consequences of the double-counting correction are most pronounced in its impact on the local energy-shifts η c,c;σ which we may write as For nickel, the cubic symmetry guarantees that i.e., the correlated and uncorrelated local densities agree with each other. The doublecounting correction (100) leads to η dc,1 c,c;σ = U(n − 1/2) + J(n σ − 1/2). It has been argued in Ref. [11] that this double-counting correction is insufficient for the investigation of Cerium and some of its compounds. Instead, the authors of that work propose two alternative double-counting corrections which, in effect, correspond to the energy shifts As we will demonstrate in the following section, these three double-counting corrections lead to noticeably different results for nickel. This is a rather unsatisfactory observation because it compromises the predictive power of the method if the results strongly depend on the particular choice of the double-counting correction. In Ref. [25] a scheme has been proposed which does not rely on the subtraction of double-counting operators and instead addresses the density functional directly. It needs to be seen if this method can provide a more general way to tackle the double-counting problem within the Gutzwiller DFT.

Implementation in DFT
We implemented our Gutzwiller scheme in the QuantumEspresso DFT code; for details on QuantumEspresso, see Ref. [26]. Due to the cubic symmetry of nickel, the single-particle density matrixC is diagonal with the local occupancies C t;σ ≡ C ξ,ξ;σ = C η,η;σ = C ζ,ζ;σ in the t 2g -orbitals and C e;σ ≡ C u,u;σ = C v,v;σ in the e g -orbitals. Likewise, the matrixη is diagonal with the corresponding entries η t;σ and η e;σ . Moreover, the q-matrix is diagonal, q b,σ a,σ = δ a,b q a,σ with identical entries for the three t 2g -orbitals, q ξ,σ = q η,σ = q ζ,σ ≡ q t,σ , and the two e g -orbitals, q u,σ = q v,σ ≡ q e,σ , respectively. Formulae for q t,σ and q e,σ as a function of the Gutzwiller parametersλ and of C e,;σ and C t;σ are given in Refs. [8,23].

Setup: DFT calculation and Wannier orbitals.
As a first step, we perform a DFT calculation that corresponds to setting U = J = 0. We use the LDA exchangecorrelation potential, This calculation provides the Kohn-Sham bandstructure ǫ n,σ (k) and the coefficients C G,n,σ (k) = k, G, σ|k, n, σ of the Kohn-Sham eigenstates ψ k,n,σ (r) in the plane-wave basis. The implemented 'poor-man Wannier' program package provides the down-folded 3d Wannier orbitals φ R,c,σ (r). In the orbital Bloch basis the coefficients k, G, σ|k, c, σ describe φ k,c,σ (r) in the plane-wave basis.

Gutzwiller-Kohn-Sham loop.
At the beginning we set q b,σ a,σ = δ a,b andη = 0. Our Gutzwiller-Kohn-Sham loop consists of the following steps.
(ii) Perform the inner minimization, i.e., minimize the energy functional E inner in eq. (93). This step provides the values for the Lagrange parameters Λ l and for the Gutzwiller variational parametersλ that determine the orbital-dependent renormalization factors q c,σ in eq. (76).
(iv) If the total energy does not decrease compared with the previous iteration, the calculation has converged and the loop terminates. If not, repeat the loop starting at step (i).
The steps (ii) and (iii) are carried out following the algorithm outlined previously [23].
In the present version of the program, step (i) requires a full DFT calculation which, however, is numerically cheap for the simple nickel system. In the future, we plan to include the Gutzwiller minimization directly in the DFT minimization cycle.
The Gutzwiller approach permits the definition of correlated orbital Bloch states, see Appendix B. Therefore, we can compare our original 3d Wannier orbitals with the Gutzwiller correlated Wannier orbitals. For nickel, we find that the deviations are negligibly small. In general, we may include the correlation-induced shape changes of the correlated Wannier orbitals in our self-consistent calculations.

Results
The electronic properties of nickel have already been investigated by means of Gutzwiller wave functions in Refs. [12,13,14]. In these works we started from a paramagnetic DFT-LDA calculation that provided the band parameters for a tight-binding model. In order to overcome the deficiencies in the underlying DFT-LDA results, we fixed the magnetic moment and other single-particle properties at their experimental values. As we will show in this section, the Gutzwiller DFT mends most of the DFT-LDA shortcomings.
As a variational approach, the Gutzwiller DFT is expected to be most suitable for the calculation of ground-state properties such as the lattice constant, the magnetic moment, or the Fermi surface of a Fermi liquid. Although more speculative than the ground-state calculations, it is also common to interpret the eigenvalues of the Gutzwiller-Kohn-Sham Hamiltonian ǫ G n,σ (k) as the dispersion of the single-particle excitations [27]. We shall discuss our results for the ground-state properties and singleparticle excitations separately. Fig. 1, we show the lattice constant and the magnetic moment as a function of U (1 eV ≤ U ≤ 14 eV) for four different values of J/U (J/U = 0, 0.05, 0.075, 0.10). In these calculations we used the full local HamiltonianV full loc and the double-counting correctionV dc;1 . As is well known, the DFT-LDA underestimates the lattice constant. We obtain a LDA 0 = 6.47a B , considerably smaller than the experimental value of a 0 = 6.66a B where a B = 0.529177Å is the Bohr radius. Fig. 1 shows that the Hubbard interaction U increases the lattice constant whereby the Hund's-rule exchange J diminishes the slope. Apparently, a good agreement with the experimental lattice constant requires substantial Hubbard interactions, U > 10 eV. Fig. 1 shows the well-known fact that DFT-LDA reproduces the experimental value for the spin-only magnetic moment m so very well, m LDA so = 0.58µ B vs. m exp so = 0.55µ B . However, when the DFT-LDA calculation is performed for the experimental value of the lattice constant, the magnetic moment is grossly overestimated. As seen in Fig. 1, the Gutzwiller DFT allows us to reconcile the experimental findings both for the lattice constant and the magnetic moment if we work in the parameter range 11 eV < U < 14 eV and 0.05 < J/U < 0.07. Note that a 'fine-tuning' of parameters is not required to obtain a reasonable agreement between theory and experiment for the lattice constant and spin-only magnetic moment.

Lattice constant, magnetic moment, and bulk modulus of nickel. In
For nickel, detailed information about the quasi-particle bands is available. The quasi-particle dispersion at various high-symmetry points in the Brillouin zone is more sensitive to the precise values of U and J. As we shall show below in more detail, we obtain a satisfactory agreement with ARPES data for the choice (U opt = 13 eV, J opt = 0.9 eV) with an uncertainty of ±1 in the last digit. For our optimal values we show in Fig. 2 the ground-state energy per particle E(a)/N as a function of the fcc lattice constant a together with a second-order polynomial fit. The minimum is obtained at a 0 = 6.63a B , in good agreement with the experimental value a exp 0 = 6.66a B . For the magnetic spin-only moment we obtain m so = 0.52µ B , in good agreement with the experimental value m exp so = 0.55µ B .  From the curvature of E(a)/N at a = a 0 we can extract the bulk modulus. The bulk modulus at zero temperature is defined as the second-derivative of the ground-state energy with respect to the volume, This implies the Taylor expansion E(V ) = E(V 0 ) + (KV 0 /2)(V /V 0 − 1) 2 + . . . for the ground-state energy as a function of the volume V = a 3 . For the ground-state energy per particle we can thus write E(a)/N = E(a 0 )/N + e 2 (a/a B − a 0 /a B ) 2 + . . . with where we took into account that the fcc unit cell hosts four atoms, V 0 = Na 3 0 /4. The fit leads to K = 169 GPa, in good agreement with the experimental value, K = 182 GPa [28]. It is a well-known fact that the DFT-LDA overestimates the bulk modulus of nickel. Indeed, our DFT-LDA gives K LDA = 245 GPa. We also calculated the lattice parameter and the magnetic spin-only moment for the density-dependent interaction V dens loc , see eq. (96), with the same double-counting correctionV dc;1 . Our results do not show significant discrepancies for the groundstate properties. Note, however, that nickel is a special case because it has an almost filled 3d-shell (n 3d ≈ 9/10) such that the terms fromV n.dens. loc in eq. (99) are more or less deactivated. Preliminary calculations for iron indicate that these terms are more important for a partially filled 3d-shell.
The combination of the full local interaction V full loc with the second double-counting correction, see eq. (103), leads to very similar results for the lattice constant, spin-only magnetic moment, and compressibility for nickel. We suspect that the almost filled 3dshell is the reason for this observation. However, we could not obtain converged results for the third double-counting correction, see eq. (104). It appears as if the assumption of an average band shift is not useful, at least for nickel.

Symmetry
ExperimentV Theoretical data show the spin average and the exchange splittings in square brackets. Λ 3;1/2 denotes the point half-way on the Λ-line that links the points Γ and L. The first column gives experimental data compiled in [8], the second, third, and fourth column give theoretical results results forV full loc withV dc;1 ,V full loc withV dc;2 , andV dens loc witĥ V dc;1 , respectively, at (U opt = 13 eV, J opt = 0.9 eV). Fig. 3 we show the quasi-particle band structure of fcc nickel for (U opt = 13 eV, J opt = 0.9 eV). The most prominent effect of the Gutzwiller correlator is the reduction of the 3d bandwidth. From a paramagnetic DFT-LDA calculation one can deduce W LDA = 4.5 eV [12,13,14]. whereas we find W = 3.3 eV, in agreement with experiment. This bandwidth reduction is due to the q-factors q t,↑ = 0.851, q t,↓ = 0.824, q e,↑ = 0.852, q e,↓ = 0.819,q = σ (3q t,σ + 2q e,σ )/10 = 0.837, so that W ≈q 2 W LDA .

Quasi-particle bands of nickel. In
A more detailed comparison of the quasi-particle bands structure with experiment is given in table 1. The overall agreement between experiment and theory forV full loc witĥ V dc;1 is quite satisfactory. In particular, only one hole ellipsoid is found at the X-point, in agreement with experiment and in contrast to the DFT-LDA result [8]. Note, however, that the second double-counting correctionV dc;2 spoils this advantage. Therefore, this form of the double-correction term is not particularly useful for nickel.
We comment on two noticeable discrepancies between theory and experiment. First, the energy of the band L 2 ′ at the L-point deviates by a factor of five. This is an artifact that occurs already at the DFT-LDA level and is not cured by the Gutzwiller approach.
Since the level has pure 3p character around the L point, the origin of the discrepancy is related to the uncertainties in the partial charge densities n 3p,3s , n 3d in the 3p and 3p/3s bands. Second, the Gutzwiller DFT prediction for the exchange splitting ∆ t 2g (X 5 ) of the t 2g bands at the X-point is a factor of two larger than in experiment. This deviation is related to the fact that, quite generally, all bands are slightly too low in energy. This can be cured by decreasing U and increasing J but this deteriorates the values for the lattice constant and the magnetic moment. We suspect that the deviations are partly due to the use of a heuristic double-counting correction and the neglect of the spin-orbit coupling. Moreover, we expect the results for the band structure to improve when we replace the 'poor-man Wannier' orbitals for the correlated 3d electrons by more sophisticated wave functions. Table 1 also shows the results for V dens loc with density-density interactions only and withV dc;1 as double-counting correction. The description provides the correct Fermi surface topology but the deviations from the experimental band energies is significantly larger. In particular, the exchange splitting ∆ eg (X 2 ) of the e g bands at the X-point becomes negative, i.e., the order of the majority and minority bands is inverted. The comparison of the band structures shows that the full atomic Hamiltonian should be used for a detailed description of the quasi-particle bands in nickel.

Summary and conclusions
In this work, we presented a detailed derivation of the Gutzwiller Density Functional Theory. Unlike previous studies, our formalism covers all conceivable cases of symmetries and Gutzwiller wave functions. Moreover, our theory is not based on the 'Gutzwiller approximation' which corresponds to an evaluation of expectation values in the limit of infinite lattice coordination number. It is only in the last step that we resort to this limit.
In particular, our derivation consists of three main steps.
1. The density functional of the full many-particle system is related to that of a reference system with Hubbard-type local Coulomb interactions in the correlated orbitals. This generalizes the widely used Kohn-Sham scheme where a singleparticle reference system is used.
2. The energy functional of the Hubbard-type reference system is (approximately) evaluated by means of Gutzwiller variational wave functions.
3. Analytical results for the energy functional are derived with the Gutzwiller approximation.
In a first application we studied the electronic properties of ferromagnetic nickel. It turned out that the Gutzwiller DFT resolves the main deficiencies of DFT in describing ground-state properties such as the lattice constant, the magnetic moment, or the bulk modulus of nickel. Note that our approach requires the relatively large value U ≈ 13 eV for the local Coulomb interaction in order to obtain a good agreement with experiments.
Our results for the quasi-particle band structure are by and large satisfactory. In fact, a perfect agreement with ARPES data would be surprising because we calculate these quantities based on Fermi-liquid assumptions that are strictly valid only in the vicinity of the Fermi surface. Moreover, the quasi-particle energies strongly depend on the orbital occupations that are influenced by the somewhat arbitrary choice of the double-counting corrections. As we have also shown in this work, different forms of the double-counting correction from the literature lead to fairly different results for nickel. We consider this as the main shortcoming of the Gutzwiller DFT in its present form that we shall address in future studies. Hence,H andρ must have the same basis of (single-particle) eigenvectors and, consequently, we find an extremum of E(ρ) if we choose |Φ as an eigenstate of where V is the crystal volume. As a consequence of the lattice periodicity, the crystal momentum k from the first Brillouin zone is a good quantum number. As seen from eq. (39), the Kohn-Sham Hamiltonian is diagonalized for the singleparticle states ψ k,n,σ (r) = r|k, n, σ that obey h KS σ (r)ψ k,n,σ (r) = ǫ n,σ (k)ψ k,n,σ (r) , (A. 13) where n is the band index. Eqs. (A.13) are the Kohn-Sham equations [1]. In its eigenbasis, the Kohn-Sham Hamiltonian takes the form H KS = k,n,σ ǫ n,σ (k)b † k,n,σbk,n,σ . (A.14) Its ground state is given by |Φ 0 = σ k,n ′b † k,n,σ |vac , (A. 15) where the N levels lowest in energy are occupied as indicated by the prime at the product, ǫ n,σ (k) ≤ E F,σ . Then, f k,n,σ = Φ 0 |b † k,n,σbk,n,σ |Φ 0 = Θ (E F,σ − ǫ n,σ (k)) (A. 16) is unity for occupied levels up to the Fermi energy E F,σ , and zero otherwise.