How to Calculate Band Structures in Quantum ESPRESSO

The electronic band structure is the single most informative quantity you can extract from a density functional theory calculation. It tells you whether a material is a metal, semiconductor, or insulator, where the band gap lives, and how carriers will move. This tutorial walks through the full Quantum ESPRESSO workflow for computing and plotting a band structure, using crystalline silicon as a worked example.

Understanding the Band Structure Workflow

A band structure calculation in Quantum ESPRESSO is never a single step. The charge density and the Kohn–Sham eigenvalues \(\varepsilon_{n\mathbf{k}}\) along a \(\mathbf{k}\)-path are computed in separate runs because they have fundamentally different sampling requirements. Getting the sequence right is the difference between a clean dispersion plot and hours of confusion.

The workflow has four stages:

flowchart LR
  A["pw.x scf<br/>uniform k-grid"] --> B[Converged density]
  B --> C["pw.x bands<br/>high-symmetry path"]
  C --> D["bands.x<br/>reorder crossings"]
  D --> E["Plot bands and extract gap"]

1. SCF calculation (pw.x, calculation='scf') on a uniform Monkhorst–Pack grid to converge the self-consistent charge density.
2. Bands (non-SCF) calculation (pw.x, calculation='bands') that reads the fixed density and diagonalizes the Hamiltonian at explicit \(\mathbf{k}\)-points along a high-symmetry path.
3. Post-processing (bands.x) to reorder eigenvalues into continuous bands and write a gnuplot-ready file.
4. Plotting and analysis, including extracting the band gap.

The key idea: the density is variationally determined once, on a dense uniform grid, and then held frozen while you sample the special path through the Brillouin zone that reveals the dispersion.

Step 1: The Self-Consistent Field Calculation

The SCF run establishes the ground-state charge density. For silicon in the diamond structure, use the experimental lattice constant of 5.43 angstrom (10.26 bohr) and a well-converged k-grid.

&CONTROL
    calculation = 'scf'
    prefix      = 'si'
    outdir      = './out'
    pseudo_dir  = './pseudo'
    verbosity   = 'high'
/
&SYSTEM
    ibrav       = 2
    celldm(1)   = 10.26
    nat         = 2
    ntyp        = 1
    ecutwfc     = 40.0
    ecutrho     = 320.0
/
&ELECTRONS
    conv_thr    = 1.0d-8
    mixing_beta = 0.7
/
ATOMIC_SPECIES
 Si  28.0855  Si.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS (alat)
 Si  0.00  0.00  0.00
 Si  0.25  0.25  0.25
K_POINTS (automatic)
 8 8 8 0 0 0

Note ibrav = 2 selects the FCC Bravais lattice, so only the two-atom primitive basis is needed. The ecutrho value is 8x ecutwfc because we are using a PAW pseudopotential; for norm-conserving potentials, ecutrho defaults to 4x. Always converge both cutoffs and the k-grid before trusting any downstream result.

Run it:

pw.x -in si.scf.in > si.scf.out

Check the final total energy and confirm convergence was reached. The charge density is now stored in ./out, ready to be reused.

Step 2: Defining the High-Symmetry k-Path

The band structure is plotted along straight-line segments connecting high-symmetry points in the Brillouin zone. For an FCC lattice the conventional path is L to Gamma to X to W to K back to Gamma. Quantum ESPRESSO offers two convenient ways to specify this.

Using crystal_b coordinates

The {crystal_b} option lets you list the endpoints of each segment in crystal (reciprocal lattice) coordinates, with an integer telling pw.x how many points to generate between that point and the next.

K_POINTS {crystal_b}
 5
   0.5000  0.5000  0.5000  20  ! L
   0.0000  0.0000  0.0000  20  ! Gamma
   0.5000  0.0000  0.5000  20  ! X
   0.7500  0.2500  0.7500  20  ! W
   0.0000  0.0000  0.0000  1   ! Gamma

The first line (5) is the number of high-symmetry points. Each subsequent line gives the fractional reciprocal coordinates and the number of interpolated points to the next vertex. The final point takes a weight of 1 because there is no segment following it.

The tpiba_b alternative

If you prefer Cartesian units of 2*pi/alat, use {tpiba_b} instead. Both produce identical physics; crystal_b is generally safer because it is independent of the cell orientation.

Step 3: The Bands (Non-SCF) Calculation

Now run pw.x in bands mode. It reads the frozen density from the SCF step and computes eigenvalues along the path. Request more bands than occupied states so you can see the conduction bands.

&CONTROL
    calculation = 'bands'
    prefix      = 'si'
    outdir      = './out'
    pseudo_dir  = './pseudo'
/
&SYSTEM
    ibrav       = 2
    celldm(1)   = 10.26
    nat         = 2
    ntyp        = 1
    ecutwfc     = 40.0
    ecutrho     = 320.0
    nbnd        = 12
/
&ELECTRONS
    conv_thr    = 1.0d-8
/
ATOMIC_SPECIES
 Si  28.0855  Si.pbe-n-kjpaw_psl.1.0.0.UPF
