Universal Thermodynamic Interatomic Potentials for Crystalline Materials

Juno Nam; Bowen Deng; Xiaochen Du; Luis Barroso-Luque; Benjamin Kurt Miller; Rafael Gómez-Bombarelli

Source-linked AI summary

Universal Thermodynamic Interatomic Potentials for Crystalline Materials

Juno Nam, Bowen Deng, Xiaochen Du, Luis Barroso-Luque, Benjamin Kurt Miller, Rafael Gómez-Bombarelli

arXiv:2608.14502v1Aug 14, 2026cond-mat.mtrl-scicond-mat.stat-mechcs.AIcs.LGphysics.chem-ph

TL;DR

Free-energy calculations are needed to resolve finite-temperature phase stability but remain difficult because they require thermodynamic ensemble averages. This paper introduces TIP[UMA], a differentiable Gibbs free-energy model that achieves 10.3 meV atom^-1 median error and predicts several entropy-driven transitions within ∼250 K of experiment.

Problem
Finite-temperature phase stability requires precise free energies over temperature and pressure, but computing them requires costly ensemble averaging that includes vibrational and configurational contributions.
Method
TIP encodes a relaxed crystal as a differentiable Gibbs free-energy surface, trained across quasi-harmonic and molecular-dynamics fidelities and calibrated to higher-resolution calculations or experiment.
Results
10.3 meV atom^-1 global median absolute Gibbs energy error is achieved on held-out validation, while three of four tested entropy-driven transitions are predicted within ∼250 K of experiment.
Takeaways & Limitations
TIP makes crystal free energies and thermodynamic response functions available from a single evaluation, supporting phase-transition and alloy-solubility analyses without further sampling for in-domain crystals.
Takeaways & Limitations
The present model omits magnetic and electronic entropy, limiting its applicability to magnetic and mixed-valence solids where these contributions can qualitatively rearrange phase diagrams.

Abstract

from arXiv · show

Free energies govern solid-state phase stability, yet computational materials discovery still relies largely on ground-state energies because free energy calculations require ensemble averages. We introduce the thermodynamic interatomic potential (TIP), which extends an interatomic potential from its static energy to a thermodynamically consistent Gibbs free energy model, with thermodynamic responses following from temperature and pressure by automatic differentiation. We implement TIP[UMA] using the universal potential UMA, train it on free energies from quasi-harmonic to molecular dynamics fidelity, and calibrate it to higher-resolution calculations or experiment. From a single evaluation, it returns the equation of state of a crystal and locates phase transitions among competing branches, including dynamically stabilized phases. Fine-tuning extends the model to alloy solubility limits and miscibility gaps. TIP makes the free energy as accessible as the potential energy, opening finite-temperature phase stability to high-throughput discovery.

I. INTRODUCTION · II. RESULTS · A. Thermodynamic Interatomic Potentials

The paper introduces thermodynamic interatomic potentials (TIPs) to model Gibbs free-energy surfaces directly across temperature and pressure, addressing the limitations of ground-state energetics and repeated ensemble sampling. TIP[UMA] combines a universal interatomic potential with a thermodynamic adapter trained across increasing physical fidelity to predict phase stability and thermodynamic responses from relaxed crystal structures.

