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	any	Γ-only	❌ (queued M5)	`pbc_bipole.py` and siblings
BIPOLE	any	multi-k	❌ (queued behind multi-k BIPOLE chat)
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 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.

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.