ATOMIC_POSITIONS (alat)
 Si  0.00  0.00  0.00
 Si  0.25  0.25  0.25
K_POINTS {crystal_b}
 5
   0.5000  0.5000  0.5000  20
   0.0000  0.0000  0.0000  20
   0.5000  0.0000  0.5000  20
   0.7500  0.2500  0.7500  20
   0.0000  0.0000  0.0000  1

Silicon has 8 valence electrons per primitive cell, so 4 doubly occupied bands. Setting nbnd = 12 gives 4 valence plus 8 conduction bands, more than enough to visualize the gap region.

pw.x -in si.bands.in > si.bands.out

Step 4: Post-Processing with bands.x

The bands run writes eigenvalues but not in a plot-friendly order; band crossings can scramble the sorting. The bands.x tool resolves crossings (optionally using symmetry) and produces a clean data file.

&BANDS
    prefix   = 'si'
    outdir   = './out'
    filband  = 'si_bands.dat'
/

bands.x -in si.bands.pp.in > si.bands.pp.out

This generates si_bands.dat and si_bands.dat.gnu. The .gnu file has two columns — a cumulative distance along the k-path and the energy in eV — ready for plotting.

Step 5: Finding the Fermi Level and Band Gap

Before plotting, you need the reference energy. For a semiconductor, extract the valence band maximum (VBM) and conduction band minimum (CBM). Grep the SCF output for the highest occupied and lowest unoccupied levels:

grep -A1 "highest occupied" si.scf.out

For silicon this reports the VBM near \(\Gamma\) and the CBM roughly 85% of the way from \(\Gamma\) to \(X\) — the hallmark of an indirect band gap:

PBE typically underestimates it: you will see about \(0.6\)–\(0.7\,\mathrm{eV}\) versus the \(1.12\,\mathrm{eV}\) experimental value. This is the well-known DFT band-gap problem, corrected only by hybrid functionals or \(GW\).

To get an accurate gap, rerun the bands step with input_dft = 'hse' or use the gw route through additional packages — at significantly higher cost.

Step 6: Plotting the Dispersion

A minimal gnuplot script produces a publication-quality figure. Shift energies so the VBM sits at zero.

gnuplot -persist <<'EOF'
set ylabel "Energy (eV)"
set xlabel "Wavevector"
set yrange [-13:6]
vbm = 6.25   # replace with your VBM in eV
plot "si_bands.dat.gnu" using 1:($2-vbm) with lines lc rgb "blue" notitle
EOF

Add vertical lines at the segment boundaries and label them L, Gamma, X, W, Gamma. The x-axis tick positions are printed at the top of si.bands.pp.out. Python users can load the same .gnu file with NumPy and Matplotlib for finer control over styling and gap annotation.

Reading band character from the plot

The dispersion itself carries physical meaning beyond the gap. Steeply dispersing bands correspond to light effective masses and high carrier mobility; flat bands signal heavy, localized carriers. Near a band extremum the effective mass is

In silicon, the light and heavy hole bands are degenerate at \(\Gamma\) at the valence-band maximum and split apart as you move toward \(X\) and \(L\) — a feature directly visible once you zoom into the top of the valence manifold. The curvature of the conduction band near its minimum sets the electron effective mass that feeds into transport and device models.

Orbital-resolved (fat) band structures

Sometimes you want to know which atomic orbitals dominate a given band. Running projwfc.x after the bands step and weighting each point by its orbital projection produces a fat band plot, where line thickness or color encodes s, p, or d character. This is the band-structure analog of the projected density of states and is invaluable for understanding hybridization, band inversion in topological materials, and the orbital origin of band edges.

Common Pitfalls and Best Practices

• Insufficient SCF k-grid: a sparse grid gives a poorly converged density, which contaminates every eigenvalue on the path. Converge the grid on total energy first.
• Forgetting to set nbnd: without extra bands you cannot see conduction states or locate the CBM.
• Wrong high-symmetry coordinates: the correct k-points depend on the Bravais lattice. Tools like SeeK-path or the QE kpoints.x utility generate the standard path for you.
• Mismatched cutoffs: ecutwfc and ecutrho must be identical in the scf and bands runs, or the frozen density is inconsistent with the new Hamiltonian.
• Reading the Fermi level for metals: for metallic systems use smearing and read the Fermi energy is from the SCF output instead of the highest-occupied line.

For a broader view of how these DFT engines fit together, see our overview of the computational engines behind the platform, and the official Quantum ESPRESSO documentation for the complete input reference.

Run it on Simatra

Band structure workflows chain several pw.x and bands.x runs, and the dense SCF and long k-paths dominate the wall-clock time. Simatra runs Quantum ESPRESSO on GPU-accelerated clusters (instance GPU-Opt-V2) that deliver up to 5x faster convergence and handle supercells up to ~2,000 atoms, so you can iterate on cutoffs, functionals, and paths without waiting overnight. Start with a free trial and $100 in credits at app.simatra.io, and explore the full set of computational engines — including Quantum ESPRESSO and our native KRONOS engine — powering the platform.