I. INTRODUCTION: Free energies, rather than ground-state energies alone, determine phase diagrams, stability fields, and metastability windows, requiring precision across relevant temperatures and pressures.Ground-state workflows remain prevalent because free-energy calculations require ensemble averages over thermally accessible configurations.
I. INTRODUCTION: Existing approaches address configurational or vibrational thermodynamics but become limited by parent-lattice assumptions, anharmonicity, or the cost of first-principles sampling.Universal MLIPs reduce individual energy and force costs but still require ensemble sampling for free energies.
I. INTRODUCTION: TIP combines MLIP chemical transferability, cluster-expansion-style branch thermodynamics, and CALPHAD-style direct free-energy modeling.The approach models an ordered branch’s Gibbs free energy directly rather than applying statistical-mechanical calculations to an energy model.
II. RESULTS: A single relaxed ordered crystal generates a differentiable Gibbs free-energy surface over temperature and pressure, eliminating further sampling for new structures within the model domain.TIP encodes the structure once, while statistical-mechanical costs are confined to training.
II. RESULTS: TIP compares candidate ordered branches to resolve stable phases and transitions, including dynamically stabilized phases, polymorphs, and distinct chemical orderings.The current model identifies transitions only among supplied branch candidates, does not discover crystal structures, and excludes liquids.
A. Thermodynamic Interatomic Potentials: TIP[UMA] uses UMA and trains in three stages of increasing physical fidelity, beginning with dense quasi-harmonic free-energy grids over temperature and pressure.Later stages incorporate higher-fidelity molecular-dynamics free energies and experimental calibration.
A. Thermodynamic Interatomic Potentials: The thermodynamic adapter predicts temperature- and pressure-independent coefficients from frozen MLIP structure features, while temperature and pressure enter through an analytic thermodynamic head.The free energy is decomposed into static 0 K energy plus a smooth thermodynamic residual, reducing the structure-dependent target range.
A. Thermodynamic Interatomic Potentials: Inference requires only a relaxed unit cell, reducing costs from hours of explicit sampling or minutes of quasi-harmonic workflows to seconds of model inference.Training labels may require large supercells, but inference cost is set by the unit cell.

B. Multi-Fidelity Gibbs Free Energy Dataset · C. Validation and Error Analysis · D. Ordered Crystalline Phases

TIP[UMA] is trained on a multi-fidelity Gibbs free energy dataset combining broad quasi-harmonic labels with accurate molecular-dynamics labels, achieving low validation error after anharmonic mid-training. From a single relaxed structure, it predicts thermodynamic response functions, calibrated phase properties, and temperature-driven transitions, while revealing limitations for sparse anharmonic and magnetic effects.

B. Multi-Fidelity Gibbs Free Energy Dataset: 123,424 metastable Materials Project structures within 0.2 eV/atom of the convex hull seed the dataset, combining abundant approximate and sparse accurate labels.QHA Gibbs free energies are computed for every structure, while absolute MD free energies are computed only for representatives.
C. Validation and Error Analysis: 10.3 meV atom−1 is TIP[UMA]’s global median absolute Gibbs energy error on the filtered held-out validation set.This error is comparable to the approximately 15 meV atom−1 median energy above the convex hull of experimentally observed inorganic phases.
C. Validation and Error Analysis: The largest element-resolved median errors occur for halogens, chalcogens, alkali and alkaline-earth metals, heavy p-block elements, and actinides.Most transition metals have errors at or below the global median.
C. Validation and Error Analysis: Mid-training on absolute MD labels lowers median Gibbs energy error about threefold at every temperature relative to quasi-harmonic pre-training.For the mid-trained model, error grows weakly with temperature and more strongly at high pressure.
D. Ordered Crystalline Phases: A single relaxed structure supplies the Gibbs free energy surface from which TIP[UMA] evaluates heat capacities, equations of state, polymorphic transitions, and phase diagrams.These equilibrium response functions follow from derivatives of the Gibbs free energy.
D. Ordered Crystalline Phases: 2.28 and 1.94 J K−1 (mol atoms)−1 are the median absolute errors for standard entropy and heat capacity across 255 NIST–JANAF stoichiometric solids.Predicted heat capacities for Al2O3 and CaO track tabulated values over the full temperature range and recover CP →0 as T →0.
D. Ordered Crystalline Phases: The pre-trained model finds no ∆G = 0 crossing for Ti or Sc and predicts Hf and CaSiO3 crossings near 3200 and 2200 K, versus experimental 2016 and 1398 K.These entropy-driven transitions require anharmonic treatment because high-temperature bcc metals are dynamically unstable at 0 K and CaSiO3 requires anharmonic entropy.
D. Ordered Crystalline Phases: Post-training reduces bulk modulus error against experiment from ∼13% to a few percent and removes the ∼5% PBE-related volume overestimate.Reaction Gibbs free energy and entropy errors also decrease, while held-out quaternary phases improve alongside fitted binary and ternary oxides.

E. Disordered Solid Solutions · F. Chemical Trends in Predicted Responses · III. DISCUSSION

