Universal Thermodynamic Interatomic Potentials for Crystalline Materials

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

cond-mat.mtrl-sciAug 14, 2026

Universal Thermodynamic Interatomic Potentials for Crystalline Materials

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

Organizations: Department of Materials Science and Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA · Fundamental AI Research, Meta, San Francisco, CA 94105, USA

Abstract

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.

Figures & tables

Fig. 1: Free energy adaptation with a thermodynamic interatomic potential (TIP). (a) A TIP combines three established routes to solid-state thermodynamics: the chemical transferability of MLIPs, the configurational phase stability formalism of the cluster expansion, and the analytic, experiment-aligned free energy representation of CALPHAD. (b) A TIP augments a pre-trained MLIP with a thermodynamic adapter: for each locally relaxed branch representative $\bm{x}^{\circ}$ , the frozen MLIP encodes $\bm{x}^{\circ}$ , and the adapter combines those features with $(T,P)$ to predict the branch Gibbs free energy $\hat{G}(T,P)$ . (c) Architecture: the MLIP backbone’s energy head returns the static reference $U^{\circ}:=U(\bm{x}^{\circ})$ and its per-atom features condition the adapter, whose analytic head outputs the thermodynamic residual $\hat{G}^{\delta}(T,P)$ . Their sum defines $\hat{G}$ , and its derivatives give the response properties. (d) Three-stage training: pre-training on quasi-harmonic $(T,P)$ grids, mid-training on anharmonic absolute free energies from nonequilibrium MD, and post-training against experimental data by response matching for ordered phases (e.g., heat capacity) and differentiable reweighting for disordered phases (e.g., composition–temperature binodals). (e) Wall-clock time per Gibbs free energy evaluation on a given $\bm{x}^{\circ}$ : nonequilibrium MD ( ${\sim}10^{4}$ s) and the quasi-harmonic approximation ( ${\sim}10^{2}$ s) compared with one TIP[UMA] evaluation ( ${\sim}10^{-1}$ s), five and three orders of magnitude faster.

Fig. 2: Multi-fidelity Gibbs free energy dataset. (a) Data-generation workflow. Structures from a materials database are filtered and relaxed with the base potential $U_{\mathrm{base}}$ (UMA [ 66 ] ); the quasi-harmonic approximation (QHA) on each yields the low-fidelity QHA Gibbs dataset $G^{\mathrm{QHA}}(T,P)$ . Clustering selects representatives, which are equilibrated under $NPT$ dynamics with a fast potential $U_{\mathrm{fast}}$ (Orb-v3 [ 47 ] ) and then processed by $NVT$ Frenkel–Ladd switching [ 24 ] (absolute free energy against a harmonic Einstein crystal reference) and $NPT$ alchemical switching (transfer from the fast to the base Hamiltonian); combining the two contributions gives the high-fidelity MD Gibbs dataset $G^{\mathrm{MD}}(T,P)$ . (b) Elemental coverage of the MD Gibbs dataset: the periodic table is colored by the number of dataset entries containing each element. (c) Distribution of the residual Gibbs free energy $G^{\delta}$ of the MD Gibbs dataset versus temperature (left) and pressure (right), shown as two-dimensional histograms with the marginal count distributions above; white lines mark the quartiles (Q1, median, Q3).

Fig. 3: Validation error analysis. (a) Element-wise validation error. Periodic-table map of the median absolute Gibbs error for each element, computed over the filtered validation data points that contain that element; elements represented by fewer than $20$ validation state points are shown in gray. The black tick on the colorbar marks the global median absolute error. (b) Thermodynamic and dissipation trends. Signed Gibbs error distributions as functions of temperature, pressure, and Frenkel–Ladd (FL) dissipation for the filtered validation set. Left panels show two-dimensional count histograms of the mid-trained model on a logarithmic color scale; right panels show the binned median absolute error and the signed median error (meV atom -1 ) for the pre-trained (light) and mid-trained (dark) models. The dotted line in the absolute-median panels marks the global median absolute error.

