Hybrid functionals: formalism

In hybrid functionals the exchange energy is a mixture of semilocal (SL) and nonlocal Hartree-Fock (HF) types. They can be categorized into different families according to the type of semilocal approximation (LDA, GGA, or MGGA) or the treatment of the short- and long-range parts of the exchange. A rather general formula that encompasses the different families of hybrid functionals is given by

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{hybrid}}=a_{\mathrm{SR}} E_{\mathrm{x,SR}}^{\mathrm{HF}}(\mu) + a_{\mathrm{LR}} E_{\mathrm{x,LR}}^{\mathrm{HF}}(\mu) + (1-a_{\mathrm{SR}})E_{\mathrm{x,SR}}^{\mathrm{SL}}(\mu) + (1-a_{\mathrm{LR}})E_{\mathrm{x,LR}}^{\mathrm{SL}}(\mu) + E_{\mathrm{c}}^{\mathrm{SL}} }[/math]

where

• [math]\displaystyle{ a_{\mathrm{SR}} }[/math] and [math]\displaystyle{ a_{\mathrm{LR}} }[/math] are the mixing parameters (fraction of HF exchange) at short and long range, respectively.
• [math]\displaystyle{ \mu }[/math] is the screening parameter that determines the separation between short range (SR) and long range (LR).

The SR and LR components of the full-range [math]\displaystyle{ E_{\mathrm{x}}^{\mathrm{SL}} }[/math] and [math]\displaystyle{ E_{\mathrm{x}}^{\mathrm{HF}} }[/math] exchange energies are constructed such that at all values of [math]\displaystyle{ \mu }[/math]

• [math]\displaystyle{ E_{\mathrm{x}}^{\mathrm{HF}}=E_{\mathrm{x,SR}}^{\mathrm{HF}}(\mu)+E_{\mathrm{x,LR}}^{\mathrm{HF}}(\mu) }[/math]
• [math]\displaystyle{ E_{\mathrm{x}}^{\mathrm{SL}}=E_{\mathrm{x,SR}}^{\mathrm{SL}}(\mu)+E_{\mathrm{x,LR}}^{\mathrm{SL}}(\mu) }[/math]

The HF exchange energy (full-range, SR, or LR) is given by

[math]\displaystyle{ E_{\mathrm{x,(SR/LR)}}^{\rm HF}(\mu)= -\frac{1}{2}\sum_{n\mathbf{k},m\mathbf{q}} f_{n\mathbf{k}} f_{m\mathbf{q}} \int \int d^3\mathbf{r} d^3\mathbf{r}' v(\mu,|\mathbf{r}-\mathbf{r}'|) \psi_{n\mathbf{k}}^{*}(\mathbf{r})\psi_{m\mathbf{q}}^{*}(\mathbf{r}') \psi_{n\mathbf{k}}(\mathbf{r}')\psi_{m\mathbf{q}}(\mathbf{r}) }[/math]

with [math]\displaystyle{ \{\psi_{n\mathbf{k}}(\mathbf{r})\} }[/math] being the set of one-electron Bloch states of the system, and [math]\displaystyle{ \{f_{n\mathbf{k}}\} }[/math] the corresponding set of (possibly fractional) occupational numbers. The sums over [math]\displaystyle{ {\bf k} }[/math] and [math]\displaystyle{ {\bf q} }[/math] run over all k-points chosen to sample the Brillouin zone, whereas the sums over [math]\displaystyle{ m }[/math] and [math]\displaystyle{ n }[/math] run over all bands at these k-points. The corresponding nonlocal HF potential is given by

