Smearing (fractional occupations for periodic SCF)¶

TL;DR. For a metallic periodic system (Fermi level cuts through bands), turn smearing on, set opts.smearing_temperature = 0.005 (Hartree, ~1580 K), or pass smearing_temperature="metal" to vq.run_periodic_job(...). For a clear-gap insulator, leave smearing off, the default of 0.0 is correct. Smearing is not a remedy for an oscillating or diverging SCF, if your SCF runs cleanly with smearing on but diverges with it off, you are masking a backend bug, not solving a convergence problem (see § When smearing is the wrong knob).

Smearing replaces the integer step occupation \(n_i \in \{0, 2\}\) at the Fermi level with a smooth distribution that lets fractional occupations relax around \(\mu\). Mathematically this regularises the Brillouin-zone integral \(\sum_{k\in\text{BZ}} \theta(\mu - \varepsilon_i(k))\) that is ill- defined when any band crosses \(\mu\), i.e., for every metal. The trade is a finite electronic temperature \(T\): the SCF converges to Mermin’s free energy \(A = E - TS\) (Mermin 1965) instead of the strict ground-state energy \(E_0\).

For metals the smearing energy \(TS\) is necessary to make the band integral well-defined at a finite k-mesh; the price is that you must either pick \(T\) small enough that \(TS\) is negligible at your accuracy target, or extrapolate \(T \to 0\) from a ladder of runs. The recommended workflow per flavor is in § Choosing T.

Quick start¶

High-level wrapper¶

import vibeqc as vq

# "metal" preset → k_B T = 0.005 Ha (~0.136 eV, ~1580 K).
vq.run_periodic_job(
    system, basis,
    method="rks", functional="pbe",
    smearing_temperature="metal",
    kpoints=(4, 4, 4),
)

Direct driver call¶

opts = vq.PeriodicSCFOptions()
opts.lattice_opts.coulomb_method = vq.CoulombMethod.EWALD_3D
opts.smearing_temperature = 0.005   # k_B T in Hartree

result = vq.run_rhf_periodic_multi_k_ewald3d(
    system, basis, kmesh=(4, 4, 4), options=opts,
)

print(f"Energy        = {result.energy:.8f} Ha")
print(f"Free energy   = {result.free_energy:.8f} Ha   (Mermin A = E - TS)")
print(f"Entropy S/k_B = {result.entropy:.6f}")
print(f"Fermi level μ = {result.fermi_level:.6f} Ha")
print(f"Occupations   = {result.occupations}")  # list[ndarray], per k-point

The SCF converges Mermin’s free energy \(A\), not \(E\). Both are on the result struct, use \(A\) for variational comparisons (e.g., forces, lattice optimisations under fixed \(T\)); use \(E\) for “what would the strict ground-state energy be at this geometry” only after T → 0 extrapolation.

Input forms¶

smearing_temperature on the high-level wrapper (vq.run_periodic_job) accepts any of:

Form	Example	Notes
Numeric width in Hartree (canonical)	`0.005`	Same units as the on-options-struct field
Named preset (string)	`"metal"`, `"small-gap"`, `"insulator"`, `"debug"`, `"off"`	See `vq.SMEARING_PRESETS`
Numeric string with unit	`"1000 K"`, `"0.1 eV"`, `"0.02 Ha"`, `"0.2 Ry"`	Resolved via `vq.resolve_smearing_temperature`
`"auto"` + hints	`smearing_temperature="auto", smearing_metallic=True`	Cautious, defaults to 0.0 if no hint
`None` / `"none"` / `"off"`	`None`	Disables smearing

The on-options-struct field PeriodicRHFOptions.smearing_temperature (and siblings on PeriodicSCFOptions / PeriodicKSOptions) is always the canonical electronic temperature \(k_B T\) in Hartree, no implicit unit conversion happens there.

run_periodic_job also accepts the v0.10.x migration surface:

vq.run_periodic_job(
    system, basis,
    method="rks", functional="pbe",
    smearing=vq.SmearingOptions.from_user("metal"),
)

Common conversions¶

Macroscopic value	\(k_B T\) in Hartree
100 K	\(3.17 \times 10^{-4}\)
300 K (room temperature)	\(9.50 \times 10^{-4}\)
1000 K	\(3.17 \times 10^{-3}\)
1580 K (the `"metal"` preset)	\(5.00 \times 10^{-3}\)
3000 K	\(9.50 \times 10^{-3}\)
0.1 eV	\(3.67 \times 10^{-3}\)
0.01 Ry	\(5.00 \times 10^{-3}\)

Helpers:

