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

Scalable long-range machine learning interatomic potentials with electrical response for atomistic material modeling
Oct 16, 2026Yicheng Chen, Xiaoyun Yu, Hyejeong Song +5Machine Learning Interatomic Potentials Computational Materials Science
Department of Materials Science and Engineering, National University of Singapore, Singapore · Department of Physics, National University of Singapore, Singapore
Polarizable atomic multipoles for learning long-range electrostatics
May 7, 2026Yoonjae Park, Dongjin Kim, Daniel S. King +5Machine Learning Interatomic Potentials Materials Property Prediction
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
Distilling latent electrostatics from foundation machine learning interatomic potentials
Jun 12, 2026Xiaoyu Wang, Bingqing ChengInfrared Spectroscopy Machine Learning Interatomic Potentials
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

Oct 16, 2026J/K move · Enter open · S save

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

Yicheng Chen, Xiaoyun Yu, Hyejeong Song, Xintian Wang+4

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

Open paper AskPDF arXiv

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.

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