TIP[UMA] extends thermodynamically consistent free-energy modeling from ordered crystals to disordered solid solutions, while its learned entropy and volume descriptors reflect established chemical trends. The discussion identifies missing magnetic and electronic entropy, reference-label quality, and liquid-state representations as limitations toward a broader differentiable thermodynamic foundation model.

E. Disordered Solid Solutions: TIP evaluates branch free energies across composition to model configurational disorder, supporting solubility limits and miscibility gaps without fitting a separate configurational Hamiltonian.Calibration can proceed bottom-up against computational references or top-down against experiment.
E. Disordered Solid Solutions: The fine-tuned Al–Sc model more closely reproduces the dedicated quasi-harmonic reference, whose high-temperature solubility exceeds experiment by a factor of three to four.This example evaluates fidelity to the computational reference rather than experimental agreement.
F. Chemical Trends in Predicted Responses: TIP[UMA] predicts per-atom entropy and volume contributions that sum to structure totals, enabling tests of learned chemical trends despite structure-level supervision.These decompositions are learned descriptors rather than unique physical observables.
F. Chemical Trends in Predicted Responses: Predicted entropy and volume generally increase with atomic mass, with large, soft, polarizable elements highest and compact, stiff refractory species lowest.Ta and W are exceptions, remaining compact with low entropy because of strong bonding.
F. Chemical Trends in Predicted Responses: Local environment separates predicted entropy and volume: heavier, more polarizable neighbors raise Ti-site values, while stiff covalent bonds fall below soft, high-coordination metallic bonding.This agrees qualitatively with the microscopic relation between bond stiffness, vibrational frequency, and entropy.
III. DISCUSSION: TIP predicts a branch-conditioned Gibbs free-energy surface whose automatic derivatives yield volume, entropy, and heat capacity across quasi-harmonic and molecular-dynamics fidelities.Calibration can target higher-resolution calculations or experiment.
III. DISCUSSION: The model omits magnetic and electronic entropy, while accuracy also depends strongly on the convergence and coverage of high-fidelity free-energy labels.Magnetic and mixed-valence solids may require additional entropy terms, and tighter switching reversibility is prioritized.
III. DISCUSSION: Extending TIP to liquids requires new absolute reference states and representations independent of relaxed ordered branches, supporting a broader foundation model for routine finite-temperature discovery.The proposed direction targets first-principles and experimental calibration at MLIP-like inference compute cost.

IV. METHODS · A. Branch-conditioned Gibbs free energy · B. Gibbs free energy model

TIP defines Gibbs free energy branch by branch, restricting NPT fluctuations to the basin of each relaxed reference structure. Its atomwise thermodynamic model separates static energy from finite-temperature and pressure contributions, with responses obtained from derivatives of the same free energy.

A. Branch-conditioned Gibbs free energy: Each ordered branch is labeled by a symmetry-constrained relaxed reference structure and restricts the NPT partition function to its corresponding basin.The branch free energy marginalizes vibrational and cell fluctuations within that basin.
A. Branch-conditioned Gibbs free energy: The model splits static branch energy U° from a thermodynamic residual ˆGδ containing finite-(T, P) contributions.Zero-point energy is excluded to match classical MD references and avoid ill-defined contributions for dynamically unstable branches.
A. Branch-conditioned Gibbs free energy: First-order equilibrium responses and second-order responses, including heat capacity, bulk modulus, and thermal expansion, follow from derivatives of G.These responses are thermodynamically linked through the same Gibbs free energy model.
B. Gibbs free energy model: The thermodynamic residual is parameterized atomwise to preserve extensivity and uses per-atom latent features h_i from an equivariant PaiNN backbone.The backbone refines frozen MLIP scalar and vector features through equivariant message passing over local geometry.
B. Gibbs free energy model: Temperature and pressure enter only at the prediction head, while the geometry encoder uses rmax = 6 Å, d = 64 latent features, and 16 radial basis functions.The encoder itself receives neither temperature nor pressure.
B. Gibbs free energy model: At zero pressure, the branch reference combines an Einstein oscillator term with a cubic polynomial in τ := T/Tmax, using Tmax = 4000 K.Pre-training predicts five quantities: c0–c3 and the Einstein temperature θE.
B. Gibbs free energy model: Mid-training adds gated anharmonic τ ln τ and 1/τ terms above θE, preserving the low-temperature quasi-harmonic branch.The gate uses p = 4 and confines corrections to temperatures above θE.
B. Gibbs free energy model: A Murnaghan-like pressure response provides closed-form atomwise volumes and pressure integrals, yielding analytic Gibbs free energies without root solving.The Murnaghan-like and Vinet forms showed comparable accuracy in an ablation.

