In recent years the interest in hybrid functionals (that is the incorporation of parts of Hartree-Fock exchange in calculations based on common approximations of exchange-correlations such as LDA, GGA and so on) has steadily increased owing to its improvement over most common functionals, especially when it comes to band-gap calculation of extended solids – see for example this work for a comparison.
Quantum Espresso now offers a variety of hybrid functionals (for a complete list see the header of funct.f90) but currently the code can only use hybrid functionals for self-consistent calculations. The reason being that:

The problem is quite fundamental, because in order to get the Fock operator at a certain k-point you need the wavefunctions on a grid that is commensurate with it, this can only be done self-consistently. (L. Paulatto)

But there is a quite elegant solution via maximally localized Wannier functions. If you worked with wannier functions before you know they required the “full grid” (conveniently generated using the tool kmesh.pl) but a recently introduced tool allows to “unfold” the SCF calculation of the reduced grid to the full grid using the executable open_grid.x. An example is provided with recent versions of QE as well.
In the following I will use this approach to improve the band-gap of MgO. Using PBE pseudopotentials the band-gap is underestimated at about 4.3 eV but using the HSE hybrid the band-gap can be improved to about 7 eV, which is within acceptable range of the experimental value of ~7.7 eV.The workflow is fairly straightforward:

run a (converged) SCF calculation with input_dft=’HSE’ and a number of empty bands. You have to ensure convergence with respect to the usual parameters (k-points, cutoff, …) AND the mesh for the Fock operator (nqx)

unfold the reduced grid onto the full grid using open_grid.x

wannierize the obtained wave-functions using wannier90.x and plot the band-structure along a desired high symmetry path

The wannierization is by far the most tricky part in this particular example but by projecting on O:p and Mg:s one can accurately describe the valence band and an additional single conduction band. For the resulting band-structure see the figure below.
So this is actually quite straight-forward. please note that in above calculation the Fock operator was calculated on a very coarse 1x1x1 grid. The calculation on such a coarse grid actually seems to over-estimate the band-gap and a converged energy can be obtained on a 6x6x6 nqx grid with ~6.7 eV.
The entire calculation can be run using the attached script – I hope you find it helpful!

thanks for the tutorial and references; it was very insightful. I was wondering how to plot the results. We essentially have to harvest the orbital eigenvalues and match with the kpts file from open_grid.x right?

Dear zww4855, could you let me know what exactly you want to plot? The bands? The easiest way would be to let wannier90 do the work by using bands_plot = .true. in the win file and then plot the prefix_band.dat

In my example I get quite nice bands, at least close to the Gamma point; if you want to improve the bands you probably have to get better MLWFs!

Thanks for sharing your experience of hybrid functional band structure calculation!

I have a problem with the script of open_grid.x , everytime the script will show the segmentation fault when I tried to unfold the reduced grids obtained by paw pseudopotentials. However, it works fine with the results obtained by the NC pseudo, do you know why? Or can the open_grid.x grid only deal with the results from NC pseudo?

I also have to admit that I rarely use PAW as I am quite satisfied with the PPs from the ONCV library which are “almost as soft” as PAW/US (well,.. “almost”)! Until this is resolved I think PAW cannot be used with open_grid.x

Honestly, I don’t think it’s only an issue with PAW PP. I’m using norm-conserving ones (pbesol-nc), and I’m still getting a segmentation fault when running open_grid.x:

————————————————————————–
mpirun noticed that process rank 41 with PID 0 on node n006.cluster.com exited on signal 11 (Segmentation fault).
————————————————————————–

ATOMIC_POSITIONS (angstrom)
C -0.014436109 0.000000000 -0.008945780
C 0.633731142 0.000000000 1.185927865
H -1.102040351 0.000000000 -0.004392100
H 1.721345317 0.000000000 1.181410015

Thanks for sharing!! Very nicely described

I would like to try if phonon calculation can be performed for more accurate epsilon results.

LikeLike

Hi Christoph,

thanks for the tutorial and references; it was very insightful. I was wondering how to plot the results. We essentially have to harvest the orbital eigenvalues and match with the kpts file from open_grid.x right?

LikeLike

Dear zww4855, could you let me know what exactly you want to plot? The bands? The easiest way would be to let wannier90 do the work by using bands_plot = .true. in the win file and then plot the prefix_band.dat

In my example I get quite nice bands, at least close to the Gamma point; if you want to improve the bands you probably have to get better MLWFs!

HTH!

Chris

LikeLike

Hi Christoph,

Thanks for sharing your experience of hybrid functional band structure calculation!

I have a problem with the script of open_grid.x , everytime the script will show the segmentation fault when I tried to unfold the reduced grids obtained by paw pseudopotentials. However, it works fine with the results obtained by the NC pseudo, do you know why? Or can the open_grid.x grid only deal with the results from NC pseudo?

Thanks!

Ruishen

LikeLike

Dear Ruishen,

This seems to be an open (?) problem: https://gitlab.com/QEF/q-e/issues/109

I also have to admit that I rarely use PAW as I am quite satisfied with the PPs from the ONCV library which are “almost as soft” as PAW/US (well,.. “almost”)! Until this is resolved I think PAW cannot be used with open_grid.x

HTH!

Chris

LikeLike

Honestly, I don’t think it’s only an issue with PAW PP. I’m using norm-conserving ones (pbesol-nc), and I’m still getting a segmentation fault when running open_grid.x:

————————————————————————–

mpirun noticed that process rank 41 with PID 0 on node n006.cluster.com exited on signal 11 (Segmentation fault).

————————————————————————–

Do you know what may be the issue?

LikeLike

My input scf file looks like this:

&control

calculation = ‘scf’,

prefix = ‘pa’,

verbosity = ‘high’,

pseudo_dir = ‘/home/ccevallos/pslibrary.1.0.0/pbesol/PSEUDOPOTENTIALS’,

outdir = ‘./tmp’,

wf_collect=.true.

/

&SYSTEM

ibrav = 0,

nat=4,

ntyp=2,

ecutwfc=100.00,

ecutrho=1000.00,

occupations = ‘fixed’,

input_DFT=’PBE0′,

nqx1=1,

nqx2=1,

nqx3=25,

/

ATOMIC_SPECIES

C 12.011 C.pbesol-nc.UPF

H 1.00784 H.pbesol-n-nc.UPF

CELL_PARAMETERS {angstrom}

10.000000000000000 0.000000000000000 0.000000000000000

0.000000000000000 10.000000000000000 0.000000000000000

0.000000000000000 0.000000000000000 2.46

ATOMIC_POSITIONS (angstrom)

C -0.014436109 0.000000000 -0.008945780

C 0.633731142 0.000000000 1.185927865

H -1.102040351 0.000000000 -0.004392100

H 1.721345317 0.000000000 1.181410015

K_POINTS automatic

1 1 25 1 1 1

LikeLike