Requests for technical support from the VASP group should be posted in the VASP-forum.

# Hartree-Fock and HF/DFT hybrid functionals

## Fock exchange

The non-local Fock exchange energy, Ex, can be written as

${\displaystyle E_{\mathrm {x} }=-{\frac {e^{2}}{2}}\sum _{n{\mathbf {k} },m{\mathbf {q} }}f_{n{\mathbf {k} }}f_{m{\mathbf {q} }}\times \int \int d^{3}{\mathbf {r} }d^{3}{\mathbf {r} }'{\frac {\psi _{n{\mathbf {k} }}^{*}({\mathbf {r} })\psi _{m{\mathbf {q} }}^{*}({\mathbf {r} }')\psi _{n{\mathbf {k} }}({\mathbf {r} }')\psi _{m{\mathbf {q} }}({\mathbf {r} })}{\vert {\mathbf {r} }-{\mathbf {r} }'\vert }}}$

with ${\displaystyle \{\psi _{n{\mathbf {k} }}({\mathbf {r} })\}}$ being the set of one-electron Bloch states of the system, and ${\displaystyle \{f_{n{\mathbf {k} }}\}}$ the corresponding set of (possibly fractional) occupational numbers. The sums over k and q run over all k-points chosen to sample the Brillouin zone (BZ), whereas the sums over m and n run over all bands at these k-points.

The corresponding non-local Fock potential is given by

${\displaystyle V_{\mathrm {x} }\left({\mathbf {r} },{\mathbf {r} }'\right)=-{\frac {e^{2}}{2}}\sum _{m{\mathbf {q} }}f_{m{\mathbf {q} }}{\frac {\psi _{m{\mathbf {q} }}^{*}({\mathbf {r} }')\psi _{m{\mathbf {q} }}({\mathbf {r} })}{\vert {\mathbf {r} }-{\mathbf {r} }'\vert }}=-{\frac {e^{2}}{2}}\sum _{m{\mathbf {q} }}f_{m{\mathbf {q} }}e^{-i{\mathbf {q} }\cdot {\mathbf {r} }'}{\frac {u_{m{\mathbf {q} }}^{*}({\mathbf {r} }')u_{m{\mathbf {q} }}({\mathbf {r} })}{\vert {\mathbf {r} }-{\mathbf {r} }'\vert }}e^{i{\mathbf {q} }\cdot {\mathbf {r} }}}$

where ${\displaystyle u_{m{\mathbf {q} }}({\mathbf {r} })}$ is the cell periodic part of the Bloch state, ${\displaystyle \psi _{n{\mathbf {q} }}({\mathbf {r} })}$, at k-point, q, with band index m.

Using the decomposition of the Bloch states, ${\displaystyle \psi _{m{\mathbf {q} }}}$, in plane waves,

${\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} }}}$

the Fock exchange potential may be written as

