In this paper, we provide a comprehensive exploration of the LES framework, detailing its theoretical foundation, possible extensions, and application to a range of test systems. Importantly, we show that, when limited to a single charge channel, the LES charges can be interpreted as physical partial charges. In ref. 24, LES was compared to other LR methods such as LODE8,9 and density-based long-range descriptor21 that do not explicitly learn charges. Here, we further compare LES to existing methods that incorporate long-range interactions via explicit charge learning and show that LES achieves superior performance.

We first briefly recap LES, and then make an explicit connection between LES and physical charges. Finally, we briefly demonstrate how different global charge states can be encoded in the LES framework.

The total potential energy of a system with N atoms is split into short-range (SR) and long-range (LR) components, . The short-range energy is the sum of atomic energies, each depending on local B features of atom i. The B features can be local atomic environment descriptors such as ACE,

or learned features in message passing neural networks (MPNNs). For the long-range part, a multilayer perceptron with parameters ϕ maps the invariant features of each atom i to a hidden variable:

In general, q can be multi-dimensional to represent the generalized long-range interactions. When q is restricted to be one-dimensional, it can be interpreted as the atomic charge as we discuss later.

Suppose that the potential-generating field by a single particle with unity latent variable is proportional to u(r) = ∣r∣, with p being a fixed exponent. Following the standard range-separation formalism, one can express short-range and long-range interactions by multiplying the interaction by a convergence function φ(r) with φ(0) = 1 decreasing rapidly to zero as r increases:

Both and are short-ranged in nature and can be described by the short-ranged MLIP based on the local features.

For p = 1, which corresponds to electrostatics, one choice for the convergence function can be expressed as the complementary error function . For isolated systems without periodic boundary conditions, one can compute the term directly in the real space based on enumerating pairwise distances between atoms. For periodic systems, the corresponding long-range electrostatics can be computed in the reciprocal space as

where the structure factor S(k) of the hidden variable is defined as

The omission of the k = 0 term in Eq. (3) means the tinfoil boundary condition is applied. The detailed derivations and the case for p = 6 which corresponds to London dispersions can be found in the Supplementary Information.

When training the MLIP, the total potential energy E, interatomic forces F = -∂E/∂r, and sometimes virial stress are fitted to the reference values from the dataset. In LES, unlike methods that explicitly learn partial charges, the hidden variables q are hypothesized to represent flexible atomic charges when the physical electrostatic constant 1/4πϵ is included. In particular, when LES is limited to a single charge channel, we find that the charge used to compute the long-range energy in Eq. (3) is physically meaningful and can be used to predict physical observables such as the dipole moment of gas phase molecules and Born effective charges. However, it is noted that because the structure factor is squared in Eq. (3), the predicted charges do not distinguish the charge parity, as the total energy stays the same if all signs of the charges are flipped. In practice, it is easy to unflip the signs of atomic charges based on the known electronegativity of elements.

We note the success of the LES method in predicting charge locally while computing energy globally. This choice reflects the generally nearsightedness of electron matter. Additionally, as LES learns the charges via the energy and forces, its learning is flexible to arbitrary charge distributions (e.g., different oxidation states) as long as they have an impact on the energy in the training set. Indeed, the LES approach proves appropriate for a wide range of systems such as electrolyte/electrodes, charged molecules, and doped surfaces, as we will show in the examples. However, it is important to note that while the local charge assumption works well empirically, it lacks theoretical guarantees and may encounter limitations in specific edge cases, such as systems involving long-range charge transfer.

Empirically, in the examples below and in previous work we have found it unnecessary to explicitly enforce charge neutrality or fixed total charge state in the training process of LES. In practice, we have found that the sum of q is usually close to the total charges for both neutral and charged systems without enforcing neutrality. Additionally, any residual difference is treated as a uniform background charge, which does not affect the total energy as the k = 0 term is omitted in the reciprocal space computation of electrostatic interactions in Eq. (3). In all the examples we have tested, we did not observe any loss of accuracy or artifacts due to the lack of charge equilibration. In contrast, for the ML models that explicitly learn charges such as 3G-HDNNP, the lack of charge equilibration may result in dramatically larger errors, and sometimes pathological behaviors were observed for systems involving charge transfer and change of charge states.

