MACE-POLAR-1: A Polarisable Electrostatic Foundation Model for Molecular Chemistry
Authors: Ilyes Batatia, William J. Baldwin, Domantas Kuryla, Joseph Hart, Elliott Kasoar, Alin M. Elena, Harry Moore, Mikołaj J. Gawkowski, Benjamin X. Shi, Venkat Kapil, Panagiotis Kourtis, Ioan-Bogdan Magdău, Gábor Csányi
Organizations: Engineering Laboratory, University of Cambridge, Trumpington St, Cambridge, UK · Yusuf Hamied Department of Chemistry, University of Cambridge, Lensfield Road, Cambridge, UK · Cavendish Laboratory, University of Cambridge, J. J. Thomson Avenue, Cambridge, UK · Scientific Computing Department, Science and Technology Facilities Council, Daresbury Laboratory, Keckwick Lane, Daresbury WA4 4AD, UK · Ångström AI, San Francisco, California, USA · Department of Physics and Astronomy, University College London, 7-19 Gordon St, London WC1H 0AH, United Kingdom · Initiative for Computational Catalysis, Flatiron Institute, 160 5th Avenue, New York, NY 10010 · School of Natural and Environmental Sciences, Newcastle University, Newcastle upon Tyne, UK · Max Planck Institute for Polymer Research, Ackermannweg 10, Mainz, Germany
Accurate modelling of electrostatic interactions and charge transfer is fundamental to computational chemistry, yet most machine learning interatomic potentials (MLIPs) rely on local atomic descriptors that cannot capture long-range electrostatic effects. We present a new electrostatic foundation model for molecular chemistry that extends the MACE architecture with explicit treatment of long-range interactions and electrostatic induction. Our approach combines local many-body geometric features with a non-self-consistent field formalism that updates learnable charge and spin densities through polarisable iterations to model induction, followed by global charge equilibration via learnable Fukui functions to control total charge and total spin. This design enables an accurate and physical description of systems with varying charge and spin states while maintaining computational efficiency. Trained on the OMol25 dataset of 100 million hybrid DFT calculations, our models achieve chemical accuracy across diverse benchmarks, with accuracy competitive with hybrid DFT on thermochemistry, reaction barriers, conformational energies, and transition metal complexes. Notably, we demonstrate that the inclusion of long-range electrostatics leads to a large improvement in the description of non-covalent interactions and supramolecular complexes over non-electrostatic models, including sub-kcal/mol prediction of molecular crystal formation energy in the X23-DMC dataset and a fourfold improvement over short-ranged models on protein-ligand interactions. The model's ability to handle variable charge and spin states, respond to external fields, provide interpretable spin-resolved charge densities, and maintain accuracy from small molecules to protein-ligand complexes positions it as a versatile tool for computational molecular chemistry and drug discovery.
Figures & tables
Figure 1: Overview of the MACE-POLAR-1 architecture and benchmarked applications. (A) Model architecture. Atomic positions and species are passed to a MACE model, which predicts local node features, a local energy contribution, and an initial set of spin-charge multipoles. The local features richly encode semi-local geometry and chemistry. Then, the spin-charge multipoles are iteratively refined through a long-range operation. First, a physically inspired global convolution maps the spin-charge multipoles into atom-centred electrostatic features vi,nlm . Then, a local operation followed by a global normalisation predicts a new set of spin-charge multipoles. The process is iterated twice, before the final set of multipoles is used to compute a Coulomb energy and an additional learned non-local energy contribution. (B) Construction of long-range electrostatic features. Information is propagated between distant atoms by constructing a smooth charge density ρ(r) from the spin-charge multipoles and convolving with the Coulomb kernel to give an electric potential v(r) . Then, one can project the potential onto atom-centred functions to give equivariant electrostatic features. (C) Physical extrapolation capabilities: charge localisation when fragmenting a cluster, response of the model charge density to an electric field, and correct prediction of the oxidation state of transition metal ions in water. (D) Application domains benchmarked: protein-ligand binding, molecular crystals, transition-metal redox potentials, and supramolecular complexes.
Model
Receptive Field
Long-Range Coulomb
MACE-OMOL
18 Å
✗
MACE-POLAR-1-M
12 Å
✓
MACE-POLAR-1-L
18 Å
✓
UMA-S-1-OMOL
30 Å
✗
UMA-M-1-OMOL
66 Å
✗
ORBMOL
30 Å
✗
Table 1: Receptive fields and explicit long-range physics in benchmarked models.
Figure 2: Thermochemistry and reaction barrier benchmarks on GSCDB138 subsets. Bar heights show the weighted total mean absolute deviation (WTMAD-2) in kcal/mol for each model, where lower values indicate better accuracy. WTMAD-2 rescales errors by the characteristic energy scale of each subset to enable fair comparison across datasets of different magnitudes (definition in SI V.1 ). For these summary bars, extreme outliers are filtered by excluding points with ∣ΔE∣>100 kcal/mol (details in SI V.2 ). (a) Thermochemistry subsets grouped by property class: ionisation potentials, electron affinities, proton affinities, and bond dissociation energies. (b) Reaction-barrier subsets: barrier heights, proton-transfer reactions, and general reaction energies. Models compared: ω B97M-V (reference hybrid DFT), g-xTB (semi-empirical), UMA-S-1P1 /UMA-M-1P1 , MACE-OMOL , ORBMOL (local MLIPs), and the electrostatic MACE-POLAR-1-M /MACE-POLAR-1-L .
Subset
Nsets
wB97M-V
g-xTB
UMA-S-1P1
UMA-M-1P1
MACE-OMOL
ORBMOL
MACE-POLAR-1-M
MACE-POLAR-1-L
SO3LR
AIMNet2-NSE
FENNIX-BIO-2
TC: Bond Energies + HAT
7
1.57
5.30
3.16
2.18
4.68
5.63
2.75
3.80
26.46
9.27
14.18
TC: Reaction Energies
6
0.94
3.61
1.06
1.14
1.77
1.31
1.51
1.47
18.66
3.37
6.76
TC: Thermochemistry
9
1.22
5.31
2.22
2.47
3.27
2.94
2.39
2.34
27.22
4.87
11.15
BH: Barrier + Proton
10
1.27
3.56
1.79
1.56
2.06
2.62
2.69
1.91
14.27
6.75
18.32
Table 2: WTMAD-2 (kcal/mol) for thermochemistry and reaction-barrier subsets with only H, C, N, O, F, P, S, and Cl elements to support SO3LR and AIMNet2-NSE. Datasets for each subset are listed in SI Sec. V.4 .
Figure 3: Comprehensive evaluation of non-covalent interaction accuracy. (a) Bar heights show mean absolute errors in kcal/mol for small-molecule non-covalent interaction datasets with CCSD(T)/CBS references: S22 and S66 (hydrogen bonding and dispersion), XB20 (halogen bonding), X40 (mixed interactions), WATER27 (water clusters), HB49 (diverse hydrogen bonds), NC11 (charge-transfer complexes), and O24x4 (potential energy curves). (b) Mean absolute errors for protein-ligand fragment benchmarks: QUID (quantum-chemistry dimers from pocket-ligand motifs) and PLF547 (protein-ligand fragments with MP2-F12 + DLPNO-CCSD(T) references). (c) Mean absolute errors on the IHB100x10 ionic hydrogen bond dataset from NCI Atlas, where electrostatic polarisation is critical. (d) Potential energy curves for gas-phase alkali halide dissociation (LiCl, NaCl, KBr), plotting energy versus interatomic distance to test long-range 1/r Coulombic behaviour. (e) Mean absolute errors for PLA15 complete protein-ligand active sites (259–584 atoms). (f) Mean absolute errors for S30L supramolecular host-guest complexes (up to 200 atoms, charge states −1 to +4 ). (g) Mean absolute errors for X23-DMC molecular crystal lattice energies with diffusion Monte Carlo references. In all bar charts, lower values indicate better accuracy.
Figure 4: Absolute lattice energy errors for CPOSS209 molecular crystals. Bar heights show mean absolute errors in kcal/mol for predicted lattice formation energies, grouped by molecular family. The dataset comprises 209 experimental and predicted polymorphs from 20 small drug molecules. Reference calculations are performed at the ω B97M-D3(BJ) level with 1-body CCSD(T) corrections. Lower values indicate better accuracy.
Figure 5: Accuracy on transition metal complexes and molecular conformers. Bar heights show mean absolute errors in kcal/mol for each model; lower values indicate better accuracy. (a) Transition metal datasets on logarithmic scale due to large error ranges: CUAGAU83 (coinage metal Cu, Ag, Au complexes), DAPD (palladium diatomics), MOBH28 (organometallic barrier heights), and TMD10 (transition metal diatomics). (b) Transition metal datasets on linear scale: 3dTMV (vertical ionisation energies, ph-AFQMC references), MME52 (metalloenzyme models, DLPNO-CCSD(T) references), ROST61 (open-shell reactions), MOR13 (closed-shell reactions), and TMB11 (barrier heights). (c) Conformational energy benchmarks: 37CONF8 (small organics), ACONFL (n-alkane conformers), DipConfS (amino acids and dipeptides), Maltose222 (carbohydrates), MPCONF196 (medicinal fragments), OpenFF-Tors (torsional profiles), and UPU46 (RNA backbone fragments). All reference values are CCSD(T) or equivalent. The GSCDB138 transition metal sets use updated references with spin-contaminated structures removed.
Figure 6: Liquid water density as a function of temperature. Data points show the equilibrium density (g/cm 3 ) from NPT molecular dynamics simulations of 333 H 2 O molecules at 1 atm pressure. Each point represents the mean density from 500 ps of production dynamics following 500 ps of equilibration. Coloured symbols correspond to different models as indicated in the legend. Experimental data (black circles) show the characteristic density maximum near 277 K arising from competition between thermal expansion and the tetrahedral hydrogen-bond network.
Figure 7: Differences between predicted and experimental radial distribution functions for O-O, H-H, and O-H pairs in liquid water at 300 K. The radial distribution function g(r) measures the probability of finding an atom pair at distance r relative to an ideal gas. Left panels show differences between predicted and experimental g(r) from a 500 ps NPT simulation; right-panel dashed lines show differences in g(r) from NVT simulations at the experimental density of 0.997 g cm -3 . The inset indicates the reference experimental neutron diffraction data. Coloured lines correspond to different models as indicated in the legend.
Figure 8: Absolute errors on experimental densities for 62 organic liquids. Each point represents the absolute error on the equilibrium density (g/cm 3 ) of one molecular liquid computed from NPT molecular dynamics at 1 atm with approximately 1000 atoms for 500 ps of equilibration followed by 500 ps of simulation. The dataset spans alkanes (n-hexane to n-dodecane), alcohols (methanol to octanol), ethers (diethyl ether, THF), aromatics (benzene, toluene, xylenes), ketones (acetone, MEK), and polar solvents (acetonitrile, DMF, DMSO). The diagonal black dashed line indicates perfect agreement; points above the line indicate overestimated density.
Model
Density MAE, g cm -3
MACE-OMOL
0.126
MACE-POLAR-1-M
0.091
MACE-POLAR-1-L
0.077
UMA-S-1P1
0.102
Table 3: Mean absolute errors (MAE) in predicted organic liquid densities.
Figure 9: Stability of solvated Cl 2 during geometry relaxation. A neutral Cl 2 molecule is solvated in water alongside two distant Cl - ions (total system charge −2 ). Because the chloride ions lie beyond the message-passing cutoff ( > 6 Å), local models cannot distinguish neutral Cl 2 from hypothetical Cl 22− based on local environment alone. Left: Cl-Cl bond length (Å) versus relaxation step. Non-electrostatic models (MACE-OMOL , UMA-S-1P1 , ORBMOL ; red/pink/gray traces) drive unphysical bond dissociation, while electrostatic MACE-POLAR-1 models (teal traces) preserve the covalent bond near its gas-phase equilibrium value ( ∼ 2.0 Å). Shaded regions indicate bonded (green, < 2.5 Å) and dissociated (pink, > 3.5 Å) regimes. Right: Molecular snapshots showing the initial solvated Cl 2 configuration (top) and final relaxed structures for each model.
Figure 10: Solvation structure and charge dynamics of aqueous iron chloride systems. Left column: Fe atomic partial charge (in charge atomic units, ∣e∣ ) as a function of simulation time over 20 ps of NVT molecular dynamics at 300 K. Results are shown for MACE-POLAR-1-M (top) and MACE-POLAR-1-L (bottom), comparing three systems: FeCl 2 with net charge 0 (blue), FeCl 3 with net charge 0 (red), and FeCl 2 with net charge +1 (teal). Dotted horizontal lines indicate time-averaged Fe charges. The similar charge values for FeCl 3 (0) and FeCl 2 (+1) indicate correct localisation of charge to Fe 3+ in both cases. Right column: Fe-O radial distribution functions g(r) showing the probability of finding water oxygen atoms at distance r from the Fe centre. All three systems are overlaid in each panel for MACE-OMOL , UMA-S-1P1 , MACE-POLAR-1-M , and MACE-POLAR-1-L . Vertical dotted lines mark first and second coordination shell peaks. The first-shell contraction from Fe 2+ to Fe 3+ reflects the smaller ionic radius at higher oxidation state.
Figure 11: Vertical ionisation energies of hydrated first-row transition metal ions. Line plots show the computed ionisation energy in eV for removing one electron from M 2+ to form M 3+ at fixed geometry, for metals Ti, V, Cr, Mn, Fe, Ni, and Cu. (a) M-W6 clusters: metal ion coordinated by six water molecules in octahedral geometry. (b) M-W18 clusters: same ions with an explicit second solvation shell of 12 additional water molecules. Black circles show DLPNO-CCSD(T) reference values from Bhattacharjee et al. [ 13 ] . Coloured lines show ML potential predictions: MACE-POLAR-1-M (cyan), MACE-POLAR-1-L (teal), MACE-OMOL (red), UMA-S-1P1 (light gray), and UMA-M-1P1 (dark gray). Mean absolute errors relative to the reference are given in the legend. The 4–5 eV reduction from M-W6 to M-W18 reflects electrostatic screening by the second solvation shell.
Self-distributed ( E , V vs SHE)
On the MACE-POLAR-1-L trajectories ( E , V vs SHE)
Ion
Spin multiplicities
Eexp
MACE-POLAR-1-M
MACE-POLAR-1-L
MACE-OMOL
UMA-S-1P1
MACE-POLAR-1-M
MACE-POLAR-1-L
MACE-OMOL
UMA-S-1P1
Co 2+ /Co 3+
(4,1)
1.92
0.309
0.663
Crashed.
Crashed.
0.383
0.664
10.412
14.394
Fe 2+ /Fe 3+
(5,6)
0.77
1.561
1.380
Crashed.
Crashed.
1.557
1.380
3.313
2.223
Mn 2+ /Mn 3+
(6,5)
1.50
1.734
1.280
Crashed.
Crashed.
1.590
1.280
2.892
-1.664
Ti 2+ /Ti 3+
(3,2)
-0.37
-0.367
-0.273
10.498
Crashed.
-0.350
-0.273
-3.248
-11.807
V 2+ /V 3+
(4,3)
-0.255
0.329
0.514
Crashed.
Crashed.
0.385
0.514
-9.804
0.419
Table 4: Redox potentials E in volts for M 2+ /M 3+ aqueous couples computed from Marcus-theory free energy. Model values have been shifted by a constant estimated from a linear regression fit to account for differences in reference potential. Experimental potentials Eexp (V vs SHE) are taken from the same dataset as in the main text. Cell colours indicate the absolute deviation ∣E−Eexp∣ , from green (small error) to red (large error).
Figure 12: External electric field response properties. Scatter plots compare model predictions (y-axis) against ω B97M-V reference values (x-axis) for molecular response properties. Left: Dip146 dataset of molecular dipole magnitudes (Debye). Right: HR46 dataset of static isotropic polarisabilities (Å 3 ). The diagonal line indicates perfect agreement; points closer to the diagonal represent more accurate predictions. Notably, the MACE-POLAR-1 models were trained only on ground-state energies and forces, yet predict these response properties through the learned electrostatic physics without explicit supervision.
Figure 13: Lanthanide isomerisation energies. Relative energies (kcal/mol) for seven Ln complexes (La, Ce, Nd, Sm, Eu, Lu; 18 isomers total) from the GFN-FF study [ 154 ] compared against r2SCAN-3c references (Table S4). Points show model ΔE versus reference; dashed line is perfect agreement. Models: MACE-POLAR-1-M /MACE-POLAR-1-L (teal), UMA-S-1P1 /UMA-M-1P1 (gray), ORBMOL (purple). MACE-OMOL is excluded due to unstable topologies on this set.
Figure 14: Charge localisation during water cluster dissociation. Line plots show the partial charge (in charge atomic units, ∣e∣ ) on each of two water cluster fragments as a function of interfragment separation (Å), for systems with total charge (a) −1 , (b) 0 , and (c) +1 (all in ∣e∣ ). Solid and dashed lines show charges on fragments 1 and 2, respectively, as predicted by MACE-POLAR-1-M (blue) and MACE-POLAR-1-L (orange). Black solid lines show ω B97M-V Hirshfeld reference charges; dotted horizontal lines indicate the expected integer charges at infinite separation. At short distances (1–2 Å), polarisation in the contact regime distributes charge across both fragments. As fragments separate beyond 4 Å, excess charge localises onto one fragment, converging to the physical −1/0 or +1/0 limits.
Model
Single-atom
Outlier >100
Total
% excluded
MACE-POLAR-1-M
86
9
95
3.08
MACE-POLAR-1-L
86
11
97
3.14
MACE-OMOL
86
96
182
6.26
UMA-S-1P1
86
66
152
4.92
UMA-M-1P1
86
59
145
4.69
ORBMOL
86
121
207
6.70
Table 5: Breakdown of excluded structures for the thermochemistry and reaction-barrier summary plots (Fig. 2 ). Columns report single-atom exclusions ( num_atoms ≤1 ) and outliers with ∣ΔE∣>100 kcal/mol (summary bars only). The final column gives the fraction of excluded structures relative to the total number of structures evaluated for each model on these subsets.
Models
hyper-parameter
MACE-POLAR-1-M
MACE-POLAR-1-L
max_ell
3
3
correlation
3
3
max_L
1
1
num_channels_edge
128
128
num_channels_node
512
512
Table 6: Hyper-parameter settings for the two MACE-POLAR-1 variants used in this work.
Figure 15: Distribution of signed interaction - energy errors (kcal mol -1 ) across non-covalent benchmarks. Boxplots compare ω B97M - V, r2SCAN, UMA-M-1P1 , UMA-S-1P1 , MACE-POLAR-1-L , and MACE-POLAR-1-M . Panels use symmetric log scaling only when outliers require it; y - axis ticks are placed at the lower bound, 0 (if within range), and the upper bound.
Figure 16: Distribution of signed errors (kcal mol -1 ) for the GSCDB138 reaction - barrier benchmarks. Each panel shows a dataset; boxplots compare six methods ( ω B97M - V, r2SCAN, UMA-M-1P1 , UMA-S-1P1 , MACE-POLAR-1-L , MACE-POLAR-1-M ). Boxes denote the interquartile range with median; whiskers extend to 1.5×IQR; outliers are small, colour - matched dots. The y - axis uses a symmetric logarithmic scale when extreme outliers are present; major ticks are placed at the lower bound, 0 (if within range), and the upper bound.
Figure 17: Distribution of signed thermochemical energy errors (kcal mol -1 ) for GSCDB138. Total atomisation energy (TAE) datasets (e.g., PlatonicTAE6 and W4 - 17 TAE subsets) are excluded by design.
Figure 18: Distribution of signed errors (kcal mol -1 ) for transition - metal datasets in GSCDB138, comparing ω B97M - V, r2SCAN, UMA-M-1P1 , UMA-S-1P1 , MACE-POLAR-1-L , and MACE-POLAR-1-M .
Liquid
Exp. T, K
Exp. density, g cm -3
MACE-OMOL
MACE-POLAR-1-M
MACE-POLAR-1-L
UMA-S-1P1
1,2-dichloroethane
303.0
1.245
1.354
1.370
1.333
1.369
1,2-dimethoxy-ethane
303.0
0.864
0.983
0.963
0.971
0.956
o-xylene
278.0
0.880
1.031
0.965
0.978
0.984
1,3-dioxolane
298.0
1.060
1.166
1.169
1.150
1.176
1,4-dioxane
298.0
1.034
1.168
1.134
1.147
1.160
1-butanol
298.0
0.809
0.916
0.901
0.891
0.894
Table 7: Calculated densities of organic liquids. Experimental temperatures and densities for each liquid are also given in the table. NPT simulations were performed on systems of around 1000 atoms, with 500 ps equilibration and 500 ps production runs.
Figure 19: Relative lattice energy errors for CPOSS209 molecular crystals. Bar heights show mean absolute errors in kcal/mol for predicted relative lattice energies (polymorph energy differences), grouped by molecular family. The dataset comprises 209 experimental and predicted polymorphs from 20 small drug molecules. Reference calculations combine ω B97M-D3 crystal-phase energies with 1-body CCSD(T) corrections. Lower values indicate better accuracy.
Long-range electrostatics and polarization remain central obstacles to extending machine learning interatomic potentials (MLIPs) to ionic, polar, and interfacial systems. Here we introduce a semi-local framework for learning electrostatics from energies and forces using polarizable atomic multipoles. Local equivariant descriptors predict environment-dependent latent monopoles, dipoles, and quadrupoles, while residual non-local charge transfer and polarization are captured by non-self-consistent linear response in induced charges and dipoles. Across four diverse benchmarks and four short-range MLIP architectures, the multipole hierarchy and response terms systematically improve potential energy surface accuracy, with the largest gains in systems where long-range effects are essential. More importantly, physically meaningful electrical responses emerge without direct supervision. The learned latent multipoles yield accurate Born effective charge tensors and infrared spectra in close agreement with experiments. The induced-dipole extension introduces new capabilities: it predicts polarizabilities and thereby enables semi-quantitative Raman spectra for bulk water and hybrid MAPbI3 perovskite, as well as the essential features of the surface-specific vibrational sum-frequency generation spectrum at the water-air interface. In ferroelectric HfO2, the predicted electrical response also captures LO-TO splitting and polarization switching. This systematically improvable, physically transparent framework enables MLIPs trained on standard energy and force labels to predict polarization-sensitive observables.
Yoonjae Park, Dongjin Kim, Daniel S. King +5
Department of Chemistry, UC Berkeley, California 94720, United States · Bakar Institute of Digital Materials for the Planet, UC Berkeley, California 94720, United States · Lawrence Livermore National Laboratory, Livermore, CA, USA +1
Machine learning interatomic potentials (MLIPs) have become widely used tools in atomistic simulations. For much of the history of this field, the most commonly employed architectures were based on short-ranged atomic energy contributions, and the assumption of locality still persists in many modern foundation models. While this approach has enabled efficient and accurate modelling for many use cases, it poses intrinsic limitations for systems where long-range electrostatics, charge transfer, or induced polarization play a central role. A growing body of work has proposed extensions that incorporate electrostatic effects, ranging from locally predicted atomic charges to self-consistent models. While these models have demonstrated success for specific examples, their underlying assumptions, and fundamental limitations are not yet well understood. In this work, we present a framework for treating electrostatics in MLIPs by viewing existing models as coarse-grained approximations to density functional theory (DFT). This perspective makes explicit the approximations involved, clarifies the physical meaning of the learned quantities, and reveals connections and equivalences between several previously proposed models. Using this formalism, we identify key design choices that define a broader design space of self-consistent electrostatic MLIPs. We implement salient points in this space using the MACE architecture and a shared representation of the charge density, enabling controlled comparisons between different approaches. Finally, we evaluate these models on two instructive test cases: metal-water interfaces, which probe the contrasting electrostatic response of conducting and insulating systems, and charged vacancies in silicon dioxide. Our results highlight the limitations of existing approaches and demonstrate how more expressive self-consistent models are needed to resolve failures.
William J. Baldwin, Ilyes Batatia, Martin Vondrák +2
Engineering Laboratory, University of Cambridge, Trumpington St, Cambridge, UK · University of Bayreuth, Bavarian Center for Battery Technology (BayBatt), Bayreuth, Germany · Fritz-Haber-Institute of the Max-Planck-Society, Berlin, Germany
Foundation machine learning interatomic potentials (MLIPs) have enabled atomistic simulations across broad regions of chemical and materials space, but many remain computationally expensive and lack explicit electrostatics, limiting their use for systems governed by long-range interactions and electrical response. Previously, we introduced Latent Ewald Summation (LES), which learns latent atomic charges and long-range electrostatics from density functional theory (DFT) energy and force labels alone. Here, we use LES to extract electrostatics that are latent in foundation models: energies and forces predicted by a teacher model are used to train a lightweight LES-augmented student MLIP, with optional fine-tuning on additional DFT data. The resulting models reduce computational cost while providing access to Born effective charge tensors, and infrared spectra. We benchmark student models distilled from a broad set of foundation MLIPs, including UMA, MACE, Orb, eSEN, GemNet-OC, PET, and EquiformerV2-based models, against experimental infrared spectra for liquid water, concentrated hydrochloric acid, and the anatase TiO2(101)-water interface. Across these systems, electrostatic response can be extracted from most foundation MLIPs. The benchmark further shows that the underlying DFT level and dataset used to train the teacher model play a larger role than architecture in determining electrostatic and spectroscopic accuracy. For the TiO2-water interface, fine-tuning with a modest amount of higher-level DFT data improves structural and infrared predictions. LES-based distillation therefore provides a practical route for converting foundation MLIPs into efficient, electrically responsive models, while also testing the physical fidelity encoded in foundation models.
Xiaoyu Wang, Bingqing Cheng
Department of Chemistry, UC Berkeley, California 94720, United States · Bakar Institute of Digital Materials for the Planet, UC Berkeley, California 94720, United States · Chemical Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California, 94720, United States