C. Multi-fidelity dataset

The multi-fidelity dataset combines broad quasi-harmonic Gibbs free-energy labels with higher-fidelity nonequilibrium molecular-dynamics labels on representative structures. Curation, transfer, and quality filters produce accepted thermodynamic surfaces while excluding unstable or non-crystalline simulations.

Source curation and splits: 123,424 metastable Materials Project structures within 0.2 eV/atom of the convex hull were UMA-relaxed, clustered into 26,226 representatives, and split into training and validation sets.The representative split contains 23,603 training and 2,623 validation structures.
Low-fidelity quasi-harmonic labels: The low-fidelity signal uses quasi-harmonic phonons, excludes imaginary-frequency modes, and fits pressure-dependent Gibbs branches with a Vinet equation of state.One (T, P) point is sampled per structure, and surfaces failing equation-of-state residual or bulk-modulus checks are rejected.
Low-fidelity quasi-harmonic labels: 122,913 accepted quasi-harmonic surfaces remain after filtering, comprising 111,402 training and 11,511 validation surfaces.Rejected surfaces have equation-of-state RMS residuals above 5 meV/atom or negative fitted bulk moduli.
High-fidelity molecular dynamics labels: High-fidelity labels use absolute Gibbs free energies from nonequilibrium molecular dynamics at several independently sampled state points on each representative.State points draw T from [100 K, Tm + 500 K] and P from [0, 40] GPa, averaging approximately four per material.
High-fidelity molecular dynamics labels: The high-fidelity free energy transfers from Orb-v3 to UMA through alchemical switching at the same (T, P), avoiding a full Frenkel–Ladd calculation on UMA.The transfer uses forward and backward 5 ps legs with 1 ps UMA re-equilibration and symmetrization.
Quality filtering and dataset: Quality filters reject simulations with cell-edge scale factors outside [0.5, 2.0] or time-averaged atomic displacements above 1.0 Å, excluding collapse, sublimation, melting, or reconstruction.The NPT cell volume must equilibrate within 8,000 of 10,000 steps.

D. Multi-fidelity training · E. Experimental alignment for ordered phases · F. Disordered solid solutions

The model is trained across thermodynamic fidelities, calibrated to experimental ordered-phase behavior, and extended to disordered solid solutions through lightweight adapter updates. These procedures enable response-consistent free-energy modeling, phase-diagram construction, dilute-solubility prediction, and phase-boundary calibration.

D. Multi-fidelity training: Pre-training supervises the thermodynamic residual and its differentiated responses, covering the full first- and second-order (T, P) structure of G.The responses include volume, entropy, heat capacity, isothermal bulk modulus, and thermal expansion, with extensive terms normalized per atom.
D. Multi-fidelity training: Mid-training applies a rank-8, scaling-16 LoRA update with all other weights frozen, adding ∼44,000 trainable parameters while switching targets to the MD Gibbs residual.Two anharmonic coefficients are initialized to zero, and the UMA-relaxed structure remains the input.
D. Multi-fidelity training: A smooth-L1 heat-capacity penalty with tolerance 1 J K−1 (mol atoms)−1 preserves the pre-trained temperature dependence while permitting an absolute anharmonic correction.Optimization uses Adam at learning rate 2 × 10−4 for approximately ∼87,000 steps.
E. Experimental alignment for ordered phases: For ordered phases, Holland–Powell calibration matches relative Gibbs stabilities and response derivatives rather than absolute G because earlier stages retain a PBE-level offset.The calibration uses the ds62 assessment for predefined Si–Al–Mg–Ca–O minerals.
E. Experimental alignment for ordered phases: The calibrated SiO2 surface constructs a pressure–temperature phase diagram by selecting the minimum-G polymorph across 300–2400 K and 0–14 GPa, with boundaries from zero pairwise Gibbs differences.Five polymorphs are evaluated on the temperature-pressure grid.
F. Disordered solid solutions: For disordered solid solutions, ordered branches are calibrated bottom-up against higher-resolution computation or top-down against measured phase boundaries using lightweight adapters with the encoder backbone frozen.Bottom-up fitting updates LoRA layers, while top-down calibration also updates analytic thermodynamic heads.
F. Disordered solid solutions: Bottom-up Al–Sc calibration fine-tunes rank-8, scaling-16 adapters on 491 substitutional supercells spanning pure Al to pure Sc across 65 Sc fractions.The target is the UMA quasi-harmonic Gibbs free energy, and vibrational entropy of dissolution sets the solubility’s order of magnitude.
F. Disordered solid solutions: Top-down Au–Pt solvus calibration differentiably reweights stored SGCMC and metadynamics samples, equates semigrand free energies of two basins, and regularizes binodal fitting with energy-drift and L2 penalties.Reweighting updates adapter parameters without rerunning the sampler.