In a standard MLIP, the atomic features B depend on the chemical elements and the coordinates of the atoms surrounding atom i, and are agnostic to the charge or oxidation state. This means that two systems with identical atomic positions but different net charges Q will have degenerate features. Although this degeneracy does not affect the training or prediction for systems with a fixed net charge, it can cause problems when handling systems with varying charge states simultaneously. To resolve this, in training sets containing multiple net charges (only one of the examples below, ) we concatenate the total charge Q of the system with the local atomic features B, B ⊕ Q, and use this combined feature as the input for predicting short-range atomic energies and local hidden variables. Note that this global charge embedding scheme can have limitations, for example, when dealing with charged systems with varying sizes.

As an initial test, a gas of point charges was constructed. As shown in Fig. 1a, each configuration consists of 128 atoms, with 64 carrying a positive charge of +1e and the remaining 64 carrying a negative charge of -1e. The atoms interact through the Coulomb potential and the repulsive component of a Lennard-Jones potential. This benchmark aims to evaluate the learning efficiency of the LES framework and assess whether the correct atomic charges can be accurately learned. Unlike in density functional theory (DFT), where the precise values of partial charges depend on the chosen definition, the charges in this system are unambiguously defined.

For the short-range component, we employed CACE with different cutoff distances of r = 4 Å, 5 Å, and 5.5 Å. For the long-range interactions, we used a one-dimensional q with σ = 1 Å in the Ewald summation, without enforcing a net charge constraint. Figure 1b presents the parity plot of the CACE-LR model with r = 4 Å, comparing the true and predicted charges (after unflipping the charge parity) for various numbers of training samples. Remarkably, even with just 10 training configurations, the predicted charges are nearly exact.

Figure 1c illustrates the learning curves for the mean absolute errors (MAEs) in energy, forces, and charges, using short-range (SR) and long-range (LR) models with different cutoffs. The SR models exhibit slow learning and significant errors for this dataset, with performance improving as r increases. In contrast, the LR models achieve errors more than an order of magnitude lower, with learning efficiency improving as r decreases. This example highlights that, unlike the typical behavior of SR MLIPs, long-range potentials achieve more efficient learning with appropriately small r values.

We constructed a dataset of potassium fluoride (KF) aqueous solutions with concentrations ranging from 0 to approximately 2 mol/L. The dataset includes both bulk electrolyte solution configurations and electrolyte-vapor interfaces, as illustrated in Fig. 2. The reference energies and forces were computed using the flexible SPC/Fw water model (with oxygen carrying a charge of -0.8476e and hydrogen carrying a charge of +0.4238e), alongside ions with fixed charges (K: +1e, F: -1e). It is worth noting, however, scaled charge ion models are typically better at capturing implicit effects of liquid-phase polarization and modeling electrolytes, although this example aims to demonstrate the learning ability of MLIPs rather than to accurately model electrolytes. This electrolyte dataset is significantly more challenging than the random charge example, as it involves multiple species with distinct atomic charges. Additionally, water acts as a dielectric medium, and the presence of interfaces introduces diverse screening effects that vary with depth from the surface.

Figure 2 b shows that the CACE-LR model with r = 4.5 Å is able to recover the true charges after a couple of hundred training samples. Figure 2c shows the learning curves for the MAEs on forces and charges, and the MAEs on energies are all pretty small for all models (< 0.3 meV/atom for ⪆ 100 samples). While a larger cutoff or a message passing layer (MP1) improves the SR model, the LR model with a smaller cutoff r = 4.5 Å achieves better learning efficiency. Adding a message-passing layer to the LR model has little effect in this case. See below in the Methods section, we also show the learning curves from just the bulk or just the interfacial configurations. This electrolyte example shows that the LR model is able to learn the charges and energetics of systems involving different species and a dielectric medium that screens electrostatics.

