4
$\begingroup$

I used QE to calculate DOS for hBN. I expect to get results as shown here. The behavior of the computed DOS is as expected but QE gives Fermi energy = 1.974 eV. So, when I plot DOS concerning E-Ef, the Fermi level isn't located at zero:

enter image description here

If we look at the figure, we can observe that the highest occupied level is about 5.126 eV (as in the DOS data file).

Why QE gives Fermi energy = 1.974 eV, not 5.126 eV?

Edit: Here is the input file I used:

&CONTROL
  calculation = 'scf'
  outdir = './out/'
  prefix = 'hBN'
  pseudo_dir = '.'
/
&SYSTEM
  ecutwfc = 40
  ibrav = 0
  nat = 4
  ntyp = 2
/
&ELECTRONS
  conv_thr =   1.0d-10
/
ATOMIC_SPECIES
B     10.0 B_ONCV_PBE_FR-1.0.upf
N     14.0 N_ONCV_PBE_FR-1.0.upf

ATOMIC_POSITIONS crystal
B            0.6666666667       0.3333333333       0.7500000000 
B            0.3333333333       0.6666666667       0.2500000000 
N            0.6666666667       0.3333333333       0.2500000000 
N            0.3333333333       0.6666666667       0.7500000000 

CELL_PARAMETERS angstrom
     -1.2562141150      -2.1758266724      -0.0000000000
     -1.2562141150       2.1758266724       0.0000000000
      0.0000000000       0.0000000000      -7.7072650000

K_POINTS automatic
9 9 9 0 0 0

For the nscf I have added

nbnd = 16
occupations = 'tetrahedra'

to the &SYSTEM section.

Here is a comparison with the results in the materialsproject

enter image description here

$\endgroup$
5
  • $\begingroup$ Could you share your input and output files please? As a side-note: never trust any DFT calculation using the method of D. Bagayoko, it is a completely unphysical approach with no valid justification. $\endgroup$ Dec 2, 2023 at 3:02
  • $\begingroup$ @PhilHasnip I Edited the question, please take a look $\endgroup$
    – Bekaso
    Dec 4, 2023 at 11:40
  • $\begingroup$ How do you know that the states upto 5.126 eV are filled? $\endgroup$ Dec 4, 2023 at 12:15
  • $\begingroup$ The results from the Materials Project look identical to yours, within the usual uncertainty in Ef. Why do you think yours are incorrect? $\endgroup$ Dec 4, 2023 at 14:33
  • $\begingroup$ I think you need to do (E-Ef). Quantum espresso gives Energy including fermi energy. $\endgroup$ Dec 12, 2023 at 12:15

1 Answer 1

3
$\begingroup$

I have repeated your SCF calculation with SSSP precision pseudopotentials without doing any convergence testing and using the suggested minimum values. The input and the output files are attached at the end. But here are the results: the band gap is found to be 4.5354 eV as expected. The highest occupied and lowest unoccupied level are 3.1546 eV and 7.7535 eV. The DOS also seems to be correct. I arbitrarily set the Fermi energy in the highest occupied level since Fermi energy can be anywhere in the band gap in an insulator. So, your result seems correct to me but I just used a more strict set of parameters. I also use the nbnd parameter manually and occupations='fixed' in order to show the band gap in the output file of the SCF run. enter image description here

Input:

&CONTROL
  calculation = 'scf'
  outdir = './out/'
  prefix = 'hBN'
  pseudo_dir = '.'
/
&SYSTEM
  ecutwfc = 80
  ecutrho = 440
  ibrav = 0
  nat = 4
  ntyp = 2
  smearing='gauss'
  occupations='fixed'
  degauss=0.01
  nbnd=16
/
&ELECTRONS
  conv_thr =   1.0d-10
  !startingpot='file'
  !startingwfc='file'
/
ATOMIC_SPECIES
B     10.0 B_pbe_v1.01.uspp.F.UPF
N     14.0 N.oncvpsp.upf
ATOMIC_POSITIONS crystal
B            0.6666666667       0.3333333333       0.7500000000 
B            0.3333333333       0.6666666667       0.2500000000 
N            0.6666666667       0.3333333333       0.2500000000 
N            0.3333333333       0.6666666667       0.7500000000 
CELL_PARAMETERS angstrom
     -1.2562141150      -2.1758266724      -0.0000000000
     -1.2562141150       2.1758266724       0.0000000000
      0.0000000000       0.0000000000      -7.7072650000