[math]\displaystyle{ V_{\mathrm{x,(SR/LR)}}^{\mathrm{HF}}\left(\mu,\mathbf{r},\mathbf{r}'\right)= -\sum_{m\mathbf{q}}f_{m\mathbf{q}}v(\mu,|\mathbf{r}-\mathbf{r}'|)\psi_{m\mathbf{q}}^{*}(\mathbf{r}')\psi_{m\mathbf{q}}(\mathbf{r}) }[/math]

The non-multiplicative form of the HF potential is such that hybrid functionals are implemented within the generalized KS scheme. Thus, the total energy is minimized with respect to the orbitals (instead of the electron density), which means that the HF exchange leads to a nonlocal operator as in the Hartree-Fock-Roothaan theory.

Expressions of Hartree-Fock energy and potential for the plane-wave basis set

Using the decomposition of the Bloch states [math]\displaystyle{ \psi_{m\mathbf{q}} }[/math] in plane waves,

[math]\displaystyle{ \psi_{m\mathbf{q}}(\mathbf{r})= \frac{1}{\sqrt{\Omega}} \sum_\mathbf{G}C_{m\mathbf{q}}(\mathbf{G})e^{i(\mathbf{q}+\mathbf{G}) \cdot \mathbf{r}} }[/math]

the HF exchange potential can be written as

[math]\displaystyle{ V_{\mathrm{x,(SR/LR)}}^{\mathrm{HF}}\left(\mu,\mathbf{r},\mathbf{r}'\right)= \sum_{\mathbf{k}}\sum_{\mathbf{G}\mathbf{G}'} e^{i(\mathbf{k}+\mathbf{G})\cdot\mathbf{r}} V_{\mathbf{k}}\left(\mu, \mathbf{G},\mathbf{G}'\right) e^{-i(\mathbf{k}+\mathbf{G}')\cdot\mathbf{r}'} }[/math]

where

[math]\displaystyle{ V_{\mathbf{k}}\left(\mu, \mathbf{G},\mathbf{G}'\right)= \langle \mathbf{k}+\mathbf{G} | V_{\mathrm{x,(SR/LR)}}^{\mathrm{HF}} | \mathbf{k}+\mathbf{G}'\rangle }[/math]

Types of potentials

For most hybrid functionals proposed in the literature, the interelectronic Coulomb potential [math]\displaystyle{ v(\mu,|\mathbf{r}-\mathbf{r}'|) }[/math] is one of these types ([math]\displaystyle{ r=|\mathbf{r}-\mathbf{r}'| }[/math]):

• Full range (bare Colomb potential):

[math]\displaystyle{ v^{\mathrm{bare}}(r)=\frac{1}{r} }[/math]

[math]\displaystyle{ V_{\mathbf{k}}^{\mathrm{bare}}\left( \mathbf{G},\mathbf{G}'\right)= -\frac{4\pi}{\Omega} \sum_{m\mathbf{q}}f_{m\mathbf{q}}\sum_{\mathbf{G}''} \frac{C^*_{m\mathbf{q}}(\mathbf{G}'-\mathbf{G}'') C_{m\mathbf{q}}(\mathbf{G}-\mathbf{G}'')} {|\mathbf{k}-\mathbf{q}+\mathbf{G}''|^2} }[/math]

• Short-range with error function:

[math]\displaystyle{ v_{\mathrm{SR}}^{\mathrm{erf}}(\mu,r)=\frac{\mathrm{erfc}(\mu r)}{r} }[/math]

[math]\displaystyle{ \begin{align} V_{\mathbf{k},\mathrm{SR}}^{\mathrm{erf}}\left(\mu, \mathbf{G},\mathbf{G}'\right)= -\frac{4\pi}{\Omega} \sum_{m\mathbf{q}}f_{m\mathbf{q}}\sum_{\mathbf{G}''} \frac{C^*_{m\mathbf{q}}(\mathbf{G}'-\mathbf{G}'') C_{m\mathbf{q}}(\mathbf{G}-\mathbf{G}'')} {|\mathbf{k}-\mathbf{q}+\mathbf{G}''|^2} \left( 1-e^{-|\mathbf{k}-\mathbf{q}+\mathbf{G}''|^2 /(4\mu^2)} \right) \end{align} }[/math]

Short- and long-range potentials using the error or exponential screening compared to the bare potential.

• Short-range with exponential function:

[math]\displaystyle{ v_{\mathrm{SR}}^{\mathrm{exp}}(\mu,r)=\frac{e^{-\mu r}}{r} }[/math]

[math]\displaystyle{ V_{\mathbf{k},\mathrm{SR}}^{\mathrm{exp}}\left(\mu, \mathbf{G},\mathbf{G}'\right)= -\frac{4\pi}{\Omega} \sum_{m\mathbf{q}}f_{m\mathbf{q}}\sum_{\mathbf{G}''} \frac{C^*_{m\mathbf{q}}(\mathbf{G}'-\mathbf{G}'') C_{m\mathbf{q}}(\mathbf{G}-\mathbf{G}'')} {|\mathbf{k}-\mathbf{q}+\mathbf{G}''|^2 + \mu^2} }[/math]

The corresponding long-range potentials are given by [math]\displaystyle{ v_{\mathrm{LR}}^{\mathrm{erf/exp}}=v^{\mathrm{bare}}-v_{\mathrm{SR}}^{\mathrm{erf/exp}} }[/math].

In VASP, these expressions are implemented within the PAW formalism.

The families of hybrid functionals implemented in VASP are listed below along with examples, whose corresponding INCAR files can be found at the page list of hybrid functionals.

Important: The screening [math]\displaystyle{ \mu }[/math] (HFSCREEN tag) can be used only when the semilocal functional is PBE, PBEsol, or LDA (GGA=PE, PS, or CA, respectively). The other GGA and METAGGA functionals have no screened version available in VASP.

Families of hybrid functionals

HF exchange at full range

There is no range separation, i.e. the same fraction of HF exchange is applied at full range, [math]\displaystyle{ a_{\mathrm{SR}}=a_{\mathrm{LR}}=a }[/math] (AEXX tag):

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{hybrid}}=a E_{\mathrm{x}}^{\mathrm{HF}} + (1-a)E_{\mathrm{x}}^{\mathrm{SL}} + E_{\mathrm{c}}^{\mathrm{SL}} }[/math]