COMPETING INTERESTS · Supplementary Information for Universal Thermodynamic Interatomic Potentials for Crystalline Materials

The supplementary information identifies the authors and their institutional affiliations, and the authors declare no competing interests.

COMPETING INTERESTS: The authors declare no competing interests.
COMPETING INTERESTS: A related-work fragment mentions K. A. Persson and commentary on The Materials Project.
Supplementary Information for Universal Thermodynamic Interatomic Potentials for Crystalline Materials: The paper lists Juno Nam, Bowen Deng, and Xiaochen Du among its authors.
Supplementary Information for Universal Thermodynamic Interatomic Potentials for Crystalline Materials: Luis Barroso-Luque is also listed as an author.
Supplementary Information for Universal Thermodynamic Interatomic Potentials for Crystalline Materials: Juno Nam, Xiaochen Du, and Rafael Gómez-Bombarelli are affiliated with MIT’s Department of Materials Science and Engineering.
Supplementary Information for Universal Thermodynamic Interatomic Potentials for Crystalline Materials: Bowen Deng, Luis Barroso-Luque, and Benjamin Kurt Miller are affiliated with Fundamental AI Research at Meta in San Francisco.
Supplementary Information for Universal Thermodynamic Interatomic Potentials for Crystalline Materials: Benjamin Kurt Miller and Rafael Gómez-Bombarelli are also listed as authors.
Supplementary Information for Universal Thermodynamic Interatomic Potentials for Crystalline Materials: The supplementary information is dated August 17, 2026.

A. THEORETICAL FORMALISM AND STATISTICAL BASIS … 3. Branch-Conditioned Gibbs Free Energy

The TIP formalism defines a Gibbs free energy conditioned on a crystal branch, with branch representatives, stability assumptions, and thermodynamic responses specified consistently. It marginalizes vibrational and cell-shape fluctuations within each branch while enabling comparison between competing branches.

1. Notation: A periodic crystal is represented by species, atomic positions modulo lattice translations, and a cell matrix in an upper-triangular positive-diagonal gauge.The cell matrix defines the Bravais lattice and positive cell volume V := det L > 0.
2. Assumptions: The TIP models the Gibbs free energy of an ordered bulk crystal branch indexed by a reference structure x◦.A deterministic, idempotent branch map assigns periodic structures to branch representatives.
2. Assumptions: Branch representatives are obtained by minimizing the base potential under a reference space-group constraint, which admits dynamically stabilized phases such as bcc Ti.For ordinary phases, the constraint is inactive; for dynamically stabilized phases, it prevents requiring the high-symmetry structure to be a mechanical minimum.
2. Assumptions: The formalism assumes branch representability, finite-(T, P) stability or metastability, and a fixed-branch scope for evaluating G(x◦; T, P).Within the queried branch, equilibration is assumed on the simulation timescale while escape to other branches is exponentially rare.
3. Branch-Conditioned Gibbs Free Energy: A TIP targets the Gibbs free energy conditioned on one branch label rather than the free energy summed over all branches.Competing branches are compared after evaluating their respective branch-conditioned free energies.
3. Branch-Conditioned Gibbs Free Energy: The branch-conditioned Helmholtz and Gibbs free energies arise from canonical and isothermal-isobaric partition functions that integrate atomic and cell degrees of freedom within the branch.At fixed cell, the canonical partition function includes indistinguishability and species-resolved thermal de Broglie factors; the NPT partition function integrates over branch-consistent cell matrices.
3. Branch-Conditioned Gibbs Free Energy: Equilibrium observables follow from thermodynamic derivatives of G(x◦; T, P), which is decomposed into static reference energy U◦ and finite-(T, P) residual Gδ.The branch-conditioned G marginalizes vibrational and cell-shape fluctuations within one branch; chemically disordered systems require configurational sampling over disjoint branches on a parent lattice.

