11. Hybrid functionals

In this section, we introduce hybrid functional calculations in RESCU. Hybrid functional calculations are so named because the exchange functional is a weighted sum of a semi-local exchange functional (e.g. PBE) and the exact exchange functional (employed in the Hartree-Fock method). Such calculations are more computationally demanding because the exact exchange functional gives rise to a non-local potential in the Kohn-Sham equation, but they are also more accurate than semilocal DFT calculations for a number of applications. Moreover, the non-locality of the functional brings about some infinity-related issues in periodic systems. RESCU implements the NOPAIN method [C18, CCMG18] to deal with both problems. Simply put, the NOPAIN method performs an analytical low-rank decomposition of the Coulomb repulsion integral, which is the atomic orbital representation of the Coulomb interaction kernel. Most importantly, basis functions retain their locality during the decomposition and the coefficients are obtained from the planewave representation of the Coulomb operator.

11.1. Example: diamond band gap

Certain hybrid functionals famously correct band gap in the Kohn-Sham spectrum. In this example, we show how to calculate the DOS, and hence band gap, of the diamond crystal. First, we perform a self-consistent calculation to obtain the ground state density and properties. Please create the input file c_lcao_scf_hse as follows

LCAO.status = true % required for hybrid
info.calculationType = 'self-consistent'
info.savepath = 'results/c_lcao_scf_hse'
atom.element = [1 1]
atom.fracxyz = [0 0 0;0.25 0.25 0.25]
domain.latvec = 6.74*[.5 .5 0;.5 0 .5;0 .5 .5]
domain.lowres = 0.3
element(1).species = 'C'
element(1).path = './C_OV_LDA.mat'
functional.libxc = true % required for hybrid
functional.list = {'XC_HYB_GGA_XC_HSE06'}
kpoint.gridn = [3 3 3] % must be gamma centered
Exx.Gcutoff = 4.0

We comment on the purpose of each keyword in a hybrid calculation context.

  • LCAO.status must be set to true. As mentioned in the section introduction, the NOPAIN method represents the Coulomb kernel using atomic orbitals, and hence the calculation ought to be done in an NAO basis;

  • element.path uses a pseudopotential drawn from the LDA-Troullier-Martins database. For certain atoms (e.g. transition metals, heavy atoms) hybrid calculations may require a pseudopotential including more electronic shells in the reference valence configuration and a more extensive atomic orbital basis. The standard pseudopotential file is sufficient for our simple example;

  • functional.libxc must be set to true. RESCU draws the semi-local part of hybrid functionals from the LibXC library;

  • functional.list includes a single functional (2006-version of the HSE functional) combining the semi-local exchange and correlation functionals in this example. At the moment, RESCU implements only the PBE0 and HSE functionals;

  • kpoint.gridn defines the grid sampling reciprocal space. This grid must be \(\Gamma\)-centered, so set kpoint.shift = 0.5*[1 1 1] for even grids. In HSE calculations, the non-local exchange potential couples the various k-points such that the workload and memory requirements quickly grow with increasing k-sampling. Moreover, RESCU will ignore any symmetry when performing hybrid calculations, which means that the full Brillouin zone is resolved;

  • Exx.Gcutoff is the planewave cutoff, in Hartree, used in the Coulomb kernel decomposition. The higher the more accurate, but also the more computationally expensive. In a research project, convergence with respect to this parameter is crucial. The default value is 4, which is sufficient for several simple materials.

You may execute RESCU as follows

rescu -i c_lcao_scf_hse

The standard output will ensue with the addition of something like

Initializing YqG-matrix                               68.365
Initializing Hartree-Fock potential                    8.282

The first operation corresponds to the NOPAIN decomposition mention previously. In small systems, its cost may dominate the entire calculation, but this is less and less the case as the system size increases. The second operation is the calculation of the exchange exchange potential matrix, which is required to build the Hamiltonian in the first iteration.

Next, we write the DOS input file c_lcao_dos_hse as follows

info.calculationType = 'dos'
info.savepath = 'results/c_lcao_dos_hse'
info.outfile = 'fc_lcao_dos_hse.out'
rho.in = './results/c_lcao_scf_hse'
kpoint.gridn = [10 10 10]
kpoint.sampling = 'tetrahedron'
dos.range = [-20,10]
units.energy = 'ev'

This is a standard DOS input file. We employ the tetrahedron method and set the range from -20 eV to +10 eV around the Fermi level. You can calculate the DOS and plot it typing

rescu -i c_lcao_dos_hse
rescu -p results/c_lcao_dos_hse

The band gap can be extracted as follows

load('results/c_lcao_dos_hse','energy')
nrg = energy.ksnrg;
bg = (min(nrg(5,:))-max(nrg(4,:)))*27.211

In the last line, we subtract the minimum of the fifth band (conduction band minimum) from the maximum of the fourth band (valence band maximum) and multiply by 27.211 to convert the units from Hartree to eV. You should get a value close to 6.2 eV, which is 0.7 eV above the experimental value of 5.5 eV. Following Ref. [C18], increasing the k-sampling to [7,7,7], the planewave cutoff to 6 Ha and using an optimized DZP basis brings the band gap in line with the experimental value. If you repeat the calculation using PBE, keeping all other parameters the save, you should get a value close to 4.7 eV.

[CCMG18]

Chen Ying-Chih, Chen Jing-Zhe, Vincent Michaud-Rioux, and Hong Guo. Efficient evaluation of nonlocal operators in density functional theory. Physical Review B, 97(7), 075139, 2018.