Mind: These functionals are set with LHFCALC=.TRUE. By default AEXX=0.25, but can be set to another value. The semilocal part can be of the LDA, GGA or MGGA type.

These are the original and most simple forms of hybrid functionals. Two examples, PBE0 and B3LYP, are given below.

• PBE0:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{PBE0}}=\frac{1}{4} E_{\mathrm{x}}^{\mathrm{HF}} + \frac{3}{4} E_{\mathrm{x}}^{\mathrm{PBE}} + E_{\mathrm{c}}^{\mathrm{PBE}} }[/math]

It is based on the PBE GGA functional and [math]\displaystyle{ a=1/4 }[/math].

• B3LYP, well known and popular amongst quantum chemists:

[math]\displaystyle{ \begin{align} E_{\mathrm{x}}^{\mathrm{B3LYP}} &=0.8 E_{\mathrm{x}}^{\mathrm{LDA}}+ 0.2 E_{\mathrm{x}}^{\mathrm{HF}} + 0.72 (E_{\mathrm{x}}^{\mathrm{B88}}-E_{\mathrm{x}}^{\mathrm{LDA}}) + 0.19 E_{\mathrm{c}}^{\mathrm{VWN3}}+ 0.81 E_{\mathrm{c}}^{\mathrm{LYP}} \end{align} }[/math]

The exchange part consists of 80% of LDA exchange plus 20% of HF exchange, and 72% of the gradient corrections of the B88 GGA functional. The correlation consists of 81% of LYP correlation energy, which contains a LDA and a GGA part, and 19% of the LDA Vosko-Wilk-Nusair correlation functional III, which was fitted to the correlation energy in the random phase approximation of the homogeneous electron gas.

HF exchange at short range (error-function screening)

