Gutzwiller density functional theory: a formal derivation and application to ferromagnetic nickel

The formal task is to find the minimum of the energy functional F in (8) with respect to $|{{\Psi }^{(n)}}\rangle$,

$$\begin{eqnarray}&&\bar{F}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]={\rm Mi}{{{\rm n}}_{\left\{ |{{\Psi }^{(n)}}\rangle \right\}}}F\left[ \{{{n}_{\sigma }}({\bf r})\},\left\{ |{{\Psi }^{(n)}}\rangle \right\} \right]\ .\end{eqnarray} \tag{ 9 }$$

Recall that the electron densities ${{n}_{\sigma }}({\bf r})$ are fixed in this step. We denote the resulting optimal many-particle state $|\Psi _{0}^{(n)}\rangle$. Thus, we may write

$$\begin{eqnarray}&&\bar{F}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]=F\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ |\Psi _{0}^{(n)}\rangle \right\} \right]=\langle \Psi _{0}^{(n)}|{{\hat{H}}_{{\rm e}}}|\Psi _{0}^{(n)}\rangle .\end{eqnarray} \tag{ 10 }$$

For later use, we define the functionals for the kinetic energy

$$\begin{eqnarray}&&K\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]=\langle \Psi _{0}^{(n)}|{{\hat{H}}_{{\rm kin}}}|\Psi _{0}^{(n)}\rangle \end{eqnarray} \tag{ 11 }$$

and the exchange-correlation energy

$$\begin{eqnarray}&&{{E}_{{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]=\langle \Psi _{0}^{(n)}|{{\hat{V}}_{{\rm xc}}}|\Psi _{0}^{(n)}\rangle \end{eqnarray} \tag{ 12 }$$

so that

$$\begin{eqnarray}&&\bar{F}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]=K\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{E}_{{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right].\end{eqnarray} \tag{ 13 }$$

After the constrained search as a first step, we are left with the density functional

$$\begin{eqnarray}\begin{array}{rcl} D\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] & = & \bar{F}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+U\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{V}_{{\rm Har}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] \\ {} & = & K\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+U\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{V}_{{\rm Har}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{E}_{{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] \\ \end{array}\end{eqnarray} \tag{ 14 }$$

with

$$\begin{eqnarray}\begin{array}{rcl} U\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] & = & \mathop{\sum }\limits_{\sigma }\int {\rm d}{\bf r}U({\bf r}){{n}_{\sigma }}({\bf r}), \\ {{V}_{{\rm Har}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] & = & \mathop{\sum }\limits_{\sigma ,\sigma ^{\prime} }\int {\rm d}{\bf r}\int {\rm d}{\bf r}^{\prime} V({\bf r}-{\bf r}^{\prime} ){{n}_{\sigma }}({\bf r}){{n}_{\sigma ^{\prime} }}({\bf r}^{\prime} ). \\ \end{array}\end{eqnarray} \tag{ 15 }$$

According to the Ritz variational principle, the ground-state energy E₀ is found from the minimization of this functional over the densities ${{n}_{\sigma }}({\bf r})$,

$$\begin{eqnarray}&&{{E}_{0}}={\rm Mi}{{{\rm n}}_{\left\{ {{n}_{\sigma }}({\bf r}) \right\}}}D\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right].\end{eqnarray} \tag{ 16 }$$

The ground-state densities $n_{\sigma }^{0}({\bf r})$ are those where the minimum of $D\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]$ is obtained.

As in section 2.2, we carry out a constrained search in the subset of states $|{{\Phi }^{(n)}}\rangle$. The task is the minimization of the kinetic-energy functional ${{F}_{{\rm sp}}}[\left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\},\left\{ |{{\Phi }^{(n)}}\rangle \right\}]$. We denote the optimized single-particle product state $|\Phi _{0}^{(n)}\rangle$ so that we find the density functional for the kinetic energy as

$$\begin{eqnarray}&&{{\bar{F}}_{{\rm sp}}}\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right]=\langle \Phi _{0}^{(n)}|{{\hat{H}}_{{\rm kin}}}|\Phi _{0}^{(n)}\rangle \equiv {{K}_{{\rm sp}}}\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right].\end{eqnarray} \tag{ 19 }$$

As the density functional ${{D}_{{\rm sp}}}\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right]$ that corresponds to the single-particle problem we define

$$\begin{eqnarray}\begin{array}{rcl} {{D}_{{\rm sp}}}\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right] & = & {{K}_{{\rm sp}}}\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right]+U\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right]+{{V}_{{\rm Har}}}\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right] \\ {} & {} & +{{E}_{{\rm sp},{\rm xc}}}\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right], \\ \end{array}\end{eqnarray} \tag{ 20 }$$

with the kinetic energy term from (19), the contributions from the external potential and the Hartree terms $U\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right]$ and ${{V}_{{\rm Har}}}\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right]$ from equation (15), and the single-particle exchange-correlation potential ${{E}_{{\rm sp},{\rm xc}}}\left[ \left\{ n_{\sigma }^{{\rm sp}}({\bf r}) \right\} \right]$ that we will specify later. The functional (20) defines our single-particle reference system.

In order to link the many-particle and single-particle approaches we make the assumption of non-interacting V-representability [1]: for any given (physical) densities ${{n}_{\sigma }}({\bf r})$ we can find a subset of normalized single-particle product states $|{{\Phi }^{(n)}}\rangle$ with N electrons such that

$$\begin{eqnarray}&&n_{\sigma }^{{\rm sp}}({\bf r})={{n}_{\sigma }}({\bf r}).\end{eqnarray} \tag{ 21 }$$

Moreover, we demand that the density functionals $D\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]$ (14) for the interacting electrons and ${{D}_{{\rm sp}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]$ (20) for the single-particle problem agree with each other⁶ ,

$$\begin{eqnarray}&&{{D}_{{\rm sp}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]=D\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right].\end{eqnarray} \tag{ 22 }$$

Then, the single-particle problem leads to the same ground-state density $n_{\sigma }^{0}({\bf r})$ and ground-state energy E₀ as the interacting-particle Hamiltonian because the density variation is done with the same density functional (Hohenberg–Kohn theorem) [1].

The condition (22) is equivalent to

$$\begin{eqnarray}&&{{K}_{{\rm sp}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{E}_{{\rm sp},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]=K\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{E}_{{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]\end{eqnarray} \tag{ 23 }$$

because the interaction with the external potential and the Hartree term only depend on the densities. Equation (23) then leads to an exact expression for the single-particle exchange-correlation energy

$$\begin{eqnarray}&&{{E}_{{\rm sp},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]=K\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]-{{K}_{{\rm sp}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{E}_{{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right].\end{eqnarray} \tag{ 24 }$$

This is our defining equation for ${{E}_{{\rm sp},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]$ in equation (20).

We introduce the single-particle density matrix $\tilde{\rho }$. Its elements in the general single-particle basis read

$$\begin{eqnarray}&&{{\rho }_{(i,\sigma ),(j,\sigma )}}=\langle \Phi |\hat{c}_{j,\sigma }^{\dagger }\hat{c}_{i,\sigma }^{{}}|\Phi \rangle \equiv {{\rho }_{i,j;\sigma }}.\end{eqnarray} \tag{ 30 }$$

Then, the densities are given by

$$\begin{eqnarray}&&{{n}_{\sigma }}({\bf r})=\mathop{\sum }\limits_{i,j}\xi _{i,\sigma }^{*}({\bf r}){{\xi }_{j,\sigma }}({\bf r}){{\rho }_{j,i;\sigma }}.\end{eqnarray} \tag{ 31 }$$

Using these definitions, we can write the energy functional in the form

$$\begin{eqnarray}\begin{array}{rcl} E\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\tilde{\rho } \right] & = & \mathop{\sum }\limits_{i,j}\mathop{\sum }\limits_{\sigma }{{T}_{i,j;\sigma }}{{\rho }_{j,i;\sigma }}+U\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{V}_{{\rm Har}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] \\ {} & {} & +{{E}_{{\rm sp},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]. \\ \end{array}\end{eqnarray} \tag{ 32 }$$

The fact that $|\Phi \rangle$ are normalized single-particle product states is encoded in the matrix relation

$$\begin{eqnarray}&&\tilde{\rho }\cdot \tilde{\rho }=\tilde{\rho }.\end{eqnarray} \tag{ 33 }$$

This is readily proven by using a unitary transformation between the operators ${{\hat{c}}_{i,\sigma }}$ and the single-particle operators ${{\hat{b}}_{k,\sigma }}$ that generate $|\Phi \rangle$, see appendix A.1.

When we minimize $E\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\tilde{\rho } \right]$ with respect to $\tilde{\rho }$ we must take the condition (33) into account using a matrix $\tilde{\Omega }$ of Lagrange multipliers ${{\Omega }_{l,m;\sigma }}$. Moreover, we use the Lagrange multipliers ${{\kappa }_{\sigma }}({\bf r})$ to ensure equation (31), i.e., altogether we address the functional ${{G}_{{\rm DFT}}}\equiv {{G}_{{\rm DFT}}}\left[ \tilde{\rho },\tilde{\Omega },\left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ {{\kappa }_{\sigma }}({\bf r}) \right\} \right]$

$$\begin{eqnarray}\begin{array}{rcl} {{G}_{{\rm DFT}}} & = & E\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\tilde{\rho } \right]-\mathop{\sum }\limits_{l,m,\sigma }{{\Omega }_{l,m;\sigma }}\left( \mathop{\sum }\limits_{p}{{\rho }_{l,p;\sigma }}{{\rho }_{p,m;\sigma }}-{{\rho }_{l,m;\sigma }} \right) \\ {} & {} & -\mathop{\sum }\limits_{\sigma }\int {\rm d}{\bf r}{{\kappa }_{\sigma }}({\bf r})\left( {{n}_{\sigma }}({\bf r})-\mathop{\sum }\limits_{i,j}\xi _{i,\sigma }^{*}({\bf r}){{\xi }_{j,\sigma }}({\bf r}){{\rho }_{j,i;\sigma }} \right)\;. \\ \end{array}\end{eqnarray} \tag{ 34 }$$

When we minimize ${{G}_{{\rm DFT}}}$ in equation (34) with respect to ${{n}_{\sigma }}({\bf r})$ we find

$$\begin{eqnarray}&&{{\kappa }_{\sigma }}({\bf r})=U({\bf r})+{{V}_{{\rm Har}}}({\bf r})+{{v}_{{\rm sp},{\rm xc},\sigma }}({\bf r}),\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{{V}_{{\rm Har}}}({\bf r})\equiv \mathop{\sum }\limits_{\sigma ^{\prime} }\int {\rm d}{\bf r}^{\prime} 2\;V({\bf r}-{\bf r}^{\prime} )n_{\sigma ^{\prime} }^{0}({\bf r}^{\prime} ),\end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray}\begin{array}{rcl} {{v}_{{\rm sp},{\rm xc},\sigma }}({\bf r}) & \equiv & {{\left. \frac{\partial {{E}_{{\rm sp},{\rm xc}}}\left[ \left\{ {{n}_{\sigma ^{\prime} }}({\bf r}^{\prime} ) \right\} \right]}{\partial {{n}_{\sigma }}({\bf r})} \right|}_{{{n}_{\sigma }}({\bf r})=n_{\sigma }^{0}({\bf r})}} \\ {} & = & {{\left. \frac{\partial \left[ K\left[ \left\{ {{n}_{\sigma ^{\prime} }}({\bf r}^{\prime} ) \right\} \right]-{{K}_{{\rm sp}}}\left[ \left\{ {{n}_{\sigma ^{\prime} }}({\bf r}^{\prime} ) \right\} \right]+{{E}_{{\rm xc}}}\left[ \left\{ {{n}_{\sigma ^{\prime} }}({\bf r}^{\prime} ) \right\} \right] \right]}{\partial {{n}_{\sigma }}({\bf r})} \right|}_{{{n}_{\sigma }}({\bf r})=n_{\sigma }^{0}({\bf r})}}, \\ \end{array}\end{eqnarray} \tag{ 37 }$$

where ${{V}_{{\rm Har}}}({\bf r})$ is the Hartree interaction and ${{v}_{{\rm sp},{\rm xc},\sigma }}({\bf r})$ is the single-particle exchange-correlation potential.

The minimization with respect to $\tilde{\rho }$ is outlined in appendix A.2 [18]. It leads to the Kohn–Sham single-particle Hamiltonian

$$\begin{eqnarray}&&{{\hat{H}}^{{\rm KS}}}=\mathop{\sum }\limits_{i,j,\sigma }T_{i,j;\sigma }^{{\rm KS}}\hat{c}_{i,\sigma }^{\dagger }\hat{c}_{j,\sigma }^{{}},\end{eqnarray} \tag{ 38 }$$

where the elements of the Kohn–Sham Hamilton matrix $\tilde{T}_{\sigma }^{{\rm KS}}$ are given by

$$\begin{eqnarray}&&T_{i,j;\sigma }^{{\rm KS}}=\frac{\partial E\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\tilde{\rho } \right]}{\partial {{\rho }_{j,i;\sigma }}}+\int {\rm d}{\bf r}{{\kappa }_{\sigma }}({\bf r})\xi _{i,\sigma }^{*}({\bf r}){{\xi }_{j,\sigma }}({\bf r}).\end{eqnarray} \tag{ 39 }$$

Explicitly,

$$\begin{eqnarray}&&T_{i,j;\sigma }^{{\rm KS}}=\int {\rm d}{\bf r}\xi _{i,\sigma }^{*}({\bf r})h_{\sigma }^{{\rm KS}}({\bf r}){{\xi }_{j,\sigma }}({\bf r})\end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray}&&h_{\sigma }^{{\rm KS}}({\bf r})\equiv -\frac{{{\Delta }_{{\bf r}}}}{2\;m}+V_{\sigma }^{{\rm KS}}({\bf r}),\end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray}&&V_{\sigma }^{{\rm KS}}({\bf r})\equiv {{\kappa }_{\sigma }}({\bf r})=U({\bf r})+{{V}_{{\rm Har}}}({\bf r})+{{v}_{{\rm sp},{\rm xc},\sigma }}({\bf r}).\end{eqnarray} \tag{ 42 }$$

Here, we defined the 'Kohn–Sham potential' $V_{\sigma }^{{\rm KS}}({\bf r})$ that, in our derivation, is identical to the Lagrange parameter ${{\kappa }_{\sigma }}({\bf r})$.

The remaining task is to find the basis in which the Kohn–Sham matrix $\tilde{T}_{\sigma }^{{\rm KS}}$ is diagonal, see appendix A.3.

A better description of transition metals and their compounds can be achieved by supplementing the single-particle reference system resulting from ${{\hat{H}}_{{\rm kin}}}$ in section 2.3 by a multi-band Hubbard interaction. Then, our multi-band reference system follows from

$$\begin{eqnarray}&&{{\hat{H}}_{{\rm H}}}={{\hat{H}}_{{\rm kin}}}+{{\hat{V}}_{{\rm loc}}}-{{\hat{V}}_{{\rm dc}}},\end{eqnarray} \tag{ 43 }$$

where ${{\hat{V}}_{{\rm loc}}}$ describes local interactions between electrons in Wannier orbitals on the same site R. The local single-particle operator ${{\hat{V}}_{{\rm dc}}}$ accounts for the double counting of their interactions in the Hubbard term ${{\hat{V}}_{{\rm loc}}}$ and in the single-particle exchange-correlation energy ${{E}_{{\rm sp},{\rm xc}}}$. We assume that ${{\hat{V}}_{{\rm loc}}}$ and ${{\hat{V}}_{{\rm dc}}}$ do not depend on the densities ${{n}_{\sigma }}({\bf r})$ explicitly.

For the local interaction we set

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{{\hat{V}}}_{{\rm loc}}} & = & \mathop{\sum }\limits_{{\bf R}}{{{\hat{V}}}_{{\rm loc}}}({\bf R}), \\ {{{\hat{V}}}_{{\rm loc}}}({\bf R}) & = & \frac{1}{2}\mathop{\sum }\limits_{({{c}_{1}},{{\sigma }_{1}}),\ldots ,({{c}_{4}},{{\sigma }_{4}})}U_{({{c}_{3}},{{\sigma }_{3}}),({{c}_{4}},{{\sigma }_{4}})}^{({{c}_{1}},{{\sigma }_{1}}),({{c}_{2}},{{\sigma }_{2}})}\;\hat{c}_{{\bf R},{{c}_{1}},{{\sigma }_{1}}}^{\dagger }\hat{c}_{{\bf R},{{c}_{2}},{{\sigma }_{2}}}^{\dagger }\hat{c}_{{\bf R},{{c}_{3}},{{\sigma }_{3}}}^{{}}\hat{c}_{{\bf R},{{c}_{4}},{{\sigma }_{4}}}^{{}} \\ \end{array}.\end{eqnarray} \tag{ 44 }$$