K_POINTS automatic
10 10 10 0 0 0

NSCF Input:

&CONTROL
  calculation = 'nscf'
  outdir = './out/'
  prefix = 'hBN'
  pseudo_dir = '.'
/
&SYSTEM
  ecutwfc = 80
  ecutrho = 440
  ibrav = 0
  nat = 4
  ntyp = 2
  smearing='gauss'
  occupations='tetrahedra'
  degauss=0.01
  nbnd=16
/
&ELECTRONS
  conv_thr =   1.0d-10
  !startingpot='file'
  !startingwfc='file'
/
ATOMIC_SPECIES
B     10.0 B_pbe_v1.01.uspp.F.UPF
N     14.0 N.oncvpsp.upf
ATOMIC_POSITIONS crystal
B            0.6666666667       0.3333333333       0.7500000000 
B            0.3333333333       0.6666666667       0.2500000000 
N            0.6666666667       0.3333333333       0.2500000000 
N            0.3333333333       0.6666666667       0.7500000000 
CELL_PARAMETERS angstrom
     -1.2562141150      -2.1758266724      -0.0000000000
     -1.2562141150       2.1758266724       0.0000000000
      0.0000000000       0.0000000000      -7.7072650000
K_POINTS automatic
20 20 6 0 0 0

DOS Input:

&DOS
  prefix='hBN'
  outdir='./out/'
/

Truncated Output:

     highest occupied, lowest unoccupied level (ev):     3.1546    7.7535

!    total energy              =     -54.87830935 Ry
     estimated scf accuracy    <          5.8E-12 Ry

     The total energy is the sum of the following terms:
     one-electron contribution =     -27.98073829 Ry
     hartree contribution      =      21.53196044 Ry
     xc contribution           =     -18.59693952 Ry
     ewald contribution        =     -29.83259198 Ry

     convergence has been achieved in  11 iterations

  PWSCF        :   1m54.99s CPU   2m 1.43s WALL

Python script to plot the DOS:

import numpy as np 
import matplotlib.pyplot as plt 

for ws in [10]:         # This is the window size 
    data = np.loadtxt('hBN.dos', skiprows=1, dtype='float') 

    mov_avg_up = [] 
    arr_up = data[:,1] 

    i = 0 
    while i < len(arr_up)-ws + 1: 
        window = arr_up[i : i + ws] 
        avg = sum(window)/ws 
        mov_avg_up.append(avg) 
        i += 1 

    fig,ax = plt.subplots(figsize=(15,10)) 

    plt.xlabel("Energy (eV)", fontsize=15) 
    plt.ylabel("DOS (states/eV)", fontsize=15) 
    plt.ylim(-.25,4)
    #plt.xlim(-50,50)

    fermi=3.235

    xval = data[int(ws/2)-1:int(-1*ws/2),0]

    ax.plot(xval-fermi, mov_avg_up, label="hBN DOS", linewidth=2, color='green')

    plt.title("DOS of hBN, PBE (moving average windows size="+str(ws)+")", fontsize=20)
    
    plt.fill_between(xval-fermi,0,mov_avg_up,where=(xval <fermi), facecolor='green', alpha=0.1)

    plt.axvline(x=0, linewidth=0.5, color='k', linestyle=(0, (8, 10)))
    plt.text(0.1, 1.7, 'Fermi energy', rotation=90, fontsize=15)

    ax.legend() 

plt.show()
$\endgroup$
2
  • $\begingroup$ Thank you for your effort. Are the pseudopotentials you used Ultrasoft? Could you please let me know the K_POINTS values you used with the nscf calculations? $\endgroup$
    – Bekaso
    Dec 12, 2023 at 16:12
  • $\begingroup$ I have edited my answer and added the nscf and dos input files as well as the python script to plot the DOS diagram from the output. I used the pseudopotentials from SSSP precision library where B has a ultrasoft pseudopotential and N has optimized norm-conservinng Vanderbilt pseudopotential. The K_POINTS value has been chosen to get a very fine mesh (i.e. 1 k-point per 0.15 Å with a smearing of 0.1 eV) as suggested in QE input generator. $\endgroup$ Dec 13, 2023 at 7:03

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .