Deep, code-cited reference. See the readable summary for the public version. This page is the full provenance. The track shipped on 2026-06-27 (archive vacancy_archive_20260627, benchmark commit 8e047ba); §2 and §9 describe the live site, and the Shipped-reality note at the bottom reconciles the earlier 2026-06-08 verification addendum (which predates the run).

Monovacancy Formation Energy Benchmark — Methodology

Scientific and numerical specification of the vacancy-formation track in mlip-elemental-benchmark. Covers the atomistics (supercell construction, site enumeration, the relaxation protocol) and the optimiser specifics. Operational details (SLURM, claim queue, drivers) live in VACANCY.md; this file is the what and why of the physics.

Reference implementation: benchmark/vacancy_compute.py (compute_vacancies).

1. Quantity computed

For each elemental crystal and each symmetry-inequivalent atomic site, the constant-volume monovacancy formation energy:

E_f = E_vac − (N−1)/N · E_bulk

The (N−1)/N · E_bulk term is the chemical potential of the removed atom taken as the per-atom bulk energy at the same supercell and cell, giving a constant-volume formation energy referenced to each MLIP's own equilibrium lattice. This is the verified normalisation from the pyiron_workflow_atomistics vacancy node; we compose the package primitives rather than calling its single-site macro (the macro does no bulk pre-relax and recomputes the supercell).