Fig. 4: Thermodynamics of ordered crystalline phases from a learned free energy surface. (a) Mid-trained mean absolute error in heat capacity $C_{P}$ and standard entropy $S(298~\mathrm{K})$ across the NIST–JANAF benchmark of $255$ stoichiometric solids, with medians of $1.94$ ( $C_{P}$ ) and $2.28$ ( $S$ ) J K -1 (mol atoms) -1 . (b) Mid-trained heat capacity $C_{P}(T)$ (solid lines) against NIST–JANAF values [ 12 ] (open circles). The larger deviations for $\text{Cr}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}{\vphantom{\text{X}}}_{\smash[t]{\text{3}}}$ and $\text{Fe}{\vphantom{\text{X}}}_{\smash[t]{\text{3}}}\text{O}{\vphantom{\text{X}}}_{\smash[t]{\text{4}}}$ occur near magnetic $\lambda$ -anomalies (Néel and Curie transitions) outside the model’s vibrational scope. (c) Isothermal compression $V(P)$ at $300$ K and (d) thermal expansion $V(T)$ at $P\approx 0$ for forsterite $\text{Mg}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{4}}}$ , comparing experiment (Holland–Powell [ 29 ] ), the UMA quasi-harmonic reference, and TIP[UMA] across pre-, mid-, and post-training. The quasi-harmonic reference and the pre- and mid-trained models share the known PBE volume overestimate, which post-training removes. (e–h) Temperature-driven polymorphic transitions, shown as the per-atom Gibbs difference $\Delta G(T)$ between the high- and low-temperature phases for the pre-trained and mid-trained models: (e) Ti , (f) Hf , (g) $\text{CaSiO}{\vphantom{\text{X}}}_{\smash[t]{\text{3}}}$ , and (h) Sc . The $\Delta G=0$ crossing (circle) is the predicted transition temperature, and the dashed line and star mark experimental transition [ 22 , 48 ] with the shaded band spanning a $\pm 250$ K tolerance. The pre-trained surface shows no crossing ( Ti , Sc ) or one far above experiment ( Hf , $\text{CaSiO}{\vphantom{\text{X}}}_{\smash[t]{\text{3}}}$ ), whereas mid-training recovers the transitions, Sc only in part. (i) Accuracy of the pre-, mid-, and post-trained models on the Holland–Powell ordered-oxide benchmark for the reaction Gibbs free energy, standard entropy, and bulk modulus $B_{0}$ , on a training partition (binary and ternary oxides) and a held-out partition (quaternary oxides). Post-training improves every metric on both, including the held-out quaternaries. (j–l) $\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$ pressure–temperature phase diagram for (j) the mid-trained model, (k) the post-trained model, and (l) the experimentally assessed topology [ 56 ] . The mid-trained model spuriously stabilizes cristobalite at ambient conditions, whereas post-training recovers the low-pressure quartz field and the quartz, coesite, stishovite ordering within the supplied candidate set. Red asterisks mark comparisons for which the experimental value was included in calibration.

Fig. 5: Calibration of TIP[UMA] for disordered phases using computational and experimental references. Lightweight fine-tuning fits the frozen free energy surrogate to a chosen reference and extends it from ordered branches to configurationally disordered solid solutions, preserving thermodynamic consistency. (a) Dilute Sc solubility in fcc Al versus temperature, obtained from the dissolution free energy relative to the adjacent Al – $\text{Al}{\vphantom{\text{X}}}_{\smash[t]{\text{3}}}\text{Sc}$ convex hull ( Section IV ). By accounting for vibrational entropy, TIP[UMA] predicts the measured order of magnitude and tracks the quasi-harmonic UMA reference (UMA QHA), whereas the $0$ K static estimate falls about an order of magnitude below the measured value; Al – Sc -focused fine-tuning (pre to post) moves the model toward that reference. The computed reference itself lies above the measured solubility, an offset set by the reference thermodynamics rather than the surrogate. (b) Workflow for a disordered phase: TIP[UMA] scores lattice configurations, well-tempered semigrand Monte Carlo samples the mixing free energy, a common-tangent construction yields the miscibility gap, and differentiable reweighting calibrates the surrogate. (c) Au – Pt miscibility gap, composition in Pt atomic percent (at.%): iterated differentiable reweighting moves the TIP[UMA] binodal from its initial prediction to the experimental solvus, the Pt-concentration error falling over the iterations (inset). Horizontal bars are the seed spread, one standard deviation of the branch Pt concentration across four independent well-tempered metadynamics seeds. Red asterisks mark comparisons for which the experimental value was included in calibration.

