How to do a comparaison study of the Total DOS of Silicon

Silicon's total density of states, Elk versus Quantum ESPRESSO -- a type 3, comparison calculation with the QELK suite: both codes start from the same cell, neither uses the other's result, and the two DOS curves are overlaid on one figure. The single-code versions are the Si_tdos_elk and Si_tdos_qe examples.

The input file

We need only the input file

Si_tdos_elk_vs_qe.qelk

version 1

cell bohr
  -5.1   0.0   5.1
   0.0   5.1   5.1
  -5.1   5.1   0.0
end

atoms crystal
  Si  0.00  0.00  0.00
  Si  0.25  0.25  0.25
end

species Si
  elk_file   Si.in
  qe_pseudo  Si.pz-vbc.UPF
  qe_mass    28.086
end

step tdos_elk elk
  from initial
  task    0 10
  rgkmax  7.0
  ngridk  12 12 12
  xctype  2
end

step scf_qe qe
  from initial
  calculation  scf
  ecutwfc      20
  kpoints      6 6 6
end

step nscf_qe qe
  from scf_qe
  calculation  nscf
  ecutwfc      20
  kpoints      12 12 12
end

step tdos_qe tdos
  from nscf_qe
  fildos  si.dos
end

Note: the Elk step and the first QE step say from initial, and the other QE steps continue from the QE step before them (nscf_qe from scf_qe, tdos_qe from nscf_qe) -- Elk never hands anything to QE or the other way round, which is what makes this a comparison and not a combined calculation. Elk's task 0 10 is one LIST in one step (task 10 reads the STATE.OUT that task 0 just wrote). On the QE side the nscf step uses the same 12x12x12 mesh as Elk, and tdos_qe runs dos.x on the nscf wavefunctions, which it reads through the shared outdir/prefix (declared with from nscf_qe; dos.x takes no geometry) -- no hand-written dos.in. xctype 2 (Elk) and Si.pz-vbc.UPF (QE) are both LDA.

Calculation

algerien1970@linux-ml0a:~/abinitio/QELK/examples/Si_tdos_elk_vs_qe> qelk Si_tdos_elk_vs_qe.qelk

A first check with --dry-run writes elk.in, pw.in and dos.in for all four steps without launching any binary:

algerien1970@linux-ml0a:~/abinitio/QELK/examples/Si_tdos_elk_vs_qe> qelk Si_tdos_elk_vs_qe.qelk --dry-run

The steps run in file order, so the directories are run/00_tdos_elk, run/01_scf_qe, run/02_nscf_qe and run/03_tdos_qe. The important files are also copied next to the .qelk file, prefixed with the step id (tdos_elk.TDOS.OUT, scf_qe.pw.out, nscf_qe.pw.out, tdos_qe.dos.out, ...). The DOS data of QE, si.dos, is not copied up: it stays in run/03_tdos_qe/si.dos.

SCF results

Elk:

grep -A1 "Convergence targets achieved" tdos_elk.INFO.OUT
grep " total energy" tdos_elk.INFO.OUT | tail -1
grep " Fermi " tdos_elk.INFO.OUT | tail -1
grep "Estimated" tdos_elk.INFO.OUT | tail -2

 total energy                :     -578.077094713     (Hartree)
 Fermi                       :     0.212309003552     (Hartree)
 Estimated indirect band gap :   0.1633150043E-01     (from k-point 1 to k-point 40)
 Estimated direct band gap   :   0.9365530730E-01     (at k-point 1)

72 irreducible k-points, converged in 18 iterations (a 19th loop is run after the targets are met).

QE:

grep -E "convergence has been achieved|^!|Fermi energy" scf_qe.pw.out
grep -E "number of k points|Fermi energy" nscf_qe.pw.out

     the Fermi energy is     6.5511 ev
!    total energy              =     -15.84588576 Ry
     convergence has been achieved in   5 iterations
     number of k points=    72  Methfessel-Paxton smearing, width (Ry)=  0.0100
     the Fermi energy is     6.4867 ev

Note: the two QE Fermi energies differ because the scf step uses 6x6x6 (16 k-points) and the nscf step 12x12x12 (72 k-points). dos.x reads the nscf wavefunctions, so the EFermi in si.dos's header should be the nscf value, 6.4867 eV. Elk's and QE's total energies (-578.0771 Hartree versus -15.8459 Ry) are not on a comparable scale -- all-electron versus pseudopotential zeros of energy -- and are not meant to agree.

Plot

Each code on its own:

plot_dos.py run/03_tdos_qe/si.dos -o qe_dos.png

plot_elk_dos.py run/00_tdos_elk/TDOS.OUT -o elk_dos.png

And the overlay:

overlay_dos.py run/03_tdos_qe/si.dos run/00_tdos_elk/TDOS.OUT -o overlay_dos.png

Note: the two codes have no shared absolute energy zero, so each curve is aligned to its own Fermi level at 0 eV. QE's value is read from si.dos's header; Elk writes TDOS.OUT already relative to its Fermi level, so it needs no shift -- the same point made in Si_tdos_elk and Si_band_elk.

Analysis

The numbers that can be compared directly, both codes on the same 12x12x12 mesh (72 irreducible k-points each):

                      Elk                    QE (nscf_qe.pw.out)
indirect gap          0.4444 eV              0.4907 eV
                      (k-point 1 -> 40)      (VBM 6.2534 eV at k-point 1,
                                              CBM 6.7441 eV at k-point 40)
direct gap            2.5485 eV (k-point 1)  2.5638 eV (k-point 1)
bottom of valence     -12.19 eV (DOS onset)  -12.30 eV (lowest band)
band, rel. E_Fermi
smearing              Fermi-Dirac,           Methfessel-Paxton (order 1),
                      0.001 Ha (~316 K)      0.01 Ry (dos.x: ngauss=1,
                                              degauss=0.010)

Both codes put the valence maximum at G and the conduction minimum away from it, so both give silicon's indirect gap; the direct gaps at G agree to 16 meV and the valence bandwidths to about 0.1 eV. The indirect gaps differ by 46 meV, about the size of the change Elk itself shows between an 8x8x8 SCF mesh (0.4916 eV in Si_tdos_elk) and the 12x12x12 mesh used here (0.4444 eV), so it is not evidence that either code is wrong. Both values sit at the usual LDA level for silicon, well below the experimental 1.17 eV.

Elk's TDOS.OUT runs from -13.61 to +13.55 eV around its Fermi level and integrates to exactly 8.000 states up to E_Fermi (2 atoms x 4 valence electrons). Its DOS is zero from about -0.16 to +0.22 eV. QE's Fermi level sits 0.23 eV above the valence maximum and 0.26 eV below the conduction minimum, so in both codes E_Fermi lies inside the gap and aligning the curves on it is reasonable to about a tenth of an eV.

Elk peaks (eV from E_Fermi, states/eV):
  valence      : -9.90 (0.87), -7.02 (1.55), -3.05 (1.54), -2.78 (1.47)
  valence dip  : -8.00 (0.006) -- the near-empty gap between the deep
                 s-derived band and the p-derived manifold
  conduction   : +2.18 (2.02), +3.32 (2.21), +6.75 (1.52), +11.16 (1.60)

With each curve aligned only on its own Fermi level, the overlay agrees on every large feature, with no further adjustment: the gap sits at 0 for both codes; the bottom of the valence band is at about -12.3 eV in both; both curves drop almost to zero at about -8 eV, between the deep s-derived band (broad hump around -10 eV) and the p-derived valence manifold; both have a sharp peak at about -7 eV; and the broad structure between -4.5 and -1.5 eV, with Elk's peak at -3.05 eV, is followed by QE's spikes in the same region. Above the gap both show the strong group of peaks at about +2 to +3.5 eV, and the peak of Elk at +6.75 eV has a counterpart in QE just below +7 eV (positions read off the figure, so approximate).

Note: two differences are about the calculation, not the physics. (1) QE's curve is much spikier: dos.x broadens only 72 k-points with a 0.01 Ry (about 0.14 eV) Methfessel-Paxton smearing, which also lets the curve dip very slightly below zero (visible near -7.7 eV and again just above +9 eV), while Elk's curve is smoother. Compare peak positions rather than heights or fine structure. (2) Above about +7.5 eV the curves separate: QE's DOS falls away and is gone by +10 eV, while Elk stays between 1.1 and 1.3 states/eV from +7.5 to +10.5 eV and continues to +13.5 eV. This is a limit of the QE input, not of silicon: the nscf step keeps pw.x's default of 8 bands (4 occupied and 4 empty), whose highest energy is 9.92 eV above E_Fermi, so the higher conduction bands are simply missing from QE's DOS. Trust the comparison up to about +7 eV, and add nbnd 16 to the nscf_qe step (the same keyword Si_band_qe uses in its bands step) to extend it.

The point of the example is the workflow: two very different methods, one input file, four steps, and a single command that puts their DOS on the same energy axis.