One structure yields multiple rows — one E_f per inequivalent site (Han's per-site refinement). DB key: (structure_id, potential, site_index).

2. Structure corpus

• 834 MP elemental structures across 89 elements (all polymorphs across the elemental phase space — matches the elastic track's scope). The campaign ran, and the site ships, the full 834; this is not a ground-states-only run.
• Each is a Materials Project conventional/primitive cell stored as a pickled ase.Atoms blob in the structure table.
• The is_ground_state flag drives a ground-state facet on the leaderboard (toggle: Ground-state only / All polymorphs / Metastable), which is the default view because the high-energy polymorphs stress the relaxers (and the consensus) hardest. The ground-state subset is 88 structures / 166 inequivalent sites / 88 elements.
• Archive snapshot (vacancy_archive_20260627): 233,388 done, 24,142 unconverged, 1,821 error rows in mlip_vacancy, across 60 potentials. 797 structures have a result for all 60 potentials (the clean rectangular set); 37 are ragged (NaN for ≥1 potential) after the 2026-06-27 scope purge (see the shipped-reality note).

3. Supercell construction

3.1 Minimum-image criterion — perpendicular inter-plane spacing

To keep a vacancy from interacting with its periodic image, each axis is repeated until its true perpendicular inter-plane spacing reaches the minimum-image target:

d_i = V / |a_j × a_k|        (i,j,k cyclic),   V = |det(cell)|
n_i = ceil(min_image / d_i)  (≥ 1)

Minimum-image target: min_image = 12.0 Å.

This uses the inter-plane spacing d_i, not the cell-vector length |a_i|. For skewed or layered cells the vector length badly overestimates isolation — e.g. a flat carbon polymorph can have |c| ≥ 12 Å yet only ~2 Å between a vacancy and its image; the naïve |a_i| criterion under-expands such cells. (_repeats in vacancy_compute.py.)

Consequence: median supercell ≈ 150 atoms; maximum ≈ 3000 atoms (a couple of low-density carbon polymorphs at 6×6×6).

3.2 Determinism across potentials

The integer repeats (n_x, n_y, n_z) and the inequivalent-site representative indices are derived from the original MP cell, i.e. they are potential-independent. The same fixed repeats are applied to each MLIP's relaxed bulk. Because the supercell uses ASE repeat ordering with the image-0 block first, primitive-cell representative indices are valid removal indices in the supercell. Every potential therefore sees the identical supercell size and the identical set of removed sites for a given structure — the analog of a shared deformation set in the elastic track.

4. Inequivalent-site enumeration

Sites are potential-independent, so they are enumerated once for all 834 structures into benchmark/site_map.json (precompute_sites.py); workers load it read-only via sites_override, so no symmetry/descriptor code runs in the GPU hot loop and the removed- site set is byte-identical on every potential.

4.1 Authoritative backend — SOAP

• Descriptor: SOAP (dscribe), computed on the supercell.
• Cutoff: r_cut = 8 Å.
• Post-processing: PCA-whiten the SOAP vectors, then group sites by cosine similarity ≥ 0.999; one representative per group = one inequivalent site.
• Run in a separate envs_soap environment (single dscribe import for the whole corpus, BLAS pinned to one thread per worker).
• r_cut = 8 (not 6) is threshold-stable across cosine 0.99–0.999 for all structures; r_cut = 6 over-splits and 0.9999 introduces wobble on a few elements.

4.2 Cross-check — spacegroup-locked symmetry

A spglib enumeration on the primitive cell, with the symprec ladder locked to the MP structure.spacegroup, is recorded as an audit cross-check (sym_n, agree_sym) and supplies the Wyckoff labels carried on each row. SOAP and symmetry disagree on a small minority of soft/low-symmetry crystals where SOAP splits one Wyckoff orbit into near- identical sub-groups; their per-site E_f come out near-identical, confirming the split is benign. SOAP is authoritative; symmetry is the labelled cross-check.

The live (non-precomputed) fallback inside compute_vacancies is site_method="symmetry", used only if no override is supplied.

5. Relaxation protocol (constant-volume, three stages)

All three stages use the same calculator (the MLIP under test) and the same convergence target. Composed from pyiron_workflow_atomistics ASEEngine + CalcInputMinimize.

Rationale for the split: the bulk pre-relax finds each potential's equilibrium cell so the formation energy is referenced to that potential's lattice (decision A); the vacancy and perfect-supercell relaxations are then done at fixed cell so the formation energy is genuinely constant-volume (the vacancy is not allowed to collapse the cell). The supercell is rebuilt from the relaxed bulk with the same fixed repeats, preserving atom indexing.

6. Optimiser specifics

The optimiser is matched to the relaxation type (_relax in vacancy_compute.py):

optimizer_class = BFGS if relax_cell else FIRE

6.1 Shared optimiser settings

6.2 Engine wiring note

ASEEngine.get_calculate_fn routes optimizer_class and max_steps into the ASE minimisation for CalcInputMinimize inputs. CalcInputStatic would force optimizer_class=None (single-point), but all three vacancy stages are minimisations, so the matched BFGS/FIRE choice applies throughout.

7. Convergence semantics

A site's E_f is trustworthy only if all three feeding relaxations reached fmax. The row carries per-stage flags and a combined converged:

converged = bulk_converged AND supercell_converged AND vacancy_converged

Row status:

Stored columns include bulk_converged, supercell_converged, converged, n_atoms_supercell, supercell_nx/ny/nz, min_image_angstrom, relaxed_*_volume, relaxed_bulk_energy_per_atom, fmax, plus the per-site formation_energy, vacancy_energy, wyckoff_symbol, site_multiplicity, site_element.

Convergence is a real benchmark signal: a potential with a rough force/stress landscape (e.g. one that cannot relax variable cells) shows a low converged fraction — that is reported honestly, not forced by inflating the step cap.

8. Error handling

Per-structure failures are caught and recorded (status='error', full traceback in error_msg); the worker continues to the next structure. Two classes are distinguished:

• Recoverable / per-structure (e.g. plain CUDA out-of-memory on a huge cell): flagged error and skipped; the worker keeps going. Accepted as graceful degradation.
• Sticky CUDA faults (illegal memory access, device-side assert, CUBLAS failure): these corrupt the process's CUDA context, so every subsequent kernel call fails identically — a cascade of bogus errors. The worker detects these (_is_cuda_poison), writes the single real error row, and exits so the SLURM resubmit chain respawns a worker with a fresh CUDA context. One real fault ⇒ one error row, not a cascade. Cascaded rows are recovered by re-queuing those structures and deleting their error rows, then rerunning the potential.

Plain "out of memory" is deliberately not treated as a sticky fault (it is usually recoverable and per-structure).

9. Reference and aggregation

No DFT or MP reference column. The reference is the cross-potential consensus of E_f, overlaid where available with a curated experimental value. The site rollup (tools/regen_vacancy_from_export.py, run off the archive parquet, not the DB-side build_vacancy_rollup.py) groups done rows by (structure_id, site_index) — a stable cross-potential site key, since the removed-site set is potential-independent — and writes four vendored CSVs under src/assets/elemental/vacancy/:

unconverged and error rows are excluded from the consensus and the leaderboard. conv_pct = done / (done + unconverged) (error rows are not in the denominator).

9.1 Leaderboard metrics

Registered in src/components/elemental-leaderboard-agg.js (the vacancy track). All are median / median-of-absolute, so they are robust to the heavy tail (see §9.3):

vacancy_consensus_mae is the implemented MAE-vs-median: the cross-potential spread is the reference model, so it ranks "most typical", not "most correct" — read it alongside the vs PBE-DFT column (dft_dev_ev = E_f − the published PBE value from Medasani et al. 2015, 24 phase-matched metals; the closest like-for-like comparand for PBE-trained MLIPs) and the experiment column. The track's default facet is ground states. conv_pct, runtime and GPU are appended as per-potential lookup columns from vacancy_potentials.csv (not metrics).

9.2 Experimental literature overlay

src/assets/elemental/vacancy/vacancy_literature.csv curates 20 common elemental metals from the Korhonen–Puska–Nieminen compilation (Phys. Rev. B 51, 9526 (1995)) plus primary positron-annihilation / dilatometry papers (Fluss–Smedskjaer, Triftshäuser, Maier–Seeger, Feder–Charbnau, …). regen_vacancy_from_export.py::load_literature keys it by element (element,e_f_ev,spacegroup, source[,note]) and phase-gates the attachment: a value is attached to an element's ground-state site only when the corpus 0 K ground state IS the experimentally-measured phase (spacegroup match). For 5 metals (Ag, Mg, In, Na, Ti) the lowest-hull DFT polymorph differs from the measured phase, so they are listed but not attached, leaving 15 phase-matched, single-site overlay elements (Al, Cu, Au, Ni, Pt, Pd, Fe, W, Mo, Ta, Nb, Cr, V, Pb, Co). This avoids a wrong-phase reference injecting spurious per-element "vs experiment" deviations (e.g. Ti's measured-hcp 1.55 eV vs the corpus sg191 polymorph), and incidentally removes every multi-site overlay case. Finding: the consensus systematically under-predicts noble-metal and Ni vacancies vs experiment (Pt −0.65, Pd −0.51, Au −0.46, Cu −0.20, Ni −0.30 eV) while matching Al / Fe / Mo / Nb / Ta within ~0.1 eV; a few early-3d metals (Cr +0.84, Co +0.43) instead over-predict. The full per-element value/phase/source table and the bibliography (with links) are on the concise methodology page.

9.3 Robustness — why everything is a median

~4 % of raw E_f are non-physical blow-ups (up to ~1e25 eV) from potentials that "converge" to a garbage PES minimum — which is exactly why every metric is median / median-of-absolute rather than a mean. The ground-state-default view is clean (median E_f range −1.76 … 6.92 eV; only the Br₂ / I₂ molecular crystals have an unphysical negative consensus median, with a few open non-metal sites — P, B — carrying comparably large MAD). For ~200 hard metastable sites >50 % of potentials blow up, so even the per-site median is unreliable — hence the ground-state default and the "most typical, not most correct" caveat.

10. Parameter summary

Implementation: benchmark/vacancy_compute.py · site map: precompute_sites.py → benchmark/site_map.json · worker / error handling: benchmark/vacancy_worker.py · DB rollup: benchmark/eval/build_vacancy_rollup.py · site rollup (shipped): tools/regen_vacancy_from_export.py → the four src/assets/elemental/vacancy/*.csv · leaderboard metrics: src/components/elemental-leaderboard-agg.js · operations: VACANCY.md.

Shipped-reality note (campaign closed 2026-06-27, archive vacancy_archive_20260627)

The 2026-06-08 verification addendum below predates the production run; two of its points are now superseded by ship:

• Corpus is the full 834 (not ground-states-only). The campaign ran all 834 structures across 89 elements; the site ships them with a ground-state-default facet. The ground-state subset is 88 structures / 166 inequivalent sites / 88 elements (the earlier "90 / 164" was the pre-run plan). Snapshot: 233,388 done / 24,142 unconverged / 1,821 error across 60 potentials. The campaign was closed by purging the pending stress-test polymorphs (P1 cells, ~100 sites × ~800-atom supercells) on 2026-06-27; existing results kept, the rest left NaN — 797 structures complete across all 60 potentials, 37 ragged.
• mae_vs_median IS implemented — as the headline vacancy_consensus_mae (medianAbsSkipNull over deviation_from_median) in elemental-leaderboard-agg.js, alongside vacancy_consensus_mbe and the experimental vacancy_lit_mae. The shipped site rollup is tools/regen_vacancy_from_export.py (keys the site consensus by (structure_id, site_index) and the literature overlay by element), not the DB-side build_vacancy_rollup.py.

Verification addendum (run wf_7fbcd118-848, 2026-06-08, pre-run)

• Production corpus is ground-states-only: 90 structures / 164 inequivalent sites — superseded above (the run is the full 834; the GS facet is 88 / 166).
• mae_vs_median is an unimplemented planned metric — superseded above (shipped as vacancy_consensus_mae). The rollup's n_converged counts converged potentials at a site (sp["n"]), not "sites a potential converged" — still true.
• error_msg stores only the last 4000 characters of the traceback (vacancy_worker.py:221).
• Sticky-CUDA poison set is five markers (vacancy_worker.py:141-147): "illegal memory access", "device-side assert", "CUBLAS_STATUS", "CUDA error:", "an illegal instruction" — not three.