We revisit an example from a molecular dimer dataset used to benchmark LODE and LES. This example consists of the binding curve between two charged molecules of (shown in Fig. 3a). The training set of this example is tiny: it consists of 10 configurations of the dimer pair, with the internal coordinates of the molecules frozen and only dimer separation distances varying between ~5 Å and 12 Å. The test set includes 3 configurations with separations between approximately 12 Å and 15 Å. The dataset includes energy and force information calculated using the HSE06 hybrid DFT with a many-body dispersion correction.

For the CACE-LR model, here we use a one-dimensional q, whereas the original LES paper used a four-dimensional hidden variable, and the model test errors are comparable. Figure 3b compares the predicted forces and dispersion curves for the LR and SR models. The SR model has one message-passing layer, but as the two molecules can have a distance beyond the cutoff of r = 5 Å, the message-passing scheme does not help. Figure 3c shows the predicted charge distribution. The total predicted charges on molecules are +0.83e/ -1.08e, and +1.01 e/ -1.01 e after removing the mean charge of each atom . The reason why the mean charge deviates from zero is due to the tiny training sizes. Nevertheless, the mean-adjusted charges are very close to the ground truth of +1 e/ -1 e molecular charges, despite the fact that the MLIP training is agnostic about these charge states. Even though the atomic charges are not quantitative due to the minimal training set, the learned charges are broadly consistent with chemical intuitions: The two under-coordinated oxygen atoms in have the same strong negative charge, while the rest of the molecule is positively charged. The undercoordinated carbon in has a positive charge, while the other atoms have smaller positive charges.

Since atomic charges in quantum mechanics are not well-defined quantities, a key question is whether the LES charges can be used to predict physical observables such as dipole and quadrupole moments. To answer this question, we turn to the SPICE dataset, which contains DFT dipole and quadrupole moments as well as minimal basis iterative stockholder (MBIS) charges for a wide array of drug-like molecules. Specifically, we fit CACE-LR on a dataset of polar dipeptides, just by learning from the energy and forces. Then we determine whether LES is able to infer the DFT dipole and quadrupole moments on a holdout test set of unseen polar dipeptides (illustrated in Fig. 4a). We compute the predicted LES dipole via and quadrupole via where q are the charges predicted by LES and r are the positions of atoms i. To make the comparison translationally invariant, we additionally subtract the trace from the calculated and DFT quadrupole moments ().

Figure 4b compares the charges predicted by LES to the MBIS charges from SPICE. As is seen, the charges predicted by LES correlate well with the MBIS charges, and agree with the usual ordering of electronegativities (O > N > C > H). However, we note that such agreement can only be qualitative (R = 0.87, MAE = 0.24). The reason behind this is that there is no rigorous definition of atomic charge. To show this, we also compare between different definitions of DFT charges (MBIS, Mulliken charges, and Hirshfeld charges), as illustrated below in the Methods section. Indeed, the extent of disagreement between the LES charges and any definition of these DFT charges is similar to that between different definitions of DFT charges.

To evaluate the quality of the LES charges quantitatively, Fig. 4c compares the dipole moments (a well-defined experimental observable) derived from LES to that from DFT. Remarkably, we find that the derived dipoles from the LES charges are in excellent agreement with those from DFT (R = 0.991), even though the LR model is not trained explicitly on any charge or dipole information. In absolute terms, the LES mean absolute error (MAE) for dipole moments is 0.089 e-Å, comparable to the 0.063 e-Å MAE of MBIS charges derived directly from DFT densities. Figure 4d compares the calculated quadrupole moments to those of DFT. Again, we see good agreement of the LES quadrupoles with the physical DFT values (R = 0.911).

