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

May 15, 2026

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

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.

Sauradeep Majumdar, Miguel Steiner, Johannes C. B. Dietschreit +4

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

Ask PDF

Apr 30, 2026

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

While machine-learned interatomic potentials (MLIPs) accelerate phonon dispersion calculations, merely identifying dynamical instabilities in computationally predicted materials is insufficient; automated pathways to resolve them are required. We introduce VibroML, an open-source Python toolkit driven by foundational MLIPs that shifts the paradigm from stability verification to automated structural remediation. VibroML employs an energy-guided genetic algorithm that vastly outperforms traditional soft-mode following, efficiently navigating the potential energy surface to uncover diverse, dynamically stable polymorphs. As 0 K harmonic stability does not guarantee macroscopic viability, an automated molecular dynamics workflow evaluates finite-temperature structural retention. VibroML also couples with ProtoCSP, our combinatorial structure prediction engine, to stabilize frustrated crystal topologies via targeted alloying, successfully rescuing functional perovskite networks like Cs

_2

KInI

_6

and KTaSe

_3

. Demonstrating broader applicability, we mined the Alexandria database -- where ~50% of quaternary and 99.5% of quinary elemental combinations lack any structural entries -- to identify thousands of abandoned, high-symmetry stoichiometries. Deploying ProtoCSP's "cold start" retrieval and VibroML's evolutionary search on a sample, we successfully identified dynamically stable low-symmetry candidates. Through integrated structural remediation, thermal validation, and systematic compositional exploration, VibroML enables a comprehensive deep-screening approach, yielding physically sound structural propositions that far surpass standard high-throughput workflows.

Rogério Almeida Gouvêa, Gian-Marco Rignanese

Institute of Condensed Matter and Nanosciences, Universit´e Catholique de Louvain, Louvain-la-Neuve,2026 Belgium. · WEL Research Institute, Avenue Pasteur 6, Wavre, Belgium.

Crystal Structure Prediction Materials Science Crystal Structure Generation

Ask PDF

Apr 15, 2026

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

We present an open Molecular Crystal (MC) database of Machine-Learned Interatomic Potentials (MLIP) called MolCryst-MLIPs. The first release comprises fine-tuned MACE models for nine molecular crystal systems---Benzamide, Benzoic acid, Coumarin, Durene, Isonicotinamide, Nicotinic acid , Niacinamide, Pyrazinamide, and Resorcinol---developed using the Automated Machine Learning Pipeline (AMLP), which streamlines the entire MLIP development workflow, from reference data generation to model training and validation, into a reproducible and user-friendly pipeline. Models are fine-tuned from the MACE-MH-1 foundation model omol head), yielding a mean energy MAE of 0.141 kJ/mol/atom and a mean force MAE of 0.648 kJ/mol/Angstrom across all systems. Benchmarked against three state-of-the-art foundation models on the DFT-labelled polymorph set, only the fine-tuned models resolve the polymorphic energy landscape. Dynamical stability and structural integrity, as assessed through energy conservation, P2 orientational order parameters, and radial distribution functions, are evaluated using molecular dynamics simulations. The released models and datasets constitute a growing open database of validated MLIPs, ready for production MD simulations of molecular crystal polymorphism across the polymorphic landscape of each target compound under different thermodynamic conditions.

Adam Lahouari, Shen Ai, Jihye Han +16

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

Machine Learning Interatomic Potentials

Ask PDF

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