PaperlayerView exact revision on arXiv

Source-linked AI summary

MACE-POLAR-1: A Polarisable Electrostatic Foundation Model for Molecular Chemistry

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

arXiv:2602.19411v1physics.chem-phcs.LG

TL;DR

Local MLIPs struggle with long-range electrostatics, charge transfer, and polarisation despite their success for short-range chemistry. MACE-POLAR-1 extends MACE with polarisable spin-charge updates and global Fukui equilibration, achieving broad accuracy across molecular benchmarks while improving long-range interaction predictions. Its scope remains limited by benchmark and reference-data challenges, including systematic liquid-density errors and difficult transition-metal redox physics.

  • Problem

    Local MLIPs cannot represent long-range electrostatics beyond their cutoff, limiting modelling of charged molecules, extended biomolecular systems, distant charge transfer, and polarisation.

  • Method

    MACE-POLAR-1 extends MACE with a non-self-consistent field-dependent induction mechanism that updates learnable spin-charge densities and equilibrates total charge and spin with Fukui functions.

  • Results

    Across diverse chemical systems, MACE-POLAR-1 shows broad quantitative accuracy, including sub-kcal/mol molecular-crystal errors and substantially improved protein-ligand and redox predictions over local models.

  • Takeaways & Limitations

    Explicit long-range electrostatics supports a transferable molecular foundation model spanning small molecules, protein-ligand complexes, molecular crystals, and solvated transition-metal redox systems.

  • Takeaways & Limitations

    Liquid densities remain systematically 5–10% from experiment, while transition-metal redox accuracy is limited by multireference character, strong correlation, and uncertain hybrid-DFT references.

Abstract

from arXiv · show

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.

I. Introduction

MACE-POLAR-1 addresses the inability of local MLIPs to represent long-range electrostatics by extending MACE with explicit, field-dependent electrostatic interactions. The model combines local geometric features with a learnable spin-charge density and polarisable updates for broader molecular applicability.

  • Motivation: Local MLIPs capture short-range bonding and repulsion but cannot represent interactions beyond their cutoff, limiting charged, ionic, and extended biomolecular applications.This also prevents reliable responses to external fields, distant charge transfer, and polarisation in extended systems.
  • Contribution: MACE-POLAR-1 introduces physics-based long-range electrostatics through a non-self-consistent field update that refines atomic multipoles in response to electrostatic potentials.The architecture retains MACE for short-range interactions while adding field-dependent induction.
  • Energy model: The total energy combines a local MLIP contribution, smeared long-range Coulomb interactions, and a learned non-local correction for residual effects such as dispersion.This decomposition separates short-range chemistry from explicit electrostatics and additional non-local contributions.
  • Representation: The model learns a smooth, coarse-grained spin-charge density whose low-frequency component represents long-range electrostatics, while rapidly varying electronic structure remains in the local energy.Broad Gaussian multipoles are sufficient for the long-range residual density after nuclear and core-electron screening.

3. Local guess to the spin-charge density

The model first predicts a local spin-charge density from MACE features, providing a geometry- and chemistry-dependent initial estimate. This local guess is then prepared for later long-range refinement and global charge constraints.

  • Local prediction: MACE node features are used to predict atom-centred spin-charge multipoles at each local-model layer.The learned weight matrix maps local node features into multipole coefficients.
  • Local prediction: These initial multipoles capture the geometry and chemical dependence of the density that can be described within the local receptive field.They serve as a local baseline before non-local electrostatic information is incorporated.
  • Initial equilibration: Learnable Fukui features equilibrate the monopoles so the spin-charge density satisfies the target total charge and total spin.The equilibration changes the monopole channels, while higher multipoles remain equal to their unequilibrated values.

4. Long-range Polarisable Field Updates and Fukui Equilibration