4. Statistical Basis of the Learned Free Energy · B. DERIVATIONS OF THE FREE ENERGY AND SAMPLING METHODS · 1. Quasi-Harmonic Approximation

TIP[UMA] combines quantum harmonic free energies with classical anharmonic corrections from molecular dynamics, using QHA as a low-fidelity training signal. The QHA derives branch-conditioned Gibbs free energies from volume-dependent phonons but neglects mode coupling, anisotropic relaxation, and unstable harmonic modes.

4. Statistical Basis of the Learned Free Energy: QHA and nonequilibrium molecular-dynamics free energies use different nuclear statistics: QHA is quantum after zero-point removal, whereas Frenkel–Ladd sampling is classical.The classical Frenkel–Ladd free energy uses a classical Einstein reference and classical molecular dynamics.
4. Statistical Basis of the Learned Free Energy: The analytic head retains nuclear quantum effects at the harmonic level while incorporating classical anharmonic corrections from molecular dynamics.This combines the quantum QHA target with classical molecular-dynamics corrections, apart from the zero-point offset.
1. Quasi-Harmonic Approximation: QHA supplies TIP[UMA] with a low-fidelity approximation to G(x◦; T, P) that retains volume-dependent vibrational physics while neglecting anharmonic mode coupling.The approximation is used as the low-fidelity training signal.
B. DERIVATIONS OF THE FREE ENERGY AND SAMPLING METHODS: At fixed cell, the potential expands to second order around cell-constrained relaxed positions, with the force-constant matrix defining harmonic vibrational modes.Lattice periodicity block-diagonalizes the mass-weighted force constants over Brillouin-zone wavevectors and phonon branches.
1. Quasi-Harmonic Approximation: The QHA vibrational free energy sums quantum harmonic-oscillator contributions over Brillouin-zone wavevectors, including zero-point and thermal-occupation terms.Cell-dependent phonon frequencies capture thermal expansion and Grüneisen-type effects; force constants come from finite displacements and Phonopy interpolation.
1. Quasi-Harmonic Approximation: Branch-conditioned QHA Gibbs free energy minimizes F(L, T) + PV over branch-consistent cells, using isotropic strain, a strain grid, and a Vinet equation of state.Subtracting the static reference produces the QHA thermodynamic residual used as the low-fidelity regression target.
1. Quasi-Harmonic Approximation: Within isotropic strain, QHA is exact for harmonic vibrations and captures leading volumetric anharmonicity, but neglects anisotropic relaxation and mode-mode coupling.It breaks down when any phonon squared frequency becomes negative along the relevant cell trajectory.
1. Quasi-Harmonic Approximation: For dynamically stabilized branches whose static references are unstable, QHA must be augmented or replaced by finite-temperature methods such as Frenkel–Ladd switching.The limitation arises because unconstrained relaxation destabilizes the static reference.

2. Frenkel–Ladd Switching

Frenkel–Ladd switching computes a crystalline branch’s absolute Helmholtz free energy by thermodynamic integration between the physical potential and an analytically tractable harmonic reference. Converting at the equilibrated cell gives a branch-conditioned Gibbs free-energy estimate, with cell-fluctuation effects treated as a finite-size approximation.

