Band structure and density of states

Band structure is to periodic systems what the MO diagram is to molecules — the ordered energy levels that tell you whether a crystal is a metal, semiconductor, or insulator, and what the gap is when it’s not metallic. vibe-qc samples the one-electron Hamiltonian on a k-path (for bands) or a k-mesh (for DOS), plus a matching matplotlib plotter for publication-style figures.

This tutorial uses the same 1D H2 molecular crystal from the periodic HF tutorial — two bands (bonding + antibonding), clear band-gap picture, tiny compute cost.

Minimum example

import vibeqc as vq

# 1D H2 chain: unit cell = one H2 molecule (bond 1.4 bohr),
# lattice vector a = 6 bohr along x.
system = vq.PeriodicSystem(
    dim=1,
    lattice=[[6.0, 0.0, 0.0], [0.0, 30.0, 0.0], [0.0, 0.0, 30.0]],
    unit_cell=[vq.Atom(1, [0.0, 0.0, 0.0]),
               vq.Atom(1, [1.4, 0.0, 0.0])],
)
basis = vq.BasisSet(system.unit_cell_molecule(), "sto-3g")

# Γ → X along the chain axis: fractional k from (0,0,0) to (½,0,0).
kpath = vq.kpath_from_segments(
    system,
    [((0.0, 0.0, 0.0), "Γ", (0.5, 0.0, 0.0), "X")],
    points_per_segment=40,
)

bands = vq.band_structure_hcore(system, basis, kpath, n_electrons_per_cell=2)
dos   = vq.density_of_states_hcore(
    system, basis, mesh=[80, 1, 1],
    sigma=0.02, n_electrons_per_cell=2,
)

You now hold two convenience objects:

  • bands.energies — (n_kpoints, n_bands) numpy array (Hartree).

  • bands.e_fermi — the HOMO level (chemical potential at T=0).

  • dos.energies / dos.dos — Gaussian-broadened DOS on an auto-chosen energy grid.

At the edges:

Γ:  bonding = -2.6660 Ha,   antibonding = -1.8335 Ha
X:  bonding = -2.6542 Ha,   antibonding = -1.9638 Ha
Fermi level: -2.6542 Ha   (filled bonding band, empty antibonding)

The plot

from vibeqc.plot import bands_dos_figure

fig = bands_dos_figure(bands, dos, title="H2 chain (STO-3G, Hcore)")
fig.savefig("h_chain_bands.png", dpi=150)

You get a two-panel figure: bands on the left with high-symmetry labels (Γ, X) along the bottom, DOS on the right sharing the same energy axis, Fermi level marked in both panels as a horizontal line.

What the picture teaches

  • Two bands, as expected. STO-3G gives each H one s orbital; one H2 per cell = two AOs = two bands.

  • Bonding band is almost flat. Dispersion along Γ → X is only ~12 mHa = 0.3 eV, because neighbouring H2 molecules (3.2 bohr apart) barely overlap their bonding orbitals.

  • Antibonding band disperses more (~130 mHa = 3.5 eV). The anti-bonding orbital has a node between H atoms of one molecule and constructive overlap with the nearest H of the next molecule, so its energy is more sensitive to the k-point.

  • There’s a gap everywhere — this is an insulator. Shrink the lattice parameter (a) and the antibonding band eventually crosses down below the bonding at Γ; at that point the system is a semimetal and HF refuses to converge, consistent with the Peierls-instability story from examples/input-h-chain-peierls.py.

A multi-segment k-path

Real 3D Brillouin zones have many high-symmetry points. Stitch segments together:

kpath = vq.kpath_from_segments(
    system,
    [
        ((0.0, 0.0, 0.0), "Γ", (0.5, 0.0, 0.0), "X"),
        ((0.5, 0.0, 0.0), "X", (0.5, 0.5, 0.0), "M"),
        ((0.5, 0.5, 0.0), "M", (0.0, 0.0, 0.0), "Γ"),
    ],
    points_per_segment=40,
)

The plotter adds labels and ticks at each segment boundary.

Custom plotting with matplotlib

If you want to overlay multiple runs or style the plot for a paper, grab the raw arrays and use matplotlib directly:

import matplotlib.pyplot as plt

fig, (ax_b, ax_d) = plt.subplots(1, 2, sharey=True, figsize=(8, 4),
                                  gridspec_kw={"width_ratios": [3, 1]})

# Bands
x = bands.kpath.distances
for i in range(bands.energies.shape[1]):
    ax_b.plot(x, bands.energies[:, i] - bands.e_fermi,
              color="C0", lw=1.2)
ax_b.axhline(0.0, color="k", lw=0.5, ls="--")
# kpath.labels is [(distance, name), ...] at each high-symmetry point.
tick_pos, tick_txt = zip(*bands.kpath.labels)
ax_b.set_xticks(tick_pos)
ax_b.set_xticklabels(tick_txt)
ax_b.set_xlim(x[0], x[-1])
ax_b.set_ylabel("E - E_F (Ha)")

# DOS
ax_d.fill_betweenx(dos.energies - bands.e_fermi, dos.dos, color="C0", alpha=0.5)
ax_d.axhline(0.0, color="k", lw=0.5, ls="--")
ax_d.set_xlabel("DOS (states / Ha)")

fig.tight_layout()
fig.savefig("bands_custom.png", dpi=200)

A note on “Hcore” bands

band_structure_hcore / density_of_states_hcore diagonalise the kinetic + nuclear-attraction part of the Fock matrix — the non-interacting electronic structure. This is:

  • Fast. No SCF required, just Bloch-summed one-electron integrals and a diagonalisation at each k-point.

  • Correct for tight-binding-like bands. The s-only H-chain has exactly the structure you’d get from a Hückel model; Hcore nails it.

  • Missing electron-electron interaction. For real bands on a real solid you want HF or KS-DFT self-consistent. That path is wired through run_rhf_periodic + band_structure_from_real_space, though the solver is still being hardened for tight-cell 2D/3D geometries.

For quick pedagogical plots Hcore is plenty.

Saving the DOS for external plotting

import numpy as np
np.savetxt("h_chain_dos.dat",
           np.column_stack([dos.energies, dos.dos]),
           header="E (Ha)    DOS (states/Ha)")

Open in gnuplot, Origin, or any plotting tool.

Periodic volumetric output

Crystal orbital and density visualisation uses XSF / BXSF format instead of cube — see volumetric data user guide for write_xsf_volume (isosurfaces in VESTA / XCrySDen) and write_bxsf (Fermi-surface plots).

Next

This is the first solid-state tutorial built on the matured periodic infrastructure; more will follow (full-SCF bands, phonons, surface cells) as the periodic SCF convergence tooling lands.