Furthermore, we compared the Born effective charge (BEC) tensor, another well-defined physical quantity that corresponds to the derivative to the dipole moment with respect to atomic positions, i.e. . Figure 4e and 4f compare the BECs predicted using the LES charges to the DFT reference values, for both the diagonal (R = 0.976) and off-diagonal (R = 0.838) BEC elements. There is again good agreement between the LES BECs and the DFT values. Overall, the agreement between DFT and LES dipoles, BECs, and quadrupoles shows that LES is able to convincingly model observables of the molecular charge density even though no charge information is explicitly input into the model training.

Again, we emphasize that DFT partial charges, such as MBIS, Mulliken charges, and Hirshfeld charges, are not physical observables - although there is significant disagreement between the such DFT charges and LES charges (Fig. 4b), there are similar disagreement between the different flavors of DFT charges. Nevertheless, they are all good predictors of the observable molecular dipole and quadrupole moments, as shown below in the Methods section. In other words, the LES charges are just as physical as any definitions of DFT partial charges. The ability of LES to infer dipole and quadrupole moments as well as BECs just from energies and forces strongly supports the thesis that it is not necessary to explicitly learn a specific definition of DFT charges or electronegativities.

Ko et al. compiled four datasets (, , , and Au on MgO(001), illustrated in Fig. 5) that specifically target systems in different charge states or where charge transfer mediated by long-range electrostatic interactions is significant. In Table 1, we compare the CACE-LR errors with the values obtained with CACE-SR, 3G-HDNNP and 4G-HDNNP, as well as a charge constraint ACE model through a local many-body expansion (χ+η(ACE)). The comparison between CACE and ACE is a rather direct one: their descriptors are mathematically equivalent. 4G-HDNNP and χ+η(ACE) both fit charges explicitly, while CACE-LR only fits to energy and forces and no total charge constraint was used. We used a 90% train and 10% test split, consistent with ref. .

The set contains carbon chains terminated with hydrogen atoms in the neutral or positively charged state. With and without the added proton on the right-hand side of Fig. 5a, the atoms in the left half of the molecule can have almost identical environments but different atomic charges, which results in high fitting errors in 3G-HDNNP due to the contradictory information.

The example illustrated in Fig. 5b contains Ag trimers in two different charge states. As the system size is small such that there are no long-range interactions, we used only a short-ranged CACE MLIP with embedded charge states. Since the energies depend on the overall charge states of the clusters, this causes the degeneracy issue between atomic structures and potential energy surfaces, leading to the poor performance of the 3G-HDNNP and the charge-agnostic ACE methods. Both the charge constraint χ+η(ACE) model and the charge-state-embedded CACE lift such degeneracies, leading to drastically improved descriptions.

The set (Fig. 5c) contains the ionic clusters and when a neutral Na atom is removed. This is also an example where global charge transfer is present. CACE-LR achieves the lowest errors in this case.

The Au - MgO(001) set (Fig. 5d) has a diatomic gold cluster supported on the MgO(001) surface with two adsorption geometries: an upright non-wetting orientation of the dimer attached to a surface oxygen, and a parallel wetting configuration on top of two Mg atoms. Moreover, three Al dopant atoms were introduced into the fifth layer below the surface (the gray atoms in the left panel of Fig. 5d). Despite having large distances of more than 10 Å, the dopant atoms have a major influence on the electronic structure and the relative stability between the wetting and the non-wetting configurations.