Method: Frenkel–Ladd switching connects the physical Hamiltonian to a harmonic Einstein-crystal reference to compute a single crystalline branch’s absolute Helmholtz free energy.The implementation follows the optimized switching variant of de Koning et al.
Method: The mixed Hamiltonian U(r; λ) := (1 − λ) U(r) + λ U_harm(r) interpolates between physical and harmonic endpoints.The switching path spans λ ∈ [0, 1].
Switching protocol: The smooth schedule dλ/d˜t = 630 ˜t^4(1 − ˜t)^4 vanishes at both endpoints, reducing dissipation relative to a linear ramp.Forward and backward trajectories are combined symmetrically so leading-order dissipation cancels.
Gibbs free-energy estimate: The branch-conditioned Gibbs estimate is G_FL(x◦; T, P) = F_harm(T) − ΔF + PV + F_COM(T), combining switching, pressure–volume conversion, and the center-of-mass correction.The reference free energy is available in closed form for the classical Einstein crystal.
Approximation: Evaluating F at one equilibrated cell rather than integrating cell fluctuations approximates the full NPT Gibbs free energy.The neglected contribution is of order k_BT per simulation cell and negligible per atom for the ~1,000-atom supercells used here, though it can grow for strongly flexible systems.

3. Alchemical Switching · 4. Semigrand Canonical Monte Carlo

Alchemical switching transfers Gibbs free energies between potential surfaces through a symmetrized NPT path, while semigrand Monte Carlo samples composition fluctuations on a shared parent lattice. Well-tempered metadynamics enables traversal of composition basins, and reweighting recovers equilibrium concentrations and phase boundaries.

3. Alchemical Switching: Alchemical switching estimates the Gibbs free energy difference between two potential energy surfaces at the same thermodynamic state without requiring an absolute reference.The method uses a mixed Hamiltonian interpolating between surfaces while holding species and temperature fixed.
3. Alchemical Switching: Both alchemical-switching legs run in the NPT ensemble, with forward-backward symmetrization yielding the Gibbs free energy difference between surfaces.Cell variables respond to temperature and pressure along the switching path.
3. Alchemical Switching: Combining an absolute Frenkel–Ladd estimate on surface A with alchemical switching transfers the absolute Gibbs free energy to target surface B through a triangle relation.For A = Orb-v3 and B = UMA, this avoids a direct Frenkel–Ladd calculation on UMA.
4. Semigrand Canonical Monte Carlo: Semigrand canonical Monte Carlo samples equilibrium species concentrations on a fixed parent lattice shared by all considered species assignments.Continuous coordinates and cell variables are relaxed onto branch representatives and contribute through marginalized branch-conditioned Gibbs free energies.
4. Semigrand Canonical Monte Carlo: The semigrand ensemble fixes N, T, P, and chemical potentials while allowing site assignments to fluctuate, with metadynamics promoting visits across composition basins.Identity swaps change species at individual sites, and only chemical-potential differences are physical.
4. Semigrand Canonical Monte Carlo: Well-tempered metadynamics biases a composition collective variable on the (|A|−1)-simplex, using an effective collective-variable temperature Teff := γT.In the long-time limit, the bias converges to the projected semigrand Gibbs free energy plus a constant.
4. Semigrand Canonical Monte Carlo: Reweighting the biased trajectory recovers the equilibrium concentration distribution, while sweeping chemical potentials at fixed temperature and pressure traces composition-space phase boundaries.Equilibrium concentrations are identified as basin-conditioned mean concentrations of the unbiased semigrand distribution.

5. Differentiable Reweighting

Differentiable reweighting fine-tunes a Gibbs free energy model to reproduce experimental binodal concentrations using stored SGCMC plus WTMetaD samples. Coexistence chemical potentials are obtained by solving equal basin free energies, enabling gradient-based calibration without additional simulations.

Calibration: Fine-tuning uses differentiable reweighting of stored SGCMC plus WTMetaD samples to match experimental binodal concentrations.The method specializes to binary alloys using concentration and chemical-potential difference as collective variables.
Reweighting: Candidate semigrand observables are estimated by importance weights that correct Gibbs-model changes and unbias the accumulated metadynamics bias.The simulator itself is not differentiated; gradients flow through reweighted observables.
Binodal calibration: Binodal concentrations are calibrated from basin-conditioned mean concentrations defined using masks around the two coexisting composition basins.The method uses basin-local weights rather than the full free energy profile.
Coexistence: Coexistence is enforced by solving equal semigrand basin free energies for the chemical potential difference, using stored samples without additional simulation.This equality is identified with the common-tangent construction for the composition free energy.
Optimization: Adam optimization uses learning rate 5 × 10−6, gradient clipping at norm 5, and 160 gradient steps while regularization limits distribution drift.The drift penalty maintains overlap with the sampled distribution, and an L2 penalty keeps parameters near the reference.

