Method & caveats

Each potential relaxes every elemental crystal to its own equilibrium (FrechetCellFilter + FIRE), then phonopy finite-displacement supercells (0.01 Å, ≥12 Å perpendicular spacing, q-mesh L = 100) give the harmonic spectrum: max frequency ω_max, an imaginary-mode flag (min ω < −0.1 THz), and the free energy / entropy / Cv. Quasi-harmonic thermodynamics come from an 11-volume Vinet-EOS fit (unstable volumes dropped), yielding 300 K thermal expansion, bulk modulus B, Cp and Gibbs energy.

There is no per-structure DFT phonon ground truth for most of the corpus, so the headline is cross-potential agreement: the median |ω_max − 58-potential consensus median| per potential (lower = more typical). This is the design's MLIP-vs-MLIP-spread reference — it ranks typicality, not correctness, so it is read alongside three other axes: the done-only imaginary-mode fraction (the 9 % ground-state vs 34 % metastable physics signal), the QHA bulk modulus vs DFT K_VRH (an order-of-magnitude anchor — 300 K thermally-softened B vs 0 K static K, hence the expected negative bias), and a Materials-Project DFT ω_max overlay for the ~half the corpus MP covers. All metrics use medians, robust to the campaign's single-structure force-constant and EOS blow-ups (the raw mean ω_max and mean QHA MAE are not). Full protocol, the as-run verification, and the MP-DFT provenance: Phonons & QHA methodology →.

Reproduce these phonons

Every harmonic row in the leaderboard is the output of one self-contained pipeline: relax the cell to the potential's own equilibrium, build phonopy finite-displacement supercells, evaluate the displaced forces with the MLIP, assemble the force constants, and read off the spectrum. The campaign packages exactly these steps as compute_harmonic() / compute_qha() in benchmark/phonon_compute.py (branch dais-port); the worker that drives them over the 834-structure queue is phonon_worker.py. The block below is that same recipe distilled to a single structure with a generic calculator — only one line is model-specific.

Dependencies

# the phonon engine — framework-agnostic, identical across every potential
pip install "phonopy==4.1.0" ase numpy

# + ONE machine-learning potential exposed as an ASE Calculator, e.g.
pip install mace-torch          # from mace.calculators import mace_mp
#   sevenn        -> from sevenn.calculator import SevenNetCalculator
#   fairchem-core -> OCPCalculator           (eSEN / eqV2 / UMA)
#   orb-models    -> from orb_models.forcefield import pretrained
#   tensorpotential / grace, mattersim, nequix (JAX), chgnet, matgl …
# exact per-family pins live in the env-build scripts (ligerzero-ai/mlip-env-setup).

A CUDA GPU is optional — pass device="cpu" to any calculator to run on CPU. The Postgres claim-queue and psycopg are campaign infrastructure for fanning the work over a GPU fleet; none of it is needed to compute a single structure.

Harmonic phonons (the eight steps, end to end)

import numpy as np
from ase import Atoms
from ase.build import bulk
from ase.filters import FrechetCellFilter
from ase.optimize import FIRE
from phonopy import Phonopy
from phonopy.structure.atoms import PhonopyAtoms
from phonopy.harmonic.dynmat_to_fc import get_commensurate_points

# 0. Pick a structure and attach ANY MLIP as an ASE calculator — the only model-specific line.
atoms = bulk("Cu", "fcc", a=3.6)          # the campaign fed in the 834 Materials-Project elemental cells
calc  = ...                               # e.g. mace_mp(model="medium-mpa-0", device="cuda", default_dtype="float32")

# 1. Relax cell + ions to the MLIP's OWN equilibrium V0  (FIRE + FrechetCellFilter, fmax = 1e-4 eV/Å).
atoms.calc = calc
converged = FIRE(FrechetCellFilter(atoms), logfile=None).run(fmax=1e-4, steps=500)
# status = "done" if converged else "unconverged" — the FC2 is still built + usable on the near-minimum cell.