In this example, CACE-LR achieves errors that are approximately an order of magnitude smaller than those of the other methods compared. As an additional test, we performed geometry optimizations of the positions of the gold atoms, with the substrate fixed, for both doped and undoped surfaces. The results were compared to reference DFT calculations and previous results using the 4G-HDNNP method. Note that the reference DFT results have been updated using tighter convergence settings of the geometry optimization, as performed by the authors of ref. . For the pure MgO substrate, the non-wetting configuration is energetically favored, whereas doping stabilizes the wetting geometry. The energy differences between the wetting and non-wetting configurations for both doped and undoped substrates are presented in Table 2. Short-range models, such as 2G-HDNNP and CACE-SR, predict nearly degenerate energy values for these configurations, as expected. In contrast, CACE-LR delivers highly accurate predictions, closely matching the reference results. Consistent with findings in ref. , we also present the potential energy surface for the non-wetting geometry on doped and undoped substrates as a function of the distance between the bottom Au atom and its neighboring oxygen atom, shown in Fig. 5e. Equilibrium bond lengths and energies derived from DFT are marked with black symbols. Notably, CACE-LR accurately resolves the distinct equilibrium bond lengths, with a slight shift in the potential energy surface likely attributable to differences in DFT convergence settings.

We rationalize why the CACE-LR method delivers significantly more accurate predictions compared to other long-range methods that explicitly fit atomic charges. In Fig. 5f, we compare the atomic charge distribution from the underlying DFT data, obtained via Hirshfeld population analysis, with the charges predicted by CACE-LR. The charges from CACE-LR are generally much smaller in magnitude and are primarily localized on the Au dimer and the dopant. In contrast, the DFT charges show sharp positive values for metal atoms and sharp negative values for oxygen atoms in the substrate. We hypothesize that explicitly modeling such DFT-derived charges for metals and oxygen is unnecessary for accurately predicting energy and forces. Short-ranged MLIPs are already well-suited to describe bulk oxides without dopants due to the screening effects that diminish the influence of these charge extremes. In Fig. 5g, we plot the changes in atomic charges resulting from doping, by taking the atomic charge difference for each atom from relaxed doped and doped structures, which shows a clear correlation between DFT and CACE-LR results. This example suggests that the charges predicted by CACE-LR can be interpreted as response charges rather than DFT partial charges, focusing on the aspects of charge redistribution relevant to energy and force predictions.

As example applications to electrolyte/solid interfaces, we selected two sets of systems. The first is the Pt(111)/KF(aq) interface dataset from ref. , which describes the Pt electrode with the (111) surface forming an interface with K and F ions in water solutions. For training the MLIP, ref. used a DPLR model: the short-ranged part is a standard Deep Potential (DP) model with a cutoff of 5.5 Å, and the long-range electrostatics is computed using spherical Gaussian charges associated with the nuclei (i.e., 6 e, 1 e, 9 e, 7 e, and 0 e for O, H, K, F, and Pt atoms, respectively) and the average positions of the MLWCs with a total charge of -8 e associated with each O, K, and F atom. Note that such MLWC schemes are not applicable to conductors, so ref. used the classical Siepmann-Sprik model to describe the Pt electrode in MD simulations.

The second dataset from ref. is for modeling the anatase TiO (101) surface in contact with NaCl-water electrolyte solutions at various pHs. This dataset comprehensively spans the configurational space of bulk anatase TiO, water, and various aqueous electrolyte solutions (NaCl, NaOH, HCl, and their mixtures), as well as anatase (101) interfaces with each of these liquids. ref. trained a standard short-ranged DP and a DPLR MLIP. The LR part in the DPLR model is also based on the electrostatics of spherical Gaussian charges associated with the ions (nuclei + core electrons) and the valence electrons. More specifically, 4 e, 1 e, 6 e, 9 e, and 7 e for Ti, H, O, Na, and Cl ions, and each O, Na, and Cl ion has four WCs each carrying -2e.

We fitted the CACE-SR and CACE-LR models, without message passing. The results are presented in Table 3. We speculate that the improved performance of the CACE models compared to the DP models can be attributed to two reasons: First, the DP descriptors are restricted to two-body and three-body terms, while the ACE framework can include higher-body-order interactions and in this case we truncate to four-body terms. The inclusion of higher-body terms makes the model more expressive and helps alleviate the degeneracy problem. Second, the LES scheme allows each atom to carry a flexible learned charge, in contrast with the fixed charge in the DPLR method.