MACE-POLAR-1 refines its spin-charge density through repeated long-range updates inspired by a self-consistent-field loop. Each update constructs electrostatic features, predicts new multipoles using local and non-local information, and globally equilibrates charge and spin.

  • Update cycle: Each long-range update consists of electrostatic feature construction, a local multipole update, and Fukui equilibration.The sequence is repeated across update iterations to model polarisation and long-range charge transfer.
  • Electrostatic features: The model convolves the spin-charge density with the Coulomb kernel to obtain a spin-resolved electrostatic potential, then projects that potential onto atom-centred Gaussian functions.An applied potential can be included alongside the density-generated potential.
  • Multipole update: Updated spin-charge multipoles are predicted by combining non-local electrostatic features with local MACE node features.The resulting operation incorporates distant information while computing each atom’s new multipoles locally from its surroundings and electrostatic features.
  • Fukui equilibration: After each prediction, learnable Fukui functions globally renormalise the monopoles to enforce total charge and total spin.The Fukui functions can incorporate non-local information through the field update, while higher multipoles remain unchanged by equilibration.

5. Non-local energy

The model evaluates electrostatic energy from the final multipole charge density and adds a learned non-local correction. Separate real-space and reciprocal-space treatments support isolated and periodic systems, while analytic Gaussian representations enable efficient evaluation and force response.

  • Non-local correction: A learned field- and charge-dependent non-local term captures energetic contributions beyond Coulomb interactions.It is intended to represent residual non-local effects such as dispersion.
  • Electrostatic energy: The electrostatic energy is computed from the final multipole charge-density coefficients using Gaussian smearing for smooth short-distance blending.The final density after U long-range updates is used in the energy calculation.
  • Non-periodic systems: For isolated molecules, the Hartree term is evaluated by direct real-space summation over Gaussian multipoles.The implementation includes self-energy and Gaussian-damped Coulomb interactions, with dipolar terms evaluated through an efficient displaced-charge construction.
  • Periodic systems: For periodic systems, the Hartree energy is computed in reciprocal space from the analytically available Fourier series of the Gaussian charge density.The formulation assumes a neutral supercell; charged-defect calculations require correction terms analogous to those used in DFT.
  • Forces and response: Differentiating the electrostatic energy with respect to atomic positions and multipole coefficients supplies electrostatic forces and field-response terms.These derivatives enter both the training loss and molecular-dynamics propagation.

7. Summary of the model architecture

MACE-POLAR-1 combines local MACE interactions with iteratively updated spin-charge densities and Fukui-based charge equilibration, then computes local, Coulomb, and non-local energy contributions. Trained on 100 million OMol25 structures, it is evaluated across thermochemistry, reaction barriers, and related benchmarks against local and electrostatic models.

  • Architecture: The architecture starts from atomic positions and species, producing local MACE features and a local energy contribution for bonding, Pauli repulsion, and short-range electrostatics.
  • Architecture: The model predicts an atom-centred Gaussian spin-charge density, converts it into per-spin electrostatic features, and uses those features to update the density.
  • Architecture: Learnable Fukui features are normalised across atoms and used to equilibrate charges, with this update cycle repeated two times in current models.
  • Energy and forces: The final density supplies Coulomb and non-local energy terms, which are summed with the local energy before forces are obtained by differentiation.
  • Training and evaluation: 100 million OMol25 structures were used for training, while MACE-POLAR-1-M and MACE-POLAR-1-L vary the number of local MACE interaction layers.
  • Benchmark results: On GSCDB138, electrostatic models roughly halve other MLIPs’ ionisation-potential error, improve electron and proton affinities, and perform well on reaction barriers.The benchmark uses WTMAD-2; barrier-category differences are described as modest, with accuracies on the order of a few kcal/mol.

3. Comparison to other pre-trained models on a subset of elements

