Silicon's band structure, computed with Elk and Quantum ESPRESSO together -- a type 3, comparison calculation with the QELK suite. Same cell, same k-path, same functional on both sides; only the method changes (all-electron LAPW against pseudopotential plane waves).
The input file
We need only the input file
Si_band_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 bands_elk elk
from initial
task 0 20
rgkmax 7.0
ngridk 6 6 6
xctype 2
npp1d 200
bandpath
0.0 0.5 0.0 L
0.0 0.0 0.0 G
0.5 0.5 0.0 X
0.5 0.75 0.25 W
0.0 0.0 0.0 G
end
end
step scf_qe qe
from initial
calculation scf
ecutwfc 20
kpoints 6 6 6
end
step bands_qe qe
from initial
calculation bands
ecutwfc 20
nbnd 8
kpath
0.35355 0.35355 0.35355 1 L
0.28284 0.28284 0.28284 2
0.21213 0.21213 0.21213 3
0.14142 0.14142 0.14142 4
0.07071 0.07071 0.07071 5
0.00000 0.00000 0.00000 6 G
0.00000 0.00000 0.07071 7
0.00000 0.00000 0.14142 8
0.00000 0.00000 0.21213 9
0.00000 0.00000 0.28284 10
0.00000 0.00000 0.35355 11
0.00000 0.00000 0.42426 12
0.00000 0.00000 0.49497 13
0.00000 0.00000 0.56569 14
0.00000 0.00000 0.63640 15
0.00000 0.00000 0.70711 16 X
0.00000 0.07071 0.70711 17
0.00000 0.14142 0.70711 18
0.00000 0.21213 0.70711 19
0.00000 0.28284 0.70711 20
0.00000 0.35355 0.70711 21 W
0.00000 0.28284 0.56569 22
0.00000 0.21213 0.42426 23
0.00000 0.14142 0.28284 24
0.00000 0.07071 0.14142 25
0.00000 0.00000 0.00000 26 G
end
end
Note: nothing is chained with "from" between Elk and QE -- both start from the SAME initial cell, independently. What is kept identical on purpose: the path (L-G-X-W-G), the SCF mesh (6x6x6) and the functional (xctype 2, LDA, matched to Si.pz-vbc.UPF), so a difference in the bands is a difference between the methods, not a setup mismatch. Elk's task 0 20 is one LIST in one step (task 20 reads the STATE.OUT that task 0 just wrote); QE needs a separate scf step, and bands_qe finds its density through the shared outdir/prefix.
Note -- the k-point coordinates are NOT the same numbers in the two blocks. The two codes read k-points in different coordinate systems, and the path is only the same if each block is written in its own:
Elk bandpath : lattice coordinates (fractions of the reciprocal
primitive vectors of this cell)
QE kpath : Cartesian, in units of 2*pi/alat, where alat is the LENGTH
OF THE FIRST CELL VECTOR (7.2125 bohr here, printed in
scf_qe.pw.out as "lattice parameter (alat)"), not the
cubic constant 10.2 bohr
point Cartesian [2pi/a, a=10.2] Elk bandpath QE kpath
L 0.5 0.5 0.5 0.0 0.5 0.0 0.35355 (x3)
G 0 0 0 0 0 0 0
X 0 0 1 0.5 0.5 0.0 0 0 0.70711
W 0 0.5 1 0.5 0.75 0.25 0 0.35355 0.70711
Q
E's numbers are the Cartesian ones times alat/a = 0.70711. Writing the familiar 0.5 0.5 0.5 / 0 0 1 / 0 0.5 1 into both blocks sends the two codes along different lines: Elk would run along the G-L line and never reach X or W, QE would overshoot X, W and L by a factor of 1.41 -- and the overlay would compare unrelated k-points. (An earlier version of this example did exactly that; checking the k-points that actually reach elk.in and pw.in showed it, and the comparison below is the corrected one.)
Calculation
algerien1970@linux-ml0a:~/abinitio/QELK/examples/Si_band_elk_vs_qe> qelk Si_band_elk_vs_qe.qelk
A first check with --dry-run writes elk.in and both pw.in files without launching any binary:
algerien1970@linux-ml0a:~/abinitio/QELK/examples/Si_band_elk_vs_qe> qelk Si_band_elk_vs_qe.qelk --dry-run
The steps run in file order: run/00_bands_elk, run/01_scf_qe, run/02_bands_qe. The important files are also copied next to the .qelk file, prefixed with the step id (bands_elk.BAND.OUT, bands_elk.INFO.OUT, scf_qe.pw.out, bands_qe.pw.out, ...).
SCF results
grep -A1 "Convergence targets achieved" bands_elk.INFO.OUT grep " total energy" bands_elk.INFO.OUT | tail -1 grep " Fermi " bands_elk.INFO.OUT | tail -1 grep "Estimated" bands_elk.INFO.OUT | tail -2
total energy : -578.076131082 Fermi : 0.214810646258 Estimated indirect band gap : 0.2103472743E-01 (from k-point 1 to k-point 13) Estimated direct band gap : 0.9350519618E-01 (at k-point 1)
grep -E "convergence has been achieved|^!|Fermi energy" scf_qe.pw.out
the Fermi energy is 6.5511 ev
! total energy = -15.84588576 Ry
convergence has been achieved in 5 iterations
Both SCF runs use 16 irreducible k-points (6x6x6). Elk's targets were met at iteration 18. The total energies are on different zeros (all-electron against pseudopotential) and are not meant to agree.
The overlay
overlay_bands.py bands_qe.pw.out bands_elk.BAND.OUT --nocc 4 -o overlay_bands.png
QE: VBM = 6.2534 eV CBM = 6.7541 eV gap = 0.5007 eV Elk: VBM = -0.2733 eV CBM = 0.1645 eV gap = 0.4378 eV wrote overlay_bands.png
Note: the two codes share no absolute energy zero, so each is shifted so its own valence-band maximum is 0 eV; --nocc 4 is the number of doubly occupied bands (8 valence electrons / 2). The x-axis of the two codes is rescaled to a common axis, which is only meaningful because the two paths are now geometrically identical (same vertices, same segment lengths).
Analysis
The two band structures lie on top of each other over the whole path. To put numbers on it, Elk's 200-point curve was read at the 26 k-points of the QE path (same geometry), each code aligned on its own VBM, and compared band by band, Elk minus QE:
band(s) rms difference largest difference 1-4 (valence) 0.022 eV 0.041 eV 5-6 (lowest conduction) 0.036 eV 0.073 eV 7-8 (upper conduction) 0.102 eV 0.223 eV (at X)
The valence bands agree to a few hundredths of an eV; the difference grows with energy, as expected, because higher states depend more on what differs between the methods -- Elk's finite angular-momentum expansion inside the muffin-tin spheres and all core electrons, against QE's plane-wave cutoff and valence-only pseudopotential. At the high-symmetry points, energies in eV relative to each code's VBM:
bands 1 ... 8
G QE -12.06 0.00 0.00 0.00 2.56 2.56 2.56 3.47
Elk -12.09 0.00 0.00 0.00 2.54 2.54 2.54 3.38
X QE -7.85 -7.85 -2.92 -2.92 0.63 0.63 10.16 10.16
Elk -7.88 -7.88 -2.91 -2.91 0.57 0.57 10.28 10.28
L QE -9.67 -7.07 -1.22 -1.22 1.56 3.34 3.34 7.58
Elk -9.71 -7.09 -1.22 -1.22 1.49 3.33 3.33 7.59
W QE -7.67 -7.67 -3.96 -3.96 4.23 4.23 5.03 5.03
Elk -7.71 -7.71 -3.96 -3.96 4.26 4.26 5.04 5.04
Both codes reproduce the same degeneracy pattern independently. At G it is 1 + 3 + 3 + 1: one low band, the threefold valence top, the threefold conduction level, and one higher band -- the same pattern found separately in the Elk-only and QE-only band tutorials. At X and W every band is doubly degenerate in both codes, as expected for the diamond structure. The valence bandwidth (bottom of band 1 at G to the VBM) is 12.06 eV in QE and 12.09 eV in Elk.
The gap, along this path:
Elk QE VBM at G at G (k-point 6) CBM 0.84 of the way G-X 0.80 of the way G-X (k-point 14) indirect gap on path 0.4378 eV 0.5007 eV direct gap at G 2.54 eV 2.56 eV
Both codes put the valence maximum at G and the conduction minimum along G-X, so both give silicon's indirect gap, with the minimum where the real conduction valley lies (about 85% of the way to X). QE's path has a point only every 0.1 of G-X, so its minimum can only be 0.8 or 0.9; Elk's 200-point path finds it at 0.84. That sampling accounts for only about 12 meV in Elk's own curve (0.4497 eV at 0.8, 0.4378 eV at 0.84). The remaining difference is a real, small method difference: at the same point, 0.8 of G-X, Elk gives 0.450 eV and QE 0.501 eV, and at X itself the lowest conduction level is 0.57 eV (Elk) against 0.63 eV (QE) -- an offset of about 50-60 meV along the whole G-X line. Both values are at the usual LDA level for silicon, well below the experimental gap of 1.17 eV.
Note: the gap printed by any SCF run depends on how it is measured. Elk's own whole-zone estimate from the 6x6x6 SCF mesh is 0.5724 eV (from k-point 1 to k-point 13), larger than the 0.4378 eV found on the path, because a coarse mesh contains no point near 0.84 of G-X. With denser meshes the same estimate falls toward the path value: 0.4916 eV at 8x8x8 (Si_dos_elk) and 0.4444 eV at 12x12x12 (Si_dos_elk_vs_qe). For a gap, use a path or mesh that samples the conduction valley, and treat a number read off a coarse mesh as an upper bound.
Note: two things in the figure are about the calculation, not the physics. QE's curves show small kinks because QE was asked for 26 k-points, against Elk's 200. And the overlay compares only the eight bands QE computed (nbnd 8): Elk also computed bands 9 to 13, which appear alone at the top of the plot with no QE partner.
The point of the example is the workflow: one input file, two very different methods, and a single command that puts their band structures on the same axes -- provided each code is given the same k-points in its own coordinates.
My name is Dr Abderrahmane Reggad. I'm a 54 year old and I am a researcher in materials' sciences. 
0 Comments