Calculating the band gap using tb-mbj and the full hybrid method YS-PBE0

We will take the example of MgO

Reference http://susi.theochem.tuwien.ac.at/events/ws2017/notes/tutorial-mBJ-hybrid-DFT.pdf

   ----------------------------PBE method -----------------------

MgO

(NaCl, a=7.96 bohr; default initialization; scf-cycle)

PBE: check the gap (:GAP from “anaylysis”),

plot a band structure in PBE

(E-range from -19 to 18 eV)

    
 We will use thr GUI interface w2web

Analysis of parameter:
:GAP
in MgO.scf (showing last 10 / 1 lines)

--- GAP -----------
in  1 files:
MgO.scf::GAP  :    0.3492 Ry =     4.749 eV   (provided you have a proper k-mesh)

   ---------------------------------------------------------------------
------------ tb-mbj method  -----------------------------------

B-mBJ:

save the PBE calculation, execute:

init_mbj_lapw (in utils) „phase 1“ of

the initialization (see also in the UG 4.5.9)

run_lapw -NI -i 1

rm *.bro*

init_mbj_lapw    „phase 2“,use original mBJ parameters

run scf cycle (note, it may not converge in 40 cycles, submit another runwith -NI option)

monitor the change of the :GAP

plot a band structure (fcc) and compare with PBE

NB. You need to use the .machines file to do parallel calculation

wadha@linux-57qi:~/WIEN2k/MgO-PBE> save_lapw -d pbe
wadha@linux-57qi:~/WIEN2k/MgO-PBE> init_mbj_lapw
prepared MgO-PBE.inm_vresp and changed to R2V in MgO-PBE.in0
Now do:  
run_lapw -i 1 -NI    #   to prepare the r2v and vresp files
save_lapw -d pbe      #   save the pbe run
init_mbj_lapw         #   rerun the init script to finish mbj-initialization
wadha@linux-57qi:~/WIEN2k/MgO-PBE> run_lapw -i 1 -NI
LAPW0 END
 LAPW1 END
 LAPW2 END
 CORE  END
 MIXER END
 MIXER END
 MIXER END
ec cc and fc_conv 0 1 1

>   energy in SCF NOT CONVERGED

