How to calculate the Band Structure of Silicon

Band structure of silicon with the QELK suite -- a combined calculation: Elk runs the SCF, then Quantum ESPRESSO computes the band structure on that same structure.

 

The three QELK calculation types

 

Every QELK example name says which of the three kinds it is:

1. Independent   <System>_<task>_<code>                     e.g. Si_band_qe
2. Combined      <System>_<task1>_<code1>_<task2>_<code2>    e.g. Si_scf_elk_band_qe
3. Comparison    <System>_<task>_elk_vs_qe                  e.g. Si_band_elk_vs_qe

Today we do type 2: Elk converges the ground state first, then QE continues from that structure to compute the bands.

 

Why three steps, not two

 

QE's band-structure calculation cannot read Elk's charge density -- an all-electron (LAPW) density and a pseudopotential plane-wave density are not interchangeable formats. So QE has to build its own density first with a quick SCF step, the same way the independent Si_band_qe example hides its own SCF step under the name “band_qe”. Here that becomes three chained steps:

scf_elk    Elk,  task 0        : converges the ground state (all-electron)
scf_qe     QE,   calculation=scf : converges QE's own density on Elk's structure
bands_qe   QE,   calculation=bands : the band structure, L - G - X - W - G

The last two steps share the suite's default outdir/prefix, so bands_qe automatically finds scf_qe's density -- no extra wiring needed.

 

The input file

 

We need only the input file

Si_scf_elk_band_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

# Elk converges the ground state at the all-electron level
step scf_elk elk
  from initial
  task    0
  rgkmax  7.0
  ngridk  4 4 4
  xctype  2
end

# QE needs its own converged density before it can compute bands
step scf_qe qe
  from scf_elk
  calculation  scf
  ecutwfc      20.0
  kpoints      4 4 4
end

# Band structure along L - Gamma - X - W - Gamma
step bands_qe qe
  from scf_qe
  calculation  bands
  ecutwfc      20.0
  nbnd         8
  kpath
    0.5 0.5 0.5 1  L
    0.4 0.4 0.4 2
    0.3 0.3 0.3 3
    0.2 0.2 0.2 4
    0.1 0.1 0.1 5
    0.0 0.0 0.0 6  G
    0.0 0.0 0.1 7
    0.0 0.0 0.2 8
    0.0 0.0 0.3 9
    0.0 0.0 0.4 10
    0.0 0.0 0.5 11
    0.0 0.0 0.6 12
    0.0 0.0 0.7 13
    0.0 0.0 0.8 14
    0.0 0.0 0.9 15
    0.0 0.0 1.0 16  X
    0.0 0.1 1.0 17
    0.0 0.2 1.0 18
    0.0 0.3 1.0 19
    0.0 0.4 1.0 20
    0.0 0.5 1.0 21  W
    0.0 0.4 0.8 22
    0.0 0.3 0.6 23
    0.0 0.2 0.4 24
    0.0 0.1 0.2 25
    0.0 0.0 0.0 26  G
  end
end

Note: xctype 2 (Elk) and Si.pz-vbc.UPF (QE) are both LDA, so the two sides of this combined calculation stay on the same exchange-correlation footing. The cell and k-path match the suite's own Si_band_elk / Si_band_qe examples.

 

Calculation

 

Run it from inside the QELK tree:

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

A first check with --dry-run writes elk.in / pw.in for all three steps without launching either binary:

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

 

Results

 

Every step's files land in run/NN_<step>/, and the important ones are also copied next to the .qelk file, prefixed with the step id.

Elk SCF convergence:

cat scf_elk.INFO.OUT

QE SCF -- convergence, total energy, Fermi energy:

grep -E "convergence has been achieved|^!|Fermi energy" scf_qe.pw.out

Band energies at every k-point of the path:

less bands_qe.pw.out

 

Plot

 

Now we can use the suite's own plotting script to get a labelled band structure:

algerien1970@linux-ml0a:~/abinitio/QELK/examples/Si_scf_elk_band_qe> plot_qe_bands.py bands_qe.pw.out --qelk Si_scf_elk_band_qe.qelk --step bands_qe -o si_bands.png

--qelk / --step read the L / G / X / W labels straight from the .qelk file's kpath block, so nothing needs to be typed by hand. QE's bands sit on QE's own absolute eV scale (unlike Elk's, which are relative to the Fermi level), so pull the Fermi energy out of scf_qe.pw.out and pass --efermi <value> to put it at 0 eV, plus --emin/--emax to set the plotted window.

 

Note: Since this is SCF only on both sides (no relaxation), the geometry never changes across the three steps -- the point of this example is the Elk to QE workflow itself, not a structural update. The "from" keyword in a .qelk file only ever carries geometry between steps; the density that bands_qe actually reuses comes from sharing scf_qe's outdir/prefix by default, not from "from" itself.