MACE-POLAR-1 consistently improves electrostatic and non-covalent interaction modelling across protein-ligand, supramolecular, and molecular-crystal benchmarks, while retaining strong accuracy on transition metals, conformers, and liquid-water structure.

  • Non-covalent interactions: MACE-POLAR-1 variants reduce ionic hydrogen-bond errors threefold versus local baselines and recover long-range 1/r Coulombic behaviour in alkali-halide dissociation.Local models flatten at long range because of their hard cutoff distance.
  • Supramolecular complexes: 3.52–4.78 kcal/mol MAEs on S30L improve 40% over MACE-OMOL’s 7.31 kcal/mol, with largest gains for charged host-guest and π-stacked complexes.S30L contains complexes of up to 200 atoms with charge states from −1 to +4.
  • Molecular crystals: 0.46 kcal/mol on X23-DMC is the best lattice-formation-energy accuracy, below 1 kcal/mol and threefold better than MACE-OMOL.This benchmark contains 23 neutral molecular crystals evaluated against diffusion Monte Carlo references.
  • Molecular crystals: 1.06 and 1.22 kcal/mol MAEs on CPOSS209 improve over MACE-OMOL’s 2.73 kcal/mol, although relative lattice-energy errors are similar across models.CPOSS209 contains 209 polymorphs from 20 small drug molecules and precursors.
  • Protein-ligand interactions: 0.37 and 0.47 kcal/mol MAEs on PLF547 improve over MACE-OMOL’s 1.08 kcal/mol for protein-fragment interaction energies within local cutoffs.The benchmark probes subtle intermolecular interactions rather than interactions beyond the local cutoff.
  • Protein-ligand interactions: 3.35–3.68 kcal/mol MAEs on PLA15 outperform MACE-OMOL’s 29.9 kcal/mol, demonstrating substantially better protein-ligand active-site interactions.The models capture interactions from individual hydrogen bonds through the collective electrostatic environment of the binding pocket.
  • Conformers: MLIPs achieve below-0.5 kcal/mol MAEs for most conformer systems, with similar accuracy across test sets and overall performance close to ωB97M hybrid DFT.Adding electrostatics does not materially change conformer accuracy because these tests primarily probe covalent or short-range interactions.

I. Organic liquid densities

MACE-POLAR-1 models were evaluated on organic-liquid densities and aqueous structural properties, with long-range electrostatics supporting charge localisation and improved redox behaviour. Liquid-density accuracy was strongest for MACE-POLAR-1-L, while water structure retained residual deviations from experiment.

  • Liquid-water structure: Water radial distribution functions remained imperfect, with deviations in peak heights and second-shell structure despite close agreement among models.The comparison used O-O, H-H, and O-H pairs from 300 K NPT simulations.
  • Charge localisation: MACE-POLAR-1 models preserved solvated Cl2 near its covalent bond length, whereas non-electrostatic models drove unphysical bond dissociation.The system placed two distant Cl− ions beyond the local message-passing cutoff, testing charge localisation.
  • Redox chemistry in solution: 0.60–0.64 V: MACE-POLAR-1-M and MACE-POLAR-1-L achieved these mean absolute errors across five redox couples, while local models exceeded 10 V on stable MACE-POLAR-1-L trajectories.The authors attribute the remaining ∼0.6 V accuracy challenge to multireference character, strong electron correlation, and subtle spin-state energetics.

K. Effect of external fields

MACE-POLAR-1 is evaluated on external-field response, lanthanide isomerisation, charge localisation, and broader chemical-system benchmarks. The models show physically meaningful response and broad quantitative accuracy, while revealing limitations for low-polarisability molecules and some liquid properties.

  • External-field response: MACE-POLAR-1 models predict external-field dipoles and polarisabilities despite training exclusively on ground-state energies and forces.The response properties are computed by finite differences of total energy under a uniform applied field.
  • External-field response: Dipole predictions achieve MAE values of 0.197 D and 0.174 D for MACE-POLAR-1-M and MACE-POLAR-1-L, respectively.These errors represent approximately threefold degradation relative to the reference while retaining quantitative accuracy for most molecules.
  • External-field response: HR46 polarisability errors are 1.76 Å3 and 1.85 Å3 for MACE-POLAR-1-M and MACE-POLAR-1-L, versus 0.187 Å3 for wB97M-V.The models show an order-of-magnitude degradation on this second-order response benchmark.
  • External-field response: The models show systematic errors for low-polarisability molecules, with predictions clustered near zero, indicating missing higher-order correlation effects.The reported limitation is attributed to the many-body polarisation tensor not being fully encoded in the short-range message-passing architecture.
  • Lanthanide isomerisation: Lanthanide MLIPs match reference isomer ordering within approximately 1 kcal/mol on average, while ORBMOL has several outliers.The tightest La and Eu splittings cause visible offsets, and local MACE-OMOL failed to converge for several topologies.
  • Broader benchmarks: The study reports quantitative accuracy across small organics, protein–ligand complexes, redox ions, transition-metal complexes, and molecular crystals, while liquid densities remain 5–10% off experiment.The authors identify reference-functional, architecture, and simulation-protocol effects as possible contributors to the liquid-density discrepancy.

