Charge self-consistent many-body corrections using optimized projected localized orbitals

To perform DFT+DMFT calculations one needs a way to transform between the basis of the Kohn–Sham (KS) states and the localized basis of the subspace used to define lattice models, e.g. a Hubbard model. The most obvious way to do this is to perform a unitary transformation of a subset of KS states to construct a set of local states. This method is at the heart of various types of Wannier function based methods, such as the maximally-localized Wannier functions [24].

In many cases, when KS bands with different characters are strongly entangled, it becomes however difficult to construct a well-defined unitary transformation. In this case it is more advantageous to use projector operators which, acting on the KS Hilbert space, project out only KS states with a desired character. This is the basis of projected localized orbitals (PLO) [8].

In the PLO formalism one starts by defining an orthonormal localized basis set $\newcommand{\ket}[1]{|{#1}\rangle} \ket{\chi_{L}}$ associated with each correlated site, which is typically indexed by local quantum numbers, e.g. $L = \{l, m, \sigma, ...\}$ with ($l, m$) and σ being orbital and spin angular-momentum quantum-numbers, respectively. $\newcommand{\ket}[1]{|{#1}\rangle} \{\ket{\chi_{L}}\}$ spans a correlated subspace $\newcommand{\corrspace}{\mathcal{C}} \corrspace$ at each correlated site. Assuming $\newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \newcommand{\bra}[1]{\langle{#1}|} \braket{\chi_{L}}{\chi_{L'}} = \delta_{LL'}$, any operator acting on this space can be constructed by projecting onto $\newcommand{\corrspace}{\mathcal{C}} \corrspace{}$ using a projector operator $\newcommand{\corrspace}{\mathcal{C}} \newcommand{\projC}{\hat{P}^{\corrspace}} \projC$, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\projC}{\hat{P}^{\corrspace}} \newcommand{\corrspace}{\mathcal{C}} \displaystyle \hat{A}^{{\rm imp}} = \projC \hat{A} \projC, \nonumber \end{align} \tag{ 1 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\projC}{\hat{P}^{\corrspace}} \newcommand{\corrspace}{\mathcal{C}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\ket}[1]{|{#1}\rangle} \displaystyle \projC = \sum_{L} \ket{\chi_{L}}\bra{\chi_{L}}. \label{eq:Pc} \nonumber \end{align} \tag{ 2 }$$

An arbitrary vector $\newcommand{\ket}[1]{|{#1}\rangle} \ket{\Psi}$ of the Hilbert space can be decomposed in terms of local states by writing $\newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\corrspace}{\mathcal{C}} \newcommand{\projC}{\hat{P}^{\corrspace}} \projC \ket{\Psi} = \sum_{L} \ket{\chi_{L}} \braket{\chi_{L}}{\Psi}$. If we now consider a complete basis $\newcommand{\ket}[1]{|{#1}\rangle} \ket{\Psi_{\mu}}$, the projector operator is completely defined by specifying PLO functions $\newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\e}{{\rm e}} P_{L, \mu} \equiv \braket{\chi_{L}}{\Psi_{\mu}}$.

As described in detail in [8, 9], the PLO projector in the PAW framework [7, 18] can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\psp}{\ps{p}} \newcommand{\eaw}{\phi} \newcommand{\psKS}{\ps{\Psi}_{\nu\kv}} \newcommand{\ps}[1]{\tilde{#1}} \newcommand{\Rv}{{\bf R}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \displaystyle P^{\Rv}_{L, \nu}(\kv) = \sum_{i} \braket{\chi_{L}^{\Rv}}{\eaw_{i}} \braket{\psp_{i}}{\psKS}, \nonumber \end{align} \tag{ 3 }$$

where $\newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\Rv}{{\bf R}} \ket{\chi_{L}^{\Rv}}$ are localized basis functions associated with each correlated site $\newcommand{\Rv}{{\bf R}} {\Rv}$, $\newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\eaw}{\phi} \ket{{\eaw_{i}}}$ are all-electron partial waves, and $\newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\ps}[1]{\tilde{#1}} \newcommand{\psp}{\ps{p}} \ket{\psp_{i}}$ are the standard PAW projectors. In the following, we will omit the site index $\newcommand{\Rv}{{\bf R}} \Rv$ unless it leads to a confusion. The index i stands for the PAW channel $\newcommand{\ich}{n} \ich$, the angular momentum quantum number l, and its magnetic quantum number m. $\newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\ps}[1]{\tilde{#1}} \newcommand{\psKS}{\ps{\Psi}_{\nu\kv}} \ket{\psKS}$ are pseudo-KS states. For the PAW potentials distributed with vasp, generally a minimum of two channels $\newcommand{\ich}{n} \ich$ for each angular quantum number are used. For instance for 3d elements, the first 3d channel is placed at the energy of the bound 3d state in the atom, and a second channel is added a few eV above the bound atomic 3d state. Inclusion of the second channel improves the transferability of the PAW potentials and the description of the scattering properties greatly. For s and p states in, e.g. 3d transition metals, the first channel usually describes 3s and 3p semi-core states and the second channel is placed somewhere in the valence regime (4s and 4p states) or even above the vacuum level. Although a description of how the projectors have been obtained is stored in more recent PAW potential files, it is not always straightforward for the user to identify the best suitable PLO projectors.

Specifically, various choices are possible for $\newcommand{\ket}[1]{|{#1}\rangle} \ket{\chi_{L}}$ [9]. Conveniently, one might simply use the all-electron partial waves for a particular PAW channel $\newcommand{\ich}{n} \ich$, $\newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\ich}{n} \newcommand{\eaw}{\phi} \ket{\eaw_{L \ich}}$, as the local basis. However, doing so in vasp can lead to a non-optimal projection. For instance, for transition metals to the left of the periodic table, the d states in the solid might be more or less contracted than those in the atom, so that the 'best' $\newcommand{\ket}[1]{|{#1}\rangle} \ket{\chi_{L}}$ is a linear combination of the two available 3d channels. Likewise, projection onto the TM valence s states is often difficult, because of the presence of semi-core s states in the PAW potential file. Ideally, the user should not need to make a manual choice of the local $\newcommand{\ket}[1]{|{#1}\rangle} \ket{\chi_{L}}$ functions.

In the following we explain a protocol that largely resolves these issues. The first step is not strictly required, but simplifies the subsequent coding somewhat. We first construct a set of PAW projectors and partial waves that are orthogonalized inside the PAW sphere. This can be achieved by diagonalizing the all-electron one-center overlap matrix $\newcommand{\ich}{n} \newcommand{\jch}{n'} O_{\ich\jch}$ (for each angular and magnetic quantum number L),

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\eaw}{\phi} \newcommand{\jch}{n'} \newcommand{\ich}{n} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \displaystyle O_{\ich\jch} = \braket{\eaw_{L \ich}}{\eaw_{L \jch}}, \label{eq:overlap_Q} \nonumber \end{align} \tag{ 4 }$$

which gives the eigenvalues $\newcommand{\ich}{n} \lambda_{\ich}$ and the eigenvector matrix U,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Lambda = {\rm diag}(\lambda_{1}, \dots, \lambda_{n}), \nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \displaystyle \Lambda = U^{\dagger} O U. \nonumber \end{align} \tag{ 6 }$$

The new partial waves and the corresponding PAW projectors can now be defined as linear combinations of the original ones:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\eaw}{\phi} \newcommand{\jch}{n'} \newcommand{\ich}{n} \newcommand{\ket}[1]{|{#1}\rangle} \displaystyle \ket{\xi_{L \ich}} = \frac{1}{\sqrt{\lambda_{\ich}}} \sum_{\jch} U_{\ich\jch} \ket{\eaw_{L \jch}}, \nonumber \end{align} \tag{ 7 }$$

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\psp}{\ps{p}} \newcommand{\jch}{n'} \newcommand{\ich}{n} \newcommand{\ps}[1]{\tilde{#1}} \newcommand{\bra}[1]{\langle{#1}|} \displaystyle \bra{\tilde{\beta}_{L \ich}} = \sqrt{\lambda_{\ich}} \sum_{\jch} U^{\dagger}_{\jch\ich} \bra{\psp_{L \jch}}. \nonumber \end{align} \tag{ 8 }$$

This new set of projectors and partial waves have the important property that the all-electron partial waves $\newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\ich}{n} \ket{\xi_{L \ich}}$ are orthogonal to each other and that they are normalized to one. Thus, we will now refer to them as 'orthonormal' partial waves. The choice above is not unique, though, for instance one might also adopt a Gram–Schmidt orthogonalization procedure.

To construct a unique projector, we seek a unitary transformation that maximizes the overlap between the projector and KS valence states within a chosen (user supplied) energy window $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} [\eps^P_{\min}, \eps^P_{\max}]$. Specifically, we diagonalize a matrix

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\jch}{n'} \newcommand{\ich}{n} \newcommand{\psKS}{\ps{\Psi}_{\nu\kv}} \newcommand{\ps}[1]{\tilde{#1}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \displaystyle M_{\ich \jch} = \sum_{\nu \kv} \braket{\tilde{\beta}_{L \ich}}{\psKS} \braket{\psKS}{\tilde{\beta}_{L \jch}}, \nonumber \end{align} \tag{ 9 }$$

where the sum is restricted to states with energies $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \eps_{\nu \kv}\in[\eps^P_{\min}, \eps^P_{\max}]$. We then pick the eigenvector $\newcommand{\up}{{\uparrow}} \newcommand{\ich}{n} \upsilon_{\ich}$ with the largest absolute eigenvalue to construct local basis functions and corresponding PAW projectors,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ich}{n} \newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\up}{{\uparrow}} \displaystyle \ket{\chi_{L}} = \sum_{\ich} \upsilon_{\ich} \ket{\xi_{L \ich}}, \nonumber \end{align} \tag{ 10 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ich}{n} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\up}{{\uparrow}} \displaystyle \bra{\tilde{\pi}_{L}} = \sum_{\ich} \upsilon^{*}_{\ich} \bra{\tilde{\beta}_{L \ich}}, \nonumber \end{align} \tag{ 11 }$$

with projector functions given by⁸

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\psKS}{\ps{\Psi}_{\nu\kv}} \newcommand{\ps}[1]{\tilde{#1}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \displaystyle P_{L,\nu}(\kv) = \braket{\tilde{\pi}_{L}}{\psKS}. \label{eq:projector-functions} \nonumber \end{align} \tag{ 12 }$$

We note in passing that both steps could be combined into a single computational step, by diagonalization of a generalized eigenvalue problem involving the overlap matrix O and the matrix M evaluated using the original set of projectors and partial waves.

Note that since the all-electron partial waves do not form an orthonormal basis set between different PAW spheres and because projection is performed for a subset of KS states, the above procedure will usually produce non-normalized local states. However, normalized local states can be constructed if we orthonormalize the projector functions as follows,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle O_{L L'} = \sum_{\nu\kv} P_{L,\nu}(\kv) P^{*}_{\nu, L'}(\kv), \label{eq:projector_overlap} \nonumber \end{align} \tag{ 13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle P_{L,\nu}(\kv) \leftarrow \sum_{L'} O_{L L'}^{-\frac{1}{2}} P_{L',\nu}(\kv). \label{eq:projector-orthonormalized} \nonumber \end{align} \tag{ 14 }$$

On a technical level, the optimal projectors are constructed internally in vasp according to (4)–(12) for each sites $\newcommand{\Rv}{{\bf R}} \Rv$ given an energy window $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} [\eps^P_{\min}, \eps^P_{\max}]$. The last orthonormalization step (14) is done externally as post-processing. This allows one to choose a different energy window for the summation in (13) in order to fine-tune the degree of localization of the impurity Wannier functions.

The projection scheme has been used quite widely in the solid state community [8, 9], even in combination with VASP. When using VASP, these calculations were often not based on a solid mathematical foundation. VASP used to project onto all available PAW projectors with a certain angular and momentum quantum number lm, summed the resulting densities, averaged the phase factors and intensities over all lm projectors, and finally wrote the results to a file (PROCAR). Such a prescription, in particular, the averaging of the intensities and phase factors was done ad hoc without a proper mathematical prescription. This can result in unsatisfactory results if the two projectors span a completely different subspace, for instance Ni 3d and 4d states, or transition metal semi-core p states and valence p states. We will demonstrate this issue for NiO in section 3.1.

A general extension of DFT Kohn–Sham equations to a many-body problem based on the Hubbard model can be formulated using the total-energy functional [10, 25]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Tr}{{\rm Tr}} \newcommand{\ks}{\mathtt{KS}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\rmi}{{\rm i}} \newcommand{\dc}{\textsc{dc}} \newcommand{\xc}{{\rm xc}} \newcommand{\rv}{{\bf r}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle E [n,& {\hat G}^\mathcal{C} ] = \frac{1}{\beta} \sum_{\iomn\kv} \Tr\{\hat G(\kv,\iomn) \hat H_{\ks}(\kv)\} \nonumber \\ & \qquad +E_{\rm H}[n] + E_{\xc}[n] - \int {\rm d}\rv\, (v_{\xc} + v_{\rm H}) n(\rv) \nonumber \\ & \qquad +E_{{\rm corr}}[{\hat G}^\mathcal{C}] - E_{\dc}[{\hat G}^\mathcal{C}], \label{eq:etot_dmft} \nonumber \end{align} \tag{ 15 }$$

where $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\rmi}{{\rm i}} \newcommand{\iomn}{\rmi\omega_{n}} {\hat G}(\kv, \iomn)$ is the Green's function containing many-body effects, and $\newcommand{\ket}[1]{|{#1}\rangle} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\ks}{\mathtt{KS}} \hat H_{\ks}(\kv)=\sum_{\nu} \ket{\Psi_{\nu\kv}} \varepsilon_{\nu \kv} \bra{\Psi_{\nu\kv}}$ is the KS Hamiltonian. The charge density, $\newcommand{\rv}{{\bf r}} n(\rv)$, and Green's function in the correlated subspace, $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\rmi}{{\rm i}} \newcommand{\iomn}{\rmi\omega_{n}} {\hat G}^\mathcal{C}(\kv, \iomn)$, are both related to $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\rmi}{{\rm i}} \newcommand{\iomn}{\rmi\omega_{n}} \hat G(\kv, \iomn)$ by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Tr}{{\rm Tr}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\rmi}{{\rm i}} \newcommand{\braopket}[3]{\langle{#1}|{#2}|{#3}\rangle} \newcommand{\rv}{{\bf r}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \displaystyle n(\rv) = \frac{1}{\beta} \Tr \sum_{\iomn\kv} \braopket{\rv}{\hat G(\kv,\iomn)}{\rv}, \nonumber \end{align} \tag{ 16 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\rmi}{{\rm i}} \newcommand{\projC}{\hat{P}^{\corrspace}} \newcommand{\corrspace}{\mathcal{C}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle {\hat G}^\mathcal{C}(\kv,\iomn) = \projC(\kv) \hat G(\kv,\iomn) \projC(\kv), \nonumber \end{align} \tag{ 17 }$$

where $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\corrspace}{\mathcal{C}} \newcommand{\projC}{\hat{P}^{\corrspace}} \projC_\kv$ projects onto correlated localized states, see (2). $E_{{\rm corr}}$ is the energy contribution of the interaction term (Hubbard U-term) and E_dc is the double-counting correction.

The first four terms of the above expression form the usual DFT total energy calculated for the density matrix and charge density of the interacting system, which is non-diagonal in this case. It is thus clear that to get the correct value of the total energy the density matrix must be calculated in a self-consistent manner.

In the framework of DFT+DMFT [4, 6] the interacting Green's function is defined as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ks}{\mathtt{KS}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\rmi}{{\rm i}} \newcommand{\idop}{\hat {1\!\!1}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle \hat G(\kv,\iomn) =\left[(\iomn + \mu) \idop - \hat H_{\ks}(\kv) - \hat \Sigma^{\ks}(\kv, \iomn) \right]^{-1}, \nonumber \end{align} \tag{ 18 }$$

where $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\rmi}{{\rm i}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\ks}{\mathtt{KS}} \hat \Sigma^{\ks}(\kv, \iomn)$ is obtained by up-folding the local self-energy,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ks}{\mathtt{KS}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\rmi}{{\rm i}} \newcommand{\KSp}{\Psi_{\nu'\kv}} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\ket}[1]{|{#1}\rangle} \displaystyle \hat \Sigma^{\ks}(\kv, \iomn) = &\sum_{\nu \nu'} \ket{\KS}\bra{\KSp} \nonumber \\ &\cdot \sum_{m m'} P^{*}_{\nu,m}(\kv) \Sigma_{mm'}(\iomn) P_{m',\nu'}(\kv) \nonumber \end{align} \tag{ 19 }$$

where the local self energy's matrix elements in the basis of localized impurity states, $\newcommand{\ket}[1]{|{#1}\rangle} \ket{\chi_{lm}}$, are

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\rmi}{{\rm i}} \newcommand{\Sigimp}{\Sigma^{{\rm imp}}} \newcommand{\dc}{\textsc{dc}} \displaystyle \Sigma_{mm'}(\iomn) = \Sigimp_{mm'}(\iomn) - \Sigma^{\dc}_{mm'}, \label{eq:self_energy_dc} \nonumber \end{align} \tag{ 20 }$$

which consist of the purely local impurity self-energy $\newcommand{\rmi}{{\rm i}} \newcommand{\iomn}{\rmi\omega_{n}} \hat \Sigma^{{\rm imp}}(\iomn)$, obtained from an impurity solver, and the double counting $\newcommand{\dc}{\textsc{dc}} \hat \Sigma^{\dc}$. Using the PLO functions we can write this Green's function in terms of its KS matrix elements

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\ks}{\mathtt{KS}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\rmi}{{\rm i}} \newcommand{\idop}{\hat {1\!\!1}} \newcommand{\eps}{\varepsilon} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle G_{\nu\nu'}(\kv, \iomn) = \Big [ (\iomn +& \mu - \eps_{\nu\kv}) \idop - \hat \Sigma^{\ks}(\kv, \iomn)\Big]^{-1}_{\nu\nu'}. \nonumber \end{align} \tag{ 21 }$$