vq.kelvin_to_hartree_temperature(300.0)         # → 9.5e-4
vq.electronvolt_to_hartree_temperature(0.1)     # → 3.67e-3
vq.rydberg_to_hartree_temperature(0.01)         # → 5.0e-3
vq.hartree_to_kelvin_temperature(0.005)         # → 1579.3

Choosing T¶

The Mermin free-energy formalism is exact in \(T\), but only the \(T \to 0\) limit gives the strict ground-state energy. Finite \(T\) adds a \(TS\) shift that for Fermi-Dirac smearing decays only linearly in \(T\) as \(T \to 0\) (Marzari et al. 1999), so the strategy is one of:

Pick \(T\) small enough that the \(TS\) correction is below your accuracy target. For ~mHa/atom accuracy on typical 3d / 4d transition metals at room temperature, the "metal" preset (\(k_B T = 5 \times 10^{-3}\) Ha) sits at the upper edge of the “small” range. Tighter targets need \(T \sim 10^{-3}\) Ha (~300 K), which slows convergence but keeps the entropy correction near µHa.
Extrapolate \(T \to 0\) from a ladder of runs. Two to three \(T\) values plus a polynomial fit in the asymptotic regime gives a defensible ground-state estimate. The order of the extrapolation depends on flavor (linear for Fermi-Dirac, quadratic for Methfessel-Paxton, cubic for Marzari-Vanderbilt cold). A user-facing helper (vq.extrapolate_t_zero) is queued for a later milestone, until then, fit manually.

Recommended starting points¶

System type	Suggested `smearing_temperature`	Why
Clear-gap insulator (gap ≳ 1 eV)	`0.0` (default)	Integer Aufbau is exact; smearing only adds noise
Small-gap semiconductor (gap ~0.1 eV)	`"small-gap"` (\(2 \times 10^{-3}\) Ha)	Avoid integer-occupation jitter when the gap is comparable to k-mesh resolution
Metal	`"metal"` (\(5 \times 10^{-3}\) Ha)	Standard for converged metallic SCF
Convergence debugging	`"debug"` (\(10^{-2}\) Ha)	Useful to confirm a hard SCF can converge at high T; not for production runs (see § When smearing is the wrong knob)

Output surface¶

When smearing_temperature > 0, the SCF result struct populates:

Field	Meaning
`result.energy`	Strict electronic energy \(E\) at the converged density (no \(TS\) subtraction)
`result.free_energy`	Mermin free energy \(A = E - TS\), the variational quantity under smearing
`result.entropy`	Dimensionless electronic entropy \(S/k_B\) per unit cell
`result.fermi_level`	Chemical potential \(\mu\) in Hartree
`result.occupations`	List of per-\(k\) occupation arrays; each entry in \([0, 2]\) for closed-shell
`result.smearing_temperature`	Echo of the temperature actually used

When smearing_temperature == 0, the smearing fields collapse: free_energy == energy, entropy == 0, occupations lists hard-Aufbau {0, 2} arrays, and fermi_level is the conventional midgap value.

Per-backend availability¶

Smearing is being rolled out backend-by-backend as the v0.10.x smearing program (see docs/design_smearing.md for the full design and milestone plan). The state as shipped on the current main:

Backend	Spin	k-mesh	Fermi-Dirac smearing	Notes
Ewald-3D	RHF / RKS	multi-k	✅	`run_rhf_periodic_multi_k_ewald3d`, `run_rks_periodic_multi_k_ewald3d`
Ewald-3D (Python Γ-dispatch)	RHF / RKS	Γ-only	✅ (inherits via `n_kpts=1`)	`run_rhf_periodic_gamma_scf`, `run_rks_periodic_gamma_scf`
Ewald-3D (C++-direct Γ)	RHF / RKS	Γ-only	❌ (queued M2)	`run_rhf_periodic_gamma_ewald3d`, `run_rks_periodic_gamma_ewald3d`
Ewald-3D	UHF / UKS	any	❌ (queued M3)	Driver raises on `smearing_temperature > 0` today
GDF	RHF / RKS	multi-k	✅	`run_krhf_periodic_gdf`, `run_krks_periodic_gdf`
GDF	RHF / RKS	Γ-only	✅	`run_rhf_periodic_gamma_gdf`
GDF	UHF / UKS	any	❌ (queued M4 after open-shell GDF lands)
BIPOLE	RKS / UHF / UKS	Γ-only + multi-k	✅	`run_pbc_bipole_rks`, `run_pbc_bipole_uhf`, `run_pbc_bipole_uks`; analytic KS gradients reject smeared results, so production forces use FD for forces
BIPOLE	RHF	any	❌	Closed-shell HF BIPOLE keeps integer occupations today
GAPW	-	-	n/a (driver not yet a thing, queued M8)

