Design Space of Self--Consistent Electrostatic Machine Learning Interatomic Potentials

William J. Baldwin; Ilyes Batatia; Martin Vondrák; Johannes T. Margraf; Gábor Csányi

physics.chem-phMar 16, 2026

Design Space of Self--Consistent Electrostatic Machine Learning Interatomic Potentials

Authors: William J. Baldwin, Ilyes Batatia, Martin Vondrák, Johannes T. Margraf, Gábor Csányi

Organizations: 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

Abstract

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.

Figures & tables

Figure 1: Required properties of the density expansion and AIM Fukui functions. (a): The potential at $r_{1}$ can be calculated by replacing all the nuclei and electrons outside of some radius $R$ by the coarse-grained density. The plots of $\tilde{n}$ and $\rho$ show charge density, with the repeating pattern coming from negative electrons and positive nuclei. The shaded regions indicate some excess charge due to, for example, defects. The coarse-grained density smooths out all the detail of the atoms, leaving only the net charge or dipole around the defects. (b): The AIM Fukui function $\tilde{f}_{k}$ must give rise to the same far-field potential as the basis function $\phi_{k}$ . A spherically symmetric basis function associated with $p_{k}$ is shown on the left, and two possible AIM Fukui functions for this variable are shown on the right. The one at the top has more structure, but will still give rise to the same electric potential outside the dashed line as the $\phi_{k}$ so long as their multipole moments agree at all orders. The one at the bottom has a significant dipole moment, meaning the potential outside the dashed line will not match that of $\phi_{k}$ .

Figure 2: (a): Derivative discontinuity of the energy of an isolated system, (b) discontinuity of the total charge as a function of chemical potential for an isolated system. (c): Infinite susceptibility of one variable, $p_{k}$ , as one would find in a classical perfect conductor, appears as zero curvature in the energy. (d): Similarly, zero susceptibility appears as zero response to a change in the potential corresponding to that variable.

Category	Model	Description of charge density	Locality	Functional form
Energy functional	QEq	Partial charges	Local	Quadratic
	4GNN [ 35 ]	Partial charges	Local	Quadratic
	CENT [ 58 ]	Partial charges	Local	Quadratic
	kQEq [ 59 ]	Partial charges	Local	Quadratic
	BpopNN [ 38 ]	Spin resolved partial charges	Semilocal	GAP+empirical baseline
	ACE‐Q [ 96 ]	Partial charges	Semilocal	Quadratic

Table 1: Comparison of methods for incorporating a model of the charge density into machine learning force fields. This table only presents models which are strictly self-consistent, and not message passing networks where the contributions with iterations get smaller but in which there is no ‘convergence criteria’. The column ‘Locality’ describes whether the fixed-point function or the energy functional is local or semilocal, in the sense discussed in section II.4 . GAP refers to the Gaussian Approximation Potential MLIP [ 2 ]

Figure 4: Overview of the implementation of each approach.

Figure 5: Training Methods for our self-consistent models.

Figure 6: Representative configurations from the metal–water dataset. (a) Aluminium slab where $z$ is the out-of-plane direction, with an applied field. (b) Liquid water slab also with an applied field. (c) Water cluster. (d) Water cluster containing a hydronium ion, with a total charge of $+1$ .

Configuration	Number
Aluminium slab with field	308
Water slab with field	438
Water slab	644
Water cluster (Hydronium)	600
Water cluster (Neutral)	600
Total	2590

Table 2: Water–metal dataset composition.

Hyperparameter	Value
Number of chemical channels	128
Maximum equivariance order $L$	1
Single layer cutoff radius (Å)	6

Table 3: Hyperparameters for the local component of all models.

Property	Loss weight $w$
Energy	100
Force	100
Atomic multipoles ( $l=1$ )	100
Total dipole	1 000 000
Total charge	10 000

Table 4: Loss weights for all water–metal fits.

Figure 7: Training characteristics of the energy functional (a) and fixed-point models (b) for the metal–water slab example. The difference between direct and direct (converged) is explained in the text.

Figure 8: Training a fixed-point model with a truncated SC-cycle. 9 different runs are shown with different choices of the number of SC-steps to perform and the mixing parameter. The legend shows this as “number-of-steps, mixing-parameter”.