Note that only electrons in the small subset of correlated orbitals (index c) experience the two-particle interaction ${{\hat{V}}_{{\rm loc}}}$: when there are two electrons in the Wannier orbitals ${{\phi }_{{\bf R},{{c}_{3}},{{\sigma }_{3}}}}({\bf r})$ and ${{\phi }_{{\bf R},{{c}_{4}},{{\sigma }_{4}}}}({\bf r})$ centered around the lattice site R, they are scattered into the orbitals ${{\phi }_{{\bf R},{{c}_{1}},{{\sigma }_{1}}}}({\bf r})$ and ${{\phi }_{{\bf R},{{c}_{2}},{{\sigma }_{2}}}}({\bf 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_{\cdots }^{\cdots }$ 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 section 5 and appendix C. Due to screening, the effective Hubbard interaction U is smaller than its bare, atomic value. In general, U and J are chosen to obtain good agreement with experiment, see section 5.

According to Levyʼs constrained search, we must find the minimum of the functional

$$\begin{eqnarray}&&{{F}_{{\rm H}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ |{{\Psi }^{(n)}}\rangle \right\} \right]=\langle {{\Psi }^{(n)}}\mid {{\hat{H}}_{{\rm H}}}\mid {{\Psi }^{(n)}}\rangle \end{eqnarray} \tag{ 45 }$$

in the subset of normalized states with given (physical) density ${{n}_{\sigma }}({\bf r})$, see equation (4). The minimum of ${{F}_{{\rm H}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ |{{\Psi }^{(n)}}\rangle \right\} \right]$ over the states $|{{\Psi }^{(n)}}\rangle$ is the ground state $|\Psi _{{\rm H},0}^{(n)}\rangle$ of the Hamiltonian ${{\hat{H}}_{{\rm H}}}$ for fixed densities ${{n}_{\sigma }}({\bf r})$. In analogy to section 2.3, we define the Hubbard density functional

$$\begin{eqnarray}\begin{array}{rcl} {{D}_{{\rm H}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] & = & {{K}_{{\rm H}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+U\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{V}_{{\rm Har}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] \\ {} & {} & +{{V}_{{\rm loc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]-{{V}_{{\rm dc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{E}_{{\rm H},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right], \\ \end{array}\end{eqnarray} \tag{ 46 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} {{K}_{{\rm H}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] & = & \langle \Psi _{{\rm H},0}^{(n)}|{{{\hat{H}}}_{{\rm kin}}}|\Psi _{{\rm H},0}^{(n)}\rangle , \\ {{V}_{{\rm loc}/{\rm dc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] & = & \langle \Psi _{{\rm H},0}^{(n)}|{{{\hat{V}}}_{{\rm loc}/{\rm dc}}}|\Psi _{{\rm H},0}^{(n)}\rangle , \\ \end{array}\end{eqnarray} \tag{ 47 }$$

and ${{E}_{{\rm H},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]$ is the exchange-correlation energy for ${{\hat{H}}_{{\rm H}}}$. As in section 2.3, the Hubbard density functional agrees with the exact density functional if we choose

$$\begin{eqnarray}\begin{array}{rcl} {{E}_{{\rm H},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] & = & K\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]-{{K}_{{\rm H}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] \\ {} & {} & +{{E}_{{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]-\left( {{V}_{{\rm loc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]-{{V}_{{\rm dc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right] \right). \\ \end{array}\end{eqnarray} \tag{ 48 }$$

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}_{\sigma }}({\bf r})$.

When we directly apply the Ritz principle, we have to minimize the energy functional $E\equiv E\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ |\Psi \rangle \right\} \right]$

$$\begin{eqnarray}&&E=\langle \Psi |{{\hat{H}}_{{\rm H}}}|\Psi \rangle +U\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{V}_{{\rm Har}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{E}_{{\rm H},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right].\end{eqnarray} \tag{ 49 }$$

We include the constraints equation (4) and the normalization condition using the Lagrange parameters ${{\kappa }_{\sigma }}({\bf r})$ and E₀ in the functional ${{G}_{{\rm H}}}\equiv {{G}_{{\rm H}}}\left[ \left\{ |\Psi \rangle \right\},\left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ {{\kappa }_{\sigma }}({\bf r}) \right\},{{E}_{0}} \right]$

$$\begin{eqnarray}\begin{array}{rcl} {{G}_{{\rm H}}} & = & E\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ |\Psi \rangle \right\} \right]-{{E}_{0}}\left( \langle \Psi |\Psi \rangle -1 \right) \\ {} & {} & -\mathop{\sum }\limits_{\sigma }\int {\rm d}{\bf r}{{\kappa }_{\sigma }}({\bf r})\left( {{n}_{\sigma }}({\bf r})-\langle \Psi |\hat{\Psi }_{\sigma }^{\dagger }({\bf r})\hat{\Psi }_{\sigma }^{{}}({\bf r})|\Psi \rangle \right). \\ \end{array}\end{eqnarray} \tag{ 50 }$$

As in section 2.4, see equation (35) and (42), the variation of ${{G}_{{\rm H}}}$ with respect to ${{n}_{\sigma }}({\bf r})$ gives the single-particle potential

$$\begin{eqnarray}\begin{array}{rcl} V_{\sigma }^{{\rm H}}({\bf r}) & \equiv & U({\bf r})+{{V}_{{\rm Har}}}({\bf r})+{{v}_{{\rm H},{\rm xc},\sigma }}({\bf r}), \\ {{v}_{{\rm H},{\rm xc},\sigma }}({\bf r}) & \equiv & {{\left. \frac{\partial {{E}_{{\rm H},{\rm xc}}}\left[ \left\{ {{n}_{\sigma ^{\prime} }}({\bf r}^{\prime} ) \right\} \right]}{\partial {{n}_{\sigma }}({\bf r})} \right|}_{{{n}_{\sigma }}({\bf r})=n_{\sigma }^{0}({\bf r})}}. \\ \end{array}\end{eqnarray} \tag{ 51 }$$

The Hubbard-model approach is based on the idea that typical approximations for the exchange-correlation energy, e.g., the LDA, are suitable for the Hubbard model,

$$\begin{eqnarray}&&{{E}_{{\rm H},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]\approx {{E}_{{\rm LDA},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]\ .\end{eqnarray} \tag{ 52 }$$

Indeed, as seen from equation (48), in the Hubbard exchange-correlation energy ${{E}_{{\rm H},{\rm xc}}}$ the exchange-correlation contributions in the exact ${{E}_{{\rm xc}}}$ are reduced by the Hubbard term ${{V}_{{\rm loc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]-{{V}_{{\rm dc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]$, reflecting a more elaborate treatment of local correlations.

The minimization of (49) with respect to $|\Psi \rangle$ constitutes an unsolvable many-particle problem. The ground state $|{{\Psi }_{0}}\rangle$ is the solution of the many-particle Schrödinger equation with energy E₀,

$$\begin{eqnarray}&&\left( {{{\hat{H}}}_{0}}+{{{\hat{V}}}_{{\rm loc}}}-{{{\hat{V}}}_{{\rm dc}}} \right)|{{\Psi }_{0}}\rangle ={{E}_{0}}|{{\Psi }_{0}}\rangle \end{eqnarray} \tag{ 53 }$$

with the single-particle Hamiltonian

$$\begin{eqnarray}&&{{\hat{H}}_{0}}=\mathop{\sum }\limits_{\sigma }\int {\rm d}{\bf r}\hat{\Psi }_{\sigma }^{\dagger }({\bf r})\left( -\frac{{{\Delta }_{{\bf r}}}}{2\;m}+U({\bf r})+{{V}_{{\rm Har}}}({\bf r})+{{v}_{{\rm H},{\rm xc},\sigma }}({\bf r}) \right)\hat{\Psi }_{\sigma }^{{}}({\bf r}).\end{eqnarray} \tag{ 54 }$$

The Schrödinger equation (53) can be used as starting point for further approximations, for example, the DMFT. In the following we will address the functional in equation (49) directly.

In order to formulate the Gutzwiller variational ground state [4, 8], we consider the local (atomic) states $|\Gamma {{\rangle }_{{\bf R}}}$ that are built from the correlated orbitals. The local Hamiltonians take the form

$$\begin{eqnarray}&&{{\hat{V}}_{{\rm loc}/{\rm dc}}}({\bf R})=\mathop{\sum }\limits_{\Gamma ,\Gamma ^{\prime} }E_{\Gamma ,\Gamma ^{\prime} }^{{\rm loc}/{\rm dc}}({\bf R})|\Gamma {{\rangle }_{{\bf R}}}_{{\bf R}}\langle \Gamma ^{\prime} |=\mathop{\sum }\limits_{\Gamma ,\Gamma ^{\prime} }E_{\Gamma ,\Gamma ^{\prime} }^{{\rm loc}/{\rm dc}}({\bf R}){{\hat{m}}_{{\bf R};\Gamma ,\Gamma ^{\prime} }},\end{eqnarray} \tag{ 55 }$$

where $|\Gamma {{\rangle }_{{\bf R}}}$ contains $|{{\Gamma }_{{\bf R}}}|$ electrons. Here, we introduced

$$\begin{eqnarray}&&E_{\Gamma ,\Gamma ^{\prime} }^{{\rm loc}/{\rm dc}}({\bf R}){{=}_{{\bf R}}}{{\langle \Gamma |{{\hat{V}}_{{\rm loc}/{\rm dc}}}({\bf R})|\Gamma ^{\prime} \rangle }_{{\bf R}}}\end{eqnarray} \tag{ 56 }$$

and the local many-particle operators ${{\hat{m}}_{{\bf R};\Gamma ,\Gamma ^{\prime} }}=|\Gamma {{\rangle }_{{\bf R}}}_{{\bf R}}\langle \Gamma ^{\prime} |$.

The Gutzwiller correlator and the Gutzwiller variational states are defined as

$$\begin{eqnarray}&&{{\hat{P}}_{{\rm G}}}=\mathop{\prod }\limits_{{\bf R}}\mathop{\sum }\limits_{\Gamma ,\Gamma ^{\prime} }{{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R}){{\hat{m}}_{{\bf R};\Gamma ,\Gamma ^{\prime} }},\qquad |{{\Psi }_{{\rm G}}}\rangle ={{\hat{P}}_{{\rm G}}}|\Phi \rangle .\end{eqnarray} \tag{ 57 }$$

Here, $|\Phi \rangle$ is a single-particle product state, and ${{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})$ defines the matrix $\tilde{\lambda }({\bf R})$ of, in general, complex variational parameters.

We evaluate the energy functional (49) in the restricted subset of Gutzwiller variational states,

$$\begin{eqnarray}\begin{array}{rcl} E\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ |{{\Psi }_{{\rm G}}}\rangle \right\} \right] & = & \mathop{\sum }\limits_{{\bf R},b,{\bf R}^{\prime} ,{{b}^{\prime }},\sigma }{{T}_{({\bf R},b),({{{\bf R}}^{\prime }},{{b}^{\prime }});\sigma }}\rho _{({\bf R}^{\prime} ,b^{\prime} ),({\bf R},b);\sigma }^{{\rm G}}+V_{{\rm loc}}^{{\rm G}}-V_{{\rm dc}}^{{\rm G}} \\ {} & {} & +U\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{V}_{{\rm Har}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]+{{E}_{{\rm H},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right], \\ V_{{\rm loc}/{\rm dc}}^{{\rm G}} & = & \mathop{\sum }\limits_{{\bf R}}\mathop{\sum }\limits_{\Gamma ,\Gamma ^{\prime} }E_{\Gamma ,\Gamma ^{\prime} }^{{\rm loc}/{\rm dc}}({\bf R})m_{{\bf R};\Gamma ,\Gamma ^{\prime} }^{{\rm G}}. \\ \end{array}\end{eqnarray} \tag{ 58 }$$

Note that we work with the orbital Wannier basis, see appendix A.3,

$$\begin{eqnarray}&&{{T}_{({\bf R},b),({{{\bf R}}^{\prime }},{{b}^{\prime }});\sigma }}=\int {\rm d}{\bf r}\phi _{{\bf R},b,\sigma }^{*}({\bf r})\left( -\frac{{{\Delta }_{{\bf r}}}}{2\;m} \right){{\phi }_{{{{\bf R}}^{\prime }},{{b}^{\prime }},\sigma }}({\bf r}).\end{eqnarray} \tag{ 59 }$$

The elements of the Gutzwiller-correlated single-particle density matrix are

$$\begin{eqnarray}&&\rho _{({\bf R}^{\prime} ,b^{\prime} ),({\bf R},b);\sigma }^{{\rm G}}=\frac{\langle {{\Psi }_{{\rm G}}}|\hat{c}_{{\bf R},b,\sigma }^{\dagger }\hat{c}_{{\bf R}^{\prime} ,b^{\prime} ,\sigma }^{{}}|{{\Psi }_{{\rm G}}}\rangle }{\langle {{\Psi }_{{\rm G}}}|{{\Psi }_{{\rm G}}}\rangle }=\frac{\langle \Phi |\hat{P}_{{\rm G}}^{\dagger }\hat{c}_{{\bf R},b,\sigma }^{\dagger }\hat{c}_{{\bf R}^{\prime} ,b^{\prime} ,\sigma }^{{}}\hat{P}_{{\rm G}}^{{}}|\Phi \rangle }{\langle \Phi |\hat{P}_{{\rm G}}^{\dagger }\hat{P}_{{\rm G}}^{{}}|\Phi \rangle },\end{eqnarray} \tag{ 60 }$$

and the densities become

$$\begin{eqnarray}&&{{n}_{\sigma }}({\bf r})=\mathop{\sum }\limits_{{\bf R},b,{\bf R}^{\prime} ,b^{\prime} }\phi _{{\bf R},b,\sigma }^{*}({\bf r}){{\phi }_{{\bf R}^{\prime} ,b^{\prime} ,\sigma }}({\bf r})\rho _{({\bf R}^{\prime} ,b^{\prime} ),({\bf R},b);\sigma }^{{\rm G}}.\end{eqnarray} \tag{ 61 }$$

The expectation values for the atomic operators are given by

$$\begin{eqnarray}&&m_{{\bf R};\Gamma ,\Gamma ^{\prime} }^{{\rm G}}=\frac{\langle {{\Psi }_{{\rm G}}}|{{{\hat{m}}}_{{\bf R};\Gamma ,\Gamma ^{\prime} }}|{{\Psi }_{{\rm G}}}\rangle }{\langle {{\Psi }_{{\rm G}}}|{{\Psi }_{{\rm G}}}\rangle }=\frac{\langle \Phi |\hat{P}_{{\rm G}}^{\dagger }{{{\hat{m}}}_{{\bf R};\Gamma ,\Gamma ^{\prime} }}\hat{P}_{{\rm G}}^{{}}|\Phi \rangle }{\langle \Phi |\hat{P}_{{\rm G}}^{\dagger }\hat{P}_{{\rm G}}^{{}}|\Phi \rangle }.\end{eqnarray} \tag{ 62 }$$

The diagrammatic evaluation of $\rho _{({\bf R}^{\prime} ,b^{\prime} ),({\bf R},b);\sigma }^{{\rm G}}$ and of $m_{{\bf R};\Gamma ,\Gamma ^{\prime} }^{{\rm G}}$ shows that these quantities are functionals of the non-interacting single-particle density matrices $\tilde{\rho }$, see equation (30), and of the variational parameters ${{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})$. Moreover, it turns out that the local, non-interacting single-particle density matrix $\tilde{C}({\bf R})$ with the elements

$$\begin{eqnarray}&&{{C}_{b,b^{\prime} ;\sigma }}({\bf R})\equiv {{\rho }_{({\bf R},b),({\bf R},{{b}^{\prime }});\sigma }}\end{eqnarray} \tag{ 63 }$$

plays a prominent role in the Gutzwiller energy functional, in particular, for infinite lattice coordination number. Therefore, we may write

$$\begin{eqnarray}&&E\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ |{{\Psi }_{{\rm G}}}\rangle \right\} \right]\equiv {{E}^{{\rm G}}}\left[ \tilde{\rho },\left\{ \tilde{\lambda }({\bf R}) \right\},\left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ \tilde{C}({\bf R}) \right\} \right].\end{eqnarray} \tag{ 64 }$$

In the Lagrange functional we shall impose the relation (63) with the help of the Hermitian Lagrange parameter matrix $\tilde{\eta }$ with entries ${{\eta }_{({\bf R},b),({\bf R},{{b}^{\prime }});\sigma }}$. Lastly, for the analytical evaluation of equation (64) it is helpful to impose a set of (real-valued) local constraints ($l=1,2,\ldots ,{{N}_{{\rm con}}}$)

$$\begin{eqnarray}&&{{g}_{l,{\bf R}}}\left[ \tilde{\lambda }({\bf R}),\tilde{C}({\bf R}) \right]=0,\end{eqnarray} \tag{ 65 }$$

which we implement with real Lagrange parameters ${{\Lambda }_{l}}({\bf R})$; for explicit expressions, see equations (74) and (75).

In the following, we abbreviate $i=({\bf R},b)$ and $j=({\bf R}^{\prime} ,b^{\prime} )$. Consequently, in analogy with section 2.4, we address

$$\begin{eqnarray}G_{{\rm DFT}}^{{\rm G}}\equiv G_{{\rm DFT}}^{{\rm G}}\left[ \begin{array}{ccccccccccccccc} \tilde{\rho }, & \left\{ {{n}_{\sigma }}({\bf r}) \right\}, & \left\{ \tilde{C}({\bf R}) \right\}, & \left\{ \tilde{\lambda }({\bf R}) \right\} \\ \tilde{\Omega }, & \left\{ {{\kappa }_{\sigma }}({\bf r}) \right\}, & \left\{ \tilde{\eta }({\bf R}) \right\}, & \left\{ {{\Lambda }_{l}}({\bf R}) \right\} \\ \end{array} \right]\end{eqnarray} \tag{ 66 }$$

as our Lagrange functional,

$$\begin{eqnarray}\begin{array}{rcl} G_{{\rm DFT}}^{{\rm G}} & = & {{E}^{{\rm G}}}\left[ \tilde{\rho },\left\{ \tilde{\lambda }({\bf R}) \right\},\left\{ {{n}_{\sigma }}({\bf r}) \right\},\left\{ \tilde{C}({\bf R}) \right\} \right]-\mathop{\sum }\limits_{l,m,\sigma }{{\Omega }_{l,m;\sigma }}{{\left( \tilde{\rho }\cdot \tilde{\rho }-\tilde{\rho } \right)}_{m,l;\sigma }} \\ {} & {} & -\mathop{\sum }\limits_{\sigma }\int {\rm d}{\bf r}{{\kappa }_{\sigma }}({\bf r})\left( {{n}_{\sigma }}({\bf r})-\mathop{\sum }\limits_{i,j}\phi _{i,\sigma }^{*}({\bf r}){{\phi }_{j,\sigma }}({\bf r})\rho _{j,i;\sigma }^{{\rm G}} \right) \\ \end{array}\end{eqnarray} \tag{ 67 }$$

$$\begin{eqnarray*}&&+\mathop{\sum }\limits_{l,{\bf R}}{{\Lambda }_{l}}({\bf R}){{g}_{l,{\bf R}}}-\mathop{\sum }\limits_{{\bf R},b,{{b}^{\prime }},\sigma }{{\eta }_{b,{{b}^{\prime }};\sigma }}({\bf R})\left( {{C}_{{{b}^{\prime }},b;\sigma }}({\bf R})-{{\rho }_{({\bf R},{{b}^{\prime }}),({\bf R},b);\sigma }} \right),\end{eqnarray*}$$

see equation (34). Here, we took the condition (61) into account using Lagrange parameters ${{\kappa }_{\sigma }}({\bf r})$ because the external potential, the Hartree term and the exchange-correlation potential in equation (58) depend on the densities.

The functional G_DFT^G in equation (67) has to be minimized with respect to ${{n}_{\sigma }}({\bf r})$, $\tilde{C}({\bf R})$, $\tilde{\lambda }({\bf R})$, and $\tilde{\rho }$. The variation with respect to the Lagrange parameters ${{\kappa }_{\sigma }}({\bf r})$, $\tilde{\eta }({\bf R})$, ${{\Lambda }_{l}}({\bf R})$, and $\tilde{\Omega }$ gives the constraints (61), (63), (65), and (33), respectively.

• (i)  
  As in the derivation of the exact Schrödinger equation (53), the variation of G_DFT^G with respect to ${{n}_{\sigma }}({\bf r})$ generates the single-particle potential,

  $$\begin{eqnarray}&&{{\kappa }_{\sigma }}({\bf r})=V_{\sigma }^{{\rm H}}({\bf r}),\end{eqnarray} \tag{ 68 }$$

  see equations (42) and (51).
• (ii)  
  The minimization with respect to $\tilde{C}({\bf R})$ gives

  $$\begin{eqnarray}\begin{array}{rcl} {{\eta }_{b,b^{\prime} ;\sigma }}({\bf R}) & = & \frac{\partial {{E}^{{\rm G}}}}{\partial {{C}_{b^{\prime} ,b;\sigma }}({\bf R})}+\mathop{\sum }\limits_{l}{{\Lambda }_{l}}({\bf R})\frac{\partial {{g}_{l,{\bf R}}}}{\partial {{C}_{b^{\prime} ,b;\sigma }}({\bf R})} \\ {} & {} & +\mathop{\sum }\limits_{i,j,\sigma ^{\prime} }\int {\rm d}{\bf r}V_{\sigma ^{\prime} }^{{\rm H}}({\bf r})\phi _{i,\sigma ^{\prime} }^{*}({\bf r}){{\phi }_{j,\sigma ^{\prime} }}({\bf r})\frac{\partial \rho _{j,i;\sigma ^{\prime} }^{{\rm G}}}{\partial {{C}_{b^{\prime} ,b;\sigma }}({\bf R})}. \\ \end{array}\end{eqnarray} \tag{ 69 }$$
• (iii)  
  The minimization with respect to the Gutzwiller correlation parameters $\tilde{\lambda }({\bf R})$ results in

  $$\begin{eqnarray}\begin{array}{rcl} 0 & = & \frac{\partial {{E}^{{\rm G}}}}{\partial {{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})}+\mathop{\sum }\limits_{l,m,\sigma }\int {\rm d}{\bf r}V_{\sigma }^{{\rm H}}({\bf r})\phi _{l,\sigma }^{*}({\bf r}){{\phi }_{m,\sigma }}({\bf r})\frac{\partial \rho _{m,l,\sigma }^{{\rm G}}}{\partial {{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})} \\ {} & {} & +\mathop{\sum }\limits_{l}{{\Lambda }_{l}}({\bf R})\frac{\partial {{g}_{l,{\bf R}}}}{\partial {{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})} \\ \end{array}\end{eqnarray} \tag{ 70 }$$

  $$\begin{eqnarray}\begin{array}{clr} {} & = & \mathop{\sum }\limits_{l,m,\sigma }h_{l,m;\sigma }^{0}\frac{\partial \rho _{m,l,\sigma }^{{\rm G}}}{\partial {{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})}+\frac{\partial \left( V_{{\rm loc}}^{{\rm G}}-V_{{\rm dc}}^{{\rm G}} \right)}{\partial {{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})}+\mathop{\sum }\limits_{l}{{\Lambda }_{l}}({\bf R})\frac{\partial {{g}_{l,{\bf R}}}}{\partial {{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})}, \\ h_{l,m;\sigma }^{0} & \equiv & \int {\rm d}{\bf r}\phi _{l,\sigma }^{*}({\bf r})\left( -\frac{{{\Delta }_{{\bf r}}}}{2\;m}+U({\bf r})+{{V}_{{\rm Har}}}({\bf r})+{{v}_{{\rm H},{\rm xc},\sigma }}({\bf r}) \right){{\phi }_{m,\sigma }}({\bf r}) \\ \end{array}\end{eqnarray} \tag{ 71 }$$

  for all ${{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})$. Note that, in the case of complex Gutzwiller parameters, we also have to minimize with respect to ${{({{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R}))}^{*}}$. Using these equations we may calculate the Lagrange parameters ${{\Lambda }_{l}}({\bf R})$ that are needed in equation (69).
• (iv)  
  The minimization of G_DFT^G with respect to $\tilde{\rho }$ generates the Landau–Gutzwiller quasi-particle Hamiltonian, see appendix A.2,

  $$\begin{eqnarray}&&\hat{H}_{{\rm qp}}^{{\rm G}}=\mathop{\sum }\limits_{i,j,\sigma }h_{i,j;\sigma }^{{\rm G}}\hat{c}_{i,\sigma }^{\dagger }\hat{c}_{j,\sigma }^{{}}\end{eqnarray} \tag{ 72 }$$

  with the entries

  $$\begin{eqnarray}\begin{array}{rcl} h_{i,j;\sigma }^{{\rm G}} & = & \frac{\partial {{E}^{{\rm G}}}}{\partial {{\rho }_{j,i;\sigma }}}+\mathop{\sum }\limits_{l,m,\sigma ^{\prime} }\int {\rm d}{\bf r}V_{\sigma ^{\prime} }^{{\rm H}}({\bf r})\phi _{l,\sigma ^{\prime} }^{*}({\bf r}){{\phi }_{m,\sigma ^{\prime} }}({\bf r})\frac{\partial \rho _{m,l,\sigma ^{\prime} }^{{\rm G}}}{\partial {{\rho }_{j,i;\sigma }}} \\ {} & {} & +\mathop{\sum }\limits_{{\bf R},b,b^{\prime} ,\sigma ^{\prime} }{{\eta }_{b,{{b}^{\prime }};\sigma ^{\prime} }}({\bf R})\frac{\partial {{\rho }_{({\bf R},{{b}^{\prime }}),({\bf R},b);\sigma ^{\prime} }}}{\partial {{\rho }_{j,i;\sigma }}} \\ {} & = & \mathop{\sum }\limits_{l,m,\sigma ^{\prime} }h_{l,m;\sigma ^{\prime} }^{0}\frac{\partial \rho _{m,l,\sigma ^{\prime} }^{{\rm G}}}{\partial {{\rho }_{j,i;\sigma }}}+\frac{\partial \left( V_{{\rm loc}}^{{\rm G}}-V_{{\rm dc}}^{{\rm G}} \right)}{\partial {{\rho }_{j,i;\sigma }}} \\ {} & {} & +\mathop{\sum }\limits_{{\bf R},b,b^{\prime} }{{\delta }_{j,({\bf R},{{b}^{\prime }})}}{{\delta }_{i,({\bf R},b)}}{{\eta }_{b,b^{\prime} ;\sigma }}({\bf R}), \\ \end{array}\end{eqnarray} \tag{ 73 }$$

  where we used equations (58) and (70). The single-particle state $|\Phi \rangle$ is the ground state of the Hamiltonian (72) from which the single-particle density matrix $\tilde{\rho }$ follows.

The minimization problem outlined in steps (i)–(iv) requires the evaluation of the energy ${{E}^{{\rm G}}}$ in equation (58). In particular, the correlated single-particle density matrix ${{\tilde{\rho }}^{{\rm G}}}$, equation (60), 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–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.

As shown in [8, 22] it is convenient for the evaluation of Gutzwiller wave functions to impose the following (local) constraints

$$\begin{eqnarray}&&\mathop{\sum }\limits_{\Gamma ,{{\Gamma }_{1}},{{\Gamma }_{2}}}\lambda _{\Gamma ,{{\Gamma }_{1}}}^{*}({\bf R}){{\lambda }_{\Gamma ,{{\Gamma }_{2}}}}({\bf R}){{\langle {{\hat{m}}_{{\bf R};{{\Gamma }_{1}},{{\Gamma }_{2}}}}\rangle }_{\Phi }}=1\end{eqnarray} \tag{ 74 }$$

and

$$\begin{eqnarray}&&\mathop{\sum }\limits_{\Gamma ,{{\Gamma }_{1}},{{\Gamma }_{2}}}\lambda _{\Gamma ,{{\Gamma }_{1}}}^{*}({\bf R}){{\lambda }_{\Gamma ,{{\Gamma }_{2}}}}({\bf R}){{\langle {{\hat{m}}_{{\bf R};{{\Gamma }_{1}},{{\Gamma }_{2}}}}\hat{c}_{{\bf R},b,\sigma }^{\dagger }\hat{c}_{{\bf R},b^{\prime} ,\sigma }^{{}}\rangle }_{\Phi }}={{\langle \hat{c}_{{\bf R},b,\sigma }^{\dagger }\hat{c}_{{\bf R},b^{\prime} ,\sigma }^{{}}\rangle }_{\Phi }},\end{eqnarray} \tag{ 75 }$$

where we abbreviated ${{\langle \hat{A}\rangle }_{\Phi }}\equiv \langle \Phi |\hat{A}|\Phi \rangle$. Note that, for complex constraints, the index l in (65) labels real and imaginary parts separately.

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 $\tilde{\lambda }({\bf R})$ and the local density matrix $\tilde{C}({\bf R})$ of the correlated bands in $|\Phi \rangle$,

$$\begin{eqnarray}&&V_{{\rm loc}/{\rm dc}}^{{\rm G}}=\mathop{\sum }\limits_{{\bf R}}\mathop{\sum }\limits_{{{\Gamma }_{1}},\ldots ,{{\Gamma }_{4}}}\lambda _{{{\Gamma }_{2}},{{\Gamma }_{1}}}^{*}({\bf R})E_{{{\Gamma }_{2}},{{\Gamma }_{3}}}^{{\rm loc}/{\rm dc}}({\bf R}){{\lambda }_{{{\Gamma }_{3}},{{\Gamma }_{4}}}}({\bf R}){{\langle {{\hat{m}}_{{\bf R};{{\Gamma }_{1}},{{\Gamma }_{4}}}}\rangle }_{\Phi }}.\end{eqnarray} \tag{ 76 }$$

The remaining expectation values ${{\langle {{\hat{m}}_{{\bf R};{{\Gamma }_{1}},{{\Gamma }_{4}}}}\rangle }_{\Phi }}$ are evaluated using Wickʼs theorem. Explicit expressions are given in [8, 23].

The local part of the correlated single-particle density matrix is given by

$$\begin{eqnarray}\begin{array}{rcl} \rho _{({\bf R},b^{\prime} ),({\bf R},b);\sigma }^{{\rm G}} & = & \mathop{\sum }\limits_{{{\Gamma }_{1}},\ldots ,{{\Gamma }_{4}}}\lambda _{{{\Gamma }_{2}},{{\Gamma }_{1}}}^{*}({\bf R}){{\lambda }_{{{\Gamma }_{3}},{{\Gamma }_{4}}}}({\bf R}){{\langle {{{\hat{m}}}_{{\bf R};{{\Gamma }_{1}},{{\Gamma }_{2}}}}\hat{c}_{{\bf R},b,\sigma }^{\dagger }\hat{c}_{{\bf R},b^{\prime} ,\sigma }^{{}}{{{\hat{m}}}_{{\bf R};{{\Gamma }_{3}},{{\Gamma }_{4}}}}\rangle }_{\Phi }} \\ {} & \equiv & C_{b^{\prime} ,b;\sigma }^{{\rm G}}({\bf R}). \\ \end{array}\end{eqnarray} \tag{ 77 }$$

It can be evaluated using Wickʼs theorem. As can be seen from equation (77), it is a function of the variational parameters ${{\lambda }_{\Gamma ,\Gamma ^{\prime} }}({\bf R})$ and of the local non-interacting single-particle density matrix $\tilde{C}({\bf R})$.

For ${\bf R}\ne {\bf R}^{\prime}$, we have for the correlated single-particle density matrix

$$\begin{eqnarray}&&\rho _{({\bf R}^{\prime} ,{{b}^{\prime }}),({\bf R},b);\sigma }^{{\rm G}}=\mathop{\sum }\limits_{a,a^{\prime} }q_{b,\sigma }^{a,\sigma }({\bf R}){{\left( q_{{{b}^{\prime }},\sigma }^{{{a}^{\prime }},\sigma }({{{\bf R}}^{\prime }}) \right)}^{*}}{{\rho }_{({{{\bf R}}^{\prime }},{{a}^{\prime }}),({\bf R},a);\sigma }}\end{eqnarray} \tag{ 78 }$$

with the orbital-dependent renormalization factors $q_{b,\sigma }^{a,\sigma }({\bf R})$ for the electron transfer between different sites. Explicit expressions in terms of the variational parameters $\tilde{\lambda }({\bf R})$ and of the local non-interacting single-particle density matrix $\tilde{C}({\bf R})$ are given in [8, 23].

When we apply the general expressions (73), the matrix elements of the quasi-particle Hamiltonian are obtained as

$$\begin{eqnarray}&&h_{b,b^{\prime} ;\sigma }^{{\rm G}}({\bf k})={{\eta }_{b,b^{\prime} ;\sigma }}+\mathop{\sum }\limits_{a,a^{\prime} ,\sigma ^{\prime} }h_{a,a^{\prime} ;\sigma ^{\prime} }^{0}({\bf k})\frac{\partial \rho _{a^{\prime} ,a;\sigma ^{\prime} }^{{\rm G}}({\bf k})}{\partial {{\rho }_{b^{\prime} ,b;\sigma }}({\bf k})},\end{eqnarray} \tag{ 80 }$$

where we have from equation (70)

$$\begin{eqnarray}&&h_{a,a^{\prime} ;\sigma }^{0}({\bf k})=\int {\rm d}{\bf r}\phi _{{\bf k},a,\sigma }^{*}({\bf r})\left( -\frac{{{\Delta }_{{\bf r}}}}{2\;m}+U({\bf r})+{{V}_{{\rm Har}}}({\bf r})+{{v}_{{\rm H},{\rm xc},\sigma }}({\bf r}) \right){{\phi }_{{\bf k},a^{\prime} ,\sigma }}({\bf r}).\end{eqnarray} \tag{ 81 }$$

Moreover,

$$\begin{eqnarray}&&{{\rho }_{b^{\prime} ,b;\sigma }}({\bf k})=\langle \Phi |\hat{c}_{{\bf k},b,\sigma }^{\dagger }\hat{c}_{{\bf k},b^{\prime} ,\sigma }^{{}}|\Phi \rangle \end{eqnarray} \tag{ 82 }$$

are the entries of the single-particle density matrix in the orbital Bloch basis, and $\rho _{b^{\prime} ,b;\sigma }^{{\rm G}}({\bf k})$ denotes the corresponding quantities in the Gutzwiller-correlated state. In the limit of infinite lattice coordination number, we may express $V_{{\rm loc}/{\rm dc}}^{{\rm G}}$ in equation (73) as a function of the Gutzwiller variational parameters ${{\lambda }_{\Gamma ,\Gamma ^{\prime} }}$ and of the local density matrix $\tilde{C}$. Therefore, $V_{{\rm loc}/{\rm dc}}^{{\rm G}}$ are formally independent of the single-particle density matrix $\tilde{\rho }$ so that they do not contribute to $h_{b,b^{\prime} ;\sigma }^{{\rm G}}({\bf k})$. Equation (80) shows that, apart from an overall shift of the orbitals through ${{\eta }_{b,b^{\prime} ;\sigma }}$, we can still work with the matrix elements $h_{a,a^{\prime} ;\sigma ^{\prime} }^{0}({\bf k})$ of the single-particle operator ${{\hat{H}}_{0}}$ that enters the many-particle Schrödinger equation (53).

In the orbital Bloch basis, equations (77) and (78) take the form

$$\begin{eqnarray}&&\rho _{b^{\prime} ,b;\sigma }^{{\rm G}}({\bf k})=\mathop{\sum }\limits_{a,a^{\prime} }q_{b,\sigma }^{a,\sigma }{{\left( q_{b^{\prime} ,\sigma }^{a^{\prime} ,\sigma } \right)}^{*}}\left( {{\rho }_{a^{\prime} ,a;\sigma }}({\bf k})-{{C}_{a^{\prime} ,a;\sigma }} \right)+C_{b^{\prime} ,b;\sigma }^{{\rm G}}.\end{eqnarray} \tag{ 83 }$$

When we insert equation (83) into equation (80) we thus find for the entries of the Gutzwiller quasi-particle Hamiltonian

$$\begin{eqnarray}&&h_{b,b^{\prime} ;\sigma }^{{\rm G}}({\bf k})={{\eta }_{b,b^{\prime} ;\sigma }}+\mathop{\sum }\limits_{a,a^{\prime} }q_{a,\sigma }^{b,\sigma }{{\left( q_{a^{\prime} ,\sigma }^{b^{\prime} ,\sigma } \right)}^{*}}h_{a,a^{\prime} ;\sigma }^{0}({\bf k}).\end{eqnarray} \tag{ 84 }$$

Recall that the local single-particle densities $\tilde{C}$ are treated as independent parameters. Since the q-factors and the correlated local single-particle density matrix $C_{b,b^{\prime} ;\sigma }^{{\rm G}}=\rho _{({\bf R},b),({\bf R},b^{\prime} );\sigma }^{{\rm G}}$ are solely functions of the variational parameters $\tilde{\lambda }$ and of $\tilde{C}$, they are treated as constants when we take the derivative of $\rho _{b^{\prime} ,b;\sigma }^{{\rm G}}({\bf k})$ with respect to ${{\rho }_{b^{\prime} ,b;\sigma }}({\bf k})$ in equation (83).

The unitary matrix $\tilde{F}_{\sigma }^{{\rm G}}({\bf k})$ diagonalizes the Gutzwiller matrix $\tilde{h}_{\sigma }^{{\rm G}}({\bf k})$,

$$\begin{eqnarray}&&\mathop{\sum }\limits_{b,b^{\prime} }{{(F_{b,{\rm n},\sigma }^{{\rm G}}({\bf k}))}^{*}}h_{b,b^{\prime} ;\sigma }^{{\rm G}}({\bf k})F_{b^{\prime} ,{\rm m},\sigma }^{{\rm G}}({\bf k})={{\delta }_{{\rm n},{\rm m}}}\epsilon _{{\rm n},\sigma }^{{\rm G}}({\bf k}),\end{eqnarray} \tag{ 85 }$$

which provides the quasi-particle dispersion $\epsilon _{{\rm n},\sigma }^{{\rm G}}({\bf k})$. We introduce the quasi-particle band operators

$$\begin{eqnarray}&&\hat{g}_{{\bf k},{\rm n},\sigma }^{\dagger }=\mathop{\sum }\limits_{b}F_{b,{\rm n},\sigma }^{{\rm G}}({\bf k})\hat{c}_{{\bf k},b,\sigma }^{\dagger },\qquad \hat{g}_{{\bf k},{\rm n},\sigma }^{{}}=\mathop{\sum }\limits_{b}{{(F_{b,{\rm n},\sigma }^{{\rm G}}({\bf k}))}^{*}}\hat{c}_{{\bf k},b,\sigma }^{{}},\end{eqnarray} \tag{ 86 }$$

in which the Gutzwiller quasi-particle Hamiltonian becomes diagonal,

$$\begin{eqnarray}&&\hat{H}_{{\rm qp}}^{{\rm G}}=\mathop{\sum }\limits_{{\bf k},{\rm n},\sigma }\epsilon _{{\rm n},\sigma }^{{\rm G}}({\bf k})\hat{g}_{{\bf k},{\rm n},\sigma }^{\dagger }\hat{g}_{{\bf k},{\rm n},\sigma }^{{}}.\end{eqnarray} \tag{ 87 }$$

In order to minimize the Gutzwiller density functional G_DFT^G, we must work with the ground state of the quasi-particle Hamiltonian $\hat{H}_{{\rm qp}}^{{\rm G}}$,

$$\begin{eqnarray}&&|{{\Phi }_{0}}\rangle ={{\mathop{\prod }\limits_{{\bf k},{\rm n},\sigma }}^{^{\prime }}}\hat{g}_{{\bf k},{\rm n},\sigma }^{\dagger }|{\rm vac}\rangle ,\end{eqnarray} \tag{ 88 }$$

where the N levels lowest in energy are occupied as indicated by the prime at the product, $\epsilon _{{\rm n},\sigma }^{{\rm G}}({\bf k})\leqslant E_{{\rm F},\sigma }^{{\rm G}}$. Using equation (82) we find

$$\begin{eqnarray}\begin{array}{rcl} \rho _{b^{\prime} ,b;\sigma }^{{\rm opt}}({\bf k}) & = & {{\langle \hat{c}_{{\bf k},b,\sigma }^{\dagger }\hat{c}_{{\bf k},b^{\prime} ,\sigma }^{{}}\rangle }_{{{\Phi }_{0}}}} \\ {} & = & \mathop{\sum }\limits_{{\rm n}}{{\left( F_{b,{\rm n},\sigma }^{{\rm G}}({\bf k}) \right)}^{*}}F_{b^{\prime} ,{\rm n},\sigma }^{{\rm G}}({\bf k}){{\langle \hat{g}_{{\bf k},{\rm n},\sigma }^{\dagger }\hat{g}_{{\bf k},{\rm n},\sigma }^{{}}\rangle }_{{{\Phi }_{0}}}} \\ {} & = & \mathop{\sum }\limits_{{\rm n}}f_{{\bf k},{\rm n},\sigma }^{{\rm G}}{{\left( F_{b,{\rm n},\sigma }^{{\rm G}}({\bf k}) \right)}^{*}}F_{b^{\prime} ,{\rm n},\sigma }^{{\rm G}}({\bf k}), \\ \end{array}\end{eqnarray} \tag{ 89 }$$

where the quasi-particle occupancies

$$\begin{eqnarray}&&f_{{\bf k},{\rm n},\sigma }^{{\rm G}}={{\langle \hat{g}_{{\bf k},{\rm n},\sigma }^{\dagger }\hat{g}_{{\bf k},{\rm n},\sigma }^{{}}\rangle }_{{{\Phi }_{0}}}}=\Theta \left( E_{{\rm F},\sigma }^{{\rm G}}-\epsilon _{{\rm n},\sigma }^{{\rm G}}({\bf k}) \right)\end{eqnarray} \tag{ 90 }$$

are unity for occupied quasi-particle levels up to the Fermi energy $E_{{\rm F},\sigma }^{{\rm G}}$, and zero otherwise. The particle densities follow from equation (89),

$$\begin{eqnarray}&&{{n}_{\sigma }}({\bf r})=\mathop{\sum }\limits_{{\bf k},b,b^{\prime} }\phi _{{\bf k},b,\sigma }^{*}({\bf r}){{\phi }_{{\bf k},b^{\prime} ,\sigma }}({\bf r})\rho _{b^{\prime} ,b;\sigma }^{{\rm G}}({\bf k}),\end{eqnarray} \tag{ 91 }$$

where $\rho _{b^{\prime} ,b;\sigma }^{{\rm G}}({\bf k})$ is given by equation (83). As in DFT, the particle densities must be calculated self-consistently.

In order to determine the band-shift parameters ${{\eta }_{b,b^{\prime} ;\sigma }}$, we must evaluate equation (69) using the Gutzwiller energy functional in the limit of infinite lattice coordination number. We define the kinetic energy of the Gutzwiller quasi-particles as

$$\begin{eqnarray}&&E_{{\rm kin}}^{{\rm G}}=\mathop{\sum }\limits_{{\bf k},{\rm n},\sigma }f_{{\bf k},{\rm n},\sigma }^{{\rm G}}\epsilon _{{\rm n},\sigma }^{{\rm G}}({\bf k})=\mathop{\sum }\limits_{{\bf k},b,b^{\prime} ,\sigma }h_{b,b^{\prime} ;\sigma }^{{\rm G}}({\bf k}){{\rho }_{b^{\prime} ,b;\sigma }}({\bf k}).\end{eqnarray} \tag{ 92 }$$

It is easy to show that

$$\begin{eqnarray}\begin{array}{rcl} \mathop{\sum }\limits_{{\bf k},b,b^{\prime} ,\sigma }h_{b,b^{\prime} ;\sigma }^{0}({\bf k})\rho _{b^{\prime} ,b;\sigma }^{{\rm G}}({\bf k}) & = & E_{{\rm kin}}^{{\rm G}}-\mathop{\sum }\limits_{{\bf k},b,b^{\prime} ,\sigma }{{\eta }_{b,b^{\prime} ;\sigma }}{{\rho }_{b^{\prime} ,b;\sigma }}({\bf k}) \\ {} & {} & +L\mathop{\sum }\limits_{b,b^{\prime} ,\sigma }h_{b,b^{\prime} ;\sigma }^{0}\left( C_{b^{\prime} ,b;\sigma }^{{\rm G}}-{{C}_{b^{\prime} ,b;\sigma }} \right), \\ h_{b,b^{\prime} ;\sigma }^{0} & = & \frac{1}{L}\mathop{\sum }\limits_{{\bf k}}h_{b,b^{\prime} ;\sigma }^{0}({\bf k}). \\ \end{array}\end{eqnarray} \tag{ 93 }$$

Then, equation (69) gives the effective local hybridizations ${{\eta }_{b,b^{\prime} ;\sigma }}$

$$\begin{eqnarray}\begin{array}{rcl} L{{\eta }_{b,b^{\prime} ;\sigma }} & = & \frac{\partial }{\partial {{C}_{b^{\prime} ,b;\sigma }}}\left( L\mathop{\sum }\limits_{l}{{\Lambda }_{l}}{{g}_{l}}+V_{{\rm loc}}^{{\rm G}}-V_{{\rm dc}}^{{\rm G}}+E_{{\rm kin}}^{{\rm G}} \right. \\ {} & {} & \left. +L\mathop{\sum }\limits_{a,a^{\prime} ,\sigma ^{\prime} }\left( C_{a^{\prime} ,a;\sigma }^{{\rm G}}-{{C}_{a^{\prime} ,a;\sigma }} \right)h_{a,a^{\prime} ;\sigma ^{\prime} }^{0} \right). \\ \end{array}\end{eqnarray} \tag{ 94 }$$

Note that the term in the second line in equation (94) often vanishes due to symmetry, e.g., in nickel, because $C_{b,b^{\prime} ;\sigma }^{{\rm G}}={{C}_{b,b^{\prime} ;\sigma }}$.

In many theoretical studies one uses a Hamiltonian with only density-density interactions,

$$\begin{eqnarray}&&\hat{V}_{{\rm loc}}^{{\rm dens}}=\mathop{\sum }\limits_{c,\sigma }U(c,c){{\hat{n}}_{c,\sigma }}{{\hat{n}}_{c,\bar{\sigma }}}+\mathop{\sum }\limits_{c(\ne )c^{\prime} }\mathop{\sum }\limits_{\sigma ,\sigma ^{\prime} }{{\tilde{U}}_{\sigma ,\sigma ^{\prime} }}(c,c^{\prime} ){{\hat{n}}_{c,\sigma }}{{\hat{n}}_{c^{\prime} ,\sigma ^{\prime} }}\;.\end{eqnarray} \tag{ 98 }$$

Here, we introduced $\bar{\uparrow }=\downarrow$ ($\bar{\downarrow }=\uparrow$) and ${{\tilde{U}}_{\sigma ,\sigma ^{\prime} }}(c,c^{\prime} )=U(c,c^{\prime} )-{{\delta }_{\sigma ,\sigma ^{\prime} }}J(c,c^{\prime} )$, where $U(c,c^{\prime} )$ and $J(c,c^{\prime} )$ are the local Hubbard and Hundʼs-rule exchange interactions. An additional and quite common approximation is the use of orbital-independent Coulomb parameters,

$$\begin{eqnarray}&&U(c,c)\equiv U,\quad {\rm and}\quad U(c,c^{\prime} )\equiv U^{\prime} ,\ J(c,c^{\prime} )\equiv J\quad {\rm for}\;c\ne c^{\prime} .\end{eqnarray} \tag{ 99 }$$

For a system of five correlated 3d orbitals in a cubic environment as in nickel, however, the Hamiltonian (98) is incomplete [24]. The full Hamiltonian reads

$$\begin{eqnarray}&&\hat{V}_{{\rm loc}}^{{\rm full}}=\hat{V}_{{\rm loc}}^{{\rm dens}}+\hat{V}_{{\rm loc}}^{{\rm n}.{\rm dens}.},\end{eqnarray} \tag{ 100 }$$

where

$$\begin{eqnarray}\begin{array}{rcl} \hat{V}_{{\rm loc}}^{{\rm n}.{\rm dens}.} & = & \mathop{\sum }\limits_{c(\ne )c^{\prime} }J(c,c^{\prime} )\left( \hat{c}_{c,\uparrow }^{\dagger }\hat{c}_{c,\downarrow }^{\dagger }{{{\hat{c}}}_{c^{\prime} ,\downarrow }}{{{\hat{c}}}_{c^{\prime} ,\uparrow }}+{\rm h}.{\rm c}. \right)+\mathop{\sum }\limits_{c(\ne )c^{\prime} ;\sigma }J(c,c^{\prime} )\hat{c}_{c,\sigma }^{\dagger }\hat{c}_{c^{\prime} ,\bar{\sigma }}^{\dagger }{{{\hat{c}}}_{c,\bar{\sigma }}}{{{\hat{c}}}_{c^{\prime} ,\sigma }} \\ {} & {} & +\left[ \mathop{\sum }\limits_{t;\sigma ,\sigma ^{\prime} }(T(t)-{{\delta }_{\sigma ,\sigma ^{\prime} }}A(t)){{{\hat{n}}}_{t,\sigma }}\hat{c}_{u,\sigma ^{\prime} }^{\dagger }{{{\hat{c}}}_{v,\sigma ^{\prime} }} \right. \\ {} & {} & +\mathop{\sum }\limits_{t,\sigma }A(t)\left( \hat{c}_{t,\sigma }^{\dagger }\hat{c}_{t,\bar{\sigma }}^{\dagger }{{{\hat{c}}}_{u,\bar{\sigma }}}{{{\hat{c}}}_{v,\sigma }}+\hat{c}_{t,\sigma }^{\dagger }\hat{c}_{u,\bar{\sigma }}^{\dagger }{{{\hat{c}}}_{t,\bar{\sigma }}}{{{\hat{c}}}_{v,\sigma }} \right) \\ {} & {} & \left. +\mathop{\sum }\limits_{t(\ne )t^{\prime} (\ne ){{t}^{^{\prime\prime} }}}\mathop{\sum }\limits_{e,\sigma ,\sigma ^{\prime} }S(t,t^{\prime} ;{{t}^{^{\prime\prime} }},e)\hat{c}_{t,\sigma }^{\dagger }\hat{c}_{t^{\prime} ,\sigma ^{\prime} }^{\dagger }{{{\hat{c}}}_{{{t}^{^{\prime\prime} }},\sigma ^{\prime} }}{{{\hat{c}}}_{e,\sigma }}+{\rm h}.{\rm c}. \right]\;. \\ \end{array}\end{eqnarray} \tag{ 101 }$$

Here, $t=\zeta ,\eta ,\xi$ and $e=u,v$ are indices for the three ${{t}_{2\,g}}$ orbitals with symmetries $\zeta =xy$, $\eta =xz$, and $\xi =yz$, and the two e_g orbitals with symmetries $u=3{{z}^{2}}-{{r}^{2}}$ and $v={{x}^{2}}-{{y}^{2}}$, respectively. The parameters A(t), T(t), $S(t,t^{\prime} ;{{t}^{\prime\prime }},e)$ in equation (101) are of the same order of magnitude as the exchange interactions $J(c,c^{\prime} )$ and, hence, there is no a-priori reason to neglect $V_{{\rm loc}}^{{\rm n}.{\rm dens}.}$. Of all the parameters $U(c,c^{\prime} )$, $J(c,c^{\prime} )$, A(t), T(t), $S(t,t^{\prime} ;{{t}^{\prime\prime }},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={{\sum }_{c}}U(c,c)/5=A+4B+3C$, the average Coulomb interaction between electrons in different orbitals, $U^{\prime} ={{\sum }_{c\lt c^{\prime} }}U(c,c^{\prime} )/10=A-B+C$, and the average Hundʼs-rule exchange interaction, $J={{\sum }_{c\lt c^{\prime} }}J(c,c^{\prime} )/10=5B/2+C$ that are related by the symmetry relation $U^{\prime} =U-2J$, see appendix C. Due to this symmetry relation, the three values of U, $U^{\prime}$, 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 in $\hat{V}_{{\rm loc}}^{{\rm full}}$ 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 for $\hat{V}_{{\rm loc}}^{{\rm dens}}$ with orbital-independent values for U, $U^{\prime}$ and $J=(U-U^{\prime} )/2$, see equation (99), with those for the full local Hamiltonian, $\hat{V}_{{\rm loc}}^{{\rm full}}$, for the same values for U and J.

There exists no systematic (let alone rigorous) derivation of the double-counting correction in equation (43). 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

$$\begin{eqnarray}&&V_{{\rm dc};1}^{{\rm G}}=\frac{U}{2}\bar{n}(\bar{n}-1)-\frac{J}{2}\mathop{\sum }\limits_{\sigma }{{\bar{n}}_{\sigma }}(1-{{\bar{n}}_{\sigma }}),\end{eqnarray} \tag{ 102 }$$

where ${{\bar{n}}_{\sigma }}\equiv \sum _{c=1}^{{{N}_{c}}}C_{c,c;\sigma }^{{\rm G}}$, $\bar{n}\equiv {{\bar{n}}_{\uparrow }}+{{\bar{n}}_{\downarrow }}$, and N_c is the number of correlated orbitals (N_c = 5 for nickel). Note that only the two mean values U and J enter this double-counting 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 ${{\eta }_{c,c;\sigma }}$ which we may write as

$$\begin{eqnarray}&&{{\eta }_{c,c;\sigma }}\equiv \eta _{c,c;\sigma }^{{\rm G}}-\eta _{c,c;\sigma }^{{\rm dc}},\quad \eta _{c,c;\sigma }^{{\rm dc}}=\frac{\partial V_{{\rm dc}}^{{\rm G}}}{\partial {{C}_{c,c;\sigma }}}\ .\end{eqnarray} \tag{ 103 }$$

For nickel, the cubic symmetry guarantees that

$$\begin{eqnarray}&&C_{c,c^{\prime} ;\sigma }^{{\rm G}}={{C}_{c,c^{\prime} ;\sigma }}={{\delta }_{c,c^{\prime} }}{{C}_{c,c;\sigma }},\end{eqnarray} \tag{ 104 }$$

i.e., the correlated and uncorrelated local densities agree with each other. The double-counting correction (102) leads to $\eta _{c,c;\sigma }^{{\rm dc},1}=U(\bar{n}-1/2)+J({{\bar{n}}_{\sigma }}-1/2)$. It has been argued in [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

$$\begin{eqnarray}&&V_{{\rm dc};2}^{{\rm G}}:\eta _{c,c;\sigma }^{{\rm dc},2}=\eta _{c,c;\sigma }^{{\rm G}}\quad {\rm so}\;{\rm that}\;\;{{\eta }_{c,c;\sigma }}\equiv 0,\end{eqnarray} \tag{ 105 }$$

$$\begin{eqnarray}&&V_{{\rm dc};3}^{{\rm G}}:\eta _{c,c;\sigma }^{{\rm dc},3}=\frac{1}{{{N}_{c}}}\mathop{\sum }\limits_{c}\eta _{c,c;\sigma }^{{\rm G}}.\end{eqnarray} \tag{ 106 }$$

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 [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.

As a first step, we perform a DFT calculation that corresponds to setting $U=J=0$. We use the LDA exchange-correlation potential,

$$\begin{eqnarray}&&{{v}_{{\rm H},{\rm xc},\sigma }}({\bf r})={{\left. \frac{\partial {{E}_{{\rm LDA},{\rm xc}}}\left[ \left\{ {{n}_{\sigma }}({\bf r}) \right\} \right]}{\partial {{n}_{\sigma }}({\bf r})} \right|}_{{{n}_{\sigma }}({\bf r})=n_{\sigma }^{0}({\bf r})}},\end{eqnarray} \tag{ 107 }$$

see equation (52), as implemented in QuantumEspresso. The Kohn–Sham equations are solved in the plane-wave basis using ultra-soft pseudo-potentials, see equation (A.23). This calculation provides the Kohn–Sham bandstructure ${{\epsilon }_{{\rm n},\sigma }}({\bf k})$ and the coefficients ${{C}_{{\bf G},{\rm n},\sigma }}({\bf k})=\langle {\bf k},{\bf G},\sigma |{\bf k},{\rm n},\sigma \rangle$ of the Kohn–Sham eigenstates ${{\psi }_{{\bf k},{\rm n},\sigma }}({\bf r})$ in the plane-wave basis. The implemented 'poor-man Wannier' program package provides the down-folded $3d$ Wannier orbitals ${{\phi }_{{\bf R},c,\sigma }}({\bf r})$. In the orbital Bloch basis the coefficients $\langle {\bf k},{\bf G},\sigma |{\bf k},c,\sigma \rangle$ describe ${{\phi }_{{\bf k},c,\sigma }}({\bf r})$ in the plane-wave basis.

At the beginning we set $q_{a,\sigma }^{b,\sigma }={{\delta }_{a,b}}$ and $\tilde{\eta }=0$. Our Gutzwiller–Kohn–Sham loop consists of the following steps.

• (i)  
  Perform a DFT calculation with the Gutzwiller Kohn–Sham Hamiltonian from equation (B.11), and the Gutzwiller Kohn–Sham densities from equation (B.5). Here, the form ${{Q}_{{\bf G},{\bf G}^{\prime} ;\sigma }}({\bf k})={{\delta }_{{\bf G},{\bf G}^{\prime} }}+{{\sum }_{c}}({{q}_{c}}-1)\langle {\bf k},{\bf G},\sigma |{\bf k},c,\sigma \rangle \langle {\bf k},c,\sigma |{\bf k},{\bf G}^{\prime} ,\sigma \rangle$ is useful where only the correlated orbitals appear explicitly. After reaching a self-consistent density ${{n}_{\sigma }}({\bf r})$, calculate the local single-particle density matrix,

  $$\begin{eqnarray}\begin{array}{rcl} {{C}_{c,c;\sigma }} & = & \frac{1}{L}\mathop{\sum }\limits_{{\bf k},{\rm n}}{{f}_{{\bf k},{\rm n},\sigma }}\mathop{\sum }\limits_{{\bf G},{\bf G}^{\prime} }\langle {\bf k},{\rm n},\sigma |{\bf k},{\bf G},\sigma \rangle \langle {\bf k},{\bf G},\sigma |{\bf k},c,\sigma \rangle \\ {} & {} & \langle {\bf k},c,\sigma |{\bf k},{\bf G}^{\prime} ,\sigma \rangle \langle {\bf k},{\bf G}^{\prime} ,\sigma |{\bf k},{\rm n},\sigma \rangle , \\ \end{array}\end{eqnarray} \tag{ 108 }$$

  and the quantities $I_{{{c}_{1}},{{c}_{3}},{{c}_{2}},{{c}_{4}}}^{\sigma }$ from equation (96) and $K_{c,c^{\prime} }^{\sigma }$ from equation (97). For a proper convergence of the Gutzwiller–Kohn–Sham loop these quantities must be calculated with a momentum-space resolution that exceeds that of an ordinary DFT calculation considerably. To achieve this goal we use a tetrahedron method with 826 ${\bf k}$-points in the symmetry-reduced Brillouin zone.
• (ii)  
  Perform the inner minimization, i.e., minimize the energy functional ${{E}^{{\rm inner}}}$ in equation (95). This step provides the values for the Lagrange parameters Λ_l and for the Gutzwiller variational parameters $\tilde{\lambda }$ that determine the orbital-dependent renormalization factors ${{q}_{c,\sigma }}$ in equation (78).
• (iii)  
  Calculate the entries of $\tilde{\eta }$ from equation (94).
• (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.

In figure 1, we show the lattice constant and the magnetic moment as a function of U ($1\;{\rm eV}\leqslant U\leqslant 14\;{\rm 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 Hamiltonian $\hat{V}_{{\rm loc}}^{{\rm full}}$ and the double-counting correction ${{\hat{V}}_{{\rm dc};1}}$.

Zoom In Zoom Out Reset image size

As is well known, the DFT-LDA underestimates the lattice constant. We obtain $a_{0}^{{\rm LDA}}=6.47{{a}_{{\rm B}}}$, considerably smaller than the experimental value of ${{a}_{0}}=6.66{{a}_{{\rm B}}}$ where ${{a}_{{\rm B}}}=0.529177\;\overset{\circ}{A}$ is the Bohr radius. figure 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\gt 10\;{\rm eV}$.

Figure 1 shows the well-known fact that DFT-LDA reproduces the experimental value for the spin-only magnetic moment ${{m}_{{\rm so}}}$ very well, $m_{{\rm so}}^{{\rm LDA}}=0.58{{\mu }_{{\rm B}}}$ vs. $m_{{\rm so}}^{{\rm exp} }=0.55{{\mu }_{{\rm 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 figure 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\;{\rm eV}\lt U\lt 14\;{\rm eV}$ and $0.05\lt J/U\lt 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.

Our effective values are chosen to fit the experimental data for the lattice constant and the magnetic moment. The size of U and J agrees with those used in previous Gutzwiller-DFT studies on nickel [10, 12–14]. For the Gutzwiller-DFT the Hubbard-U lies between the bare, atomic value ${{U}_{{\rm bare}}}\approx 25\;{\rm eV}$ [28] and the low-frequency value for the screened on-site interaction ${{U}_{{\rm LDA}+{\rm RPA}}}(\omega \to 0)\approx 4\;{\rm eV}$, as obtained from LDA+Random-Phase Approximation [29] and used in LDA+DMFT [30]. This comparison shows that the Gutzwiller-DFT works with a partly screened value for the Hubbard interaction.

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}^{{\rm opt}}}=13\;{\rm eV},{{J}^{{\rm opt}}}=0.9\;{\rm eV})$ with an uncertainty of $\pm 1$ in the last digit. For our optimal values we show in figure 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.63{{a}_{{\rm B}}}$, in good agreement with the experimental value $a_{0}^{{\rm exp} }=6.66{{a}_{{\rm B}}}$. For the magnetic spin-only moment we obtain ${{m}_{{\rm so}}}=0.52{{\mu }_{{\rm B}}}$, in good agreement with the experimental value $m_{{\rm so}}^{{\rm exp} }=0.55{{\mu }_{{\rm B}}}$.

Zoom In Zoom Out Reset image size

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,

$$\begin{eqnarray}&&K={{V}_{0}}{{\left. \frac{{{{\rm d}}^{2}}E(V)}{{\rm d}{{V}^{2}}} \right|}_{V={{V}_{0}}}}\ .\end{eqnarray} \tag{ 109 }$$

This implies the Taylor expansion $E(V)=E({{V}_{0}})+(K{{V}_{0}}/2){{(V/{{V}_{0}}-1)}^{2}}+\ldots$ 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}_{{\rm B}}}-{{a}_{0}}/{{a}_{{\rm B}}})}^{2}}+\ldots$ with

$$\begin{eqnarray}&&{{e}_{2}}=\frac{9}{8}Ka_{{\rm B}}^{3}({{a}_{0}}/{{a}_{{\rm B}}}),\end{eqnarray} \tag{ 110 }$$

where we took into account that the fcc unit cell hosts four atoms, ${{V}_{0}}=Na_{0}^{3}/4$. The fit leads to $K=169\;{\rm GPa}$, in good agreement with the experimental value, $K=182\;{\rm GPa}$ [31]. It is a well-known fact that the DFT-LDA overestimates the bulk modulus of nickel. Indeed, our DFT-LDA gives ${{K}^{{\rm LDA}}}=245\;{\rm GPa}$.

We also calculated the lattice parameter and the magnetic spin-only moment for the density-dependent interaction V^dens_loc, see equation (98), with the same double-counting correction ${{\hat{V}}_{{\rm dc};1}}$. Our results do not show significant discrepancies for the ground-state properties. Note, however, that nickel is a special case because it has an almost filled 3d-shell (${{n}_{3d}}\approx 9/10$) such that the terms from $\hat{V}_{{\rm loc}}^{{\rm n}.{\rm dens}.}$ in equation (101) are more or less deactivated. Preliminary calculations for iron indicate that the missing interaction terms are more important for a partially filled $3d$-shell.

The combination of the full local interaction V_loc^full with the second and third double-counting correction, see equation (105) and (106), does not lead to reasonable values for the lattice constant, spin-only magnetic moment, and compressibility for nickel. If we fix the lattice constant to its experimental value, the Gutzwiller–Kohn–Sham equations lead to converged results for the second (but not for the third) double-counting correction; for the third double-counting correction, the $3d$ levels are discharged. In the next section, we use these converged results for ${{\hat{V}}_{{\rm dc};2}}$ for comparison with those for the standard double-counting correction ${{\hat{V}}_{{\rm dc};1}}$.

In figure 3 we show the quasi-particle band structure of fcc nickel for (${{U}^{{\rm opt}}}=13\;{\rm eV},{{J}^{{\rm opt}}}=0.9\;{\rm 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}^{{\rm LDA}}}=4.5\;{\rm eV}$ [12–14]. whereas we find $W=3.3\;{\rm eV}$, in agreement with experiment. This bandwidth reduction is due to the q-factors ${{q}_{t,\uparrow }}=0.851$, ${{q}_{t,\downarrow }}=0.824$, ${{q}_{e,\uparrow }}=0.852$, ${{q}_{e,\downarrow }}=0.819$, $\bar{q}={{\sum }_{\sigma }}(3{{q}_{t,\sigma }}+2{{q}_{e,\sigma }})/10=0.837$, so that $W\approx {{\bar{q}}^{2}}{{W}^{{\rm LDA}}}$.

Zoom In Zoom Out Reset image size

A more detailed comparison of the quasi-particle band structure with experiment is given in table 1. The overall agreement between experiment and theory for $\hat{V}_{{\rm loc}}^{{\rm full}}$ with ${{\hat{V}}_{{\rm 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 correction ${{\hat{V}}_{{\rm dc};2}}$ spoils this advantage. Therefore, this form of the double-correction term is not particularly useful for nickel.

Table 1.  Quasi-particle band energies with respect to the Fermi energy in eV at various high-symmetry points (counted positive for occupied states). $\langle \ldots \rangle$ indicates the spin average, errors bars in the experiments without spin resolution are given as ±. Theoretical data show the spin average and the exchange splittings in square brackets. ${{\Lambda }_{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 for $\hat{V}_{{\rm loc}}^{{\rm full}}$ with ${{\hat{V}}_{{\rm dc};1}}$, $\hat{V}_{{\rm loc}}^{{\rm full}}$ with ${{\hat{V}}_{{\rm dc};2}}$, and $\hat{V}_{{\rm loc}}^{{\rm dens}}$ with ${{\hat{V}}_{{\rm dc};1}}$, respectively, at (${{U}^{{\rm opt}}}=13\;{\rm eV},{{J}^{{\rm opt}}}=0.9\;{\rm eV}$).

Symmetry	Experiment	$\hat{V}_{{\rm loc}}^{{\rm full}}$ & ${{\hat{V}}_{{\rm dc};1}}$	$\hat{V}_{{\rm loc}}^{{\rm full}}$ & ${{\hat{V}}_{{\rm dc};2}}$	$\hat{V}_{{\rm loc}}^{{\rm dens}}$ & ${{\hat{V}}_{{\rm dc};1}}$
$\langle {{\Gamma }_{1}}\rangle$	$8.90\pm 0.30$	8.95[0.08]	8.99[0.08]	8.93[0.08]
$\langle {{\Gamma }_{{{25}^{\prime }}}}\rangle$	$1.30\pm 0.06$	1.51[0.65]	1.52[0.57]	1.56[0.80]
$\langle {{\Gamma }_{12}}\rangle$	$0.48\pm 0.08$	0.73[0.15]	0.66[0.43]	0.71[0.10]
$\langle {{X}_{1}}\rangle$	$3.30\pm 0.20$	3.37[0.27]	3.26[0.56]	3.42[0.10]
$\langle {{X}_{3}}\rangle$	$2.63\pm 0.10$	2.87[0.68]	2.87[0.61]	2.87[0.77]
${{X}_{2\uparrow }}$	$0.21\pm 0.03$	0.26	0.33	0.13
${{X}_{2\downarrow }}$	$0.04\pm 0.03$	0.14	−0.06	0.21
${{X}_{5\uparrow }}$	$0.15\pm 0.03$	0.32	0.29	0.41
${{\Delta }_{{{e}_{g}}}}({{X}_{2}})$	$0.17\pm 0.05$	0.12	0.39	−0.08
${{\Delta }_{{{t}_{2\,g}}}}({{X}_{5}})$	$0.33\pm 0.04$	0.60	0.51	0.70
$\langle {{L}_{1}}\rangle$	$3.66\pm 0.10$	3.49[0.61]	3.49[0.56]	3.55[0.83]
$\langle {{L}_{3}}\rangle$	$1.43\pm 0.07$	1.58[0.38]	1.52[0.52]	1.61[0.26]
${{L}_{3\uparrow }}$	$0.18\pm 0.03$	0.37	0.38	0.34
$\langle {{L}_{{{2}^{\prime }}}}\rangle$	$1.00\pm 0.20$	0.14[0.06]	0.17[0.06]	0.12[0.06]
$\langle {{\Lambda }_{3;1/2}}\rangle$	$0.50[0.21\pm 0.02]$	0.64[0.30]	0.61[0.45]	0.60[0.16]

We comment on two noticeable discrepancies between theory and experiment. First, the energy of the band ${{L}_{2^{\prime} }}$ 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_3d, ${{n}_{3p,3s}}$ in the $3d$ and $3p/3s$ bands. Second, the Gutzwiller DFT prediction for the exchange splitting ${{\Delta }_{{{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 with ${{\hat{V}}_{{\rm 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 ${{\Delta }_{{{e}_{g}}}}({{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.

In the following we assume that the potential is lattice periodic,

$$\begin{eqnarray}&&V_{\sigma }^{{\rm KS}}({\bf r})=U({\bf r})+{{V}_{{\rm Har}}}({\bf r})+{{v}_{{\rm sp},{\rm xc},\sigma }}({\bf r})=V_{\sigma }^{{\rm KS}}({\bf r}+{\bf R}),\end{eqnarray} \tag{ A.11 }$$

where R is a lattice vector. The Fourier components are finite only for reciprocal lattice vectors G,

$$\begin{eqnarray}&&V_{{\bf G},\sigma }^{{\rm KS}}=\frac{1}{V}\int {\rm d}{\bf r}V_{\sigma }^{{\rm KS}}({\bf r}){{{\rm e}}^{-{\rm i}{\bf G}\cdot {\bf r}}},\end{eqnarray} \tag{ A.12 }$$

where V is the crystal volume. As a consequence of the lattice periodicity, the crystal momentum ${\bf k}$ from the first Brillouin zone is a good quantum number.

As seen from equation (41), the Kohn–Sham Hamiltonian is diagonalized for the single-particle states ${{\psi }_{{\bf k},{\rm n},\sigma }}({\bf r})=\langle {\bf r}|{\bf k},{\rm n},\sigma \rangle$ that obey

$$\begin{eqnarray}&&h_{\sigma }^{{\rm KS}}({\bf r}){{\psi }_{{\bf k},{\rm n},\sigma }}({\bf r})={{\epsilon }_{{\rm n},\sigma }}({\bf k}){{\psi }_{{\bf k},{\rm n},\sigma }}({\bf r}),\end{eqnarray} \tag{ A.13 }$$

where n is the band index. Equations (A.13) are the Kohn–Sham equations [1].

In its eigenbasis, the Kohn–Sham Hamiltonian takes the form

$$\begin{eqnarray}&&{{\hat{H}}^{{\rm KS}}}=\mathop{\sum }\limits_{{\bf k},{\rm n},\sigma }{{\epsilon }_{{\rm n},\sigma }}({\bf k})\hat{b}_{{\bf k},{\rm n},\sigma }^{\dagger }\hat{b}_{{\bf k},{\rm n},\sigma }^{{}}.\end{eqnarray} \tag{ A.14 }$$

Its ground state is given by

$$\begin{eqnarray}&&|{{\Phi }_{0}}\rangle =\mathop{\prod }\limits_{\sigma }{{\mathop{\prod }\limits_{{\bf k},{\rm n}}}^{^{\prime }}}\hat{b}_{{\bf k},{\rm n},\sigma }^{\dagger }|{\rm vac}\rangle ,\end{eqnarray} \tag{ A.15 }$$

where the N levels lowest in energy are occupied as indicated by the prime at the product, ${{\epsilon }_{{\rm n},\sigma }}({\bf k})\leqslant {{E}_{{\rm F},\sigma }}$. Then,

$$\begin{eqnarray}&&{{f}_{{\bf k},{\rm n},\sigma }}=\langle {{\Phi }_{0}}|\hat{b}_{{\bf k},{\rm n},\sigma }^{\dagger }\hat{b}_{{\bf k},{\rm n},\sigma }^{{}}|{{\Phi }_{0}}\rangle =\Theta \left( {{E}_{{\rm F},\sigma }}-{{\epsilon }_{{\rm n},\sigma }}({\bf k}) \right)\end{eqnarray} \tag{ A.16 }$$

is unity for occupied levels up to the Fermi energy ${{E}_{{\rm F},\sigma }}$, and zero otherwise.

From equation (26), the field operators read

$$\begin{eqnarray}&&\hat{\Psi }_{\sigma }^{{}}({\bf r})=\mathop{\sum }\limits_{{\bf k},{\rm n}}{{\psi }_{{\bf k},{\rm n},\sigma }}({\bf r})\hat{b}_{{\bf k},{\rm n},\sigma }^{{}},\qquad \hat{\Psi }_{\sigma }^{\dagger }({\bf r})=\mathop{\sum }\limits_{{\bf k},{\rm n}}\psi _{{\bf k},{\rm n},\sigma }^{*}({\bf r})\hat{b}_{{\bf k},{\rm n},\sigma }^{\dagger }.\end{eqnarray} \tag{ A.17 }$$

Therefore, the ground-state density is readily obtained as

$$\begin{eqnarray}\begin{array}{rcl} n_{\sigma }^{0}({\bf r}) & = & \langle {{\Phi }_{0}}|\hat{\Psi }_{\sigma }^{\dagger }({\bf r})\hat{\Psi }_{\sigma }^{{}}({\bf r})|{{\Phi }_{0}}\rangle =\mathop{\sum }\limits_{{\bf k},{\rm n}}{{f}_{{\bf k},{\rm n},\sigma }}{{\left| {{\psi }_{{\bf k},{\rm n},\sigma }}({\bf r}) \right|}^{2}}=\langle {\bf r}|\mathop{\sum }\limits_{{\bf k}}\hat{\rho }_{\sigma }^{(0)}({\bf k})|{\bf r}\rangle , \\ \hat{\rho }_{\sigma }^{(0)}({\bf k}) & = & \mathop{\sum }\limits_{n}{{f}_{{\bf k},{\rm n},\sigma }}|{\bf k},{\rm n},\sigma \rangle \langle {\bf k},{\rm n},\sigma |, \\ \end{array}\end{eqnarray} \tag{ A.18 }$$

see also equation (31). Since this quantity enters the Kohn–Sham Hamiltonian, its solution must be achieved self-consistently.

In many codes, the Kohn–Sham Hamiltonian is formulated in the plane-wave basis $|{\bf k},{\bf G},\sigma \rangle$ with real-space representation

$$\begin{eqnarray}&&\langle {\bf r}|{\bf k},{\bf G},\sigma \rangle =\sqrt{\frac{1}{V}}{{{\rm e}}^{{\rm i}({\bf k}+{\bf G})\cdot {\bf r}}}.\end{eqnarray} \tag{ A.19 }$$

In this basis, the field operators are given by

$$\begin{eqnarray}&&\hat{\Psi }_{\sigma }^{{}}({\bf r})=\sqrt{\frac{1}{V}}\mathop{\sum }\limits_{{\bf k},{\bf G}}{{{\rm e}}^{-{\rm i}({\bf k}+{\bf G})\cdot {\bf r}}}\hat{p}_{{\bf k},{\bf G},\sigma }^{{}},\ \hat{\Psi }_{\sigma }^{\dagger }({\bf r})=\sqrt{\frac{1}{V}}\mathop{\sum }\limits_{{\bf k},{\bf G}}{{{\rm e}}^{{\rm i}({\bf k}+{\bf G})\cdot {\bf r}}}\hat{p}_{{\bf k},{\bf G},\sigma }^{\dagger },\end{eqnarray} \tag{ A.20 }$$

and the Kohn–Sham Hamiltonian reads

$$\begin{eqnarray}&&{{\hat{H}}^{{\rm KS}}}=\mathop{\sum }\limits_{{\bf k},\sigma }\mathop{\sum }\limits_{{\bf G},{\bf G}^{\prime} }T_{{\bf G},{\bf G}^{\prime} ;\sigma }^{{\rm KS}}({\bf k})\hat{p}_{{\bf k},{\bf G},\sigma }^{\dagger }\hat{p}_{{\bf k},{\bf G}^{\prime} ,\sigma }^{{}}.\end{eqnarray} \tag{ A.21 }$$

Equation (29) shows that the entries of the Kohn–Sham Hamiltonian in reciprocal space are given by

$$\begin{eqnarray}&&T_{{\bf G},{\bf G}^{\prime} ;\sigma }^{{\rm KS}}({\bf k})={{\delta }_{{\bf G},{\bf G}^{\prime} }}\frac{1}{2\;m}{{\left( {\bf k}+{\bf G} \right)}^{2}}+V_{{\bf G}-{\bf G}^{\prime} ;\sigma }^{{\rm KS}}\end{eqnarray} \tag{ A.22 }$$

for each ${\bf k}$ from the first Brillouin zone. The eigenvalues of the Kohn–Sham matrix in reciprocal space are ${{\epsilon }_{{\rm n},\sigma }}({\bf k})$, and the solution of the eigenvalue equation [1]

$$\begin{eqnarray}&&\frac{1}{2\;m}{{\left( {\bf k}+{\bf G} \right)}^{2}}{{C}_{{\bf G},{\rm n},\sigma }}({\bf k})+\mathop{\sum }\limits_{{\bf G}^{\prime} }V_{{\bf G}-{\bf G}^{\prime} ;\sigma }^{{\rm KS}}{{C}_{{\bf G}^{\prime} ,{\rm n},\sigma }}({\bf k})={{\epsilon }_{{\rm n},\sigma }}({\bf k}){{C}_{{\bf G},{\rm n},\sigma }}({\bf k})\end{eqnarray} \tag{ A.23 }$$

for given $({\bf k},{\rm n})$ gives the entries of the eigenvectors, ${{C}_{{\bf G},{\rm n},\sigma }}({\bf k})=\langle {\bf k},{\bf G},\sigma |{\bf k},{\rm n},\sigma \rangle$. Implemented plane-wave codes provide the band energies ${{\epsilon }_{{\rm n},\sigma }}({\bf k})$ and the coefficients ${{C}_{{\bf G},{\rm n},\sigma }}({\bf k})$ so that the Kohn–Sham eigenstates are obtained as

$$\begin{eqnarray}\begin{array}{rcl} |{\bf k},{\rm n},\sigma \rangle & = & \mathop{\sum }\limits_{{\bf G}}{{C}_{{\bf G},{\rm n},\sigma }}({\bf k})|{\bf k},{\bf G},\sigma \rangle , \\ {{\psi }_{{\bf k},{\rm n},\sigma }}({\bf r}) & = & \sqrt{\frac{1}{V}}\mathop{\sum }\limits_{{\bf G}}{{C}_{{\bf G},{\rm n},\sigma }}({\bf k}){{e}^{{\rm i}({\bf k}+{\bf G})\cdot {\bf r}}}. \\ \end{array}\end{eqnarray} \tag{ A.24 }$$

In order to make contact with many-particle approaches based on Hubbard-type models, we need to identify orbitals that enter the local two-particle interaction. Implemented plane-wave codes provide the transformation coefficients ${{F}_{({\bf k},{\rm n}),({\bf R},b);\sigma }}$ from Bloch eigenstates $|{\bf k},{\rm n},\sigma \rangle$ to orbital Wannier states $|{\bf R},b,\sigma \rangle$,

$$\begin{eqnarray}&&|{\bf R},b,\sigma \rangle =\mathop{\sum }\limits_{{\bf k},n}{{F}_{({\bf k},{\rm n}),({\bf R},b);\sigma }}|{\bf k},{\rm n},\sigma \rangle ,\ {{F}_{({\bf k},{\rm n}),({\bf R},b);\sigma }}=\langle {\bf k},{\rm n},\sigma |{\bf R},b,\sigma \rangle .\end{eqnarray} \tag{ A.25 }$$

The Wannier orbitals

$$\begin{eqnarray}&&{{\phi }_{{\bf R},b,\sigma }}({\bf r})=\langle {\bf r}|{\bf R},b,\sigma \rangle \end{eqnarray} \tag{ A.26 }$$

are maximal around a lattice site R and the orbital index b resembles atomic quantum numbers, e.g., $b=3\;s,3p,3d$. In the orbital Wannier basis the field operators are given by

$$\begin{eqnarray}&&\hat{\Psi }_{\sigma }^{\dagger }({\bf r})=\mathop{\sum }\limits_{{\bf R},b}\phi _{{\bf R},b,\sigma }^{*}({\bf r})\hat{c}_{{\bf R},b,\sigma }^{\dagger },\quad \hat{\Psi }_{\sigma }^{{}}({\bf r})=\mathop{\sum }\limits_{{\bf R},b}{{\phi }_{{\bf R},b,\sigma }}({\bf r})\hat{c}_{{\bf R},b,\sigma }^{{}},\end{eqnarray} \tag{ A.27 }$$

and the Kohn–Sham Hamiltonian in the orbital Wannier basis becomes

$$\begin{eqnarray}&&{{\hat{H}}^{{\rm KS}}}=\mathop{\sum }\limits_{{\bf R},b,{\bf R}^{\prime} ,b^{\prime} ,\sigma }T_{({\bf R},b),({\bf R}^{\prime} ,b^{\prime} );\sigma }^{{\rm KS}}\hat{c}_{{\bf R},b,\sigma }^{\dagger }\hat{c}_{{\bf R}^{\prime} ,b^{\prime} ,\sigma }^{{}}\end{eqnarray} \tag{ A.28 }$$

with the overlap matrix elements

$$\begin{eqnarray}&&T_{({\bf R},b),({\bf R}^{\prime} ,b^{\prime} );\sigma }^{{\rm KS}}=\int {\rm d}{\bf r}\phi _{{\bf R},b,\sigma }^{*}({\bf r})h_{\sigma }^{{\rm KS}}({\bf r}){{\phi }_{{\bf R}^{\prime} ,b^{\prime} ,\sigma }}({\bf r}),\end{eqnarray} \tag{ A.29 }$$

see equation (41). These matrix elements appear in a tight-binding representation of the kinetic energy in Hubbard-type models.

For later use we also define the orbital Bloch basis,

$$\begin{eqnarray}&&{{\phi }_{{\bf k},b,\sigma }}({\bf r})=\sqrt{\frac{1}{L}}\mathop{\sum }\limits_{{\bf R}}{{{\rm e}}^{{\rm i}{\bf k}\cdot {\bf R}}}{{\phi }_{{\bf R},b,\sigma }}({\bf r})\;,\;{{\phi }_{{\bf R},b,\sigma }}({\bf r})=\sqrt{\frac{1}{L}}\mathop{\sum }\limits_{{\bf k}}{{{\rm e}}^{-{\rm i}{\bf k}\cdot {\bf R}}}{{\phi }_{{\bf k},b,\sigma }}({\bf r})\;,\end{eqnarray} \tag{ A.30 }$$

where ${\bf k}$ is from the first Brillouin zone and L is the number of lattice sites. The field operators are given by

$$\begin{eqnarray}&&\hat{\Psi }_{\sigma }^{\dagger }({\bf r})=\mathop{\sum }\limits_{{\bf k},b}\phi _{{\bf k},b,\sigma }^{*}({\bf r})\hat{c}_{{\bf k},b,\sigma }^{\dagger },\quad \hat{\Psi }_{\sigma }^{{}}({\bf r})=\mathop{\sum }\limits_{{\bf k},b}{{\phi }_{{\bf k},b,\sigma }}({\bf r})\hat{c}_{{\bf k},b,\sigma }^{{}}.\end{eqnarray} \tag{ A.31 }$$

In the orbital Wannier basis, the Kohn–Sham single-particle Hamiltonian reads

$$\begin{eqnarray}\begin{array}{rcl} {{{\hat{H}}}^{{\rm KS}}} & = & \mathop{\sum }\limits_{{\bf k},b,b^{\prime} ,\sigma }T_{b,b^{\prime} ;\sigma }^{{\rm KS}}({\bf k})\hat{c}_{{\bf k},b,\sigma }^{\dagger }\hat{c}_{{\bf k},b^{\prime} ,\sigma }^{{}}, \\ T_{b,b^{\prime} ;\sigma }^{{\rm KS}}({\bf k}) & = & \int {\rm d}{\bf r}\phi _{{\bf k},b,\sigma }^{*}({\bf r})h_{\sigma }^{{\rm KS}}({\bf r}){{\phi }_{{\bf k},b^{\prime} ,\sigma }}({\bf r}). \\ \end{array}\end{eqnarray} \tag{ A.32 }$$