Fig. 6: Chemical trends in predicted entropy and volume. Predicted per-atom entropy $\hat{S}$ and volume $\hat{V}$ from TIP[UMA] at $T=1000$ K and $P=0$ GPa, across the filtered high-fidelity dataset. (a) Periodic table colored by each element’s median predicted entropy. (b) Element-level median entropy versus median volume, colored by atomic mass, with faint guide lines connecting elements within the same group; the two are strongly correlated. (c–f) Site-resolved entropy–volume density maps, each grouping sites by one aspect of local chemistry over all sites of the central element, with diamonds marking the group medians. (c) Neighbor identity, Ti –X, for neighbors O , S , Se , and Te . (d,e) Bond character: covalent B – N and P – S / P – Se sites versus metallic Fe , Co , and Ni borides and phosphides. (f) Formal valence, Pb – O : compact $\text{Pb}{\vphantom{\text{X}}}^{\text{4+}}$ versus expanded $\text{Pb}{\vphantom{\text{X}}}^{\text{2+}}$ , softened by a stereochemically active $6s^{2}$ lone pair.

Appendix figures & tables4 assets

Supplementary material from the paper’s appendix.

Appendix

Symbol	Description
$N$	Number of atoms in a crystal
$\bm{x}$	Crystal structure $(\bm{L},\bm{z},[\bm{r}])$
$\bm{x}^{\circ}$	Branch representative (relaxed structure)
$\bm{L}$	Cell matrix ( $\mathrm{GL}^{+}(3,\mathbb{R})$ )
$V$	Cell volume ( $\det\bm{L}>0$ )
$\bm{z}$	Atomic species in a crystal ( $\mathcal{A}^{N}$ )

Appendix

Table S1: Summary of the notation.

Term	Units	Per atom	$\lambda$	Cap
$\mathcal{L}_{G}$	eV atom -1	yes	$1$	—
$\mathcal{L}_{V}$	Å 3 atom -1	yes	$0.25$	—
$\mathcal{L}_{S}$	J K -1 (mol atoms) -1	yes	$0.005$	$200$
$\mathcal{L}_{C_{P}}$	J K -1 (mol atoms) -1	yes	$0.003$	$200$
$\mathcal{L}_{B}$	GPa	no	$0.003$	$2000$
$\mathcal{L}_{\alpha}$	$10^{-5}$ K -1	no	$0.05$	$200$

Appendix

Table S2: Pre-training loss terms. Each is a

L_{1}

penalty on the residual values. The extensive quantities are normalized per atom; the intensive responses are not. The cap (in the units shown) bounds each per-sample contribution, and a global gradient-norm clip of

10

applies to the summed loss.

Abbrev.	Composition	Mineral
q	$\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$	quartz
trd	$\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$	tridymite
crst	$\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$	cristobalite
coe	$\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$	coesite
stv	$\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$	stishovite
per	MgO	periclase

Appendix

Table S3: Holland–Powell benchmark phases. Binary and ternary phases used for training are listed above the divider; quaternary phases used only for evaluation are listed below.

Term	Target	$\lambda$
$\mathcal{L}_{\mathrm{rxn}}$	reaction $\Delta G$	$1$
$\mathcal{L}_{\mathrm{poly}}$	polymorph $\Delta G$	$0.5$
$\mathcal{L}_{B_{0}}$	$B_{0}$	$0.003$
$\mathcal{L}_{S_{0}}$	$S(298~\mathrm{K})$	$0.003$
$\mathcal{L}_{V_{0}}$	$V_{0}$	$0.015$
$\mathcal{L}_{V}$	$V(T,P)$	$0.05$

Appendix

Table S4: Post-training loss terms for ordered phases. Per-term weights

\lambda

of the experimental calibration objective (Methods).

Computational Materials Science