To showcase the effect of long-range interactions on the structures of the electrolyte and the electric double layer (EDL), we performed MD simulations at 600 K for 5 ns on a large system of anatase TiO surface and NaCl in water solution (illustrated in Fig. 6). This is also a test that was performed in ref. . Figure 6b shows the ion distributions obtained from the MD simulations using the CACE-SR and CACE-LR models. In reality, the solution should recover its bulk properties in the central region that is away from the interface and have equal densities of Na and Cl ions. However, the SR model, lacking long-range electrostatic interactions, imposes no energy penalty for unphysical charge imbalances. Consequently, the MD simulation predicts an excess Cl density of approximately 0.05 mol/L in the center of the box. In contrast, incorporating long-range interactions with the CACE-LR model eliminates this artifact and alters the ion distributions within the EDL. These effects, including the correction of charge imbalance and modified EDL structures, were also reported in ref. . Notably, the CACE-LR model predicts a significantly lower second Na density peak near the interface compared to ref. . Figure 6c shows the predicted LES charges on atoms at different positions. Mostly notably for oxygen (red symbols) and titanium (gray symbols), the magnitude of these charges are dependent on whether the atoms are in the bulk region or at the interface. Such variance can be understood as coming from the difference of polarization environments. and may help capture the complex electrostatic interactions in interfacial systems.

Atomistic modeling of solid-solid interfaces is essential in understanding material synthesizability. The heterogeneous nature of these interfaces requires long-period structures, particularly in cases involving charge transfer, which necessitates long-range descriptions beyond standard MLIPs. To evaluate the predictive accuracy of our models, we conducted a benchmark study comparing CACE-SR and CACE-LR using the LiCl(001)/GaF(001) interfacial system. The training dataset includes bulk and interfacial configurations in the LiCl-GaF chemical space with corresponding DFT-calculated energies and interatomic forces. To assess model uncertainty, we trained an ensemble of four SR/LR models and used their predictions to estimate force uncertainties (see Methods). For in-distribution (ID) test set performance, CACE-SR and CACE-LR models achieve RMSEs of 78.8 meV/Å and 67.8 meV/Å, respectively.

To evaluate model transferability, we constructed an out-of-distribution (OOD) test set using a large solid-solid heterostructure relaxed with DFT calculations (~30 Å in the z-direction, Fig. 7a). This extended structure, containing eight Ga layers and four Li layers, represents a more realistic interface with much reduced finite-size effects compared to the training configurations. On this OOD set, the LR model demonstrates improved predictive accuracy with a force component error of 40.5 meV/Å compared to 116.3 meV/Å for the SR model. The atomic-resolved force errors are visualized in Fig. 7c, d, which were computed from the square root of the sum of force component errors in x, y, z-directions.

Force uncertainties were quantified using ensemble variance from the four trained models. The SR model exhibits lower uncertainties (Fig. 7e), indicating a good parametrization on the ID training set. In contrast, the LR model shows elevated uncertainties (Fig. 7f), effectively identifying OOD atomic environments in the heterostructure. The correlation between the absolute force errors (RMSE against DFT) and uncertainties is shown in Fig. 7b, where green dots specifically highlight the relationship between SR model errors (poor prediction) and LR model uncertainties (OOD detection). Interestingly, the LR model identifies regions of SR model failure (green dashed circle in Fig. 7b), which are further evidenced by the spatial correspondence in Fig. 7c, f. These results suggest that despite the SR MLIPs achieving adequate ID performance for this system, they lack the mathematical framework to capture long-period structure features that are essential for electrostatic interactions. In contrast, the LR models with LES overcome this limitation with improved transferability. More generally, the enhanced OOD detection capabilities are essential for robust uncertainty quantification in broader applications such as materials property predictions and generation. While our current implementation relies on computationally intensive ensemble variance, the LES framework is compatible with various uncertainty quantification methods, including Gaussian mixture models, Monte Carlo dropout, and deep evidential regression.