If we define an energy window $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} [\eps^C_{\min}, \eps^C_{\max}]$ which selects a subset $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\KSsubset}{\mathcal{W}^C} \KSsubset$ of KS states affected by correlations, the interacting charge density can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KSsubset}{\mathcal{W}^C} \newcommand{\KSp}{\Psi_{\nu'\kv}} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\rv}{{\bf r}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \displaystyle n(\rv)& = \sum_{\kv} \sum_{\nu \notin \KSsubset} f_{\nu\kv} \braket{\rv}{\KS} \braket{\KS}{\rv} \nonumber \\ &+ \sum_{\kv} \sum_{\nu\nu' \in \KSsubset} \braket{\rv}{\KS} N_{\nu\nu'}(\kv) \braket{\KSp}{\rv}, \label{eq:chden} \nonumber \end{align} \tag{ 22 }$$

where we use a non-diagonal density matrix

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\rmi}{{\rm i}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle N_{\nu\nu'}(\kv) = \frac{1}{\beta} \sum_{\iomn} G_{\nu\nu'}(\kv, \iomn), \nonumber \end{align} \tag{ 23 }$$

and a diagonal density matrix formed by the KS occupation numbers $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} f_{\nu\kv}$. Note that the chemical potential μ here is generally different from the KS chemical potential $\newcommand{\ks}{\mathtt{KS}} \mu_{\ks}$. The subset of KS states for the optimization of the projectors $\mathcal{W}^P$ and KS states affected by correlations $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\KSsubset}{\mathcal{W}^C} \KSsubset$ can be chosen to be the same but one does necessarily not have to do so.