# 2. Diagonal supercell matrix so every PERPENDICULAR width >= 12 Å (use perp spacing, never cell.lengths()).
cell   = np.asarray(atoms.get_cell())
vol    = abs(np.linalg.det(cell))
widths = [vol / np.linalg.norm(np.cross(cell[(i + 1) % 3], cell[(i + 2) % 3])) for i in range(3)]
scm    = np.diag(np.maximum(1, np.ceil(12.0 / np.array(widths))).astype(int))

# 3. Finite-displacement supercells: 0.01 Å, symmetric (+/-) displacements.
ph = Phonopy(
    PhonopyAtoms(symbols=atoms.get_chemical_symbols(), cell=atoms.get_cell().array,
                 masses=atoms.get_masses(),              # pass ASE masses — phonopy's table lacks actinides (Pu)
                 scaled_positions=atoms.get_scaled_positions()),
    supercell_matrix=scm, primitive_matrix="auto")
ph.generate_displacements(distance=0.01, is_plusminus="auto")

# 4. Single-point forces on each displaced supercell with the SAME MLIP (no relaxation of displaced cells).
forces = []
for sc in ph.supercells_with_displacements:
    a = Atoms(symbols=sc.symbols, cell=sc.cell, scaled_positions=sc.scaled_positions, pbc=True)
    a.calc = calc
    forces.append(a.get_forces())
ph.forces = forces

# 5. Force constants + acoustic sum rule (symmetrise away residual-force drift).
ph.produce_force_constants()
ph.symmetrize_force_constants()

# 6. Harmonic observables: length-based q-mesh L = 100 -> omega_max, DOS, and F/S/Cv(T); ZPE = F(0 K).
ph.run_mesh(100.0)
omega_max_thz = float(np.asarray(ph.get_mesh_dict()["frequencies"]).max())
ph.run_total_dos()
ph.run_thermal_properties(t_min=0, t_max=1000, t_step=10)

# 7. Imaginary-mode flag on the EXACT supercell-commensurate q-points (the 3 Gamma acoustic modes exempt).
qpts  = get_commensurate_points(scm)
ph.run_qpoints(qpts)
freqs = ph.get_qpoints_dict()["frequencies"]             # (n_q, n_band), THz
mins  = [(np.sort(freqs[i])[3:] if np.linalg.norm(q) < 1e-8 else np.sort(freqs[i])).min()
         for i, q in enumerate(qpts)]
has_imaginary = min(mins) < -0.1                         # min omega < -0.1 THz == dynamically unstable

# 8. (optional) Persist the FC2 — recompute any band structure / DOS / thermodynamics later with NO MLIP rerun.
ph.save("Cu.phonopy.yaml", settings={"force_constants": True})

Quasi-harmonic thermodynamics (QHA)

The QHA branch repeats steps 1–6 across an 11-point linear strain grid (ε ∈ np.linspace(-0.04, 0.06, 11)) around V₀. Each strained cell gets a constant-volume shape + ion relax (FrechetCellFilter(constant_volume=True), fmax = 1e-3) before its force constants — a free relax would collapse every cell back to V₀ and flatten the equation of state. Dynamically unstable volumes (any imaginary mode) are dropped, then phonopy.qha.QHA fits a Vinet EOS over the survivors to yield V(T), the thermal-expansion coefficient α(T), the bulk modulus B(T), Cp, and the Gibbs energy G(P, T) on a {0, 2, 5, 10} GPa pressure grid. The campaign wraps this as compute_qha(atoms, calc).

The locked parameters ("as I did here")

QHA fit quality (done vs partial). A full EOS uses all 11 strained volumes; volumes whose phonons go dynamically unstable are dropped, so many potentials yield a partial fit from 4–10 volumes. Both done and partial produce valid temperature-dependent curves and are shown — only failed fits (≤ 3 usable volumes) are excluded. The Min EOS volumes control trades coverage for confidence: the default (4) shows every usable fit, while 11 keeps only full-EOS done fits. This matters most for soft/heavy elements (Au, Re, Te), where few potentials converge all 11 volumes. There is no per-structure DFT QHA ground truth; the only anchor is the static 0 K DFT K_VRH vs the 300 K thermal B(T) — an order-of-magnitude check that expects a small negative bias from thermal softening.

Full protocol, the as-run verification, and the MP-DFT provenance: Phonons & QHA methodology → · Thermodynamics methodology → · as-run compute pipeline (deep) →.