Figure 9: Comparing the performance of charge-aware models, with and without including multipoles in the loss.

Figure 10: Comparing the energy functional and fixed-point models using different architectures.

Figure 11: Physicality of the trained models. (a): fragmentation of small water clusters. The top panel shows separation of a cluster including a hydronium ion, whereas the lower plot shows fragmentation of a neutral cluster. (b): finite field response of a combined metal and water system.

Figure 12: Quantifying the smoothness of the landscape of the self-consistent functions of each model; details are in the main text. In (a) singular values of the Jacobian of the update function of the water-metal slab fixed-point model. In (b) the eigenvalues of the Hessian of the water-metal energy model. In (c): each panel corresponds to a single eigenfunction of the Hessian of $G_{\text{ML}}[\mathbf{p}]$ . The eigenfunctions are drawn by visualising atomic dipoles as arrows, and colouring the atoms according to the charge. Since large eigenvalues in the Hessian of the energy model correspond to small singular values of fixed-point model’s update function, we have plotted the eigenvalues in increasing order and the singular values in decreasing order.

Figure 13: Comparing the effect of including atomic dipoles and vector features of the electric potential.

Figure 14: Oxygen vacancies in $\alpha$ -quartz silicon dioxide with the atoms coloured according to the change in Hirshfeld charge relative to atoms in a unit cell of the material.

Figure 15: (a): Validation errors for the different models on the $\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$ training set. The global embedding model is a local model and therefore does not have a prediction for the atomic multipoles. (b): Si–Si bond length after relaxing each of the four defect types in a supercell. (c): Supercell containing one neutral oxygen vacancy and one negative oxygen vacancy, to illustrate the double defect relaxation test. Atoms are coloured according to the difference in Hirshfeld charge relative to pristine $\alpha$ -quartz.

Model	Total	Resulting	Expected?
Model	Charge	Defects
Global Embedding	$-2$	relaxation fails	✗
	$-1$	$0,\,-1$	✓
	$0$	$0,\,\ \ 0$	✓
	$1$	$1,\,\ \ 1$	✗
	$2$	relaxation fails	✗

Table 5: Relaxing pairs of defects in different charge states. In all cases, the initial geometry included one copy of the neutral defect and one copy of the negative defect. The ‘expected’ entry is omitted when models relax to the

(0,-1)

structure as discussed in the text and table 6 .

Model	Initial	Resulting	Relative
Model	Defects	Defects	Energy (eV)
DFT ( $N_{\uparrow}-N_{\downarrow}=0$ )	$(0,-1)$	$(0,-1)$	$-0.363$
DFT ( $N_{\uparrow}-N_{\downarrow}=0$ )	$(-1,-1)$	$(-1,-1)$	$0.000$
DFT ( $N_{\uparrow}-N_{\downarrow}=2$ )	$(0,-1)$	$(0,-1)$	$0.107$
DFT ( $N_{\uparrow}-N_{\downarrow}=2$ )	$(-1,-1)$	$(-1,-1)$	$0.000$
Fixed Point	$(0,-1)$	$(0,-1)$	$-0.084$

Table 6: Defect pairs with -2 total charge do not always relax to the expected structure. Two defects are initialised in a supercell where either (i) both defects have the geometry of the isolated negative vacancy, denoted (-1, -1), or (ii) one defect has the geometry of the negative vacancy while the other has the geometry of the neutral vacancy, denoted (0, -1). For each model in the table, a defect pair with each of these initializations is relaxed to the nearest local minima with a total charge of -2 and the final energy is reported relative to the (-1,-1) case.

N_{\uparrow}-N_{\downarrow}

refers to the constraint on total projected spin.

Appendix figures & tables20 assets

Supplementary material from the paper’s appendix.

Appendix

Figure 16: Two localised, spherically symmetric charge distributions

\widetilde{f}

and

\phi

, together with their electrostatic potentials. If the two distributions have the same total charge, their far-field potentials agree and the difference

\zeta

is confined to the region in which their detailed shapes differ. For nonspherical functions, the same statement applies when the multipole moments agree at all orders.