Computational materials science uses computation and data-driven models to predict, explain, and design materials and their properties.

See all 143 papers

Reweighting free energy profiles between universal machine learning interatomic potentials for fast consensus building
May 15, 2026Sauradeep Majumdar, Miguel Steiner, Johannes C. B. Dietschreit +4Free Energy Calculation Machine Learning Interatomic Potentials
Department of Material Science and Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, 02139, United States. · Institute of Theoretical Chemistry, University of Vienna, Währinger Str. 17, 1090, Vienna, Austria. · Department of Physical and Macromolecular Chemistry, Charles University, Hlavova 8, Praha 2, 12800, Czech Republic.
VibroML: an automated toolkit for high-throughput vibrational analysis and dynamic instability remediation of crystalline materials using machine-learned potentials
Apr 30, 2026Rogério Almeida Gouvêa, Gian-Marco RignaneseCrystal Structure Prediction Materials Science
Institute of Condensed Matter and Nanosciences, Universit´e Catholique de Louvain, Louvain-la-Neuve,2026 Belgium. · WEL Research Institute, Avenue Pasteur 6, Wavre, Belgium.
MolCryst-MLIPs: A Machine-Learned Interatomic Potentials Database for Molecular Crystals
Apr 15, 2026Adam Lahouari, Shen Ai, Jihye Han +16Machine Learning Interatomic Potentials
Department of Chemistry, New York University, New York, NY 10003, USA · Department of Physics, New York University, New York, NY 10003, USA · Simons Center for Computational Physical Chemistry, New York University, New York, NY 10003, USA +5

May 15, 2026J/K move · Enter open · S save

Reweighting free energy profiles between universal machine learning interatomic potentials for fast consensus building

Sauradeep Majumdar, Miguel Steiner, Johannes C. B. Dietschreit, Swagata Roy+3

Department of Material Science and Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts, 02139, United States. · Institute of Theoretical Chemistry, University of Vienna, Währinger Str. 17, 1090, Vienna, Austria. · Department of Physical and Macromolecular Chemistry, Charles University, Hlavova 8, Praha 2, 12800, Czech Republic.

Free Energy Calculation Machine Learning Interatomic Potentials Importance Sampling Computational Materials Science

Open paper AskPDF arXiv

Free energy profiles serve as a fundamental bridge between microscopic atomic fluctuations and macroscopic thermodynamic observables. Estimating the free energy profile along a reaction coordinate, referred to as the potential of mean force (PMF), with density functional theory (DFT) accuracy is computationally expensive. Universal machine learning interatomic potentials (MLIPs) drastically reduce this cost, but their accuracy is strongly determined by their training data and hence can be uncertain for a given system. In this work, we present a systematic and scalable framework for reweighting PMFs, initially sampled with a single 'source' MLIP, across a representative suite of target MLIPs. Because traditional direct exponential reweighting fails for large system sizes due to low phase-space overlap between potentials, we deploy robust analytical corrections. Applying this to a complex 601-atom system of Li

^+

transport in a nanoconfined electrolyte, we demonstrate that a mean energy-gap approximation effectively bypasses statistical collapse, producing a highly stable PMF matching the target PMF. Using this approach, we recover high-fidelity target thermodynamics across multiple DFT reference levels (PBE+D3, PBE-sol, r

^2

SCAN,r

^2

SCAN-D4) at a fraction of the computational cost of full simulations. Furthermore, thermodynamic analysis reveals that the studied MLIPs partition into two distinct clusters driven by their training data. Our reweighting framework successfully recovers target thermodynamic properties--specifically, reaction and activation free energies--even when the phase-space overlap between potentials is critically low. Ultimately, this approach establishes a vital diagnostic protocol to achieve affordable cross-model consensus on materials chemistry properties without redundant, resource-intensive simulations.

Abstract

Computational Materials Science

Explore similar work

Reweighting free energy profiles between universal machine learning interatomic potentials for fast consensus building

VibroML: an automated toolkit for high-throughput vibrational analysis and dynamic instability remediation of crystalline materials using machine-learned potentials

MolCryst-MLIPs: A Machine-Learned Interatomic Potentials Database for Molecular Crystals