Fermi-Dirac smearing for metals and small-gap solids

The molecular-limit SCF assumes a clean separation between the occupied and virtual manifolds — every orbital below the HOMO is filled, every orbital above the LUMO is empty. Metals don’t work that way. Bands cross the Fermi level; at finite k-mesh resolution different k-points integer-occupy slightly different band counts, and the Aufbau iteration oscillates between near- degenerate occupation patterns. The classic fix, due to Mermin 1965 in the DFT context, is to give every orbital a fractional Fermi-Dirac occupation

\[ f_i(\varepsilon_F, T) \;=\; \frac{1}{1 + \exp\bigl[(\varepsilon_i - \varepsilon_F)/k_B T\bigr]} \]

at a small electronic temperature \(T\) and integrate to the right total electron count by adjusting \(\varepsilon_F\) each iteration. This smears the band-crossing into a soft step that the SCF can descend without oscillation. The trade-off is an electronic-entropy contribution to the free energy that the SCF minimises in place of the bare total energy — physically harmless at small \(T\), and the relevant thermodynamic potential for any metallic system anyway.

vibe-qc v0.8.0 wires Fermi-Dirac smearing into the periodic multi-k Ewald RHF/RKS drivers and into the Γ-only native GDF RHF/RKS path. This tutorial walks through how to enable it, what the output looks like, and how to choose the temperature for a representative case (a 1D hydrogen chain at half-filling — the simplest metallic-band testbed in the bundled examples).

The recipe

PeriodicSCFOptions.smearing_temperature accepts four forms:

opts.smearing_temperature = 0.0                # default — no smearing
opts.smearing_temperature = 0.005              # explicit Hartree
opts.smearing_temperature = "0.005 Ha"         # explicit, with unit
opts.smearing_temperature = "1580 K"           # convert from K
opts.smearing_temperature = "0.136 eV"         # convert from eV
opts.smearing_temperature = "0.01 Ry"          # convert from Rydberg
opts.smearing_temperature = "metal"            # preset = 0.005 Ha
opts.smearing_temperature = "small-gap"        # preset = 0.002 Ha
opts.smearing_temperature = "auto"             # conservative auto-guess

The same field is exposed on run_periodic_job via smearing_temperature=. Internally everything is stored as electronic \(k_B T\) in Hartree (the SCF driver consumes that directly); the unit-aware parser is just for ergonomics.

The "auto" mode is deliberately conservative — without a metallic hint it returns zero smearing rather than silently changing the physics of a likely insulator. To opt into the metallic guess, pass metallic=True or n_electrons=<odd> to vq.resolve_smearing_temperature (or just use the "metal" preset directly).

Worked example: 1D hydrogen chain at half filling

The infinite 1D H chain at uniform bond length is the simplest metallic-band testbed in vibe-qc. With one electron per H atom and two H atoms per cell at uniform spacing, the band is exactly half-filled — Mott would dimerise it via Peierls (see tutorial 17), but at the RKS/PBE level with no symmetry breaking, you get a 1D metal.

Working directory:

~/vibeqc-runs/h-chain-metallic/
    input-h-chain-smear.py

input-h-chain-smear.py:

from pathlib import Path
import numpy as np
import vibeqc as vq

# Two H atoms per cell, uniformly spaced at 1.4 bohr (compressed
# regime — well inside the metallic phase, well away from the
# Peierls dimerisation that opens a gap at larger spacing).
a = 2.8  # bohr — cell length so atoms sit at 0 and a/2
sysp = vq.PeriodicSystem(
    dim=1,
    lattice=np.diag([a, 30.0, 30.0]),
    unit_cell=[vq.Atom(1, [0.0, 0.0, 0.0]),
               vq.Atom(1, [a / 2, 0.0, 0.0])],
)
basis = vq.BasisSet(sysp.unit_cell_molecule(), "pob-tzvp")
kmesh = vq.monkhorst_pack(sysp, [16, 1, 1])

# Vanilla RKS/PBE without smearing — expect oscillation /
# convergence failure on a 16-point mesh.
opts_no_smear = vq.PeriodicRKSOptions()
opts_no_smear.lattice_opts.cutoff_bohr = 12.0
opts_no_smear.max_iter = 80
opts_no_smear.functional = "PBE"

# Same input, "metal" preset for k_B T.
opts_smear = vq.PeriodicRKSOptions()
opts_smear.lattice_opts.cutoff_bohr = 12.0
opts_smear.max_iter = 80
opts_smear.functional = "PBE"
opts_smear.smearing_temperature = "metal"     # 0.005 Ha ≈ 0.136 eV

res_no = vq.run_rks_periodic(sysp, basis, kmesh, opts_no_smear)
res_sm = vq.run_rks_periodic(sysp, basis, kmesh, opts_smear)

print(f"No smearing : converged={res_no.converged}, "
      f"iters={res_no.n_iter}, E={res_no.energy:.8f}")
print(f"Smearing    : converged={res_sm.converged}, "
      f"iters={res_sm.n_iter}, E={res_sm.energy:.8f}")