Methfessel-Paxton (M6) and Marzari-Vanderbilt cold (M7) plug into the same vibeqc.smearing.apply_smearing dispatcher and will land on every cell that has Fermi-Dirac after M2-M5.

When smearing is the wrong knob¶

Smearing fixes integer-occupation discontinuities at metallic Fermi surfaces. It does NOT fix:

a wrong Hartree term (e.g., a Madelung self-image leak, a gauge bug in the lattice sum);
a wrong exchange term (a wrong sign in the Coulomb integral, a divergent omega in the Ewald split);
a wrong nuclear-attraction term;
an over-tight DIIS subspace that’s amplifying noise.

If your periodic SCF oscillates with smearing off but converges cleanly with smearing on, the SCF is probably converging to a non-physical stationary point that happens to be smooth under fractional occupations, the underlying bug is masked, not fixed. File it against the backend chat (CLAUDE.md § 7).

The high-level run_periodic_job wrapper now emits a warning before it applies a positive smearing temperature to a system identified as an ionic or covalent insulator, or to a run with an explicitly supplied positive band gap. The automatic convergence strategy never enables smearing for those profiles. The warning is also written into the .out convergence-strategy block so a batch calculation cannot hide the choice.

The canonical cautionary tale: the v0.7.0 Madelung self-image leak, which over-bound H₂/STO-3G in a 30-bohr box by ~0.587 Ha (a Madelung-scale shift). Smearing made the symptom look like a convergence issue. The actual fix was in the gauge of the periodic Hartree term, once that landed (in Löwdin’s Compass), the SCF converged without any smearing tweaks.

Parameter-free alternative: the Gilat-Raubenheimer net¶

Smearing needs a width \(T\) that you converge / extrapolate to \(T \to 0\). The Gilat-Raubenheimer net (bz_integration="gilat" on the multi-k Ewald drivers) instead integrates the Brillouin zone analytically per microcell, the CRYSTAL SHRINK IS ISP second-net analogue, with no width to converge. It lives in vibeqc.bz_integration.

Insulators / gapped systems, GR runs self-consistently (RHF / RKS / UKS, full or symmetry-reduced meshes). For a gap its occupations are integer, so it reproduces the Aufbau result, the value is the parameter-free machinery, not a different energy.

Metals, post-SCF only. GR cannot drive a metallic SCF: the sharp \(T=0\) Fermi surface makes the occupation-versus-density map discontinuous, so the SCF oscillates and no density mixer (linear damping, DIIS, or Anderson) converges it. Drive the metal SCF with smearing, then evaluate the GR Fermi level / occupations / DOS on the converged eigenvalues, the standard tetrahedron/Gilat usage:

import numpy as np
from vibeqc.bz_integration import (
    gilat_occupations_on_kmesh, fractional_kpoints,
)

# 1. Converge the metal with smearing (a smooth map - robust).
opts.smearing_temperature = 0.005
r = vq.run_rks_periodic_multi_k_ewald3d(sysp, basis, kmesh, opts)

# 2. Parameter-free GR Fermi level + occupations on the converged spectrum.
#    (Full mesh; for an IBZ mesh expand first via expand_ibz_eigenvalues.)
frac = fractional_kpoints(sysp.reciprocal_lattice(),
                          np.asarray(kmesh.kpoints))
mesh = tuple(int(x) for x in kmesh.mesh)
occ, e_fermi = gilat_occupations_on_kmesh(
    frac, mesh, r.mo_energies, sysp.n_electrons())
# DOS: scatter the spectrum with eigenvalues_to_full_grid, then gilat_dos.

The metal total energy still comes from the smearing run (its \(T \to 0\) limit is well-defined); GR supplies the parameter-free spectral quantities. Reference: Gilat & Raubenheimer, Phys. Rev. 144, 390 (1966).

Theory references¶

Mermin, N. D. Thermal properties of the inhomogeneous electron gas. Phys. Rev. 137, A1441 (1965)., The free-energy formalism \(A = E - TS\) vibe-qc uses for all flavors.
Methfessel, M. & Paxton, A. T. High-precision sampling for Brillouin-zone integration in metals. Phys. Rev. B 40, 3616 (1989)., MP smearing (M6).
Marzari, N., Vanderbilt, D., De Vita, A. & Payne, M. C. Thermal contraction and disordering of the Al(110) surface. Phys. Rev. Lett. 82, 3296 (1999)., Cold smearing (M7) and the original \(T \to 0\) extrapolation analysis.
dos Santos, F. P. & Marzari, N. Fermi-surface effects and the electronic entropy of crystals at finite temperature. Phys. Rev. B 107, 195122 (2023)., Modern review of smearing-induced error and the practical \(T \to 0\) extrapolation ladder.