wadha@linux-57qi:~/WIEN2k/MgO-PBE> rm *.bro*
wadha@linux-57qi:~/WIEN2k/MgO-PBE> init_mbj_lapw
prepared MgO-PBE.in0_grr and changed to XC_MBJ in MgO-PBE.in0
You can use different parametrizations of mBJ:
0:    Original mBJ values (Tran,Blaha PRL102,226401) (default)
1:    New parameterization (Koller etal, PRB85, 155109)
2:    New parameterization for semiconductors (gaps up to 7 eV)
3:    Unmodified BJ potential (Becke,Johnson J.Chem.Phys 124,221101

 enter 0/1/2/3 or press enter:0
Now do the mBJ calculation:   run_lapw -i 80 ...

wadha@linux-57qi:~/WIEN2k/MgO-PBE> run_lapw -p -i 80


wadha@linux-57qi:~/WIEN2k/MgO-PBE> grep :GAP MgO.scf

:GAP  :    0.5288 Ry =     7.192 eV   (provided you have a proper k-mesh)

  ------------------------- Full hybrid   YS-PBE0 method --------------------

Perform a hybrid-DFT calculation using YS-PBE0

create a new case, perform a PBE calculation and save the results.

the setup for hybrid-calculations can be made in w2web

(Utils/init_hf_lapw), or in a terminal-window using „init_hf_lapw“.

(More details are given in the UG 4.5.8)

Select NBAND=12 (case.inhf)

and a 4x4x4 / 4x4x4 k-point mesh (no reduction)

scf cycle with -hf -p -scratch ./

(insert 4 lines with 1:localhost into .machines)

we do this in k-parallel since it will take more time, alternatively we could alsouse a „reduced“ hf-k-mesh, see UG

monitor the change of the :GAP and compare it with mBJ and exp.gaps (only every 2nd value is from HF !)

plot a band structure:

only the k-mesh selection can be done in w2web, then

open a terminal and change into the proper directory

run_bandplothf_lapw -p

cp $WIENROOT/SRC_templates/case.insp case.insp

(insert E

F

and increase the plotting energy range).

x spaghetti -hf -p

After creating a new case and doing a PBE calculation and verifying the gap value as follows

wadha@linux-57qi:~/WIEN2k/MgO-YS-PBE0> grep :GAP MgO-YS-PBE0.scf

:GAP  :    0.3492 Ry =     4.749 eV   (provided you have a proper k-mesh)

and we will save it as follows

wadha@linux-57qi:~/WIEN2k/MgO-YS-PBE0> save_lapw -d pbe
wadha@linux-57qi:~/WIEN2k/MgO-YS-PBE0> init_hf_lapw

Insulator, EF-inconsistency corrected
:GAP  :    0.3492 Ry =     4.749 eV   (provided you have a proper k-mesh)
         Bandranges (emin - emax) and occupancy:
:BAN00001:   1   -5.045696   -5.044741  2.00000000
:BAN00002:   2   -2.606216   -2.599983  2.00000000
:BAN00003:   3   -2.603726   -2.599930  2.00000000
:BAN00004:   4   -2.602120   -2.599930  2.00000000
:BAN00005:   5   -0.998725   -0.881694  2.00000000
:BAN00006:   6   -0.069527    0.269038  2.00000000
:BAN00007:   7    0.047767    0.269038  2.00000000
:BAN00008:   8    0.142198    0.269038  2.00000000
:BAN00009:   9    0.618208    1.097278  0.00000000
:BAN00010:  10    0.994285    1.432520  0.00000000
:BAN00011:  11    1.173394    1.503849  0.00000000
:BAN00012:  12    1.432520    1.686549  0.00000000
:BAN00013:  13    1.524931    1.764080  0.00000000
        Energy to separate low and high energystates:   -0.11953


:NOE  : NUMBER OF ELECTRONS          =  16.000

:FER  : F E R M I - ENERGY(TETRAH.M.)=   0.2690382780
:GMA  : POTENTIAL AND CHARGE CUT-OFF  12.00 Ry**.5

You must set NBAND to at least  NB_occ + 1  and you have 16.00 electrons
edit MgO-YS-PBE0.inhf ...

when the file opens modify the nband to 12

PS: For very accurate calc. and large NBAND you may have to increase EMAX in MgO-YS-PBE0.in1 by hand

Prepare k-mesh for HF. Use identical mesh and shift for IBZ and FBZ
Do you want to use a REDUCED k-mesh for HF (saving cpu-time) (Y/n) ?
n
 This script runs   x kgen   twice and generates equivalent meshes in the
 IBZ and FBZ.



KGEN ENDS
KGEN ENDS
How many k-points in full BZ?
If you type 0 you can give 3 numbers for nx,ny,nz
0
How many in x?
4
How many in y?
4
How many in z?
4
  NUMBER OF K-POINTS IN WHOLE CELL: (0 allows to specify 3 divisions of G)
 length of reciprocal lattice vectors:   1.367   1.367   1.367   0.000   0.000   0.000
  Specify 3 mesh-divisions (n1,n2,n3):
           8  k-points generated, ndiv=           4           4           4
KGEN ENDS
0.013u 0.001s 0:00.01 100.0%    0+0k 0+696io 0pf+0w
           1  symmetry operations without inversion
  NUMBER OF K-POINTS IN WHOLE CELL: (0 allows to specify 3 divisions of G)
 length of reciprocal lattice vectors:   1.367   1.367   1.367   0.000   0.000   0.000
  Specify 3 mesh-divisions (n1,n2,n3):
  Shift of k-mesh allowed. Do you want to shift: (0=no, 1=shift)
          64  k-points generated, ndiv=           4           4           4
KGEN ENDS
0.001u 0.001s 0:00.00 0.0%      0+0k 0+104io 0pf+0w
MgO-YS-PBE0.in0_grr and MgO-YS-PBE0.inhf and hf-kmesh prepared
Now do the hybrid calculation:   run_lapw -hf  ... 

You do the initialization in the GUI interface as follows

You must set NBAND to at least NB_occ + 1     ----> 12 nband

set larger EMAX in MgO-YS-PBE0.in1 (only in case you want NBAND > available bands)

Prepare k-mesh for HF. You may want to use a (commensurate) REDUCED k-mesh for the HF-potential to save cpu-time.                   with k-mesh (eg. 4x4x4):  nx= 4 ny= 4 nz= 4       with commensurate reduced k-mesh (eg. 2x2x2):  nx= ny= nz=

 and click on run_kgenrf_lapw

wadha@linux-57qi:~/WIEN2k/MgO-YS-PBE0> run_lapw -hf -p 

wadha@linux-57qi:~/WIEN2k/MgO-YS-PBE0> grep :GAP MgO-YS-PBE0.scf

:GAP  :    0.4749 Ry =     6.459 eV   (provided you have a proper k-mesh)

GOOD LUCK