Fit Name	Training Schedule	Architecture	Atomic multipole order	Remove from loss
Energy Models
energy_A	implicit	quadratic	1	–
energy_B	implicit	nonlinear	1	–
energy_C	implicit	manybody	1	–
energy_D	implicit	quadratic	1	all but $E$ , $F$
energy_E	implicit	quadratic	0	$q_{lm}$ weight $=1$

Appendix

Table S.I: Summary of all trained models. The ‘Architecture’ Column describes the complexity of either the fixed point update or the energy functional layer.

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$Q_{\text{tot}}$	$P$	Size	Rate
Stage 1	Implicit Differentiation	0–200	100	100	100	0	1,000,000	1	0.01

Appendix

Table S.II: Energy Functional Models Implicit Training Schedule for Al/

\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$Q_{\text{tot}}$	$P$	Size	Rate
Stage 1	Direct	0–200	100	100	100	10,000	1,000,000	2	0.01
Stage 2	Implicit Differentiation	200–230	1000	100	100	10,000	1,000,000	2	0.001

Appendix

Table S.III: Fixed Point Models Direct + Implicit Training Schedule for Al/

\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$Q_{\text{tot}}$	$P$	Size	Rate
Warm-up (Stage 0)	Direct	0–10	100	100	100	10,000	1,000,000	2	0.01
Stage 1	Shortcut SCF	10–200	100					2	0.01
Stage 2	Unrolled SCF	200–230	1000					2	0.001

Appendix

Table S.IV: Fixed Point Models Shortcut SCF Training Schedule for Al/

\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$Q_{\text{tot}}$	$P$	Size	Rate
Stage 1	–	0–200	100	100	100	10,000	1,000,000	2	0.01
Stage 2	–	200–230	1000	100	100	10,000	1,000,000	2	0.001

Appendix

Table S.V: LocalCharges Training Schedule for Al/

\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$\nabla_{q}E$	$P$	Size	Rate
Stage 1	–	0–200	100	100	100	1,000	1,000,000	1	0.01

Appendix

Table S.VI: Training Schedule for Energy Functional models when testing training dynamics on the metal water system (main text figure 7). The loss on the gradient of the functional (

\nabla_{q}E

) is omitted when training with implicit differentiation since it is zero by assumption.

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$Q_{\text{tot}}$	$P$	Size	Rate
Warn-up	direct	0–10	100	100	100	10,000	1,000,000	4	0.01
Stage 1	–	10–200	100	100	100	10,000	1,000,000	4	0.01

Appendix

Table S.VII: Training Schedule for Fixed Point models when testing training dynamics (main text figure 7).

Figure S1: Difference in force and dipole error for fixed-point models when trained with either shortcut-SCF or direct+implicit training.

Figure S2: Difference in charge localization for fixed-point models when trained with either shortcut-SCF or direct+implicit training. One can see that shortcut-SCF training is marginally worse in both cases.

Figure S3: Difference in finite field response for fixed-point models when trained with either shortcut-SCF or direct+implicit training. The legend indicates the type of update function and the training schedule. One can see that the shortcut-SCF trained models seem to have more exaggerated differences in charge density, and slightly worse screening.

Figure S4: Training a fixed-point model with a truncated SC-cycle. The legend indicates the number of SC-steps, and the mixing parameter. For example ‘8, 0.5’ means do exactly 8 SC steps with

\lambda=0.5

Figure S5: fixed-point model accuracy as a function of update function complexity and the number of electrostatic features. The column ‘No field features’ is the prediction of the MACE-LocalCharges model since this is analogous to the purely local part of the fixed-point model.

Description	$\sigma$	$l$
one scalar	2.0	0
two scalar	1.5, 2.5	0
two scalar + vector	1.5, 2.5	0, 1

Appendix

Table S.VIII: Choices of

\sigma

l

used in Figure S5

Figure S6: Histogram of Si–Si distances in the training set for the silicon dioxide example.

Defect	Si-Si Distance (Å)
$V_{\text{O}}$	2.44
$V_{\text{m}}$	2.53
$E_{\text{di}}$	3.02
$E_{\text{puck}}$	4.53

Appendix

Table S.IX: Si-Si distances for various defects after DFT relaxations with our DFT settings in a 242 atom supercell.

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$Q_{\text{tot}}$	$\mu$	Size	Rate
Stage 1	–	0–250	100	100	0	0	0	4	0.005
Stage 2	–	250–300	100	100	0	0	0	4	0.002

