Gallium arsenide's partial density of states (projected on species and orbital), 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 sets of curves are overlaid. This is the first example with two different elements and with GGA (PBE) instead of LDA.
The input file
We need only the input file
GaAs_pdos_elk_vs_qe.qelk
version 1 cell bohr -5.3415 0.0000 5.3415 0.0000 5.3415 5.3415 -5.3415 5.3415 0.0000 end atoms crystal Ga 0.00 0.00 0.00 As 0.25 0.25 0.25 end species Ga elk_file Ga.in qe_pseudo Ga.pbe-dn-kjpaw_psl.0.2.upf qe_mass 69.723 end species As elk_file As.in qe_pseudo As.pbe-n-kjpaw_psl.0.2.upf qe_mass 74.922 end step pdos elk from initial task 0 10 rgkmax 7.0 ngridk 12 12 12 xctype 20 dosmsum 1 wplot -0.85 0.5 end step scf_qe qe from initial calculation scf ecutwfc 40 kpoints 12 12 12 nbnd 13 occupations tetrahedra_opt end step nscf_qe qe from scf_qe calculation nscf ecutwfc 40 kpoints 12 12 12 nbnd 13 occupations tetrahedra_opt end step pdos_qe pdos from nscf_qe filpdos gaas end
Note: the cell is the zincblende primitive cell with a = 10.683 bohr (5.653 Å), Ga at (0,0,0) and As at (1/4,1/4,1/4); xctype 20 (Elk, PBE) matches the two PBE PAW pseudopotentials. On the Elk side, task 0 10 is one LIST in one step, dosmsum 1 sums the m sub-levels so Elk writes one curve per orbital (s, p, d) of each atom, and wplot -0.85 0.5 (Hartree) widens the energy window so that the Ga 3d semicore level is not clipped.
Note: on the QE side the pdos needs scf, then nscf on a denser mesh, then projwfc.x (the pdos code). Each step declares where it continues from: nscf_qe from scf_qe, pdos_qe from nscf_qe. The pdos step takes no geometry: projwfc.x reads the energies and wavefunctions that the nscf step left in the shared outdir/prefix. nbnd 13 gives four empty bands above the nine occupied ones, and tetrahedra_opt suits a DOS.
Calculation
algerien1970@linux-ml0a:~/abinitio/QELK/examples/GaAs_pdos_elk_vs_qe> qelk GaAs_pdos_elk_vs_qe.qelk
A first check with --dry-run writes elk.in, both pw.in files and projwfc.in without launching any binary:
algerien1970@linux-ml0a:~/abinitio/QELK/examples/GaAs_pdos_elk_vs_qe> qelk GaAs_pdos_elk_vs_qe.qelk --dry-run
The steps run in file order: run/00_pdos, run/01_scf_qe, run/02_nscf_qe, run/03_pdos_qe. Important files are also copied next to the .qelk file with the step id as prefix. Elk's PDOS_S*_A*.OUT files and QE's gaas.pdos_* files stay in their run folders.
SCF results
Elk:
grep -A1 "Convergence targets achieved" pdos.INFO.OUT grep " total energy" pdos.INFO.OUT | tail -1 grep " Fermi " pdos.INFO.OUT | tail -1 grep "Estimated" pdos.INFO.OUT | tail -2
total energy : -4205.24110722 Fermi : 0.186696425500 Estimated indirect band gap : 0.1890510335E-01 (from k-point 1 to k-point 1) Estimated direct band gap : 0.1890510335E-01 (at k-point 1)
Converged in 19 loops. Elk's valence charge is 28 (Ga 13 + As 15, both 3d shells in the valence); the plotted window starts at -23 eV, so the As 3d shell (near -35 eV) is outside it.
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
! total energy = -453.96449587 Ry
the Fermi energy is 7.8988 ev
convergence has been achieved in 8 iterations
number of k points= 72 (tetrahedron method)
the Fermi energy is 7.8990 ev
Note: the scf and nscf Fermi energies agree (7.8988 and 7.8990 eV) because both use the same 12x12x12 mesh; the nscf value is the one used below, since projwfc.x works on the nscf wavefunctions. Elk's and QE's total energies (-4205.2411 Ha against -453.9645 Ry) are on different zeros -- all-electron against pseudopotential -- and are not meant to agree.
Plot
Each code on its own (E_Fermi at 0 eV, window -12 to +12 eV on both):
plot_elk_pdos.py run/00_pdos --l-resolved -o elk_pdos_GaAs.png --xlim -12 12
plot_qe_pdos.py run/03_pdos_qe --filpdos gaas --l-resolved --efermi $(grep "Fermi energy" nscf_qe.pw.out | tail -1 | awk '{print $5}') --xlim -12 12 -o qe_pdos_GaAs.png
The overlay
overlay_pdos.py run/00_pdos run/03_pdos_qe --qe-filpdos gaas --qe-efermi $(grep "Fermi energy" nscf_qe.pw.out | tail -1 | awk '{print $5}') -o overlay_pdos_GaAs_valence.png
By default the overlay keeps the s and p orbitals; add --lmax 2 to include the d channel and see the Ga 3d semicore peak:
overlay_pdos.py run/00_pdos run/03_pdos_qe --qe-filpdos gaas --qe-efermi $(grep "Fermi energy" nscf_qe.pw.out | tail -1 | awk '{print $5}') --lmax 2 -o overlay_pdos_GaAs_all.png
Note: QE's energies are not Fermi-referenced, so the Fermi energy is required; Elk's are already relative to E_Fermi. Each curve is aligned on its own Fermi level.
Analysis
The gap. The two codes give the same direct gap at G on the same 12x12x12 mesh: Elk 0.5144 eV (0.01890510 Ha, k-point 1), QE 0.5146 eV (valence maximum 7.6149 eV, conduction minimum 8.1295 eV, both at k-point 1 of nscf_qe.pw.out). QE's Fermi level sits 0.284 eV above the valence maximum and 0.231 eV below the conduction minimum, so aligning the curves on E_Fermi puts the gap at 0 in both. Both values are the usual PBE level, well below the experimental 1.42 eV.
Shape. The overlay agrees on every main feature: the gap at 0 eV; the narrow As-s-derived band near -10.5 eV and its shoulder up to -13 eV; the sharp Ga-s peak near -6.9 eV followed by the broad p-derived valence band between about -4 and 0 eV; and the strong group of conduction peaks between +2 and +7 eV. Above about +10 eV the QE curves end, while Elk continues to 13.5 eV: this is the limit of nbnd 13 (the highest QE band lies 9.9 eV above E_Fermi). Add nbnd 16 to the scf and nscf steps to extend it.
The Ga 3d semicore peak. With --lmax 2 the Ga 3d peak dwarfs everything else. Its position agrees to 0.1 eV (-14.97 eV Elk, -15.07 eV QE from E_Fermi). Its HEIGHT does not: about 41 states/eV in Elk and 152 states/eV in QE. This is not a normalisation problem. Elk writes 500 energy points across the 1.35 Hartree window (0.073 eV per point) and smooths them, QE bins at 0.01 eV, and the 3d peak is only about 0.2 eV wide, so the height depends on the grid. What is comparable is the area (next).
Integrating the curves. Integrating each curve up to E_Fermi checks both the normalisation and the electron count:
Elk (inside muffin-tin spheres) QE (Lowdin, projwfc.out) Ga s 0.812 1.084 Ga p 0.731 1.669 Ga d 9.933 9.995 As s 1.338 1.437 As p 1.907 2.835 total projected 14.79 17.02 total DOS to E_Fermi 18.000 18.000 (see below)
Elk's total DOS integrates to exactly 18.000 electrons up to E_Fermi, the number of valence electrons with this window (Ga 13 + As 5; the As 3d shell lies below it). QE's total DOS integrates to 17.263 over its energy grid, and the full grid, which should hold 13 bands x 2 = 26 states, gives 25.263: the same 0.737 is missing in both, which is consistent with the Ga 3d bands starting right at the lowest energy of the grid (-7.30 eV), so that the first part of the peak falls below it. The Lowdin charges in projwfc.out give the Ga 3d population directly: 9.995, against 9.933 in Elk -- the 3d shell holds its ten electrons in both codes. Elk's integrals are also consistent: its PDOS total inside the spheres, 14.79, is less than 18 because the remainder lies in the interstitial region, which no sphere-projected curve counts.
So the PDOS of the two codes carry the same electrons. They differ in how much of each orbital they count: Elk integrates the charge inside a muffin-tin sphere of radius 2.29 bohr around each atom; QE projects the wavefunctions on the atomic pseudo-orbitals and, with a spilling parameter of 0.054, accounts for 94.6% of the 18 electrons. The ratio QE/Elk is 1.01 for Ga d, 1.07 for As s and 1.33 for Ga s, but 1.49 for As p and 2.28 for Ga p: the s+p curves of the overlay lie about 1.3 to 1.8 times higher in QE, while the d populations coincide. It is a difference between projection methods and has no single factor, so compare peak positions and integrals, not heights.
Note: comparing two PDOS curves point by point only makes sense on the same energy grid and with the same projection. For a like-for-like comparison of the total DOS, use the total-DOS examples (Si_tdos_elk_vs_qe); the PDOS comparison is about orbital character and peak positions.
The point of the example is the workflow: two very different methods, one input file, four steps, and a single command that puts their partial DOS on the same energy axis.
My name is Dr Abderrahmane Reggad. I'm a 54 year old and I am a researcher in materials' sciences. 
0 Comments