C. TRAINING DETAILS · 1. Multi-Fidelity Training · 2. Analytic Head Constants and Numerical Evaluation

Training combines capped, normalized multi-fidelity response losses with a compact thermodynamic adapter, while fixed normalizations, units, and numerical floors make the analytic head well-defined across structures and sampled thermodynamic conditions.

1. Multi-Fidelity Training: The pre-training objective applies L1 penalties to residual responses, normalizes extensive quantities per atom, leaves intensive responses unnormalized, caps per-sample contributions, and clips the global gradient norm at 10.The caps use the units shown in Table S2.
1. Multi-Fidelity Training: The frozen UMA-S-1.1 encoder carries ∼1.5×10^8 parameters, with ∼6 × 10^6 active per structure through mixture-of-experts routing.The active parameter count reflects the encoder’s mixture-of-experts routing.
1. Multi-Fidelity Training: The trainable thermodynamic adapter has ∼2.1 × 10^5 parameters, including ∼9 × 10^3 analytic coefficient-head parameters.These counts describe the adapter relative to the frozen encoder.
1. Multi-Fidelity Training: Mid-training adds a ∼4.4 × 10^4-parameter LoRA update, bringing the adapter to ∼2.6 × 10^5 parameters, with LoRA comprising about 17%.The LoRA update is introduced mid-training.
2. Analytic Head Constants and Numerical Evaluation: The analytic head uses fixed global scales Tmax = 4000 K and Pmax = 40 GPa, so reduced temperature and normalized pressure map sampled ranges into (0, 1] without material-specific normalization.Consequently, the predicted surface is a function of structure alone.
2. Analytic Head Constants and Numerical Evaluation: Head quantities use per-atom energies in eV, temperatures in K, per-atom volumes in ˚A3, and pressures in GPa, while reference coefficients are dimensionless.The fixed unit system applies throughout the analytic head.
2. Analytic Head Constants and Numerical Evaluation: Regularization fixes θmin = 10−4, τmin = 10−8, and gate exponent p = 4 for normalization and vibrational terms, with V0,min = 10−6 ˚A3, B0,min = 10−3 GPa, and xmin = 10−8 keeping the equation of state well-defined.The pressure-volume work is converted from ˚A3 GPa to eV using the stated conversion constant.
2. Analytic Head Constants and Numerical Evaluation: An initialization bias Boff = 100 GPa is added inside the B0,i softplus so the untrained bulk modulus starts near 100 GPa.The pressure dependence in the Murnaghan volume integral uses the dimensionless ratio B′0,iP/B0,i.

D. HOLLAND–POWELL CALIBRATION AND BENCHMARK DETAILS

The Holland–Powell ds62 dataset anchors ordered-phase calibration with experimental single-phase properties for Si–Al–Mg–Ca–O crystalline phases. Quaternary phases and their reactions are held out to evaluate transfer from binary and ternary training phases.

Dataset and benchmark: The ds62 benchmark covers pre-defined Si–Al–Mg–Ca–O crystalline phases and reactions among them.Experimental single-phase properties are provided at 298.15 K and 1 bar.
Dataset and benchmark: 298.15 K and 1 bar define ds62 measurements of B0, α, S(298 K), and V0 for each phase.These are the isothermal bulk modulus, thermal expansivity, standard entropy, and molar volume.
Dataset and benchmark: Six quaternary oxide phases and their forming reactions are excluded from fine-tuning and reserved for evaluation.Held-out accuracy measures transfer from binary and ternary training phases to quaternary phases.
Calibration objective: The post-training objective uses smooth-L1 loss with an L2 weight regularizer, normalized reaction, polymorph, and shape terms, and ds62 single-phase anchors.Reaction and polymorph losses cover 15 training reactions and 16 polymorph pairs.