${\displaystyle V_{\mathrm {x} }\left({\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({\mathbf {G} },{\mathbf {G} }'\right)e^{-i({\mathbf {k} }+{\mathbf {G} }')\cdot {\mathbf {r} }'}}$

where

${\displaystyle V_{\mathbf {k} }\left({\mathbf {G} },{\mathbf {G} }'\right)=\langle {\mathbf {k} }+{\mathbf {G} }|V_{\mathrm {x} }|{\mathbf {k} }+{\mathbf {G} }'\rangle =-{\frac {4\pi e^{2}}{\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}}}}$

is the representation of the Fock potential in reciprocal space.

In VASP these expressions are implemented within the PAW formalism.[1]

## Hybrid functionals

Hartree-Fock/DFT hybrid functionals admix a certain amount of Fock exchange to (a part of) a local or semilocal density functional.

### PBE0

The PBE0 functional is defined as:[2]

${\displaystyle E_{\mathrm {xc} }^{\mathrm {PBE0} }={\frac {1}{4}}E_{\mathrm {x} }+{\frac {3}{4}}E_{\mathrm {x} }^{\mathrm {PBE} }+E_{\mathrm {c} }^{\mathrm {PBE} }}$

where ExPBE and EcPBE denote the exchange and correlation parts of the PBE density functional, respectively.

### B3LYP

A well known hybrid functional popular amongst quantum chemists is the B3LYP functional:

{\displaystyle {\begin{aligned}E_{\mathrm {x} }^{\mathrm {B3LYP} }&=0.8E_{\mathrm {x} }^{\mathrm {LDA} }+0.2E_{\mathrm {x} }+0.72\Delta E_{\mathrm {x} }^{\mathrm {B88} }\\E_{\mathrm {c} }^{\mathrm {B3LYP} }&=0.19E_{\mathrm {c} }^{\mathrm {VWN3} }+0.81E_{\mathrm {c} }^{\mathrm {LYP} }\end{aligned}}}

where ExB3LYP and EcB3LYP are the B3LYP exchange and correlation energy contributions, respectively. ExB3LYP consists of 80% of LDA exchange plus 20% of non-local Fock exchange, and 72% of the gradient corrections of the Becke88 exchange functional. EcB3LYP consists of 81% of LYP correlation energy, which contains a local and a semilocal (gradient dependent) part, and 19% of the (local) Vosko-Wilk-Nusair correlation functional III, which is fitted to the correlation energy in the random phase approximation RPA of the homogeneous electron gas.

## Range separated hybrid functionals

### HSE

In the range separated HSE03[3][4][5] and HSE06[6] hybrid functionals the slowly decaying long-ranged part of the Fock exchange interaction is replaced by the corresponding part of the PBE density functional counterpart. The resulting expression for the exchange-correlation energy is given by:

${\displaystyle E_{\mathrm {xc} }^{\mathrm {HSE} }={\frac {1}{4}}~E_{\mathrm {x} }^{\mathrm {SR} }(\mu )+{\frac {3}{4}}E_{\mathrm {x} }^{\mathrm {PBE,SR} }(\mu )+E_{\mathrm {x} }^{\mathrm {PBE,LR} }(\mu )+E_{\mathrm {c} }^{\mathrm {PBE} }.}$

As can be seen above, the separation of the electron-electron interaction into a short- and long-ranged part, labeled SR and LR respectively, is realized only in the exchange interactions. Electronic correlation is represented by the corresponding part of the PBE density functional.

The decomposition of the Coulomb kernel is obtained using the following construction:

${\displaystyle {\frac {1}{r}}=S_{\mu }(r)+L_{\mu }(r)={\frac {{\mathrm {erfc} }(\mu r)}{r}}+{\frac {{\mathrm {erf} }(\mu r)}{r}}}$

where r =|r-r'|, and ${\displaystyle \mu }$ (=HFSCREEN) is the parameter that defines the range-separation, and is related to a characteristic distance, ${\displaystyle 2/\mu }$, at which the short-range interactions become negligible.

Note: It has been shown that the optimum μ, controlling the range separation is approximately 0.2-0.3 Å-1.[3][4][5][6] To select the HSE06 functional you need to select (HFSCREEN=0.2).

Using the decomposed Coulomb kernel one may straightforwardly rewrite the non-local Fock exhange energy:

${\displaystyle E_{\mathrm {x} }^{\rm {SR}}(\mu )=-{\frac {e^{2}}{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} }'{\frac {{\mathrm {erfc} }(\mu |{\mathbf {r} }-{\mathbf {r} }'|)}{|{\mathbf {r} }-{\mathbf {r} }'|}}\times \psi _{n{\mathbf {k} }}^{*}({\mathbf {r} })\psi _{m{\mathbf {q} }}^{*}({\mathbf {r} }')\psi _{n{\mathbf {k} }}({\mathbf {r} }')\psi _{m{\mathbf {q} }}({\mathbf {r} }).}$

The representation of the corresponding short-ranged Fock potential in reciprocal space is given by

{\displaystyle {\begin{aligned}V_{\mathbf {k} }^{\mathrm {SR} }\left({\mathbf {G} },{\mathbf {G} }'\right)&=\langle {\mathbf {k} }+{\mathbf {G} }|V_{x}^{\rm {SR}}[\mu ]|{\mathbf {k} }+{\mathbf {G} }'\rangle \\&=-{\frac {4\pi e^{2}}{\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}}}\times \left(1-e^{-|{\mathbf {k} }-{\mathbf {q} }+{\mathbf {G} }''|^{2}/4\mu ^{2}}\right).\end{aligned}}}

Clearly, the only difference to the reciprocal space representation of the complete Fock exchange potential is the second factor in the summand above, that represents the complementary error function in reciprocal space.

The short-ranged PBE exchange energy and potential, and their long-ranged counterparts, are arrived at using the same decomposition, in accordance with Heyd et al.[3] It is easily seen that the long-range term in the decomposed Coulomb kernel becomes zero for μ=0, and the short-range contribution then equals the full Coulomb operator, whereas for μ→∞ it is the other way around. Consequently, the two limiting cases of the HSE functional are a true PBE0 functional for μ=0, and a pure PBE calculation for μ→∞.

### Thomas-Fermi screening

In the case of Thomas-Fermi screening the Coulomb kernel is again decomposed in a short-range and a long-range part.[7][8][9] This decomposition can be conveniently written in reciprocal space:

${\displaystyle {\frac {4\pi e^{2}}{|{\mathbf {G} }|^{2}}}=S_{\mu }(|{\mathbf {G} }|)+L_{\mu }(|{\mathbf {G} }|)={\frac {4\pi e^{2}}{|{\mathbf {G} }|^{2}+k_{\mathrm {TF} }^{2}}}+\left({\frac {4\pi e^{2}}{|{\mathbf {G} }|^{2}}}-{\frac {4\pi e^{2}}{|{\mathbf {G} }|^{2}+k_{\mathrm {TF} }^{2}}}\right)}$

where kTF (=HFSCREEN) is the Thomas-Fermi screening length. For typical semiconductors, a Thomas-Fermi screening length 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 this parameter from the number of valence electrons (read from the POTCAR file) and the volume and writes the corresponding value to the OUTCAR file:

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

Another important detail concerns the implementation of the density functional part in the screened exchange case. Literature suggests that a global enhancement factor z (see Eq. 3.15)[8] should be used, whereas VASP implements a local density dependent enhancement factor z=kTF/k , where k is the Fermi wave vector corresponding to the local density (and not the average density as suggested Seidl et al.[8]). The VASP implementation is in the spirit of the local density approximation.

## Downsampling

Consider the description of a certain bulk system, using a supercell made up of N primitive cells, in such a way that, {A i' }, the lattice vectors of the supercell are given by A i' =n iA i (i=1,2,3), where {A i} are the lattice vectors of the primitive cell. Let Rmax=2/μ be the distance for which

${\displaystyle {\frac {{\mathrm {erfc} }\left(\mu |{\mathbf {r} }-{\mathbf {r} }'|\right)}{|{\mathbf {r} }-{\mathbf {r} }'|}}\approx 0\quad {\mathrm {for} }\quad {|{\mathbf {r} }-{\mathbf {r} }'|}>R_{\mathrm {max} }.}$

When the nearest neighbour distance between the periodically repeated images of the supercell RNN > 2Rmax (i.e. RNN > 4/μ), the short-ranged Fock potential can be represented exactly, sampling the BZ at the Γ-point only, i.e.,

${\displaystyle V_{x}^{\mathrm {SR} }\left({\mathbf {r} },{\mathbf {r} }'\right)=-{\frac {e^{2}}{2}}\sum _{m}f_{m{\Gamma }}u_{m{\Gamma }}^{*}({\mathbf {r} }')u_{m{\Gamma }}({\mathbf {r} }){\frac {{\mathrm {erfc} }(\mu |{\mathbf {r} }-{\mathbf {r} }'|)}{\vert {\mathbf {r} }-{\mathbf {r} }'\vert }}.}$

This is equivalent to a representation of the bulk system using the primitive cell and a n1×n2×n3 sampling of the BZ,

${\displaystyle V_{x}^{\mathrm {SR} }\left({\mathbf {r} },{\mathbf {r} }'\right)=-{\frac {e^{2}}{2}}\sum _{m'{\mathbf {q} }}f_{m'{\mathbf {q} }}e^{-i{\mathbf {q} }\cdot {\mathbf {r} }'}u_{m'{\mathbf {q} }}^{*}({\mathbf {r} }')u_{m'{\mathbf {q} }}({\mathbf {r} })e^{i{\mathbf {q} }\cdot {\mathbf {r} }}\times {\frac {{\mathrm {erfc} }(\mu |{\mathbf {r} }-{\mathbf {r} }'|)}{\vert {\mathbf {r} }-{\mathbf {r} }'\vert }}}$

where the set of q vectors is given by

${\displaystyle \{{\mathbf {q} }\}=\{i{\mathbf {G} }_{1}+j{\mathbf {G} }_{2}+k{\mathbf {G} }_{3}\},}$

for i =1,..,n 1, j =1,..,n 2, and k =1,..,n 3, with G1,2,3 being the reciprocal lattice vectors of the supercell.

In light of the above it is clear that the number of q-points needed to represent the short-ranged Fock potential decreases with decreasing Rmax (i.e., with increasing μ). Furthermore, one should realize that the maximal range of the exchange interactions is not only limited by the erfc(μ|r-r'|)/|r-r'| kernel, but depends on the extend of the spatial overlap of the orbitals as well (this can easily be shown for the Fock exchange energy when one adopts a Wannier representation of the orbitals); Rmax=2/μ (as defined above), therefore, provides an upper limit for the range of the exchange interactions, consistent with maximal spatial overlap of the orbitals.

It is thus well conceivable that the situation arises where the short-ranged Fock potential may be represented on a considerably coarser mesh of points in the BZ than the other contributions to the Hamiltonian. To take advantage of this situation one may, for instance, restrict the sum over q in the short range exchange potential to a subset, {qk}, of the full (N 1×N 2×N 3) k-point set, {k}, for which the following holds

${\displaystyle {\mathbf {q_{k}} }={\mathbf {b} }_{1}{\frac {n_{1}C_{1}}{N_{1}}}+{\mathbf {b} }_{2}{\frac {n_{2}C_{2}}{N_{2}}}+{\mathbf {b} }_{3}{\frac {n_{3}C_{3}}{N_{3}}},\quad (n_{i}=0,..,N_{i}-1)}$

where b1,2,3 are the reciprocal lattice vectors of the primitive cell, and C i is the integer grid reduction factor along reciprocal lattice direction b i. This leads to a reduction in the computational workload by a factor:

${\displaystyle {\frac {1}{C_{1}C_{2}C_{3}}}}$

Note: From the above one should not get the impression that the grid reduction can only be used (or is useful) only in conjunction with range separated functionals (e.g. HSE03/HSE06). It can be applied, for instance, in the PBE0 and pure HF cases as well, although from the above it might be clear that the range separated functionals, in general, will allow for a larger reduction of the grid.[10]

### Caveat: when one should not use downsampling

In metallic systems, downsampling the exact exchange potential (NKRED, NKREDX, NKREDY, and/or NKREDZ ≠ 1) must be used with great care, and results might be wrong, if downsampling is applied. Problematic cases include electron or hole doped semiconductors or insulators. If two electrons are added to a bulk TiO2 cell containing 72 atoms, and calculations are performed using 2×2×2 k-points, the following results are obtained for the one-electron energies and occupancies with and without NKRED=2 (LHFCALC=.TRUE. ; AEXX=0.2 ; HFSCREEN = 0.2):

k-point   1: 0.0000    0.0000    0.0000
DOPED NKRED = 2           DOPED NKRED = 1              UNDOPED CASE
band No.  band energies occupation   band energies occupation   band energies occupation
valence bands
262       2.4107      2.00000        2.4339      2.00000        2.4082      2.00000
263       2.4107      2.00000        2.4339      2.00000        2.4082      2.00000
264       2.8522      2.00000        2.8597      2.00000        2.8566      2.00000
conduction bands
265       5.4046      2.00000        5.8240      1.87262        5.8126      0.00000
266       5.4908      2.00000        5.8695      1.62151        5.8424      0.00000
267       5.4894      2.00000        5.8695      1.62192        5.8424      0.00000

k-point   2: 0.5000    0.0000    0.0000
DOPED NKRED = 2           DOPED NKRED = 1              UNDOPED CASE
band No.  band energies occupation   band energies occupation  band energies occupation
valence bands
262       2.0015      2.00000        2.0144      2.00000       2.0160      2.00000
263       2.5961      2.00000        2.6072      2.00000       2.6046      2.00000
264       2.5961      2.00000        2.6072      2.00000       2.6045      2.00000
conduction bands
265       6.1904      0.00000        6.1335      0.00435       6.0300      0.00000
266       6.1904      0.00000        6.1335      0.00435       6.0300      0.00000
267       6.1907      0.00000        6.1340      0.00426       6.0305      0.00000

k-point   3 :  0.5000    0.5000    0.0000
DOPED NKRED = 2           DOPED NKRED = 1              UNDOPED CASE
band No.  band energies occupation   band energies occupation  band energies occupation
valence bands
262       2.4237      2.00000        2.4433      2.00000       2.4287      2.00000
263       2.4238      2.00000        2.4432      2.00000       2.4287      2.00000
264       2.4239      2.00000        2.4433      2.00000       2.4287      2.00000
conduction bands
265       5.8966      0.42674        5.9100      1.24121       5.8817      0.00000
266       5.8780      0.54128        5.9100      1.24143       5.8817      0.00000
267       5.8826      0.50661        5.9100      1.24261       5.8817      0.00000


Without NKRED, the one electron energies are pretty similar to the one electron energies in the undoped system (last two columns), whereas using NKRED a strong reduction of the "gap" between the valence and conduction band is observed, in particular, close to the conduction band minimum (in this case the point). This result is an artefact of the approximation used for NKRED=2. The non-local exchange operator cancels the self-interaction present in the Hartree-potential. For NKRED=2 and 2×2×2 k-points, the non-local exchange operator at each k-point is evaluated using the one-electron orbitals at this k-point only, e.g.:

${\displaystyle V_{\mathbf {k} }\left({\mathbf {G} },{\mathbf {G} }'\right)=\langle {\mathbf {k} }+{\mathbf {G} }|V_{x}|{\mathbf {k} }+{\mathbf {G} }'\rangle =-{\frac {4\pi e^{2}}{\Omega }}f_{m{\mathbf {k} }}\sum _{{\mathbf {G} }''}{\frac {C_{m{\mathbf {k} }}^{*}({\mathbf {G} }'-{\mathbf {G} }'')C_{m{\mathbf {k} }}({\mathbf {G} }-{\mathbf {G} }'')}{|{\mathbf {G} }''|^{2}}}}$

The sum over q in the Fock exchange potential reduces to a single k-point. This reduces the self-interaction for states that have originally an occupancy larger one, concomitantly pulling those states to lower energies. Initially empty states (occupancy smaller one) are pushed up slightly. Since this is clearly an artefact, NKRED must be used with utmost care for large supercells with coarse k-point sampling. Please always check whether occupancies are similar at all k-points, if this is not the case, the calculations should be double checked without downsampling.

Since HF type calculations using 2×2×2 k-points without NKRED, are roughly 64 times more expensive than those using the Γ-point only, it might seem impossible to do anything but Γ-point only calculations. However, VASP allows to generate special k-points using generating lattices.

Particularly useful for HF type calculations, are the following k-point sets

k-point set generating a bcc like lattice in the BZ ->  2 k-points in BZ
0
direct
0.5 0.5 0.5
-.5 -.5 0.5
0.5 -.5 -.5
0 0 0


This KPOINTS file generates two 2 k-points, one at the Γ-point and one along the space diagonal at the BZ boundary (R-point).

The following KPOINTS file generates 4 k-points, one at the Γ-point and three at the S-points (the latter ones might be symmetry equivalent for cubic cells).

k-point set generating an fcc lattice ->  4 k-points in BZ
0
direct
0.5 0.5 0.0
0.0 0.5 0.5
0.5 0.0 0.5
0 0 0


Using such grids, sensible and fairly rapidly converging results are obtained, e.g. for electron and hole doped materials, even if the conduction or valence band is partially occupied or depleted. For instance for TiO2 the following energies are obained:

Gamma only     TOTEN  =      -837.759900 eV
2 k-points     TOTEN  =      -838.039157 eV
4 k-points     TOTEN  =      -838.129712 eV
2x2x2          TOTEN  =      -838.104787 eV
2x2x2 NKRED=2  TOTEN  =      -838.418681 eV


Note that results using NKRED not improved compared to Γ-only calculations.