The HF exchange is used only at short-range (the long-range part is fully semilocal, [math]\displaystyle{ a_{\mathrm{LR}}=0 }[/math]) and the screening is done with the error function:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{hybrid}}=a_{\mathrm{SR}} E_{\mathrm{x,SR}}^{\mathrm{HF}}(\mu) + (1-a_{\mathrm{SR}})E_{\mathrm{x,SR}}^{\mathrm{SL}}(\mu) + E_{\mathrm{x,LR}}^{\mathrm{SL}}(\mu) + E_{\mathrm{c}}^{\mathrm{SL}} }[/math]

The mixing [math]\displaystyle{ a_{\mathrm{SR}} }[/math] and screening [math]\displaystyle{ \mu }[/math] are controlled by the AEXX and HFSCREEN tags, respectively.

The most popular range-separated functional, HSE, is given below.

• HSE03 and HSE06:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{HSE}}= \frac{1}{4}E_{\mathrm{x,SR}}^{\mathrm{HF}}(\mu) + \frac{3}{4} E_{\mathrm{x,SR}}^{\mathrm{PBE}}(\mu) + E_{\mathrm{x,LR}}^{\mathrm{PBE}}(\mu) + E_{\mathrm{c}}^{\mathrm{PBE}} }[/math]

They are based on the PBE GGA functional and [math]\displaystyle{ a_{\mathrm{SR}}=1/4 }[/math]. It has been shown that the optimum [math]\displaystyle{ \mu }[/math], controlling the range separation is approximately 0.2-0.3 Å^-1. HSE03 and HSE06 correspond to HFSCREEN=0.3 and 0.2, respectively. Note that the two limit cases of HSE are PBE0 at [math]\displaystyle{ \mu=0 }[/math] and PBE at [math]\displaystyle{ \mu\rightarrow\infty }[/math].

HF exchange at short range and long range (error-function screening) with different mixings

The fractions of HF exchange at short and long range ([math]\displaystyle{ a_{\mathrm{SR}} }[/math] and [math]\displaystyle{ a_{\mathrm{LR}} }[/math], respectively) can be different:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{hybrid}}=a_{\mathrm{SR}} E_{\mathrm{x,SR}}^{\mathrm{HF}}(\mu) + a_{\mathrm{LR}} E_{\mathrm{x,LR}}^{\mathrm{HF}}(\mu) + (1-a_{\mathrm{SR}})E_{\mathrm{x,SR}}^{\mathrm{SL}}(\mu) + (1-a_{\mathrm{LR}})E_{\mathrm{x,LR}}^{\mathrm{SL}}(\mu) + E_{\mathrm{c}}^{\mathrm{SL}} }[/math]

These functionals are selected with LMODELHF=.TRUE. The mixings [math]\displaystyle{ a_{\mathrm{LR}} }[/math] and [math]\displaystyle{ a_{\mathrm{SR}} }[/math] are controlled by the AEXX and BEXX tags, respectively, and the screening [math]\displaystyle{ \mu }[/math] by the HFSCREEN tag.

Important: The possibility to set [math]\displaystyle{ a_{\mathrm{SR}} }[/math] with BEXX within the LMODELHF=.TRUE. method was introduced in VASP.6.6.0. Until VASP.6.5.1 [math]\displaystyle{ a_{\mathrm{SR}} }[/math] was fixed to 1 and could not be changed.

This functional form has been used in the context of dielectric-dependent hybrids. Examples are provided below.

• RS-DDH:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{RS-DDH}}=a_{\mathrm{SR}} E_{\mathrm{x,SR}}^{\mathrm{HF}}(\mu) + a_{\mathrm{LR}} E_{\mathrm{x,LR}}^{\mathrm{HF}}(\mu) + (1-a_{\mathrm{SR}})E_{\mathrm{x,SR}}^{\mathrm{PBE}}(\mu) + (1-a_{\mathrm{LR}})E_{\mathrm{x,LR}}^{\mathrm{PBE}}(\mu) + E_{\mathrm{c}}^{\mathrm{PBE}} }[/math]

