Construction of the exchange matrix

LDA and GGA exchange-correlation functionals typically underestimate the bandgap of a material and may even predict a metal even though the material is an insulator. One way to overcome this issue is the exployment of the LDA+U method. But this approach needs material-specific parameters. Another approach is to mix LDA and GGA with the Hartree-Fock approach which typically overestimates the bandgap. The resulting hybrid functionals require the incorporation of nonlocal (NL), spin $\sigma$ dependent Hartree-Fock-like exchange

$$\begin{equation} \renewcommand{\vec}[1]{\mathbf{#1}} V_\text{x}^{\text{NL},\sigma}(\vec{r},\vec{r'}) = -\sum_\nu^\text{occ} \sum_\vec{q}^\text{BZ} \frac{\psi_{\nu,\vec{q}}^\sigma(\vec{r}) \psi_{\nu,\vec{q}}^{\sigma \ast}(\vec{r'})}{\left| \vec{r} - \vec{r'} \right|} \end{equation}$$

in the setup of the Hamiltonian. The sums hereby go over all occupied states $\psi_{\nu,\vec{q}}^\sigma$ within the Brillouin zone (BZ).

In hybrid functionals calculations with Fleur this is finally added together with a scaling factor in terms of a nonlocal exchange matrix $V_{\text{x},\vec{G},\vec{G'}}^{\text{NL},\sigma}(\vec{k})$ to the Hamiltonian matrix. However, for reasons of efficiency the computation of this matrix is first performed in terms of an auxilliary basis of Kohn-Sham orbitals determined in a previous iteration of the SCF loop. It is then brought into the final form by a basis set transformation. In detail, the matrix in terms of the auxilliary basis is

$$\begin{equation} V_{\text{x},\nu,\nu'}^{\text{NL},\sigma}(\vec{k}) = -\sum_{\nu''}^\text{occ} \sum_\vec{q}^\text{BZ} \int \int \frac{\psi_{\nu,\vec{k}}^{\sigma \ast}(\vec{r}) \psi_{\nu'',\vec{q}}^\sigma(\vec{r}) \psi_{\nu'',\vec{q}}^{\sigma \ast}(\vec{r'}) \psi_{\nu',\vec{k}}^\sigma(\vec{r'}) }{\left| \vec{r} - \vec{r'} \right|} d^3 r ~ d^3 r'. \end{equation}$$

Of course, the size of the auxilliary basis is a cutoff parameter that has to be considered when ensuring parameter convergence.

To construct this matrix another basis set is introduced: The mixed product basis (MPB) which is capable of efficiently representing wavefunction products. It consists of plane waves projected onto the interstitial region

$$\begin{equation} M_\vec{G}^\vec{k} = \frac{1}{\sqrt{\Omega}}e^{i(\vec{k}+\vec{G}) \cdot \vec{r}} \Theta(\vec{r}) \end{equation}$$

and additional functions with Bloch function character to represent the products in each MT sphere $M_{l,m,p}^{\alpha,\vec{k}}$. These functions are directly constructed by considering the products of all possible radial functions contributing in a product to the respective $(l,m)$ channel. Since the resulting set of functions is typically numerically linearly dependent, they are subsequently used to construct a basis (indexed by $p$) spanning approximately the same space. Finally these functions are used to construct the respective Bloch functions.

Indexed by $I$ overall the mixed product basis is given as $\left\lbrace M_I^\vec{k}(\vec{r}) \right\rbrace = \left\lbrace M_{l,m,p}^{\alpha,\vec{k}}, M_\vec{G}^\vec{k} \right\rbrace$.

With this basis the expression for the exchange matrix becomes

$$\begin{equation} V_{\text{x},\nu,\nu'}^{\text{NL},\sigma}(\vec{k}) = -\sum_{\nu''}^\text{occ} \sum_\vec{q}^\text{BZ} \sum_{I,J} \left\langle \psi_{\nu,\vec{k}}^{\sigma} \middle| \psi_{\nu'',\vec{k}-\vec{q}}^\sigma M_I^\vec{q} \right\rangle v_{I,J}(\vec{q}) \left\langle M_J^\vec{q} \psi_{\nu'',\vec{k}-\vec{q}}^\sigma M_I^\vec{q} \middle| \psi_{\nu',\vec{k}}^\sigma \right\rangle, \end{equation}$$

where

$$\begin{equation} v_{I,J}(\vec{q}) = \int \int \frac{\tilde{M}_I^{\vec{q} \ast}(\vec{r}) \tilde{M}_J^{\vec{q}}(\vec{r'}) }{\left| \vec{r} - \vec{r'} \right|} d^3 r ~ d^3 r' \end{equation}$$

is the Coulomb matrix expressed in terms of the mixed product basis.

The mixed product basis comes with several cutoff parameters determining its size and composition. These have to be considered when optimizing the parameter convergence.

A detailed discussion on the construction of the exchange matrix in the context of the PBE0 hybrid functional is available in M. Betzinger, C. Friedrich, and S. Blügel, Hybrid functionals within the all-electron FLAPW method: implementation and applications of PBE0, Phys. Rev. B 81, 195117 (2010).