From a practical point of view it is convenient to split the charge density into a DFT part and a correlation-induced part,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dft}{\mathtt{DFT}} \newcommand{\rv}{{\bf r}} \displaystyle n(\rv) = n_{\dft}(\rv) + \Delta n(\rv), \nonumber \end{align} \tag{ 24 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dft}{\mathtt{DFT}} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\rv}{{\bf r}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \displaystyle n_\dft(\rv) = \sum_{\nu\kv} f_{\nu\kv} \braket{\rv}{\KS} \braket{\KS}{\rv} \nonumber \end{align} \tag{ 25 }$$

where $\newcommand{\rv}{{\bf r}} \Delta n(\rv)$ is the correlation correction,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KSsubset}{\mathcal{W}^C} \newcommand{\KSp}{\Psi_{\nu'\kv}} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\rv}{{\bf r}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \displaystyle \Delta n(\rv) = \sum_{\stackrel{\kv}{\nu\nu' \in \KSsubset} }\braket{\rv}{\KS}\Delta N_{\nu\nu'}(\kv) \braket{\KSp}{\rv}, \nonumber \end{align} \tag{ 26 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle \Delta N_{\nu\nu'}(\kv) = N_{\nu\nu'}(\kv) - f_{\nu\kv} \delta_{\nu\nu'}, \nonumber \end{align} \tag{ 27 }$$