print(f"  electronic entropy contribution: {res_sm.entropy_energy:.6f} Ha")
print(f"  Fermi level:                      {res_sm.fermi_level:.6f} Ha")

Expected behaviour:

  • Without smearing the SCF oscillates between integer-occupation patterns at the band crossing, often never converging within 80 iterations (or converging to a higher-energy stationary point that depends on the k-mesh).

  • With smearing the SCF settles in 12-20 iterations to the free-energy minimum. The reported entropy_energy is a few mHa, the fermi_level sits inside the bath of bands at the crossing.

What appears in the output

The .out log carries a smearing block right after the SCF trace:

Finite-temperature occupations (Fermi-Dirac)
----------------------------------------------------
  source            "metal" preset
  k_B T (Hartree)   5.000e-03      (1577.8 K, 0.1361 eV)
  E_free            -1.143217  Ha
  T * S             -0.001094  Ha    ← electronic entropy
  E_extrapolated    -1.144311  Ha    ← 0 K limit (E_free + ½ T S)
  E_internal        -1.145404  Ha    ← E_free - T*S, reported in `result.energy`
  Fermi level (ε_F) -0.137820  Ha

The result object carries the same fields:

res.energy           # E_internal — the usual SCF total
res.free_energy      # E_internal - T*S  (the *minimised* potential)
res.entropy_energy   # T * S  (always non-negative)
res.fermi_level      # ε_F at the converged occupations
res.occupations      # per-k-point Fermi-Dirac fractions, shape (nk, nbf)

For metallic systems the right number to report in a paper is E_extrapolated — the \(T \to 0\) extrapolation of the free energy, which is the same total energy you would have gotten from an infinitely fine k-mesh with no smearing. The other two energies are intermediates the SCF needed.

Choosing the temperature

The smearing temperature is a methodological knob, not a physical temperature — the system isn’t actually at 1580 K. Two rules of thumb:

  1. Pick the smallest \(k_B T\) that gives stable convergence. Too small and the SCF oscillates the same way it did without smearing; too large and the bandwidth gets smeared into meaninglessness. For most metals, 0.005 Ha (\(k_B T\) ≈ 0.136 eV ≈ 1580 K) is comfortably in the right regime — this is what "metal" resolves to.

  2. Verify the \(T \to 0\) extrapolation is well-behaved. Run the same calculation at 0.001, 0.002, 0.005, 0.01 Ha and check that E_extrapolated converges with \(k_B T \to 0\) faster than E_free or E_internal. If the extrapolation is unstable, you need a denser k-mesh, not a different temperature.

The "small-gap" preset (0.002 Ha) suits systems with a small but finite gap — small-bandgap semiconductors at large unit cells, or defective insulators where the defect state sits just below the conduction band. The "debug" preset (0.010 Ha) is intentionally on the high side for code-correctness tests; don’t use it for production.

The "auto" mode is conservative on purpose: without a metallic hint or recorded band gap it returns zero smearing. To opt into the metallic guess automatically pass metallic=True or band_gap_hartree=0.0:

from vibeqc import resolve_smearing_temperature
guess = resolve_smearing_temperature("auto", metallic=True)
print(guess.temperature, guess.reason)
# 0.005, 'metallic=True'

Compatibility matrix

Driver

Smearing supported?

run_rks_periodic_gamma_gdf

✓ Γ-only RHF/RKS

run_rhf_periodic_multi_k_ewald3d

✓ multi-k RHF

run_rks_periodic_multi_k_ewald3d

✓ multi-k RKS

run_uhf_periodic_* / run_uks_periodic_*

roadmap (post-v0.8.0)

Molecular run_rhf / run_rks / run_uhf / run_uks

not yet

Molecular smearing is rarely needed (open-shell molecules use UHF with explicit α/β occupations; small-gap molecular systems use the stiff-convergence recipe). Multi-k UHF / UKS smearing arrives once the multi-k open-shell drivers themselves land — see multi_k_scf.

Resources

The 1D H chain example above with a 16-point k-mesh converges in ~5 s on a modern laptop. Wider chains (32-64 k-points) or 2D graphene at a [12, 12, 1] mesh push the wall-time into the 30 s - 2 min range.

References

  • Mermin’s variational principle — N. D. Mermin, “Thermal Properties of the Inhomogeneous Electron Gas,” Phys. Rev. 137, A1441 (1965). The foundational statement that finite-temperature DFT minimises the Mermin free energy \(F = E - T S\), not the bare total energy.

  • Fermi-Dirac smearing in plane-wave DFT — M. Methfessel, A. T. Paxton, “High-precision sampling for Brillouin-zone integration in metals,” Phys. Rev. B 40, 3616 (1989). Compares Fermi-Dirac with Methfessel-Paxton smearing; the half-T-S extrapolation rule comes from this paper.

  • Smearing in Gaussian-basis periodic codes — R. Dovesi et al., “CRYSTAL14: A program for the ab initio investigation of crystalline solids,” Int. J. Quantum Chem. 114, 1287 (2014). The CRYSTAL SMEAR directive — the lineage vibe-qc’s periodic smearing follows.

Next