V. Supplementary Information

The supplementary information defines benchmark aggregation, exclusion rules, model hyperparameters, and benchmark subset composition. WTMAD-2 provides a scale-balanced summary across heterogeneous energetic datasets, while tables document exclusions and implementation settings.

  • Benchmark aggregation: WTMAD-2 rescales each benchmark subset by the inverse mean absolute reference magnitude to balance datasets with different energetic scales.The aggregate score is reported in kcal/mol and enables one-number comparison across broad benchmark suites.
  • Exclusion rules: Summary aggregation excludes single-atom configurations and outliers with |∆E| > 100 kcal/mol, while per-dataset distributions retain the full outlier range.Table V reports model-specific exclusions for thermochemistry and reaction-barrier subsets.
  • Hyperparameters: The hyperparameters include angular-resolution, correlation, channel-size, radial-basis, cutoff, charge-smearing, field-feature, and non-local-update settings.The supplementary definitions specify roles for max_ell, correlation, max_L, channel counts, num_radial_basis, r_max, smearing widths, field features, and num_update.
  • Hyperparameters: spin_charge_density_l_max sets the maximum multipole order predicted for spin and charge density, while num_update sets the number of non-local field-update iterations.These settings control the angular complexity and iterative depth of the electrostatic update.
  • Benchmark composition: The benchmark suite covers thermochemistry, reaction barriers, and disjoint WTMAD-2 subset rows for bond energies, reaction energies, thermochemistry, and barriers.The listed subsets include ionisation potentials, electron affinities, proton affinities, bond dissociation, hydrogen-atom transfer, barrier heights, proton transfer, and reaction energies.

E. Real Space Electrostatic Energy Computation

The real-space electrostatic computation represents charge and multipole densities with Gaussian basis functions and evaluates their interactions through potentials and contractions. Dipole terms are implemented using displaced Gaussian charges, including for overlapping atomic basis functions.

  • Electrostatic energy: The open-boundary electrostatic energy is expressed as a sum over charge and multipole interactions.The same interaction formulation also supplies projected potential features entering the model's non-local update.
  • Gaussian basis: For monopole interactions, the electrostatic potential and energy use Gaussian source and receiver basis widths, including a special l = l′ = 0 form.The electrostatic energy is recovered from the same contraction when source and receiver widths are identical.
  • Dipole interactions: Dipole interactions approximate an l = 1 Gaussian orbital as a linear combination of slightly displaced l = 0 Gaussians.The construction remains valid when atomic Gaussian functions overlap, not only for distant interacting atoms.
  • Dipole interactions: Using a displacement a = 0.02 Å gives sub-meV total-energy convergence for the orbital approximation, whereas a = 0.1 Å is used for projected-potential feature evaluation.The larger feature-evaluation displacement improves numerical stability while preserving long-range feature-map accuracy.

2. Self-Interaction Energy

The self-interaction and reciprocal-space procedures correct Gaussian multipole self-energy and efficiently represent periodic charge densities. Adaptive Fourier cutoffs resolve the most localised basis functions while exploiting conjugate symmetry.

  • Self-interaction correction: The implementation adds the self-interaction energy of each Gaussian charge distribution to the electrostatic energy.The self-energy is defined for Gaussian multipoles and contributes a total self-interaction correction through the overlap matrix.
  • Reciprocal-space representation: The reciprocal lattice vectors are constructed from L∗ = 2π(LT)^−1, and the periodic k-space grid retains vectors within an adaptive cutoff.The grid uses integer reciprocal-lattice indices and stores only a half-grid because scalar density and potential Fourier coefficients obey conjugate symmetry.
  • Reciprocal-space cutoff: The adaptive cutoff is determined from κ, the minimum Gaussian smearing width σmin, and the maximum angular momentum lmax.This heuristic ensures that the Fourier transform of the most localised basis function is adequately resolved; κ is typically 1.5.
  • Fourier construction: Fourier transforms of Gaussian spherical-harmonic basis functions separate angular and radial components using spherical harmonics and spherical Bessel functions.The resulting radial integrals are used to construct the Fourier series of the full charge density.
  • Fourier construction: The Fourier-series implementation computes real and imaginary parts separately, with spherical-harmonic phases producing real parts for even l and imaginary parts for odd l.This phase structure determines the trigonometric combinations used for the transformed density components.