with the last quantity having the important property

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KSsubset}{\mathcal{W}^C} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle \sum_{\nu \in \KSsubset} \sum_{\kv} \Delta N_{\nu\nu}(\kv) = 0. \nonumber \end{align} \tag{ 28 }$$

Such a definition of the new charge density is convenient because it ensures charge neutrality between DFT iterations.

In a similar way, the total energy can be split into a DFT part calculated using the correlated charge density given by (22) and a correlation part,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KSsubset}{\mathcal{W}^C} \newcommand{\dft}{\mathtt{DFT}} \newcommand{\dc}{\textsc{dc}} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\corrspace}{\mathcal{C}} \newcommand{\eps}{\varepsilon} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle E = E_{\dft}[n] + & \sum_{\kv} \sum_{\nu \in \KSsubset} \Delta N_{\nu\nu}(\kv) \eps_{\nu\kv} \nonumber \\ &+ E_{{\rm corr}}[\hat G^\corrspace] - E_{\dc}[\hat G^\corrspace]. \label{eq:vasp_toten} \nonumber \end{align} \tag{ 29 }$$

Thus, in order to calculate the total energy and to obtain the charge density for the next KS iteration, only two quantities need to be calculated after a DMFT iteration: the density-matrix correction $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \Delta N_{\nu\nu}(\kv)$ and the interaction energy (including the double-counting term) $\newcommand{\dc}{\textsc{dc}} E_{{\rm corr}} - E_{\dc}$.