where [math]\displaystyle{ a_{\mathrm{SR}}=1/4 }[/math] and [math]\displaystyle{ a_{\mathrm{LR}}=\varepsilon^{-1} }[/math] is chosen as the inverse of the dielectric constant [math]\displaystyle{ \varepsilon^{-1} }[/math].

• DD-RSH-CAM, DSH:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{DD-RSH-CAM}}=E_{\mathrm{x,SR}}^{\mathrm{HF}}(\mu) + a_{\mathrm{LR}} E_{\mathrm{x,LR}}^{\mathrm{HF}}(\mu) + (1-a_{\mathrm{LR}})E_{\mathrm{x,LR}}^{\mathrm{PBE}}(\mu) + E_{\mathrm{c}}^{\mathrm{PBE}} }[/math]

where [math]\displaystyle{ a_{\mathrm{SR}}=1 }[/math] and [math]\displaystyle{ a_{\mathrm{LR}}=\varepsilon^{-1} }[/math] is chosen as the inverse of the dielectric constant [math]\displaystyle{ \varepsilon^{-1} }[/math].

HF exchange at long range (error-function screening)

The HF exchange is used only at long-range (the short-range part is fully semilocal, [math]\displaystyle{ a_{\mathrm{SR}}=0 }[/math]) and the screening is done with the error function:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{hybrid}}=a_{\mathrm{LR}} E_{\mathrm{x,LR}}^{\mathrm{HF}}(\mu) + E_{\mathrm{x,SR}}^{\mathrm{SL}}(\mu) + (1-a_{\mathrm{LR}})E_{\mathrm{x,LR}}^{\mathrm{SL}}(\mu) + E_{\mathrm{c}}^{\mathrm{SL}} }[/math]

These functionals are selected with LRHFCALC=.TRUE. The mixing [math]\displaystyle{ a_{\mathrm{LR}} }[/math] and screening [math]\displaystyle{ \mu }[/math] are controlled by the AEXX and HFSCREEN tags, respectively.

Mind: LRHFCALC=.TRUE. automatically sets AEXX=1. However, AEXX can be set to another value.

Important: When AEXX=1 (the default for LRHFCALC=.TRUE.), the correlation [math]\displaystyle{ E_{\mathrm{c}}^{\mathrm{SL}} }[/math] is not included. However, it can be included by setting ALDAC=1.0 and AGGAC=1.0.

Long-range hybrid functionals are more popular in molecular chemistry, where a proper decay of the exchange-correlation potential at long range far from the nuclei may be important, and thus less useful for bulk solids. Examples belonging to this class of functionals are:

• RSHXLDA and RSHXPBE:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{RSHXLDA}} = E_{\mathrm{x,LR}}^{\mathrm{HF}}(\mu) + E_{\mathrm{x,SR}}^{\mathrm{LDA}}(\mu) + E_{\mathrm{c}}^{\mathrm{LDA}} }[/math]

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{RSHXPBE}} = E_{\mathrm{x,LR}}^{\mathrm{HF}}(\mu) + E_{\mathrm{x,SR}}^{\mathrm{PBE}}(\mu) + E_{\mathrm{c}}^{\mathrm{PBE}} }[/math]

When LDA is chosen, a value of [math]\displaystyle{ \mu=0.75 }[/math] Å^-1 is recommended for solids.

HF exchange at short range (exponential screening)

The HF exchange is used only at short-range (the long-range part is fully semilocal, [math]\displaystyle{ a_{\mathrm{LR}}=0 }[/math]) and the screening is done with the exponential function:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{hybrid}}=a_{\mathrm{SR}} E_{\mathrm{x,SR}}^{\mathrm{HF}}(\mu) + (1-a_{\mathrm{SR}})E_{\mathrm{x,SR}}^{\mathrm{SL}}(\mu) + E_{\mathrm{x,LR}}^{\mathrm{SL}}(\mu) + E_{\mathrm{c}}^{\mathrm{SL}} }[/math]