Appendix

Table S.X: Global Embedding MACE Training Schedule for

\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$Q_{\text{tot}}$	$\mu$	Size	Rate
Stage 1	–	0–200	10	1000	1000	0	0	4	0.005
Stage 2	–	200–250	1000	100	1000	0	0	4	0.002

Appendix

Table S.XI: MACE-QEq Training Schedule for

\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$Q_{\text{tot}}$	$\mu$	Size	Rate
Warm-up (Stage 0)	Direct	0–10	100	100	1000	10,000	1000	4	0.01
Stage 1	Shortcut SCF	10–250	100					4	0.005
Stage 2	Unrolled SCF	250–300	1000					4	0.002

Appendix

Table S.XII: Fixed Point Models Training Schedule for

\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}

	Method	Epochs	Loss weights					Batch	Learning
			$E$	$F$	$q_{ilm}$	$Q_{\text{tot}}$	$\mu$	Size	Rate
Stage 1	Unrolled SCF	0–400	100	100	1000	0	0	4	0.01
Stage 2	Unrolled SCF	400–500	1000	100	1000	0	0	4	0.01

Appendix

Table S.XIII: Energy Functional Models Training Schedule for

\text{SiO}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}

Density Functional Theory

See all 25 papers

Oct 16, 2026

Scalable long-range machine learning interatomic potentials with electrical response for atomistic material modeling

Accurate atomistic simulations of materials under electric fields require long-range electrostatics and electrical-response modeling at length scales beyond the reach of conventional first-principles methods. Yet these capabilities remain difficult to reconcile with the scalability of machine-learning interatomic potentials (MLIPs). Strictly local models enable efficient spatial parallelization but neglect long-range interactions, whereas existing electrostatic extensions introduce global dependencies that hinder distributed computation. Here, we overcome this limitation by analytically separating the long-range electrostatic forces and Born effective charges (BECs) into a global electrostatics calculation and strictly local derivatives of learned latent charges. This formulation preserves the spatial decomposability of equivariant MLIPs while enabling distributed evaluation of energies, forces, stresses, and BECs across different boundary conditions. Implemented as DPA4C-LR on the nvalchemi platform, the framework scales to 5.2 million atoms on four GPUs at approximately twice the computational cost of its short-range backbone. Trained on energies and forces, DPA4C-LR accurately infers BECs and reproduces the field-dependent infrared spectra of water. We further demonstrate its capabilities by recovering the characteristic slim polarization hysteresis of relaxor \ce{Pb(Mg_{1/3}Nb_{2/3})O3}, which is obscured by finite-size effects in smaller simulations; predicting a room-temperature field-driven ionic conductivity of

6.8 \pm 0.7

~mS/cm in amorphous \ce{Li_{0.9}NbO_{0.9}Cl_{4.1}}; and modeling interfacial thermal transport across Si/diamond interfaces with slab-corrected electrostatics. By combining long-range physics, electrical-response inference, and distributed computation, this framework extends predictive atomistic modeling toward device-relevant length scales.

Yicheng Chen, Xiaoyun Yu, Hyejeong Song +5

Department of Materials Science and Engineering, National University of Singapore, Singapore · Department of Physics, National University of Singapore, Singapore

Machine Learning Interatomic Potentials Computational Materials Science Materials Property Prediction High-Performance Computing

Ask PDF

May 7, 2026

Polarizable atomic multipoles for learning long-range electrostatics

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 MAPbI

_3

perovskite, as well as the essential features of the surface-specific vibrational sum-frequency generation spectrum at the water-air interface. In ferroelectric HfO

_2

, 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 Materials Property Prediction

Ask PDF

Jun 12, 2026

Distilling latent electrostatics from foundation machine learning interatomic potentials

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

Infrared Spectroscopy Machine Learning Interatomic Potentials Materials Property Prediction Knowledge Distillation

Ask PDF

Abstract

Density Functional Theory

Explore similar work

Scalable long-range machine learning interatomic potentials with electrical response for atomistic material modeling

Polarizable atomic multipoles for learning long-range electrostatics

Distilling latent electrostatics from foundation machine learning interatomic potentials