In this particular vasp implementation, we use $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} \Delta N_{\nu\nu}(\kv)$ to obtain the natural orbitals by a transformation $V$ given by diagonalizing the total correlated density matrix,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \displaystyle f'_{\nu\kv} \delta_{\nu\nu'} = \sum_{\mu\mu'}V_{\nu \mu} \left[\,f_{\mu\kv} \delta_{\mu\mu'} + \Delta N_{\mu\mu'}(\kv)\right] V^*_{\mu'\nu'}, \nonumber \end{align} \tag{ 30 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\ket}[1]{|{#1}\rangle} \displaystyle \ket{\KS'} = \sum_{\mu} V_{\nu\mu} \ket{\Psi_{\mu\kv}}, \nonumber \end{align} \tag{ 31 }$$

and the charge density and the one-electron energy are, then, obtained in the same way as in the normal KS cycle,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\rv}{{\bf r}} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \newcommand{\braket}[2]{\langle{#1}|{#2}\rangle} \displaystyle n(\rv) = \sum_{\nu\kv} f'_{\nu\kv} \braket{\rv}{\KS'} \braket{\KS'}{\rv}, \nonumber \end{align} \tag{ 32 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Tr}{{\rm Tr}} \newcommand{\ks}{\mathtt{KS}} \newcommand{\iomn}{\rmi\omega_{n}} \newcommand{\rmi}{{\rm i}} \newcommand{\KS}{\Psi_{\nu\kv}} \newcommand{\braopket}[3]{\langle{#1}|{#2}|{#3}\rangle} \newcommand{\kv}{\kvec} \newcommand{\kvec}{{\bf k}} \newcommand{\bra}[1]{\langle{#1}|} \displaystyle \frac{1}{\beta} \sum_{\iomn\kv} \Tr\{G(\kv, \iomn)& H_{\ks}(\kv)\} \nonumber \\ &=\sum_{\nu\kv} f'_{\nu\kv} \braopket{\KS'}{H_{\ks}}{\KS'}. \label{eq:tr_gh} \nonumber \end{align} \tag{ 33 }$$

Note that the DFT band energy is calculated using the original occupancies $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} f_{\nu\kv}$, and the second term in (29) represents essentially a band-energy correction which takes into account the change in the density matrix induced by correlations. The occupancies $\newcommand{\kvec}{{\bf k}} \newcommand{\kv}{\kvec} f'_{\nu\kv}$ will deviate from the usual Fermi–Dirac statistics as a result of particle fluctuations from the occupied non-interacting KS states into unoccupied states.

For the example of NiO, we first investigate the properties of the optimized projectors described in section 2. To this end, we analyze the Kohn–Sham Hamiltonian transformed to a localized basis, which will later serve as a starting point for including local correlation effects on a Hartree–Fock level.

In order to treat the Ni d states and ligand O p-states explicitly, we project the paramagnetic Kohn–Sham Hamiltonian of a two atomic unit cell using $48\times 48 \times 48$ k points onto Ni d and O p states according to (12). Therefore, we optimize the projectors in an energy window around the Fermi level given by $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \eps^P_{\min}=-3.3$ eV and $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \eps^P_{\max}=1.7$ eV (i.e. we mainly optimize the Ni d states) and orthonormalize according to (14) for all available energies. We perform a Gram–Schmidt procedure to orthogonally complement the Ni d and O p states with additional states to end up with the same number of localized states as Kohn–Sham states.

In figure 1 we analyze the orbital properties of selected eigenstates of the resulting Hamiltonian. The left panel shows the Ni d weight of a single band (displayed as red dashed line in the right panel) across a k-path from the Γ to the X point. For a small number of unoccupied bands, (total number of bands N_b  =  12), the non-optimized and optimized projectors give nearly identical results. However, by including additional unoccupied bands (N_b  =  24) the two schemes display huge differences. While the optimized projectors lead to only minor differences to the case of N_b  =  12 close to the Γ point, the non-optimized projectors lead to considerably less overall Ni d weight and a strongly discontinuous behavior for part of the k pathway. This underlines the advantage of the optimized projectors.

Zoom In Zoom Out Reset image size

The erroneous behavior of the non-optimized projectors roots in additional high energy bands having Ni-4d character. Using the standard VASP 'projectors', these high lying valence orbitals show appreciable overlap with the 4d partial waves (and 4d projectors). VASP used to lump the 4d and 3d contributions into single values, which causes the issues we have just observed. Without going into mathematical details one can easily understand this behavior. The 3d or 4d states are automatically orthogonal, because the 4d states possess an additional node inside the PAW sphere. By simply adding up the contributions from the 3d and 4d projectors pretending that they correspond to a single main quantum number, and then performing an orthogonalization, the total charge in each d spin channel is one instead of two, explaining the large reduction of the d character using the old scheme.

The failure of the unoptimized projectors could in principle be avoided by simply not including additional empty bands in the construction of the Wannier functions. However, there are important cases where this is not possible. First of all, the situation of 4d states close to the Fermi energy is realized in many late transition metals. Second, one might be interested in quantities that include transitions to unoccupied states over large energy scales, such as optical spectroscopy. In these cases, the optimized projectors give reliable results. Furthermore, from a physical point of view, the projectors should not depend strongly on the inclusion of unoccupied states.

The DFT+U approach improves the LDA and GGA description of NiO by supplementing the Ni d states by local Coulomb interaction which leads to a static local self energy in (20).

There are two ways to formulate DFT+U: first, the traditional one, where the Kohn–Sham potential is directly augmented with a potential arising from the Hartree–Fock decoupling of the local Hubbard interaction. Second, in terms of a charge self-consistent DFT+DMFT scheme, where the DMFT impurity problem is solved in the Hartree–Fock approximation called DFT+DMFT(HF) in the following. Both schemes are fully equivalent if the augmentation spaces are chosen to be strictly the same and the way spin-polarization is accounted for is the same.

However, many DFT+U implementations, including the one in vasp, work with angular-momentum-decomposed charge densities, which in general do not relate to proper projectors in the definitions of the '+U' potentials [16]. Additionally, DFT+DMFT and DFT+U for magnetic materials can work with or without spin polarization in the DFT part.

In the following, we compare DFT+U as implemented in vasp to DFT+DMFT(HF) with and without charge self-consistency for the testbed material NiO, and investigate the dependence of the electronic spectra on the double-counting.

We fix the double counting in the vasp DFT+U approach according to a prescription known as 'fully localized limit' (FLL) [32] which leads to an orbital independent and diagonal version of $\newcommand{\dc}{\textsc{dc}} \Sigma_{mm'}^\dc$ in (20), in short μ^dc. For a meaningful comparison of spectral features we choose the double counting potential in the DFT+DMFT(HF) approach such that we reproduce the size of the gap in vasp DFT+U, leading to $\newcommand{\dc}{\textsc{dc}} \mu^\dc = 62.5$ eV.

We sample the Brillouin zone of the unit cell of twice the size of the primitive one using $8\times8\times8$ k-points and parametrize the Coulomb interaction using Slater integrals given by effective Coulomb parameters U  =  8 eV and J  =  1 eV as obtained from constrained LDA calculations [28]. We perform spin polarized calculations considering antiferromagnetic ordering of the two Ni atoms in the unit cell, where for both, DFT+DMFT(HF) and vasp DFT+U, we only consider spin polarization in the Hartree–Fock part but not in the DFT part.

The respective density of states for spin up are presented in figure 2. In general, the DOS from DFT+U and DFT+DMFT(HF) are very similar: a gap of about 4 eV is opened between heavily spin-polarized Ni states in the conduction band and strongly hybridized Ni-O states in the valence band. Only details differ for the two approaches, such as a slightly larger spin polarization of Ni states and a larger O contribution at the valence band edge for DFT+DMFT(HF).

Zoom In Zoom Out Reset image size

The double-counting problem poses a serious problem when using many-body corrections in DFT+DMFT approaches predictively, since fundamental properties such as the single particle gap depend on the double counting. Here, we perform fully charge self-consistent (fcsc) as well as one-shot DFT+DMFT(HF) calculations for various fixed double-counting potentials μ^dc to investigate its influence on the spectra and on the gap.

Strictly speaking, this approach does not represent a proper double-counting scheme on first sight, since there seems to be no functional relation between the double-counting potential μ^dc and energy E_dc. However, one can motivate this approach following arguments in [33]. Instead of using the standard FLL form for the double-counting correction, one introduces an interaction parameter $U'$ in the double counting formulas, which is allowed to be slightly different from U used otherwise in the calculation. In that way, one can again relate μ^dc to a proper energy correction E_dc.

However, here we are not interested in total energies but in a systematic investigation of the influence of the double counting on the spectral properties. That is why we refrain from an explicit introduction of this parameter $U'$. Furthermore, we want to circumvent here further complicating ambiguities in usual double-counting approaches (e.g. if one should use the formal or self-consistent occupation in formulas for E_dc) [30].

The resulting total DOS for the fcsc and one-shot calculations are presented color coded in figure 3. In both cases the Ni states in the conduction band (dark blue around 5 eV) are shifted linearly towards the Fermi energy as μ^dc increases. The O/Ni states at the valence band edge are not affected strongly by the double counting. This is because, in general, the higher the Ni character of a band, the more its energy is shifted by the double-counting correction.

Zoom In Zoom Out Reset image size

In the conduction band we also find states with low spectral weight in the PAW spheres (i.e. they are neither located around the Ni nor the O atoms but in the interstitial region); their energies shift to lower values with decreasing μ^dc, i.e. in the opposite direction than the Ni state in the conduction band. For $\newcommand{\dc}{\textsc{dc}} \mu^\dc \lesssim 58$ eV in case of fcsc and $\newcommand{\dc}{\textsc{dc}} \mu^\dc \lesssim 57$ eV in the case of one-shot calculations, these states cross the Ni states. Then, the conduction band edge does not consist of Ni states, which is not in agreement with experiment. Note that the 'around-mean-field' AMF prescription lies in this regime both for one-shot and fcsc calculations and is thus not suitable for NiO. Similarly, the FLL prescription for fcsc calculations is very close to this regime.

For values of the double counting that give the correct conduction band character, the gap shrinks with increasing μ^dc. This dependence differs strongly for the one-shot and charge self-consistent calculations: the slope of the gap $\newcommand{\dc}{\textsc{dc}} {\rm d} E_g / {\rm d}\mu^\dc$ in the physical regime differs by a factor of nearly 2. Thus, the charge self-consistency does not cure but alleviates the double-counting problem by decreasing the influence of μ^dc on the spectrum.