Calculating HSE06 band structures on Quantum ESPRESSO

Question

I am currently trying to figure out how to compute band structures for my system, using the hybrid functional HSE06. I'm doing this on Quantum Espresso. As I understand, there are a handful of ways to do this:

1) Generate maximally localized wannier functions (MLWFs) with the HSE functional turned on

2) The 'fake scf' procedure: This is similar to what is listed in Vaspwiki(here).

3) Generate scf data for a coarse 'q' mesh (of the Fock operator) and then interpolate these bands using softwares like BoltztraP2.

I am not familiar with the Wannier module on Quantum ESPRESSO, hence option (1) is sort of last resort for me, since Wannier takes some time to learn.

The 'fake scf' procedure doesn't seem to work. The k-point mesh for my original calculation (my system is a bilayer TMD - Transition metal dichalcogenide) is 12x9x1. I tried using a commensurate 'q' grid and tried the 'fake scf' procedure, but I couldn't get it to work.

I would be grateful if anyone could help me figure out how to compute these bandstructures - either through one of the methods I mentioned, or any easy method as you see fit. I would like to add that I have access to considerable computational resources, so doing expensive calculations should not be much of an issue.

Also, I would like to say that I am open to methods that use other softwares - including VASP or CASTEP. I have both licenses (CASTEP through Materials Studio). But things can get a little tricky here because VASP for example, uses a different implementation of the Hubbard 'U' and it would corrupt my results if I port my QE input to VASP to calculate the band-structures.

Answer 1

I have seen all the methods you mentioned but have only done one myself; I'll explain here how to use Wannier90 in conjunction with Quantum Espresso to get band structures for hybrid functional calculations. It does take a bit of time to learn, but not very long! You can learn the basics and do some first calculations in an afternoon. There are subtleties of course that I won't get into, and if you have specific issues you should probably ask another question or send a message to the Wannier90 mailing list.

Initial Checklist

1. Compile and install Wannier90

2. If you are using QE <= 6.0, you will need to replace the qe/PP/src/pw2wannier90.f90 file with one provided in the Wannier90 tarball, and recompile QE. Otherwise, make sure you have a working compiled version of QE. Newer versions will typically be better for hybrid functional calculations (pairwise band parallelization, ACE algorithm, support for ultrasoft & PAW...)

3. This method requires commensurate k-point grids between QE and Wannier90. Certain k-point meshes may cause issues in Wannier90 when identifiying nearest neighbors, etc. It's recommended that you prepare a manual k-point grid with the utility located in wannier90/utility/kmesh.pl and use it in your QE calculation. Wannier90 does not recognize k-point symmetry, so it is a good idea to disable symmetry in your hybrid calculation with nosym = .true and noinv = .true. in the &SYSTEM block, along with the aforementioned manual k-points at the highest density you can afford, with corresponding q-point grid (as an aside, this type of grid is always needed for Wannier90, but normally you use this full grid with an nscf calculation after your automatic grid scf calculation when not using hybrid functionals, to reduce computational cost). I will leave it to you to determine what you need in terms of k- and q-point densities for your specific calculation.

4. The input files you need at the beginning are: DFT input for QE, Wannier90 input file (.win), and a pw2wan input file. Here are some I used for a very quick (probably unconverged) calculation I did on my desktop, just using bulk silicon. In the interest of time I used automatic projections with the SCDM method and default values, you should read up on this and decide how you should proceed in your own calculation (manual projections require the begin projections block in Wannier90 input--read the user guide for more info). It can affect the band structure significantly.

QE input file: silicon.in

 &control
    calculation = 'scf'
    restart_mode='from_scratch',
    prefix='Si-HSE',
    pseudo_dir = './',
    outdir='./TMP_DIR/'
 /
 &system    
    ibrav=  2, celldm(1) =10.20, nat=  2, ntyp= 1,
    ecutwfc =30.0,  nbnd = 8,
    input_dft='hse', nqx1 = 1, nqx2 = 1, nqx3 = 1, 
    x_gamma_extrapolation = .true.,
    exxdiv_treatment = 'gygi-baldereschi',
    nosym = .true., noinv = .true
 /
 &electrons
    mixing_beta = 0.7
 /
ATOMIC_SPECIES
 Si  28.086  Si_ONCV_PBE-1.1.upf
ATOMIC_POSITIONS alat
 Si 0.00 0.00 0.00 
 Si 0.25 0.25 0.25 
K_POINTS crystal
8
  0.00000000  0.00000000  0.00000000  1.250000e-01
  0.00000000  0.00000000  0.50000000  1.250000e-01
  0.00000000  0.50000000  0.00000000  1.250000e-01
  0.00000000  0.50000000  0.50000000  1.250000e-01
  0.50000000  0.00000000  0.00000000  1.250000e-01
  0.50000000  0.00000000  0.50000000  1.250000e-01
  0.50000000  0.50000000  0.00000000  1.250000e-01
  0.50000000  0.50000000  0.50000000  1.250000e-01