3. Coulomb Operator in Reciprocal Space

The reciprocal-space Coulomb operator converts charge density into electrostatic potential and energy using the bare Coulomb kernel. The implementation also applies finite-size corrections selectively for isolated, periodic, and slab systems.

  • The Coulomb potential is obtained from the charge density by applying the bare Coulomb kernel in reciprocal space.
  • Electrostatic energy is computed as the reciprocal-space inner product of charge density and potential.
  • Discretisation over the k-point grid expresses the energy through real and imaginary components of the Fourier-transformed density.
  • The k = 0 term is excluded for neutral systems, while sufficiently large cells yield exponential convergence with kcutoff.
  • Projected potential features use the same reciprocal-space representation and receiver-basis projection.
  • Makov–Payne–Dabo, slab-dipole, and related finite-size fields are added to reciprocal-space features before combination.
  • Corrections depend on periodicity: all three apply to non-periodic systems, none to fully periodic systems, and a slab-specific dipole correction to slab geometries.

G. Fukui functions and conceptual DFT

The Fukui-function formulation links atomic charge response to conceptual DFT and uses learnable softness weights to enforce global spin-resolved charge constraints. The surrounding benchmark suite spans non-covalent, thermochemical, reaction, transition-metal, liquid, and crystal properties.

  • Conceptual DFT: The Fukui function describes how electron density responds locally to changes in total electron number at fixed external potential.
  • Charge and spin constraints: Total charge and spin are represented through spin-resolved charges, with Q↑↓ constructed from the total charge Q and total spin S.
  • Conceptual DFT: In the model, Fukui functions are neural-network predictions that can depend on local geometry and non-local electrostatic environment.
  • Charge and spin constraints: Charge equilibration redistributes deviations from target spin-resolved charges according to Fukui coefficients.
  • Benchmark coverage: The benchmark suite includes non-covalent, reaction-barrier, thermochemical, transition-metal, and molecular-crystal datasets.
  • Benchmark coverage: The datasets use reference methods including CCSD(T)/CBS, CCSD(T), DLPNO-CCSD(T), and ωB97M-V with a 1-body CCSD(T) correction.

Molecular Crystals

The molecular-crystal evaluation includes CPOSS209 lattice-energy data and liquid-phase density and dynamics benchmarks. These tests cover polymorph energetics, organic liquids, water, and solvated iron/chloride systems.

  • Molecular crystals: CPOSS209 contains 209 molecular-crystal lattice-energy cases referenced to ωB97M-V plus a 1-body CCSD(T) correction.
  • Liquids and water: The liquid benchmarks include water density and radial distribution functions, experimental liquid densities, and organic-liquid densities.
  • Solvation and charge localisation: Solvated Fe/Cl chloride trajectories and water-cluster charge-localisation tests probe condensed-phase and charge-transfer behavior.

Solvation / Redox

The solvation and redox benchmarks test charge localisation, hydrated transition-metal ionisation, aqueous redox potentials, and lanthanide complex isomers.

  • Charge transfer: Water-cluster dissociation benchmarks evaluate charge transfer and fragment-charge localisation during dissociation.
  • Solvated transition metals: Hydrated transition-metal benchmarks cover ionisation in water for 14 systems and solvated Fe/Co chloride trajectories.
  • Redox: The redox evaluation includes redox potentials in solution for solvated Fe/Co chloride systems.
  • Lanthanides: Lanthanide evaluation uses 18 lanthanide complex isomers with r2SCAN-3c references.
Loading 2602.19411v1…