The exponential screening, also called Thomas-Fermi (TF) screening, is activated by setting LTHOMAS=.TRUE.. The mixing [math]\displaystyle{ a_{\mathrm{SR}} }[/math] and screening [math]\displaystyle{ \mu=k_{\rm TF} }[/math] are controlled by the AEXX and HFSCREEN tags, respectively.

Mind: LTHOMAS=.TRUE. automatically sets AEXX=1. However, AEXX can be set to another value.

Important: When AEXX=1 (the default for LTHOMAS=.TRUE.), the correlation [math]\displaystyle{ E_{\mathrm{c}}^{\mathrm{SL}} }[/math] is not included. However, it can be included by setting ALDAC=1.0 and AGGAC=1.0. This functional should be used only with LDA (GGA=CA).

The sX-LDA functional, which uses [math]\displaystyle{ a_{\mathrm{SR}}=1 }[/math], is probably the first hybrid using an exponential screening:

• sX-LDA:

[math]\displaystyle{ E_{\mathrm{xc}}^{\mathrm{sX-LDA}} = E_{\mathrm{x,SR}}^{\mathrm{HF}}(\mu) + E_{\mathrm{x,LR}}^{\mathrm{LDA}}(\mu) + E_{\mathrm{c}}^{\mathrm{LDA}} }[/math]

For typical semiconductors, a Thomas-Fermi screening length [math]\displaystyle{ \mu=k_{\rm TF} }[/math] of about 1.8 Å^-1 yields reasonable band gaps. In principle, however, the Thomas-Fermi screening length depends on the valence-electron density. VASP determines [math]\displaystyle{ k_{\rm TF} }[/math] from the number of valence electrons (read from the POTCAR file) and the volume (leading to an average density [math]\displaystyle{ \bar{n} }[/math]) and writes the corresponding value of [math]\displaystyle{ k_{\rm TF}=\sqrt{4\bar{k}_{\rm F}/\pi} }[/math], where [math]\displaystyle{ \bar{k}_{\rm F}=(3\pi^2\bar{n})^{1/3} }[/math] to the OUTCAR file (note that this value is only printed for information and is not used during the calculation):

 Thomas-Fermi vector in A             =   2.00000

Since VASP counts the semi-core states and d-states as valence electrons, although these states do not contribute to the screening, the values reported by VASP are often not recommended.

Another important detail concerns the implementation of the local LDA part in VASP. Literature [see Eqs. (3.10), (3.14), and (3.15) in Ref. ] suggests to use in the enhancement factor [math]\displaystyle{ F(z) }[/math] a position-independent variable [math]\displaystyle{ z=k_{\rm TF}/\bar{k}_{\rm F} }[/math] where [math]\displaystyle{ \bar{k}_{\rm F} }[/math] is as defined above but using the average density [math]\displaystyle{ \bar{n} }[/math] in the unit cell. However, implemented in VASP is a position-dependent variable [math]\displaystyle{ z({\bf r})=k_{\rm TF}/k_{\rm F}({\bf r}) }[/math], where [math]\displaystyle{ k_{\rm F}({\bf r})=(3\pi^2 n({\bf r}))^{1/3} }[/math] is the Fermi wave vector calculated with the local density [math]\displaystyle{ n({\bf r}) }[/math], while the constant [math]\displaystyle{ k_{\rm TF} }[/math] is set by HFSCREEN.

Related tags and articles

AEXX, BEXX, ALDAX, ALDAC, AGGAX, AGGAC, AMGGAX, AMGGAC, HFSCREEN, LHFCALC, LMODELHF, LTHOMAS, LRHFCALC, List of hybrid functionals, Downsampling of the Hartree-Fock operator, Coulomb singularity