Wannier90 input file: silicon.win

! Silicon HSE 1
 num_wann    = 8
 num_bands   = 8
 num_iter    = 20
 kmesh_tol   = 0.0000001
 auto_projections = .true.

! Use as much precision as you can (at least 6 decimals) to prevent issues with matching to QE output
Begin Unit_Cell_Cart
-2.698804 0.0000 2.698804
 0.0000 2.698804 2.698804
-2.698804 2.698804 0.0000
End Unit_Cell_Cart

begin atoms_frac
Si 0.00  0.00  0.00
Si 0.25  0.25  0.25
end atoms_frac

!begin projections
!Si:sp3
!end projections

! To plot the WF interpolated bandstructure
bands_plot       = .true.
bands_num_points = 200

begin kpoint_path
L 0.50000  0.50000 0.5000 G 0.00000  0.00000 0.0000
G 0.00000  0.00000 0.0000 X 0.50000  0.00000 0.5000
X 0.50000 -0.50000 0.0000 K 0.37500 -0.37500 0.0000
K 0.37500 -0.37500 0.0000 G 0.00000  0.00000 0.0000
end kpoint_path

! KPOINTS

mp_grid : 2 2 2

begin kpoints
0.0 0.0 0.0
0.0 0.0 0.5
0.0 0.5 0.0
0.0 0.5 0.5
0.5 0.0 0.0
0.5 0.0 0.5
0.5 0.5 0.0
0.5 0.5 0.5
end kpoints

pw2wannier90 input file: silicon.pw2wan

&INPUTPP
 prefix = 'Si-HSE'
 outdir = './TMP_DIR'
 seedname = 'silicon'
 scdm_proj = .true.
 scdm_entanglement = 'isolated'
 scdm_mu = 0.0
 scdm_sigma = 1.0
/

Calculation Procedure

1. Run Wannier90 to generate some needed files, with the postprocessing option: $ wannier90 -pp silicon
2. Run your QE hybrid functional calculation (obviously use the correct mpi commands for your system when running in parallel): $ pw.x -inp silicon.in > silicon.out
3. Run pw2wannier: $ pathtoqe/PP/src/pw2wannier.x -inp silicon.pw2wan
4. Run Wannier90 again to do the calculation and get band structure: $ wannier90 silicon

Doing the above I got this band structure. It obviously has some problems (some bands jumping up above the VBM) but the gap is larger than PBE and I did this with a very sparse grid, and no optimization of any options.

Answer 2

Apart from methods (1) and (2) mentioned in the question, the third method I found to be pretty viable under certain conditions. The procedure involves mathematical interpolation of the bands, so I have to caution the reader that this is not completely rigorous like Wannier interpolation, which retains the nature of the wave-functions generated from the SCF calculation. That being said, I'm pleasantly surprised by how well the BoltzTrap2 code handles interpolation, along with degeneracies.

The procedure is as follows:

i) Perform a self-consistent calculation with HSE turned on. The 'q' mesh (or Fock grid) that you choose for the calculation will affect the final band-structure. Typically, for big systems, one would not be able to choose a Fock grid that is commensurate with your K-grid. For instance, I worked on a 2D semiconductor with around 100 electrons in the system. My K-mesh was 12x9x1. In this case, I experimented with 'q' meshes of 1x1x1, 4x3x1 and 12x9x1. As one might imagine, the ideal choice would be 12x9x1 for the 'q' mesh since its commensurate with the K-grid. But this is very computationally heavy. My band-structure calculations, surprisingly, were somewhat converged at the 4x3x1 'q' mesh itself the band-gap difference with the commensurate mesh was around 0.05 eV. The bottom-line here is, the density of the 'q' mesh will affect the precision of the band-structure calculation since the basis is mathematical (for the Wannier method, this is not the case).

ii) Install BoltzTrap2. Rename the SCF input and output files as prefix.nscf.in and prefix.nscf.out . This needs to be done because the code takes in file names which are specified in the above format.

iii) Run the 'qe2boltz.py' file to convert the QE input to a format that's readable by BoltzTrap2:

python qe2boltz.py prefix pw fermi f

Replace 'f' here with the value of the Fermi level printed in your SCF output.

iv) Run the bands interpolation using

 btp2 -vv interpolate -m 'n'

where 'n' decides how fine your interpolation needs to be for the band-structure. Typically, values of 1000-2000 should be good enough.

v) Plot the bands using the plotbands command.

For any logistical issues, the BoltzTrap code documentation and also this tutorial might be helpful.