A foundation machine learning potential with polarizable long-range interactions for materials modelling

Domain: physics
Published: 2025-11-25
Editorial outcome: major_revision

Abstract
Long-range interactions are essential determinants of chemical system behavior across diverse environments. We present a foundation framework that integrates explicit polarizable long-range physics with an equivariant graph neural network potential. It employs a physically motivated polarizable charge equilibration scheme that directly optimizes electrostatic interaction energies rather than partial charges. The foundation model, trained across the periodic table up to Pu, demonstrates strong performance across key materials modeling challenges. It effectively captures long-range interactions that are challenging for traditional message-passing mechanisms and accurately reproduces polarization effects under external electric fields. We have applied the model to mechanical properties, ionic diffusivity in solid-state electrolytes, ferroelectric phase transitions, and reactive dynamics at electrode-electrolyte interfaces, highlighting the model’s capacity to balance accuracy and computational efficiency. Furthermore, we show that as a foundation model, it can be efficiently finetuned to achieve high-level accuracy for specific challenging systems.

Introduction
Molecular dynamics (MD) simulations are indispensable in describing phenomena at the atomic level in disciplines such as chemistry, materials science, and biology1–4. While ab initio MD (AIMD)5 simulation approaches offer unparalleled accuracy, their computational demands constrain their application to small timescales and length scales. On the other hand, classical MD (CMD)6,7 simulations offer computational efficiency, but their accuracy is contingent upon empirical parameters. Machine learning interatomic potentials (MLIPs)8–11 provide a solution that balances accuracy and computational efficiency. In particular, MLIPs based on equivariant graph neural networks, such as NequIP12, DimeNet13, and MACE14, achieve excellent performance by introducing equivariant and invariant symmetries. Moreover, universal foundation models trained on the periodic table, such as M3GNet15, CHGNet16, GNoME17, MACE-MP-018, MatterSim19, SevenNet20, and ORB21 have emerged, showing remarkable prospects for materials discovery.

While existing local MLIPs with a cut-off of around 5Å perform well in simulating interactions within localized chemical environments, they may fail to capture long-range phenomena. This hinders their ability to understand and elucidate the behaviors of complex materials22. One approach is to implicitly incorporate the effects of long-range interactions through message-passing mechanisms23, which can extend the cut-off through layer-wise propagation. However, if certain parts of the system are disconnected on the graph, such as two molecules separated by a distance beyond the cut-off value, the message-passing scheme becomes ineffective24. A promising direction involves investigating the potential advantages of explicitly including long-range interactions in the model formulation.

These behaviors encompass electrostatic interactions, involving forces between charged particles, and dispersion terms, which arise from dynamic fluctuations in electron distribution within molecules or atoms.

Various explicit models have been proposed to address the challenges associated with long-range interactions by including electrostatic interactions. In these methods, charges are usually obtained through charge equilibration (QEq)25 principles, including the Charge Equilibration via Neural Network Technique (CENT) method26, Fourth-generation High-Dimensional Neural Network Potential (4G-HDNNP)27,28, and Charge Equilibration Layer for Long-range Interactions (CELLI)29 method. Generally, electrostatic parameters were trained to minimize the difference between QEq charges and reference partial charges derived from quantum mechanics (QM) calculations. Notably, partial charges obtained from QM calculations are either derived from the direct partitioning of the molecular wave function into atomic contributions or calculated based on the analysis of physical observables. However, due to the incompleteness of basis sets and variances in the partitioning methods, properties calculated based on partial charges may not always be reliable16,30. The ambiguity in DFT-assigned charges suggests that directly learning of them may be inessential for constructing accurate interatomic potentials26,31,32. On the contrary, the polarizable charge equilibration (PQEq)30,33–36 method proposed by Naserifar et al. enhances the QEq method by using QM electrostatic interaction energies instead of partial charges as targets for evaluating interatomic potentials. Additionally, PQEq has proven to be effective in capturing polarization effects by explicitly introducing core-shell displacement. Another way to solve the partition-dependent issue is to take the dipole moment into the loss function, as done by SpookyNet37. However, it still uses point-like charges, making it behaves like QEq-based methods when dealing with polarization effects.

In this work, we introduce a framework that integrates the equivariant message-passing neural network potentials with the explicit polarizable long-range potential. We initially conducted benchmark tests on datasets containing systems with various charge states proposed by references28,38, showcasing the efficacy of our framework in handling both long-range interactions and varying charge environments. Building on these results, we extended the framework to develop a foundation model encompassing the periodic table up to Pu, while maintaining computational efficiency. We have applied this foundation model across different areas of applications, including the predictions of mechanical properties of materials, phase transitions in ferroelectric materials, and reactive molecular dynamics of solid-state batteries interface, thus leading to significant advancements in the fields of materials science and chemical simulations. The model also demonstrates transferability, accurately predicting interactions between clusters in various charge states and molecular polarizability. Furthermore, we show that the foundation model can be efficiently finetuned to achieve ab initio accuracy for specific challenging systems.

Results
Theoretical framework
Our framework, which integrates neural network potential and polarizable long-range interactions, is depicted in Fig. 1. For a given chemical system and boundary conditions, our goal is to construct a mapping from atomic coordinates r and atomic types Z to the total potential energy E . The potential energy is expressed as the sum of the second-order expansion30 with respect to charge fluctuations and the density functional theory39 D3 (DFT-D3) van der Waals dispersion energies correction term40,41 E , partial charges (q) and shell positions (r ). The partial charges can be obtained self-consistently in each MD step. The forces on atoms and virial stress on the cell can be calculated via automatic differentiation after partial charges are determined, as detailed in Supplementary Note 2.

Validation on diverse charge-state datasets
To evaluate the capability of our framework in capturing different charge states and long-range interactions, we validated our framework against the dataset developed by Ko et al.28 and Maruf et al.38, which encompass various charge states and charge transfer systems. They contain six distinct subsets: Ag cluster with positive and negative total charge, (Ag +/−), Na-Cl ionic cluster with one neutral Na removed (Na Cl +), hydrogenated carbon chains in both neutral and cationic states (C H /C H +), a periodic system consisting of Au clusters adsorbed on an MgO-(001) surface, and benzotriazole interactions with Cu-(111) surfaces in dry (BTA-Cu) and aqueous (BTA(HO)-Cu) environments. We demonstrate the advantage of integrating equivariant message-passing networks with explicitly polarizable long-range interaction models by benchmarking against both 4G-HDNNP, a local MLIP with explicit long-range interactions, and NequIP, which implicitly captures long-range effects through its message-passing mechanism. 4G-HDNNP differentiates charged systems by explicitly training DFT partial charges and incorporating them as descriptors into neural networks for short-range interactions. We employed the NequIP models trained in the reference38 (Maruf’s NequIP) as our baseline equivariant model. Our framework explicitly incorporates long-range interactions while adopts Maruf’s NequIP as the E component without touching its settings, thereby ensuring a fair comparative analysis. The results are presented in Table 1, with detailed neural network architectural specifications provided in the Supplementary Table 3.

For the carbon chains, the introduction of additional H+ causes some C atoms to exhibit opposite charge states43, as shown in Supplementary Fig. 1a. Although we did not fit partial charges, our model qualitatively captures the distribution of DFT partial charges, as shown also in Supplementary Fig. 1a. In the linear C H molecule, our model demonstrates symmetric charge distribution around the center which shows agreement with DFT. After explicitly incorporating physical long-range interactions, our framework achieved additional improvements in error metrics compared to NequIP and 4G-HDNNP. For positively charged NaCl clusters, the predicted energy and forces by our model are still more accurate compared to 4G-HDNNP, as well as the potential energy surfaces shown in Supplementary Fig. 1b. For Au adsorption on MgO surfaces and benzotriazole interactions with Cu-(111) surfaces in dry and aqueous environments, our model outperforms NequIP significantly. Additionally, our model correctly predicted the preferential adsorption geometries for both doped and undoped MgO, as listed in Supplementary Table 1. As for the charged Ag clusters, due to their small system size and lacking long-range effects, the entire structure falls within the typical cutoff radius used by the MLIP. The energy dependence on the overall charge state of clusters lead to degeneracy between atomic structures and potential energy surfaces, resulting in poor performance of the charge-independent NequIP model. We acknowledge that for such small, highly charged systems, the 4G-HDNNP approach to learning charge equilibration parameters is more convenient. Nonetheless, we also achieved reduced force error with the fixed physical parameters compared to 4G-HDNNP. Although PQEq partial charges are not directly comparable to DFT Hirshfeld ones, we compare them to show their relevance. Varying degrees of agreement across different chemical systems are observed, as shown in Supplementary Table 2.

Our comprehensive evaluation demonstrates that explicit incorporation of physical long-range interactions significantly enhances the performance of MLIPs across diverse charge-state systems. Notably, our proposed framework outperforms 4G-HDNNP in most cases for energies and forces without fitting partial charges. This advantage likely stems from a more physically meaningful partition of the total energy.

Foundation model benchmark
We trained a foundation model for all the periodic table elements up to Pu using the MPtrj dataset16 following our framework (our model), as described in the Methods section. We also trained a model maintaining the same MLIPs framework and dataset, but without the polarizable long-range interactions (w/o-lr model). The former one demonstrated mean absolute errors (MAE) of 18meV per atom for energies, 0.065eV Å−1 for forces, and 0.301GPa for stresses across the test set, outperforming the w/o-lr model, as shown in Supplementary Table 4. These proved the benefit of incorporating long-range interactions in the enhancement of prediction accuracy.

Message-passing mechanisms will fail when graph connections break in the system. For instance, in dimers, when the interatomic distance exceeds the graph cut-off, the atoms no longer interact with each other. We analysed interaction energies for Na-Na and Na-Cl dimers across interatomic distances ranging from 1.8Å to 12.0Å. As shown in Fig. 2a, both models (our model and the w/o-lr model) trained on the MPtrj dataset accurately reproduce the equilibrium bond lengths of the two dimers. However, the w/o-lr model failed to differentiate interaction energies beyond the MLIPs cutoff distance, i.e., 5.0Å, while our model successfully captured the DFT potential energy surface throughout the entire distance range, demonstrating its capability in describing extended ionic interactions. Moreover, the message-passing mechanism exhibits deficiencies when simulating layered materials. Taking lithium iron phosphate (LiFePO ) as a case study, we observed that the w/o-lr model predicted an irreversible phase transition, as shown in Supplementary Fig. 3. In contrast, our model correctly reproduced the experimentally44 and AIMD45 thermally stable olivine structure. These results collectively demonstrate that the explicit long-range interactions are crucial for accurately describing both ionic systems and layered materials, particularly in predicting structural stability and thermal behavior.

Polarization plays a fundamental role in determining molecular responses to external electric fields. Conventional QEq-based models have limitations in their response to an external electric field applied orthogonally to the molecular dipole. We validated this capability through a benchmark study examining a water molecule’s response to external electrostatic fields46. The water molecule was fixed in the yz-plane, and we applied varying external electric fields along the x-axis to derive energy curves, with the results presented in Fig. 2b. The model with polarizable long-range interactions demonstrates a trend consistent with the reference values obtained from DFT. Due to the external electrostatic field direction consistently being orthogonal to the molecule, the electrostatic energy in the QEq-based model (like 4G-HDNNP) is 0 and is completely independent of field strength (i.e., μ·ϵ=0, where μ, ϵ are dipole moment and electric field, respectively). This is because the QEq model treats atoms in the molecule as point charges or Gaussian charges, without distinguishing between core and shell components. We also compared the response of O atom partial charges across DFT Hirshfeld calculations, our model, and the conventional QEq model, as shown in Supplementary Fig. 4. The agreement between the PQEq model and DFT is observed, while the QEq model completely lacks response. This comparative analysis establishes that the explicit polarizable interactions are crucial for modeling molecular systems subject to external electric fields.

Besides the failure of response in static conditions, the QEq method may exhibit non-physical behavior in dynamics47. As illustrated in Supplementary Fig. 5, we conducted additional investigations of two water molecules under an external electric field of 0.25V Å−1. In charge distribution models, it is essential to allow systems to form an internal electric field opposing the external field. The QEq model assumes the system behaves as a perfect conductor without penalizing charge transfer as a function of distance. This limitation of the QEq method leads to non-physical charge transfer between molecules separated by large distances47. In contrast, within the PQEq model, the individual water molecules can polarize, generating internal electric fields that counteract the external field, thereby capturing more realistic electrostatic responses in molecular systems. Under the QEq scheme, the water molecules accumulate non-zero net charges and migrate to the top and bottom of the simulation box. In contrast, our model maintains charge neutrality of the water molecules, which remain stationary despite the applied electric field, aligning with DFT results.

We conducted benchmarking using properties that were not labeled during the training process. We applied our foundation model to predict the mechanical properties of 10,154 materials from the Materials Project48. Figure 2c–e illustrates the comparison of the bulk modulus B determined by the foundation model and DFT. Our model demonstrated impressive performance, achieving an R2 of 0.94 with the Voigt approach B and Hill average method B , and 0.93 using the Reuss method B . In comparison, as shown in Supplementary Fig. 2, the w/o-lr model achieved lower R² values of 0.89, 0.86, and 0.88 for B , B , and B , respectively. These not only highlight the robust performance of the model but also establish a solid foundation for the high-throughput screening of materials with exceptional mechanical properties.

Transferability of the foundation model
In terms of foundation model transferability, evaluations were conducted across diverse systems. The model accurately reproduces the potential energy surfaces for OH-OH systems and successfully distinguishes between neutral and ionic states in OH-OH systems, capturing the long-range Coulombic repulsion of ionic state as shown in Supplementary Fig. 8. Further validation using a periodic water system demonstrates the model’s capability to predict polarization effects under varying electric fields, achieving comparable performance to specialized models as depicted in Supplementary Fig. 10. Additionally, as shown in Supplementary Fig. 11, when tested on 7211 molecules from the QM-7b dataset49, the model shows remarkable accuracy in predicting molecular polarizability with a mean absolute error of 4.57 atomic units (a.u.), highlighting its robust transferability in capturing polarization.

Computational efficiency
Despite the incorporation of additional long-range interaction calculations, our foundation model maintains high computational efficiency. Performance benchmarks were conducted on a single NVIDIA H100 GPU to quantitatively assess the computational cost. As illustrated in Supplementary Fig. 6, for a periodic system comprising 2160 atoms, our model requires approximately 0.07s per molecular dynamics time step, representing a significant improvement over conventional universal MLIPs, which require approximately 0.91s. This computational advantage extends to larger systems, as demonstrated by simulations of a 24000-atoms system, where our model maintains efficiency with only 1.49s per time step. Such computational performance enables nanosecond-scale molecular dynamics simulations of systems containing tens of thousands of atoms, making it practical for more realistic modeling of materials.

Materials modeling applications
To demonstrate the capability of our model to reproduce results from AIMD simulations, we simulate the kinetic transport properties of a solid-state electrolyte. Specifically, we investigated lithium-ion diffusivity within the cubic phase of Li La Zr O (c-LLZO)50, a crystalline superionic conductor known for its remarkable stability as a lithium conductor. Owing to the efficiency of our foundation model, we are able to conduct MD simulations on the c-LLZO comprising 64 formula units, with a duration of 2ns for each temperature range from 800K to 1800K. Figure 3a presents the Arrhenius plots derived from our model, w/o-lr model, and AIMD51 simulations. The diffusion coefficients were determined from the slope of the logarithmic mean squared displacements (MSD) versus logarithmic time within the Fickian regime52. MSD curves of our model are depicted in Fig. 3b. Compared to previous work using AIMD simulations, our model clearly reproduces the diffusion coefficients and activation energy. The performance surpasses that of the w/o-lr model, which lacks explicit long-range interactions. With the ability to reach larger model sizes and longer simulation times, our model achieves these results with significantly reduced uncertainties. This capability not only delivers comparable results to AIMD simulations but also enables larger-scale simulations that cannot be assessed by AIMD, enhancing the statistical significance of MD studies. Consequently, the framework developed in this work provides more reliable and comprehensive data for diffusion analysis in solid-state electrolytes.

To show the application of polarizable potentials, we performed MD simulations on a typical perovskite ferroelectric material, BaTiO . As temperature increases, BaTiO undergoes a sequential phase transition process, transforming from rhombohedral to orthorhombic, then to tetragonal, and finally to a cubic structure53. Accurately simulating phase changes in ferroelectric substances requires precise potential energy functions that can respond to the small atomic shifts and structural changes, as well as account for the free energy landscape under the conditions of finite-temperature thermodynamics54. The phase sequence of BaTiO has been extensively studied using effective Hamiltonians55–58, CMD59,60, reactive force fields (ReaxFF)61, specialized MLIPs models trained on the BaTiO system54,62, and experiments63. To detect lattice distortions that distinguish different phases, a 10×10×10 supercell of BaTiO was simulated. The use of larger simulation cells (10×10×10 in our work versus 4×4×4 in Gigli et al.54) may provide some advantages for phase transition studies. Larger supercells can effectively reduce the impact of temperature fluctuations, enabling better temperature control for sampling the free energy landscape, thereby yielding results with greater statistical significance. Figure 4 illustrates the simulation results, which clearly identify four distinct phases and the three first-order phase transitions within the simulated temperature range. Below 145K, the overall average polarization of the supercell aligns with the (111) direction, indicative of the rhombohedral phase. At 145K, a phase transition to orthorhombic is suggested as the y-component of polarization approaches zero. With further increase to 205K, the polarization predominantly orients along the x-direction, indicating the tetragonal phase. This phase persists until 285K, where the cubic paraelectric phase is observed. The obtained phase transition temperatures align closely with those predicted by various models and experiments, as listed in Table 2. It is worth noting that the supercell size can affect phase transition temperatures, so that direct comparison between different simulations is not expected. Utilizing first principles methods or fitting potential energy surfaces to first principles data often underestimates these temperatures, which is attributed to the approximated exchange-correlation functional59 employed in DFT. Although the predicted phase transition temperatures are lower than the experimentally observed values, the model effectively captures the sequence of phase transitions as measured in experiments, demonstrating the efficacy of the foundation model in modeling highly polarized phase transitions in ferroelectric materials.

The final application studied in this work involves all-solid-state batteries, which represent a breakthrough in the evolution of next-generation energy solutions, owing to their superior energy density and inherent safety features64. Lithium thiophosphate, known for its exceptional ionic conductivity (~10−3Scm−1), is deemed the most promising candidate for solid electrolytes and has been widely studied through experiments65–67 and computer modeling68–70. Notably, the nanoporous β-Li PS has been validated to exhibit outstanding cycling stability71, presumably attributed to the formation of Li S and Li P solid-electrolyte interphases (SEI) during the initial battery cycles, which serves to passivate further degradation of the electrolyte69,72–74. To investigate the formation of SEI, we utilized the foundation model to simulate the interfacial reactive MD of the solid electrolyte β-Li PS in contact with the lithium metal anode, encompassing a system of 13,760 atoms. The initial β-Li PS /Li interfacial structure, measuring 33.15nm×3.1nm×2.6nm, was utilized in the MD simulations, as depicted in Fig. 5a. The highly reactive Li metal began to interact with the PS tetrahedron at the interface, triggering the formation of the SEI layer as illustrated in Fig. 5b. With the growth of the SEI layer, electrons were transferred from the Li anode to the Li PS electrolyte, which causes a gradual increase of the partial charges of Li as they migrate from the anode, across the SEI, and into the electrolyte. This results in a transition from metal lithium to ions. The dynamic behavior of partial charges at the anode-electrolyte interface is a critical aspect showcased in Fig. 5c. Initially, the partial charges of lithium changed almost linearly from the anode to the electrolyte. Over time, this distribution evolved and showed a distinct plateau. This phenomenon could be attributed to the ordered structuring of the SEI, signifying the nucleation of crystalline structures. Ultimately, a stable lithium partial charge plateau was formed between the anode and the electrolyte, spanning a range of 15 to 21.5nm. This is further supported by the visualized structure depicted in Fig. 5b, where the plateau within the SEI aligns with the Li S crystal and amorphous regions. Chen et al.68 employed AIMD simulations to investigate the radial distribution functions (RDF) of Li-P and Li-S bonds, which concluded that the decompositions of electrolytes in lithium thiophosphates are primarily due to the decomposition of PS tetrahedron by the active lithium metal, leading to the formation of new Li-P and Li-S bonds. Our RDF analysis confirmed the formation of Li-S and Li-P bonds within the SEI and electrolyte as depicted in Fig. 5d. In agreement with the AIMD simulations68, we observed the emergence of a new Li-P peak within the SEI at approximately 2.5Å. Furthermore, within the SEI, the Li-S bond displayed RDF characteristics consistent with crystalline Li₂S75. For the Li-P bonds, the peak shape closely resembles that of amorphous Li P AIMD simulations76.

As illustrated in Fig. 6a,b, the ultimate decomposition products of P form only short-range ordered structures, in contrast to S, which forms long-range ordered structures. This observation agrees with the measurements where Li S crystals were detected but Li P crystals were absent in cryogenic transmission electron microscopy (cryo-TEM) experiments, as noted in the reference66. To enhance the understanding of the formation mechanism, the atomic compositions of the SEI in both crystalline and amorphous Li S regions were examined. As depicted in Fig. 6c,d, within the first 2ns, the emergence of amorphous Li S and Li P regions was attributed to the swift diffusion of P and S atoms. Conversely, the formation of crystalline Li S regions was predominantly influenced by the swift diffusion of Li atoms and the slow diffusion of P and S atoms, which also confirmed that the formation of crystalline Li S region hindered further atomic diffusion and thus slowed the growth of the SEI. As illustrated in Fig. 5b, during the initial phase of decomposition (0–100ps), the swift diffusion of P, S, and Li led to the formation of an interface of about 3nm thick. However, in a later stage between 100 and 2000ps, the formation of Li S nuclei impeded the diffusion of P and S atoms, resulting in a period of sluggish SEI growth. Consequently, the interface experienced only a slight expansion in this period, increasing by only 1nm. Furthermore, in the following 2ns, there was essentially no SEI growth. The final interface structure (~8.5nm) could be characterized by a 4.5nm crystalline Li S region sandwiched between two 2-nm transitional layers. The cessation of growth at this interface suggests that the crystalline Li S and amorphous Li P within the crystal region contribute to the stabilization of both the lithium metal anode and the electrolyte, which agrees with the theoretical predictions made by DFT calculations72. Additionally, we studied the interface reactions between Li PSCl and lithium metal, as detailed in Supplementary Note 13. As shown in Supplementary Fig. 13, at a pressure of 1atm and a temperature of 300K, Li PSCl formed an interface of about 5-nm thick after 2ns. This indicates that the Li PSCl electrolyte exhibits a higher stability against lithium metal at room temperature compared to Li PS . Upon elevating the temperature to 500K, we detected Li S at the grain boundaries, which aligns with the cryo-TEM observations66.

Finetuning for enhanced accuracy
While our foundation model demonstrates general performance, achieving ab initio accuracy for specific systems requires finetuning. Using Na Cl + clusters and Au dimers on MgO surfaces from Ko et al.28 as test cases, we performed targeted finetuning using 20% subsets of configurations and mean squared force error as the sole loss function. After finetuning, our model achieved good agreement with DFT references, accurately distinguishing potential energy surfaces for positive charged clusters and surface conditions, as shown in Supplementary Fig. 12 and Supplementary Table 6. This success demonstrates the effectiveness of the two-step approach: foundation pretraining followed by efficient finetuning, enabling quantum-accurate predictions while maintaining the model’s fundamental physical insights and transferability.

Methods
Model training details
To train the foundation equivariance neural network potential, we used the MPtrj dataset16 sourced from Materials Projects48 as the training dataset. All configurations were calculated using DFT with the PBE82/PBE+U83 exchange-correlation functional and pseudopotential basis. The PBE/PBE+U mixing compatibility correction84 was applied to ensure energy consistency. Due to transferability issues with the Yb element85, all data containing Yb in the MPtrj dataset were excluded. The higher-order energy terms and the electrostatic Coulomb interaction energies due to charge fluctuations should be deducted from the total potential energy. Appropriately setting the PQEq parameters is crucial for training a foundation framework. Detailed information on the fitting process and the PQEq parameters used in this work can be found in the Supplementary Note 10 and Supplementary Table 7.

For equivariant neural network training, the NequIP12 model was utilized, incorporating six equivariant layers with node features set at 256×0e+256×1e. Spherical harmonic basis functions were used to represent interatomic directions, denoted as 1×0e+1×1o. Interatomic distances were described using eight Bessel basis functions, and a cut-off value of 5.0Å was chosen for building the MLIP neighbor list. The dataset is divided into a training set, validation set, and test set in a ratio of 8:1:1. The total loss function (L) of prediction values xpred and reference first principle values xref can be expressed as,

L(xref,xpred) = 1/N * sum(xref - xpred)^2

L=ωE*LE(Eref,Epred) + ωF*LF(Fref,Fpred) + ωS*LS(Sref,Spred)

where energies E, forces F, and stresses S weights were set as ωE=1, ωF=1, and ωS=10 with quantity units of eV, eV Å−1, and eV Å−3. The batch size was selected as 128, and the Adam optimizer was employed with an initial learning rate of 10–3. All the trainable parameters were initialized with the random seed 3407. The implementation of the entire neural networks and long-range interactions was based on JAX 0.4.2086. The training was conducted on a single NVIDIA H100-PCIe GPU, utilizing CUDA version 12.3 and Driver version 545.23.06.

Molecular dynamics simulations
The Atomic Simulation Environment (ASE)87 was utilized as the interface for geometry relaxations and MD simulations. Structure relaxations were conducted using the limited-memory BFGS88 method. Both atomic coordinates and cell vectors were optimized concurrently until the forces fell below the convergence threshold of 0.05 eV/Å.

The c-LLZO with a space group of Ia3d originated from experimental structures detailed in the reference50. The structure was expanded to a 2×2×2 supercell. Subsequently, 180 lithium atoms were randomly placed on the 24(d) sites, and an additional 268 lithium atoms on the 96(h) sites. After relaxation, the structure served as the starting point for MD simulations. The isothermal isobaric (NpT) ensemble simulations were performed with the temperature control via Nosé-Hoover thermostat89,90 and pressure maintained at 1atm using Parrinello-Rahman barostat91,92 to determine the lattice constants at various temperatures. Then, the canonical (NVT) ensemble simulations with the Nosé-Hoover thermostat were utilized to ascertain the lithium diffusion properties within c-LLZO. A timestep of 2 fs was used for all NpT and NVT simulations. In the case of NVT simulations, following a 100ps equilibration period, 2-ns trajectories’ mean squared displacements of Li ions were employed to calculate the self-diffusion coefficients at various target temperatures. Additionally, we conducted an uncertainty analysis of the diffusion coefficients, applying the empirical error estimation method as detailed in the reference51.

To construct the phase diagram of the BaTiO crystal structure, a supercell consisting of a 10×10×10 R3m BaTiO lattice with 5000 atoms was created. The crystal structure data for BaTiO was sourced from the Materials Project48. NpT simulations with the Nosé-Hoover thermostat and Parrinello-Rahman barostat were performed with a 1-fs timestep. Throughout the simulations, the pressure was consistently maintained at 1 atmosphere. To regulate temperature, the procedure began at 5K and involved incrementing the temperature by 10K every 20ps. The average lattice constants and the local polarization of the unit cell, denoted as P, were determined using data from the final 10ps of trajectories at each temperature step. For the comprehensive calculation method of the local polarization, please refer to the Supplementary Note 12.

We utilized a 13760-atom structure consisting of a ca. 13.6-nm thick Im3m lithium metal anode with (100) facet, paired with a ca. 19.1nm thick Pnma β-Li PS electrolyte with (010) facet to illustrate the SEI formation across realistic spatial and temporal scales. The system was structured with periodic lateral dimensions measuring 3.1nm by 2.6nm. The unit cell of anode and electrolyte structures originated from the Materials Project48. The entire cell was subject to periodic boundary conditions, with a 1nm atomic layer fixed on each side along the x-direction to ensure a singular interface reaction. NpT MD simulations were conducted at a standard temperature of 300K and a pressure of 1atm by Nosé-Hoover thermostat and Parrinello-Rahman barostat with a timestep of 2fs. Initial velocities for the particles were assigned randomly following the Maxwell-Boltzman distribution, and the simulations were carried out for a duration of 4ns.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors present a Machine Learning (ML) potential where the short-range graph neural network prediction is combined with the long-range electrostatic and dispersion interactions. While similar approaches were proposed previously, the novelty of this manuscript is in using the polarizable charge equilibration (PQEq) method instead of the charge equilibration (QEq) method. In addition, the authors do not fit the partial charges and thus avoid problems with charge partitioning schemes. Tackling long-range interactions with ML potentials is an active research area, and this manuscript presents an important contribution to the field. Nevertheless, the authors should address the following points before publication:\n\n1. In message-passing graph neural networks, the graph cut-off should not be confused with the cut-off of classical force fields since each message passing/interaction layer increases the receptive field. The authors use a NequIP model with six equivariant layers. Thus, the effective cut-off is not 5 Å but can extend up to 6 * 5 Å. This point should be made clear in the text, i.e., in the introduction where the authors discuss the “MLIPs with a cut-off of around 5 Å” and in the results section when the authors are comparing their results with strictly local ML potentials such as HDNNP, where the employed cut-off for the symmetry functions is the same as the effective cut-off.\n\n2. Related to the above, the test cases developed by Ko et al. involve very small systems and make sense for strictly local ML potentials but make little sense with a 6-layer graph neural network (GNN), where the receptive field likely encompasses the whole system. The identify the PQEq contribution, the authors should extend Table 1 to include results using only GNN and using GNN+DFT-D3. On the other hand, a comparison with 2G-HDNNP and 3G-HDNNP is not necessary since these results were already published.\n\n3. It is unclear if the resulting partial charges should be interpreted as physical partial charges, or a quantity correlated with the partial charge. The authors should report the errors for partial charge prediction for Ko et al. cases. Of course, the agreement cannot be as good as when training on a partial charge, and given the differences between the charge partitioning schemes, this is also not the goal. Nevertheless, the authors should compare the partial charge prediction errors with differences between charge partitioning schemes to clarify the meaning of predicted partial charges.\n\n4. For the universal machine learning potential model with and without polarizable long-range interactions, the energy/force/stress errors on the test set are very similar. The authors should include the results of the universal model without polarizable long-range interactions also for the reported applications, i.e., Mechanical properties, Ionic diffusivity etc.\n\n5. Figure 2 b-c: the authors should increase the x-axis range to 12 Å to actually see the long-range scaling behaviour.\n\n6. Are the model training details described in the Methods section the same for the universal equivariance neural network potential and the test cases developed by Ko et al.? If not, please add training details for Ko et al. cases.\n\n7. “... MPtrj Dataset sourced from Materials Projects as the pretraining dataset ...” Why do you call this pretraining? which dataset was used in the subsequent training?\n\n8. Why is the shell charge of all atoms set to 1?\n\n9. Missing directly related reference: 10.48550/arXiv.2501.19179\n\n10. The authors should relate their work with alternative long-range approaches not employing the charge equilibration schemes, e.g., Ewald summation, Latent Ewald Summation, LSR-MP, NeuralP3M, and So3krates", "reviewer_2": "Overall, a manuscript which should be very heavily revised - but that might decrease any novelty reported such that it would fall outside of the scope of the journal. Maybe best to reject.\n- What are the noteworthy results?\nThere are some interesting simulations using the MLIP technology developed by the authors - but there is no reliable benchmarking carried out to compare to other state-of-the-art methods. Since the focus of the manuscript is a supposedly superior method, this is a claim that should have been proven by the authors.\n- Will the work be of significance to the field and related fields? How does it compare to the established literature? If the work is not original, please provide relevant references.\nThis work is mostly combining other, existing pieces of technology in a slightly different way. In short, a short-range MLIP is combined with Goddard's QEq scheme and D3 dispersion corrections. The closest similar is Behler's 4G-HDNNPs, but these don't use the equivariant message passing short range technology. Results presented here are marginally better, but perhaps only an incremental improvement. Presumably due to the better and more expressive short-range part. Behler's scheme is superior as it uses environment dependent parameters for the QEq component, unlike here, where these are fixed. The CENT method by Goedecker (which uses a similar approach to fit the energies directly) is not even mentioned.\n- Does the work support the conclusions and claims, or is additional evidence needed?\nA lot of data which the results are compared to here are cherrypicked (there is hardly or no any mention of significant universal force fields, such as MACE-MP and its later improvements). This is a quickly moving field.\n- Are there any flaws in the data analysis, interpretation and conclusions? - Do these prohibit publication or require revision?\nPlenty of misinterpretation/misrepresentation, listed below. With very careful revision maybe these can be rectified, but most likely it will further reduce any claims of superiority and/or novelty.\n- Is the methodology sound? Does the work meet the expected standards in your field?\nMostly.\n- Is there enough detail provided in the methods for the work to be reproduced?\nYes, and the code is publicly available (I have not tested, though).\nDo not antropomorphise mathematical models (e.g. p2/l39)\n\"subtle fluctuations\" - these are not subtle (it is meaningless in this context) - they are dynamic fluctuations\nTerms in Eq. 1 should be explained immediately after in the text.\nRef. 24 is incomplete\nEqs. 2 and 3 do not follow. Since the energies also depend on the charges, which are determined self-consistently, it is far from obvious that eqs. 2 and 3 are true. There may be an expression analogous to the Hellmann-Feynman theorem, but not shown here.\nPanels 2a, 2b and 2c must be united.\nThe QEq curve on 2d seems completely flat - this cannot be true, or the QEq must be extremely rigid.\nComputational cost being linear: I find that hard to believe. Long range interactions in periodic systems are typically treated using some version of the Ewald summation method or an equivalent approach. While there exists linear scaling implementations, there is no mention of these here.\nThe sentence \"as well as account for the free energy surfaces under the conditions of finite-temperature thermodynamics\" does not make sense. The free energy is a state function and not a surface.\n\"The QEq-based machine learning models inherently struggle with atomic polarization, leading to non-physical results under external electric fields.\" there is no evidence for that. In the results quoted by the authors on 4G-HDNNPs they may be less accurate but not unphysical.\n\"outperforming conventional equivariant GNN-based MLIPs\" no evidence is showed. (Ref. 48. is not an MLIP that fits this description).", "reviewer_3": "### Summary\nThe manuscript presents a semi-heuristic approach to incorporating long-range interactions, such as polarizability, into machine learning potentials. The proposed framework integrates a physically motivated polarizable charge equilibration scheme within an equivariant graph neural network-based short-range potential. While the idea is clever, it is not entirely new—similar approaches have been explored before, such as in *SpookyNet*, although the current work extends these ideas further. The authors claim their model is universal and exhibits linear scaling in system size while accurately capturing long-range interactions.\nWhile the proposed approach is an interesting idea and appears to perform well on specific test datasets, there are significant concerns regarding the validity of the claims made in the title and throughout the manuscript. The manuscript lacks rigorous benchmarking against existing universal models, and the claim of linear scaling contradicts fundamental physical/mathematical limitations. As such, I cannot recommend publication in its current form as outlined below in more details\n---\n### Major Concerns\n#### 1. **Lack of Rigorous Benchmarking for a \"Universal\" Model**\nThe authors make a very strong claim that their model is universal, but the provided evidence does not support this assertion. There are two critical issues with this claim:\n1. **Absence of Comparison with Established Models:** The manuscript does not compare the model’s performance against other well-known universal machine learning potentials (e.g., MACE-MP-0, MATTERSim, SevenNet, ORB). Not only do the authors omit numerical comparisons, but these models are not even cited. Instead, the error analysis is restricted to their own training dataset, which is insufficient to validate universality.\n2. **Inadequate Testing of Universality:** The authors evaluate their model on the common datasets published by Ko et.al. But, the author perform their benchmarks bespoke models rather than testing their claimed universal model. If the model is indeed universal, it should be evaluated directly on these datasets without additional training. A fair comparison would involve evaluating both the universal and bespoke models to determine whether the universal model effectively learns long-range physics.\n#### 2. **Misleading Scaling Claims**\nOne of the most problematic claims in the manuscript is that the model achieves linear scaling for long-range interactions:\n- **Fundamental Computational Limits:** Long-range interactions require solving the Poisson equation, which fundamentally cannot be computed in linear time. The best possible scaling for such a computation is **O(N log N)** using Ewald summation on a grid, while traditional solvers scale as **O(N³/²)** or worse.\n- **Questionable Interpretation of Figure S4:** The claim of O(N) scaling is based on Figure S4, but the presented data does not conclusively support this claim. If the authors indeed achieve linear scaling, they must be using a cutoff solver, which is known to introduce severe errors in energy and force computations.\nGiven these issues, the scaling behavior of the proposed model needs significant clarification.\n#### 3. **Concerns Regarding Specific Applications**\nThe manuscript presents results on well-known test cases, but the methodology used for comparison raises concerns:\n- **Barium Titanate (BTO):** The authors compare phase transition temperatures in Table 2, but this comparison is misleading because the transition temperature is highly dependent on system size. The authors use a 10×10×10 cell, whereas the reference MLIP from Gigli et al. models use a 4×4×4 cell. The comparison should be made at equal system sizes or, at the very least, the system size used should be explicitly stated in the table with a comprehensive discussion.\n- **Lithium Phosphorous Sulfide (LPS):** The claim that their model captures long-range interactions in LPS is questionable. In this material, interactions are largely screened, and no significant long-range charge transfer is expected to occur. This system has already been simulated successfully with non-long-range potentials (e.g., *Phys. Rev. Materials 8, 065403*), raising doubts about the necessity and impact of the proposed approach.\n- **Lack of Fine-Tuning Capabilities:** Universal models are only widely useful if they can be fine-tuned to specific systems. This is common practice in other domains (e.g., stable diffusion for images or MATTERSim, MACE-MP0 for materials), yet the authors provide no evidence that their model supports such fine-tuning. Examination of the code suggests that this feature is not even planned.\n---\n### Minor Issues and Suggested Revisions\n1. **Dataset Clarifications:** Clearly state that MPTrj is the target training dataset in the main text (Page 6, Line 169).\n2. **Consistency in Units:** Ensure uniform usage of units, likely V/Å, in all relevant sections (e.g., Page 7, Line 221).\n3. **Charge Values for Water Molecule:** Provide explicit charge values from Figure 2d for easier comparison, similar to what is done in Figure S1a.\n4. **Thermostat Choice**: The authors use a Nosé-Hoover thermostat in their molecular dynamics simulations, which does not properly sample the canonical ensemble for small systems and can artificially suppress fluctuations. A stochastic thermostat (e.g., Langevin or Bussi-Donadio-Parrinello) would be more appropriate to ensure correct statistical sampling."}, "author_response": "Response to the reviewers\nWe sincerely thank the three reviewers for their time and invaluable comments. We have\nrevised our manuscript in accordance with their advice and criticisms and have clarified the\nissues raised by the reviewers. We have provided point-by-point responses to the questions\nraised by the reviewers in this document. The modifications to the manuscript are highlighted\nin dark red.\nTable of Contents\nResponse to reviewer #1 ............................................................................................................ 2\nResponse to reviewer #2 .......................................................................................................... 13\nResponse to reviewer #3 .......................................................................................................... 30\n1/43\nResponse to reviewer #1\nQ1. “The authors present a Machine Learning (ML) potential where the short-range graph\nneural network prediction is combined with the long-range electrostatic and dispersion\ninteractions. While similar approaches were proposed previously, the novelty of this\nmanuscript is in using the polarizable charge equilibration (PQEq) method instead of the charge\nequilibration (QEq) method. In addition, the authors do not fit the partial charges and thus avoid\nproblems with charge partitioning schemes. Tackling long-range interactions with ML\npotentials is an active research area, and this manuscript presents an important contribution to\nthe field. Nevertheless, the authors should address the following points before publication.”\nAuthors’ reply:\nWe would like to thank the reviewer for the thorough review and kind comments. We will\ncarefully address all the points raised and incorporate the necessary improvements in our\nrevised manuscript. Your expert insights are invaluable in enhancing the quality of our\nmanuscript.\nQ2. “In message-passing graph neural networks, the graph cut-off should not be confused with\nthe cut-off of classical force fields since each message passing/interaction layer increases the\nreceptive field. The authors use a NequIP model with six equivariant layers. Thus, the effective\ncut-off is not 5 Å but can extend up to 6 * 5 Å. This point should be made clear in the text, i.e.,\nin the introduction where the authors discuss the “MLIPs with a cut-off of around 5 Å” and in\nthe results section when the authors are comparing their results with strictly local ML potentials\nsuch as HDNNP, where the employed cut-off for the symmetry functions is the same as the\neffective cut-off.”\nAuthors’ reply:\nThank you for your valuable comment regarding the interpretation of cut-off distances in\nmessage-passing graph neural networks. We will revise our manuscript to clarify this\ndistinction in the following ways.\nIn the Introduction section, we now explicitly explain how message-passing mechanisms\ncan extend the cut-off through layer-wise propagation. This addresses the potential confusion\nbetween our nominal cut-off and the actual effective receptive field of message-passing graph\nneural network potential. Subsequently, we further demonstrate that implicit message-passing\nmechanisms may also face challenges in two molecules separated by a distance beyond the cut-\noff value,\nWhile existing local MLIPs with a cut-off of around 5 Å perform well in simulating\ninteractions within localized chemical environments, they may fail to capture long-range\nphenomena. This hinders their ability to understand and elucidate the behaviours of complex\nmaterials [22]. One approach is to implicitly incorporate the effects of long-range interactions\nthrough message-passing mechanisms [23] which can extend the cut-off through layer-wise\n2/43\npropagation. However, if certain parts of the system are disconnected on the graph, such as\ntwo molecules separated by a distance beyond the cut-off value, the message-passing scheme\nbecomes ineffective [24]. A promising direction involves investigating the potential\nadvantages of explicitly including long-range interactions in the model formulation.\nAnd in the Results section,\n… We demonstrate the advantage of integrating equivariant message-passing networks with\nexplicitly polarizable long-range interaction models by benchmarking against both 4G-\nHDNNP, a local MLIP with explicit long-range interactions, and NequIP, which implicitly\ncaptures long-range effects through its message-passing mechanism…\nQ3. “Related to the above, the test cases developed by Ko et al. involve very small systems\nand make sense for strictly local ML potentials but make little sense with a 6-layer graph neural\nnetwork (GNN), where the receptive field likely encompasses the whole system. The identify\nthe PQEq contribution, the authors should extend Table 1 to include results using only GNN\nand using GNN+DFT-D3. On the other hand, a comparison with 2G-HDNNP and 3G-HDNNP\nis not necessary since these results were already published.”\nAuthors’ reply:\nThank you for this insightful comment which raises an important methodological\nconsideration about our test systems and evaluation approach. We agree that clearer isolation\nof the PQEq contribution would strengthen our analysis. Following your suggestion, we have\nextended Table 1 to include results using only the GNN component (Maruf’s NequIP,\narXiv:2503.17949 https://arxiv.org/html/2503.17949v1).We do not include D3 correction to\nany of the methods because the dataset does not incorporate it. This additional analysis will\nprovide a more comprehensive picture of how each component contributes to the model’s\nperformance and will help readers better understand the specific benefits of the PQEq\ncomponent in our approach.\nRegarding the comparison with 2G-HDNNP and 3G-HDNNP, we appreciate your guidance,\nand we have removed these comparisons from Table 1 to focus on the novel aspects of our\nwork.\nIn the Results section,\nTable 1. The comparison of root means square test error metrics on the different charge-\nstate datasets.\nDataset 4G-HDNNP Maruf’s NequIP This work\nEnergy (meV/atom) 1.19 1.33 0.44\nC H /C H +\n10 2 10 3\nForce (eV/Å) 0.078 0.071 0.023\nEnergy (meV/atom) 0.48 N. A. 0.23\nNa Cl +\n8/9 8\nForce (eV/Å) 0.033 N. A. 0.006\n3/43\nDataset 4G-HDNNP Maruf’s NequIP This work\nEnergy (meV/atom) 1.32 498.39 4.87\nAg +/-\nEnergy (meV/atom) 0.22 1.03 0.14\nAu -MgO\nEnergy (meV/atom) N. A. 0.48 0.30\nBTA-Cu\nForce (eV/Å) N. A. 0.008 0.004\nEnergy (meV/atom) N. A. 0.71 0.19\nBTA (H O)-Cu\n…\nAs for the charged Ag clusters, due to their small system size and lacking long-range effects,\nthe entire structure falls within the typical cutoff radius used by the MLIP. The energy\ndependence on the overall charge state of clusters leads to degeneracy between atomic\nstructures and potential energy surfaces, resulting in poor performance of charge-independent\nNequIP model. We acknowledge that for such small, highly charged systems, 4G-HDNNP\napproach to learning charge equilibration parameters is more convenient. Nonetheless, we also\nachieved reduced force error with the fixed physical parameters compared to 4G-HDNNP. We\nre-optimized the geometries of Ag +/− using our model, with results shown in Fig. S1 c. Our\n3\nmodel accurately predicts the structures, with root mean squared displacements of only 0.002\nÅ and 0.006 Å compared to DFT-optimized structures.\n…\nOur comprehensive evaluation demonstrates that explicit incorporation of physical long-\nrange interactions significantly enhances the performance of MLIPs across diverse charge-state\nsystems. Notably, our proposed framework outperforms 4G-HDNNP in most cases for energies\nand forces without fitting partial charges. This advantage likely stems from a more physically\nmeaningful partition of the total energy.\nQ3. “It is unclear if the resulting partial charges should be interpreted as physical partial\ncharges, or a quantity correlated with the partial charge. The authors should report the errors\nfor partial charge prediction for Ko et al. cases. Of course, the agreement cannot be as good as\nwhen training on a partial charge, and given the differences between the charge partitioning\nschemes, this is also not the goal. Nevertheless, the authors should compare the partial charge\nprediction errors with differences between charge partitioning schemes to clarify the meaning\nof predicted partial charges.”\nAuthors’ reply:\nThank you for your insightful comment regarding the interpretation of partial charges in our\nmodel. The partial charges in our model should be interpreted primarily as parameters that\n4/43\neffectively reproduce the electrostatic contribution to the total energy while maintaining\nphysically reasonable characteristics.\nIn the Results section,\n… Although PQEq partial charges are not directly comparable to DFT Hirshfeld ones, we\ncompare them to show their relevance. Varying degrees of agreement across different\nchemical systems are observed as shown in Table S2…\nAnd in SI,\nWhen comparing our PQEq partial charges with reference DFT Hirshfeld charges, we\nobserve varying degrees of agreement across different chemical systems as shown in Table S2.\nWhile some systems show remarkably agreement, particularly carbon-based systems and\nmetal-oxide interfaces, it is important to acknowledge that direct numerical comparisons\nbetween different charge partitioning schemes have inherent limitations. The different\npartitioning philosophies mean that perfect agreement is neither expected nor necessarily\ndesirable as a validation metric.\nTable S2. The mean absolute difference of DFT Hirshfeld charges and PQEq partial charges\nfor the different charge state dataset, expressed in e.\nC H /C H + Na Cl + Ag +/- Au -MgO\n10 2 10 3 8/9 8 3 2\nDifference 0.0182 0.396 0.210 0.0364\nFurthermore, we would like to clarify that comparing the absolute values of partial charges\nhas limited significance. Instead, the relative trends and patterns in charge distribution can\nmeaningfully reflect the consistency of physical models. When applying an electric field\nperpendicular to the molecular plane of water molecules, we observed that the changes in PQEq\npartial charges closely track those of DFT Hirshfeld charges, despite differences in their\nabsolute values. Rather, it is the response of these charges to perturbations and their collective\ncontribution to system properties that provides more meaningful validation of our model’s\nphysical consistency.\nIn the Results section,\n… We also compared the response of O atom partial charges across DFT Hirshfeld\ncalculations, our model (PQEq model), and the conventional QEq model as shown in Fig. S4.\nExcellent agreement between PQEq model and DFT is observed, while QEq model\ncompletely lacks response…\nAnd in SI,\nIn the polarization of static water molecules under an electrostatic field shown in Fig. 2b,\nthe charge equilibration (QEq) model demonstrates limitations in handling intramolecular\npolarization when electric field strength varies in the x-axis. The QEq model represents atoms\nas either point charges or Gaussian charges, which restricts its ability to accurately capture\nmolecular polarization effects. In contrast, the PQEq model separates atomic cores and shells\n5/43\n(Eq. S6), enabling it to effectively account for intramolecular polarization phenomena. To\nillustrate this difference, we present in Fig. S4 the variations in O atom partial charge\ndifferences (𝑞 −𝑞 ) across different models as a function of the applied electric field\n!\"#\nstrengths. It is important to note that our analysis does not directly compare the absolute values\nof partial charges, but rather employs relative values as a reference point. Under different\ncharge partitioning schemes, comparing absolute values lacks meaningful interpretation,\nwhereas the relative changes more accurately reflect how various methods respond to the\napplied electrostatic fields. Fig. S4 demonstrates that our model shows consistency with first\nprinciple calculations in terms of charge variations. This comparison validates that our\napproach accurately captures the electronic response behaviour under various conditions.\nFig. S4 The variations in O atom partial charge differences (𝑞 −𝑞 ) between DFT\n!\"#\nHirshfeld calculations, our model (PQEq model), and the conventional QEq model.\nQ5. “For the universal machine learning potential model with and without polarizable long-\nrange interactions, the energy/force/stress errors on the test set are very similar. The authors\nshould include the results of the universal model without polarizable long-range interactions\nalso for the reported applications, i.e., Mechanical properties, Ionic diffusivity etc.”\nAuthors’ reply:\nThank you for this insightful comment regarding the comparative performance of our\nmodels. While our universal models with and without polarizable long-range interactions show\nsimilar error metrics on the test set, the model incorporating long-range interactions\nconsistently achieves lower errors. Following your suggestion, we extend our analysis to\ninclude results from the universal model without the interactions (w/o-lr model) for the reported\napplications, where the superiority of our model with the interactions still holds.\nFor the mechanical properties, in the Results section,\n… Fig. 2 e-g illustrate the comparison of the bulk modulus 𝐵 determined by the pretrained\nmodel and DFT. Our model demonstrated impressive performance, achieving an R2 of 0.94\nwith the Voigt approach [49] 𝐵 and Hill average method [50] 𝐵 , and 0.93 using the Reuss\n$ $%&\n6/43\nmethod [51] 𝐵 . In comparison, as shown in Fig. S2, the w/o-lr model achieved lower R²\n%\nvalues of 0.89, 0.86, and 0.88 for 𝑩 , 𝑩 , and 𝑩 , respectively.\n𝑽 𝑹 𝑽𝑹𝑯\nc d e\nFig. 2 c-e, Comparison of bulk modulus benchmarking from first principles calculations and\nour model: c, Voigt approach, d, Reuss method, and e, Hill average.\nAnd SI detailed the bulk modulus predicted by the w/o-lr model,\na b c\nFig. S2 Comparison of bulk modulus from first principles calculations and the model\nwithout polarizable long-range interactions (w/o-lr model): a, Voigt approach, b, Reuss method,\nand c, Hill average.\nFor the ionic diffusivity, in the Results section,\n… Owing to the efficiency of our pretrained model, we are able to conduct MD simulations\non the c-LLZO comprising 64 formula units, with a duration of 2 ns for each temperature range\nfrom 800 K to 1800K. Fig. 3 a presents the Arrhenius plots derived from our model, w/o-lr\nmodel, and AIMD [52] simulations. The diffusion coefficients were determined from the slope\nof the logarithmic mean square displacements (MSD) versus logarithmic time within the\nFickian regime [54]. MSD curves of our model are depicted in Fig. 3 b. Compared to previous\nwork using AIMD simulations, our model clearly reproduces the diffusion coefficients and\nactivation energy. The performance surpasses that of the w/o-lr model which lacks explicit\nlong-range interactions.\n7/43\na b\nFig. 3 a, Crystal structure of Ia3(d Li La Zr O and Arrhenius plots depicting the lithium-\n7 3 2 12\nion diffusion coefficients across varying temperatures. The dark blue polyhedron signifies La\nlocated at the 24(c) site and the light brown polyhedron indicates Zr at the 16(a) site. Li fraction\noccupies the 24(d) and 96(h) sites. Predicted diffusion coefficients of AIMD [51], our model,\nand w/o-lr model complete with error bars, are presented to calculate activation energies. b, 2-\nns mean square displacements using our model of lithium-ion in c-LLZO with different\ntemperatures ranging from 800 K to 1800 K in an increment of 200K. The linear dashed grey\nlines, with a slope of 1, are also plotted.\nQ6. “Figure 2 b-c: the authors should increase the x-axis range to 12 Å to actually see the long-\nrange scaling behaviour.”\nAuthors’ reply:\nWe appreciate the reviewer’s valuable suggestion to extend the x-axis range. However, we\nwould like to clarify two important points regarding our current presentation. First, our primary\nfocus in these figures is to demonstrate the behaviour beyond the graph cut-off distance of 5\nÅ, specifically highlighting how the model without long-range interactions (w/o-lr model) fail\nto distinguish between Na-Na and Na-Cl dimer systems. Second, even DFT, our target, fails to\ncorrectly describe molecular dissociation limit [Chem Rev 112, 289 (2012)\nhttps://doi.org/10.1021/cr200107z].\nWe note that our choice of upper limit is consistent with other recent publications in the\nfield. For example, a recent study by Shaidu et al. [npj Comput Mater 10, 47 (2024)\nhttps://doi.org/10.1038/s41524-024-01225-6] similarly presents DFT calculations only up to\n5.3 Å when demonstrating their QEq-based model’s capability to differentiate dimer\ninteractions (Fig. 5 of that reference).\nGiven these constraints, we hope the reviewer acknowledge that the current range can\nadequately demonstrate the critical difference between the models while maintaining scientific\nintegrity regarding our reference data.\n8/43\nQ7. “Are the model training details described in the Methods section the same for the universal\nequivariance neural network potential and the test cases developed by Ko et al.? If not, please\nadd training details for Ko et al. cases.”\nAuthors’ reply:\nWe thank the reviewer for pointing out this important issue. Indeed, we used different\nstructures in our initial manuscript. We now provide all the details in Table S3, SI.\nAs in Q3, we have included results on Ko’s cases from a model containing only the GNN\ncomponent, taken from the referenced work [arXiv:2503.17949\nhttps://arxiv.org/html/2503.17949v1]. For comparison, we retrained our long-range interaction\nmodel using the same number of GNN layers as specified in that reference to ensure a fair\ncomparison.\nIn the Results section,\n… The results are presented in Table 1, with detailed neural network architectural\nspecifications provided in the Table S3.\nAnd in SI,\nAs for training the different charge state dataset, the neural network architectures are\nprovided in Table S3. For feature multiplicity, we consistently employed 32 across all models\nto maintain all parameters in accordance with the reference [6].\nTable S3. The architectures for training different charge state dataset. l is the maximum\nmax\nrotation order.\nDataset Settings Maruf’s NequIP This work\nCutoff (Å) 5.00 5.00\nC H /C H + Number of layers 8 8\n10 2 10 3\nl 1 1\nmax\nCutoff (Å) N.A. 5.00\nNa Cl + Number of layers N.A. 5\n8/9 8\nl N.A. 2\nmax\nCutoff (Å) 5.29 5.29\nAg +/- Number of layers 8 8\n3\nl 2 2\nmax\nCutoff (Å) 5.50 5.50\nAu -MgO Number of layers 6 6\n2\nl 2 2\nmax\nCutoff (Å) 5.00 5.00\nBTA-Cu Number of layers 6 6\nl 1 1\nmax\nCutoff (Å) 5.00 5.00\nBTA (H O)-Cu Number of layers 6 6\n2\nl 1 1\nmax\n9/43\nWe appreciate the reviewer’s attention to this important methodological detail, which allows\nfor a more rigorous comparison between the different model architectures. This expanded\ninformation provides readers with a clearer understanding of the modelling approaches and\nstrengthens the validity of our comparative analysis.\nQ8. “‘... MPtrj Dataset sourced from Materials Projects as the pretraining dataset ...’ Why do\nyou call this pretraining? which dataset was used in the subsequent training?”\nAuthors’ reply:\nThank you for raising this terminology issue. Our use of the term “pretraining” was\nintended to convey that we have developed a universal model that can be directly applied for\nmolecular dynamics simulations and other applications without further training. To clarify, the\nmodel trained on MPtrj serves as a general-purpose, ready-to-use potential. In the revised\nmanuscript we all use “training”.\nIn the Methods section,\nTo train the universal equivariance neural network potential, we used the MPtrj dataset [16]\nsourced from Materials Projects [48] as the training dataset…\nAdditionally, in the revised manuscript, we include examples of fine-tuning, which show\nhow this pretrained model can be further refined for more specialized tasks requiring higher\nprecision. This illustrates the versatility of our approach: the model can be used as-is for many\napplications but can also serve as a starting point for more specific applications through fine-\ntuning.\nIn the Results section,\n… Our model also exhibits remarkable adaptability through fine-tuning on high-\ntemperature datasets, with performance metrics detailed in Table S6.\nAnd in SI,\n… To evaluate our model’s capacity to learn high-temperature data representations, we fine-\ntuned our pretrained universal model on the MPF-TP dataset. During this fine-tuning process,\nonly the final two convolutional layers and the multilayer perceptron layers were set as\ntrainable parameters, while all remaining layers were frozen.\nOur fine-tuning results demonstrate remarkable improvements on high-temperature, non-\nequilibrium systems as shown in Table S6. The fine-tuned model achieves state-of-the-art\naccuracy across all prediction targets on the MPF-TP dataset…\nQ9. “Why is the shell charge of all atoms set to 1?”\nAuthors’ reply:\n10/43\nThe shell charge value of 1 for all atoms was chosen following established precedent in the\nfield of original PQEq paper [J. Chem. Phys. 146, 124117 (2017)\nhttps://doi.org/10.1063/1.4978891]. In the Method section of this paper, subsection A\n“Polarizable charge equilibration (PQEq) model”, they explicitly state: “…A shell charge\nof Z = 1 is used for all atoms.”.\ni\nWe add the reference in SI,\n… Here 𝑧 is set to 1 for all atoms [1].\nThis standardized value simplifies the model while still adequately capturing polarization\neffects across different chemical environments. The physical interpretation is that the model\nadjusts the relative displacement of the shell rather than varying the shell charge itself,\neffectively modulating the induced dipole moment through positional shifts rather than charge\nmagnitude changes.\nQ10. “Missing directly related reference: 10.48550/arXiv.2501.19179”\nAuthors’ reply:\nThank you for bringing this reference to our attention. This paper combines the QEq method\nwith equivariant message-passing graph neural networks, which is indeed related to our\napproach. We have added this reference as suggested to provide readers with a more complete\nview of the current landscape in this field. We appreciate the opportunity to acknowledge this\nrelevant work and position our research within the broader context of recent developments in\ncombining charge equilibration methods with machine learning frameworks.\nIn the Introduction section,\n… In these methods, charges are usually obtained through charge equilibration (QEq) [25]\nprinciples, including the Charge Equilibration via Neural Network Technique (CENT) method\n[26], Fourth-Gneration [27, 28] High-Dimensional Neural Network Potential (4G-HDNNP),\nand Charge Equilibration Layer for Long-range Interactions (CELLI) [29] method.\n29. Fuchs, P., M. Sanocki, and J. Zavadlav, Learning Non-Local Molecular Interactions via\nEquivariant Local Representations and Charge Equilibration. arXiv preprint\narXiv:2501.19179, 2025.\nQ11. “The authors should relate their work with alternative long-range approaches not\nemploying the charge equilibration schemes, e.g., Ewald summation, Latent Ewald Summation,\nLSR-MP, NeuralP3M, and So3krates”\nAuthors’ reply:\n11/43\nWe appreciate the suggestion to relate our work with alternative long-range approaches that\ndo not employ charge equilibration schemes. We have added them in the Discussion section to\naddress this importance.\n… Additionally, several approaches that do not require charge equilibration process provide\nvaluable insights. Latent Ewald Summation [24] captures long-range interactions by learning\nhidden variables from local descriptors and applying Ewald summation to these variable. LSR-\nMP [79] employs a fragmentation-based approach with hierarchical message passing between\natoms and fragments to learn the long-range interactions. NeuralP³M [80] introduces mesh\npoints alongside atoms to discretize space and capture long-range effects. SO3krates [81]\nutilizes spherical harmonic coordinates and equivariant attention for interactions at arbitrary\nlength scales. In contrast to these approaches, our universal model explicitly simulates charge\nredistribution through a physics-driven framework that maintains a transparent connection to\nthe underlying physical principles of charge polarization and electrostatic interactions. This\ndesign choice enhances interpretability while preserving computational efficiency. A\npromising direction for future work lies in exploring strategic integrations of our approach with\nthese complementary methodologies to develop even more comprehensive and powerful\nuniversal models.\nQ12. “Remarks on code availability:\nI did not run the code but it looks OK.”\nAuthors’ reply:\nWe appreciate the reviewer’s assessment of our code. We have open-sourced the complete\nmodel codebase in the new repository and supplemented it with additional tutorials and fine-\ntuning examples to facilitate implementation for users.\nCode availability\nThe JAX implementation code is available through https://github.com/reaxnet/reaxnet.\n12/43\nResponse to reviewer #2\nQ1. “Overall, a manuscript which should be very heavily revised - but that might decrease any\nnovelty reported such that it would fall outside of the scope of the journal. Maybe best to reject.”\nAuthors’ reply:\nThank you for your feedback on our manuscript. We deeply respect your thoughtful review\nand have taken your comments very seriously. We have made substantial revisions to address\nthe concerns you raised, and we believe the revised manuscript still demonstrates adequate\nnovelty.\nQ2. “- What are the noteworthy results?\nThere are some interesting simulations using the MLIP technology developed by the authors -\nbut there is no reliable benchmarking carried out to compare to other state-of-the-art methods.\nSince the focus of the manuscript is a supposedly superior method, this is a claim that should\nhave been proven by the authors.”\nAuthors’ reply:\nWe appreciate the reviewer’s feedback regarding the need for comprehensive benchmarking\nagainst state-of-the-art methods. In response to this important point, we have significantly\nexpanded our comparative analysis in the revised manuscript.\nWe have conducted additional analyses, comparing our performance with CHGNet and\nMACE-MP-0, both of which underwent similar training on the zero-temperature MPtrj dataset.\nThese comparative evaluations were performed across the benchmark datasets proposed in\nMatterSim. Our findings demonstrate that our model consistently achieves lower force errors\nacross all datasets compared to these two universal models, while maintaining reasonable\nenergies and stresses errors.\nIn the Results section,\nWe also perform benchmarking of our model with other universal models on the test sets\nprovided by the reference [19]. Table S5 demonstrates the exceptional performance of our\nmodel across diverse test sets in force prediction capabilities, consistently outperforming both\nCHGNet [16] and MACE-MP-0 [18], while maintaining reasonable energy and stress\npredictions. We excluded MatterSim [19] results from our comparative analysis as its training\ndataset incorporates high-temperature data.\nAnd in SI,\nWe also conducted benchmarking with other universal models on the test sets provided by\nthe reference [11]. For near-equilibrium evaluations, two randomly sampled collections are\nused: 1,000 structures from the MPtrj dataset [7] (MPtrj-random-1k) and 1,000 from\n13/43\nAlexandria dataset [12] (Alexandria-1k). These assess model accuracy on energetically\nfavourable configurations, which is crucial for predicting material stability. Three increasingly\nchallenging thermodynamic perturbation (TP) test sets evaluate performance under varied\ntemperature and pressure conditions with non-equilibrium atomic configurations. MPF-Alkali-\nTP focuses on potential ionic conductors, comprising 50 randomly selected compounds\ncontaining at least one alkali metal paired with electronegative elements (N, O, P, S, Se, or\nhalogens) from MPF2021 dataset [13]. This dataset specifically tests model accuracy for\nmaterials with significant ionic character. MPF-TP broadens the scope with 50 randomly\nselected compounds from the MPF2021 dataset without elemental restrictions, presenting a\nmore diverse chemical challenge. Random-TP represents the most demanding test, constructed\nby arbitrarily positioning 20 atoms of random elements within simulation boxes, creating\nhypothetical structures that test extreme generalization capability. Table S5 demonstrates the\nperformance of our model across diverse test sets, particularly highlighting its superior force\nprediction capabilities. Our model achieves remarkable force prediction accuracy, consistently\noutperforming both CHGNet [7] and MACE-MP-0 [14] (all three models were trained on MPtrj\ndataset) across all benchmark datasets with significantly lower MAEs – nearly half the error of\ncompeting models in many cases. Impressively, our model demonstrated robust generalization\ncapability across increasingly challenging thermodynamic perturbation scenarios. Although we\nobserve slightly higher energy prediction MAEs compared to the two models, its performance\nacross the entire spectrum of materials conditions is improved particularly for force predictions\ncritical for molecular dynamics simulations. It is important to note that our comparison was\nlimited to universal models trained exclusively on zero-temperature dataset (MPtrj) and did not\ninclude MatterSim [11], which incorporates high-temperature datasets…\nTable S5. Comparison of mean absolute errors (MAE) in the test set generated by the reference\n[11].\nTest set MAE CHGNet MACE-MP-0 Our model\nEnergy (eV/atom) 0.027 0.015 0.019\nMPtrj-random-1k Force (eV/Å) 0.120 0.117 0.066\nStress (GPa) 0.290 0.468 0.271\nEnergy (eV/atom) 0.150 0.092 0.115\nAlexandria-1k Force (eV/Å) 0.108 0.095 0.064\nStress (GPa) 1.643 0.160 1.055\nEnergy (eV/atom) 0.250 1.351 0.280\nMPF-Alkali-TP Force (eV/Å) 1.636 15.819 0.717\nStress (GPa) 12.625 25.723 13.841\nEnergy (eV/atom) 0.254 256.340 0.521\nMPF-TP Force (eV/Å) 3.313 1506.854 0.768\nStress (GPa) 25.208 202.093 15.445\nEnergy (eV/atom) 0.506 9.184 3.338\nRandom-TP Force (eV/Å) 3.950 88.327 1.199\nStress (GPa) 7.230 19.224 33.282\n14/43\nFurthermore, after fine-tuning, our model can achieve even more accurate predictions on\nhigh-temperature datasets, demonstrating its adaptability to diverse thermodynamic conditions.\nIn the Results section,\n… Our model also exhibits remarkable adaptability through fine-tuning on high-temperature\ndatasets, with performance metrics detailed in Table S6.\nAnd in SI,\nOur fine-tuning results demonstrate remarkable improvements on high-temperature, non-\nequilibrium systems as shown in Table S6. The fine-tuned model achieves state-of-the-art\naccuracy across all prediction targets on the MPF-TP dataset. This suggests that while our\nuniversal model can maintain reasonable force accuracy under a wide range of conditions,\ntargeted fine-tuning is more likely to balance broad applicability with the superior accuracy\nrequired for specific applications.\nTable S6. Comparison of mean absolute errors in the MPF-TP test set of the universal\nmodels and finetuned model.\nEnergy (eV/atom) Force (eV/Å) Stress (GPa)\nCHGNet 0.254 3.313 25.208\nMACE-MP-0 256.340 1506.854 202.093\nMatterSim (M3GNet) 0.036 0.431 1.318\nMatterSim (Graphormer) 0.040 0.421 1.917\nOur universal model 0.521 0.768 15.445\nFine-tuned model 0.027 0.203 0.976\nWe believe these additional benchmarking results now provide substantial evidence to\nsupport our claims regarding the method’s capabilities.\nQ3. “- Will the work be of significance to the field and related fields? How does it compare to\nthe established literature? If the work is not original, please provide relevant references.\nThis work is mostly combining other, existing pieces of technology in a slightly different way.\nIn short, a short-range MLIP is combined with Goddard’s QEq scheme and D3 dispersion\ncorrections. The closest similar is Behler’s 4G-HDNNPs, but these don’t use the equivariant\nmessage passing short range technology. Results presented here are marginally better, but\nperhaps only an incremental improvement. Presumably due to the better and more expressive\nshort-range part. Behler’s scheme is superior as it uses environment dependent parameters for\nthe QEq component, unlike here, where these are fixed. The CENT method by Goedecker\n(which uses a similar approach to fit the energies directly) is not even mentioned.”\n15/43\nAuthors’ reply:\nWe appreciate the reviewer’s insightful comparison of our work with existing technologies\nin the field. To address the points raised, we have expanded our benchmarking analysis to better\ndemonstrate the novel contributions of our approach.\n(1) For the improvement over existing works, to show that it is not just because of a more\nexpressive short-range part, we compared to a very recent arXiv work [arXiv: 2503.17949v1\nhttps://arxiv.org/html/2503.17949v1], which used the same NequIP neural network as us, yet\nsimilar long-range interaction scheme like 4G-HDNNP (marked as arXiv-lr work). We\nextracted their results and compared to ours as shown in Table R1.\nTable R1 Energy (meV/atom) and force (eV/Å) root mean squared errors of the arXiv-lr work\nand our work.\nDataset arXiv-lr work This work\nEnergy (meV/atom) 0.84 0.44\nC H /C H +\n10 2 10 3\nForce (eV/Å) 0.055 0.023\nEnergy (meV/atom) 0.82 4.87\nAg +/-\nEnergy (meV/atom) 0.17 0.14\nAu -MgO\nEnergy (meV/atom) 0.43 0.30\nBTA-Cu\nForce (eV/Å) 0.002 0.004\nEnergy (meV/atom) 0.28 0.19\nBTA (H O)-Cu\nWe achieved some improvements with the fixed PQEq parameters over that work,\nindicating that our combined approach can indeed provide insights into this field.\n(2) Regarding the 4G-HDNNP/CENT methods, which learn environment-dependent QEq\nparameters, we have provided a detailed discussion in the Discussion section. While employing\nan additional neural network to learn environmental parameters could potentially achieve\nhigher accuracy, efficiency remains critically important for universal models designed for\nlarge-scale molecular dynamics simulations. The literature also indicates that QEq parameters\ndo not change significantly during molecular dynamics simulations. A key focus of our future\nwork will be combining efficiency and accuracy to develop even more generalizable models\nthat maintain computational performance while further enhancing predictive capabilities.\nIn the Discussion section,\nWe acknowledge that our model aims to achieve a balance between computational\nefficiency and accuracy. Several advanced machine learning approaches for long-range\ninteractions may achieve higher accuracy for specific, complex systems, including higher-order\ncharge expansion methods and simultaneous optimization of QEq parameters [32, 79] in each\n16/43\nstep. However, the more accurate methods incur higher computational complexity. The method\nthat dynamically solve charge equilibration parameters at each step require additional cost for\nMD simulations. It is also observed that the charge equilibration parameters do not change\nsignificantly during MD simulations for that method [79], suggesting that expanding\npolarization charges to higher orders while using fixed, a priori parameters may achieve higher\naccuracy with tolerable computational overheads.\nFor reference, the dynamic changes of QEq parameters (𝜒 and 𝜂) during molecular\ndynamics simulations in the reference [arXiv:2411.04062 https://arxiv.org/html/2411.04062v1]\nare reproduced in Fig. R1,\nFig. R1 Variation with time of 𝜒 and 𝜂 at the C, H and O sites through molecular dynamics.\nThe figure is extracted from arXiv:2411.04062.\nFor these parameters in different systems, fixed parameters work generally well (as Q2\ndetailed and shown in Table S5). For extreme cases (small, charged structures like Ag clusters),\nwe admit adjusting parameters would give better results, as in the revised manuscript.\nIn the Results section,\n…As for the charged Ag clusters, due to their small system size and lacking long-range\neffects, the entire structure falls within the typical cutoff radius used by the MLIP. The energy\ndependence on the overall charge state of clusters leads to degeneracy between atomic\n\nstructures and potential energy surfaces, resulting in poor performance of charge-independent\nNequIP model. We acknowledge that for such small, highly charged systems, 4G-HDNNP\napproach to learning charge equilibration parameters is more convenient. Nonetheless, we\nalso achieved reduced force error with the fixed physical parameters compared to 4G-\nHDNNP. We re-optimized the geometries of Ag +/− using our model, with results shown in Fig.\nonly 0.002 Å and 0.006 Å compared to DFT-optimized structures…\n(3) In the Introduction section, the CENT method is mentioned as it points out that learning\npartial charges is not necessary,\n17/43\n… In these methods, charges are usually obtained through charge equilibration (QEq) [25]\nprinciples, including the Charge Equilibration via Neural Network Technique (CENT)\nmethod [26], Fourth-Generation [27, 28] High-Dimensional Neural Network Potential (4G-\nHDNNP), and Charge Equilibration Layer for Long-range Interactions (CELLI) [29]\nmethod. … The ambiguity in DFT-assigned charges suggests that directly learning of them\nmay be inessential for constructing accurate interatomic potentials [26, 31-32] …\nQ4. “- Does the work support the conclusions and claims, or is additional evidence needed?\nA lot of data which the results are compared to here are cherrypicked (there is hardly or no any\nmention of significant universal force fields, such as MACE-MP and its later improvements).\nThis is a quickly moving field.”\nAuthors’ reply:\nWe appreciate your advice. As mentioned above, we have supplemented our benchmarking\nwith CHGNet and MACE-MP-0, which were also trained on zero-temperature datasets, using\nthe challenging test sets proposed by MatterSim. We did not compare with MatterSim itself\nbecause it includes high-temperature datasets in its training process, which would not allow for\na fair comparison with our model. Notably, in these additional benchmark tests, our model\ndemonstrates lower force errors than both CHGNet and MACE-MP-0 while maintaining\nreasonable energy and stress errors, further validating the effectiveness and advantages of our\napproach.\nIn SI,\nWe also conducted benchmarking with other universal models on the test sets provided by\nthe reference [11]. For near-equilibrium evaluations, two randomly sampled collections are\nused: 1,000 structures from the MPtrj dataset [7] (MPtrj-random-1k) and 1,000 from\nAlexandria dataset [12] (Alexandria-1k). These assess model accuracy on energetically\nfavourable configurations, which is crucial for predicting material stability. Three increasingly\nchallenging thermodynamic perturbation (TP) test sets evaluate performance under varied\ntemperature and pressure conditions with non-equilibrium atomic configurations. MPF-Alkali-\nTP focuses on potential ionic conductors, comprising 50 randomly selected compounds\ncontaining at least one alkali metal paired with electronegative elements (N, O, P, S, Se, or\nhalogens) from MPF2021 dataset [13]. This dataset specifically tests model accuracy for\nmaterials with significant ionic character. MPF-TP broadens the scope with 50 randomly\nselected compounds from the MPF2021 dataset without elemental restrictions, presenting a\nmore diverse chemical challenge. Random-TP represents the most demanding test, constructed\nby arbitrarily positioning 20 atoms of random elements within simulation boxes, creating\nhypothetical structures that test extreme generalization capability. Table S5 demonstrates the\nperformance of our model across diverse test sets, particularly highlighting its superior force\nprediction capabilities. Our model achieves remarkable force prediction accuracy, consistently\noutperforming both CHGNet [7] and MACE-MP-0 [14] (all three models were trained on MPtrj\ndataset) across all benchmark datasets with significantly lower MAEs – nearly half the error of\n18/43\ncompeting models in many cases. Impressively, our model demonstrated robust\ngeneralization capability across increasingly challenging thermodynamic perturbation\nscenarios. Although we observe slightly higher energy prediction MAEs compared to the two\nmodels, its performance across the entire spectrum of materials conditions is improved\nparticularly for force predictions critical for molecular dynamics simulations. It is important to\nnote that our comparison was limited to universal models trained exclusively on zero-\ntemperature dataset (MPtrj) and did not include MatterSim [11], which incorporates high-\ntemperature datasets …\nTable S5. Comparison of mean absolute errors (MAE) in the test set generated by the reference\n[11].\nTest set MAE CHGNet MACE-MP-0 Our model\nEnergy (eV/atom) 0.027 0.015 0.019\nMPtrj-random-1k Force (eV/Å) 0.120 0.117 0.066\nStress (GPa) 0.290 0.468 0.271\nEnergy (eV/atom) 0.150 0.092 0.115\nAlexandria-1k Force (eV/Å) 0.108 0.095 0.064\nStress (GPa) 1.643 0.160 1.055\nEnergy (eV/atom) 0.250 1.351 0.280\nMPF-Alkali-TP Force (eV/Å) 1.636 15.819 0.717\nStress (GPa) 12.625 25.723 13.841\nEnergy (eV/atom) 0.254 256.340 0.521\nMPF-TP Force (eV/Å) 3.313 1506.854 0.768\nStress (GPa) 25.208 202.093 15.445\nEnergy (eV/atom) 0.506 9.184 3.338\nRandom-TP Force (eV/Å) 3.950 88.327 1.199\nStress (GPa) 7.230 19.224 33.282\nQ5. “- Are there any flaws in the data analysis, interpretation and conclusions? - Do these\nprohibit publication or require revision?\nPlenty of misinterpretation/misrepresentation, listed below. With very careful revision maybe\nthese can be rectified, but most likely it will further reduce any claims of superiority and/or\nnovelty.”\nAuthors’ reply:\nWe sincerely thank the reviewer for their thorough assessment and for highlighting potential\nissues with our data analysis, interpretation, and conclusions. We take these concerns very\nseriously and have carefully addressed each point raised in the detailed comments below.\nUniversal machine learning potentials with polarizable long-range interactions are indeed still\nin their nascent stage, which underscores the significant contribution of our work. Through our\nmeticulous revisions, we believe that this novelty is not diminished.\n19/43\nQ6. “- Is the methodology sound? Does the work meet the expected standards in your field?\nMostly.”\nAuthors’ reply:\nWe appreciate the reviewer’s acknowledgment that our methodology is generally sound.\nWe have further revised our manuscript to ensure that it aligns with your expectations.\nQ7. “- Is there enough detail provided in the methods for the work to be reproduced?\nYes, and the code is publicly available (I have not tested, though).”\nAuthors’ reply:\nWe appreciate the reviewer’s assessment of our code. We have open-sourced the complete\nmodel codebase in the new repository and supplemented it with additional tutorials and fine-\ntuning examples to facilitate implementation for users.\nCode availability\nThe JAX implementation code is available through https://github.com/reaxnet/reaxnet.\nQ8. “Do not anthropomorphize mathematical models (e.g. p2/l39)\n‘subtle fluctuations’ - these are not subtle (it is meaningless in this context) - they are dynamic\nfluctuations”\nAuthors’ reply:\nWe thank the reviewer for this feedback regarding our language. We have adopted this\nprecise scientific terminology as in the Introduction section,\n… These behaviours encompass electrostatic interactions, involving forces between charged\nparticles, and dispersion terms, which arise from dynamic fluctuations in electron distribution\nwithin molecules or atoms.\nQ9. “Terms in Eq. 1 should be explained immediately after in the text.”\nAuthors’ reply:\nWe thank the reviewer for this suggestion regarding the clarity of our presentation. We have\nreorganized the structure of our manuscript to ensure that all terms in Eq. 1 are explained in\nthe text immediately following the equation.\nIn the Results section,\n20/43\n…The potential energy is expressed as the sum of the second-order expansion [30] with\nrespect to charge fluctuations and the density functional theory [39] D3 (DFT-D3) van der\nWaals dispersion energies correction term [40, 41] 𝐸 ,\n*+\n1 1\n𝐸 = /0 𝐸0(𝒓 ,𝑍 )+𝜒/𝑞 + 𝜂/𝑞1 + 𝐾0|𝒓 −𝒓 |1<\n,-. / 0 0 0 0 2 0 0 2 2 03 02\n0\n+ / 𝐶 >𝒓 ?𝑞 𝑞 +𝐸\n04,67 04,67 04 67 *+\n04867\n= /𝐸0 +𝐸 +𝐸 .\n/ 9:;< *+\n0\n(1)\nThe zeroth-order atomic energy 𝐸0 corresponds the final layer scalar output of the\n/\nequivariant graph neural networks, while the higher-order self- and interatomic Coulomb\ninteractions represent the polarizable long-range electrostatic interactions 𝐸 .\n9:;<\nTo account\nfor charge transfer and polarization effects, partial atomic charge 𝑞 of atom i is the sum of\n0\nthe nuclear core charge 𝑞 and shell charge 𝑞 , both of which assume a Gaussian charge\n03 02\ndistribution form. The first-order coefficients 𝜒/ are electronegativities, commonly defined as\nhalf of the sum of ionization potential (IP) and electron affinities (EA). The second-order\ncoefficients 𝜂/ signifies idempotential or chemical hardness, defined as 𝐼𝑃−𝐸𝐴. The spring\nconstant 𝐾0 denotes the isotropic harmonic connectivity between the shell position 𝒓 and\n2 02\ncore position 𝒓 (equal to 𝒓 ) of atom i. The Gaussian electrostatic energy is given by\n03 0\n𝐶>𝒓 ?𝑞 𝑞 , where i and j are the atomic indices, and k and l represent the core (c) or shell\n04,67 04 67\n(s), respectively. The displacement vector is given by 𝒓 = 𝒓 −𝒓 . The detailed\n04,67 04 67\nderivation of the atomic force and its response to the external electric field can be found in the\nSI…\nQ10. “Ref. 24 is incomplete”\nAuthors’ reply:\nWe sincerely thank the reviewer for bringing this bibliographic issue to our attention. We\nhave added the missing page number to the reference,\nNaserifar, S., et al., Polarizable charge equilibration model for predicting accurate\nelectrostatic interactions in molecules and solids. The Journal of Chemical Physics, 2017.\n146(12): p. 124117.\nQ11. “Eqs. 2 and 3 do not follow. Since the energies also depend on the charges, which are\ndetermined self-consistently, it is far from obvious that eqs. 2 and 3 are true. There may be an\nexpression analogous to the Hellmann-Feynman theorem but not shown here.”\nAuthors’ reply:\n21/43\nThank you for raising this important point. Indeed, the energy does depend on the self-\nconsistent determination of charges, whereby atomic types and positions determine both the\ncore charges and shell displacements. We understand your concern that when calculating forces,\nwe shall compute the negative gradient of energy with respect to position, while automatic\ndifferentiation might calculate the total derivative, causing the forces to be affected by the\nself-consistency process.\nTo be specific, in our implementation,\n𝐸 = 𝐸 (𝒓 ;𝒓 ,𝑞 ) (R1)\n,-. ,-. 03 02 0\n𝑭 = −𝜕 𝐸 , (R2)\n0 =!\" >?@\nwhere 𝒓 ,𝑞 both depend on 𝒓 and are solved self consistently. In our code, we explicitly stop\n02 0 03\nthe gradient updates with respect to 𝒓 during this process. This approach ensures that the\n03\nforces are calculated correctly after the charge distribution has been converged. For the\nstress (Eq. 3 in the main text), similar argues hold.\nIn Results section, we have refined our language,\n… . After partial charges are obtained self-consistently, the forces on atoms 𝑭 and stress\n𝒊\n𝑺 on the cell can be calculated via automatic differentiation of potential energy (Eq. 1) with\nBC\nrespect to atomic coordinates 𝒓 and strain 𝒔. Note that the gradient updates are stopped in the\n!\ncharge equilibration process to ensure partial derivatives are taken …\nWe also provided an example as in examples/basic.ipynb, section 5. Define the\ntotal energy function of our github repository (https://github.com/reaxnet/reaxnet),\ndef energy_fn(positions,\nnn_nbr,\nnb_nbr,\n):\n\"\"\"\nTotal potential energy function\nArgs:\n----\npositions: jnp.array\nAtomic positions in Angstrom\n\"\"\"\n# Update the neighbor lists first\nnn_nbr = nn_nbr.update(positions)\nnb_nbr = nb_nbr.update(positions)\npe_nn = energy_fn_nn(positions, nn_nbr)\ncharges, r_shell = charges_fn(jax.lax.stop_gradient(positions),\nnb_nbr)\npe_nb = energy_fn_nb(positions, nb_nbr, r_shell, charges)\nreturn pe_nn + pe_nb, (charges, r_shell)\n22/43\nThe ```jax.lax.stop_gradient()``` function will stop gradient updates when calculating the force\nvia automatic differentiation.\nTo ensure that our equations and implementation are correct, we test the PQEq part of forces\nof ours and LAMMPS’ (https://github.com/sabernaserifar/RexPoN) on the Li PO N (mp-\n2 2\n1020019) system. The results are listed in Table R2 and Table R3.\nTable R2. Forces comparison of LAMMPS and our code for Li PO N (mp-1020019).\n2 2\nForces (eV/Å)\n𝐹EFGG9H 𝐹-IJK 𝐹EFGG9H 𝐹-IJK 𝐹EFGG9H 𝐹-IJK\nD D L L M M\nLi 0.036 0.036 -0.044 -0.044 0.001 0.002\nLi 0.036 0.036 0.044 0.044 0.001 0.001\nLi -0.036 -0.036 0.044 0.044 0.001 0.002\nLi -0.036 -0.036 -0.044 -0.044 0.001 0.001\nLi 0.036 0.036 -0.044 -0.044 0.001 0.001\nLi 0.036 0.036 0.044 0.044 0.001 0.002\nLi -0.036 -0.036 0.044 0.044 0.001 0.001\nLi -0.036 -0.036 -0.044 -0.044 0.001 0.002\nP -0.002 -0.002 0.000 0.000 0.000 0.000\nP 0.002 0.002 0.000 0.000 0.000 0.000\nP -0.002 -0.003 0.000 0.000 0.000 0.000\nP 0.002 0.003 0.000 0.000 0.000 0.000\nN -0.057 -0.057 0.000 0.000 0.012 0.012\nN 0.057 0.057 0.000 0.000 0.012 0.012\nN -0.057 -0.057 0.000 0.000 0.012 0.012\nN 0.057 0.057 0.000 0.000 0.012 0.012\nO 0.050 0.050 0.069 0.069 -0.008 -0.008\nO 0.050 0.050 -0.069 -0.069 -0.008 -0.008\nO -0.050 -0.050 -0.069 -0.069 -0.008 -0.008\nO -0.050 -0.050 0.069 0.069 -0.008 -0.008\nO 0.050 0.050 0.069 0.069 -0.008 -0.008\nO 0.050 0.050 -0.069 -0.069 -0.008 -0.008\nO -0.050 -0.050 -0.069 -0.069 -0.008 -0.008\nO -0.050 -0.050 0.069 0.069 -0.008 -0.008\nTable R3. Stresses comparison of LAMMPS and our code for Li PO N (mp-1020019).\n2 2\n𝑆 (GPa) 𝑆 (GPa) 𝑆 (GPa) 𝑆 (GPa) 𝑆 (GPa) 𝑆 (GPa)\nDD LL MM LM DM DL\nOur code 0.3586 0.8228 0.0164 0.0000 0.0000 0.0000\nLAMMPS 0.3586 0.8228 0.0164 0.0000 0.0000 0.0000\nQ12. “Panels 2a, 2b and 2c must be united.”\nAuthors’ reply:\n23/43\nWe thank the reviewer for this suggestion on improving the presentation of our results. We\nhave revised Fig. 2 by uniting panels 2a, 2b, and 2c into a single cohesive panel with a shared\naxis system.\nIn the Results section,\n… As shown in Fig. 2a, both models (our model and w/o-lr model) trained on MPtrj dataset\naccurately reproduce the equilibrium bond lengths of the two dimers. However, the w/o-lr\nmodel failed to differentiate interaction energies beyond the MLIPs cutoff distance, i.e., 5.0 Å,\nwhile our model successfully captured the DFT potential energy surface throughout the entire\ndistance range, demonstrating its capability in describing extended ionic interactions…\nFig. 2 a, The interaction energies for Na-Na (green) and Na-Cl (purple) dimers predicted by\nour model, w/o-lr model, and reference calculations using DFT.\nQ13. “The QEq curve on 2d seems completely flat - this cannot be true, or the QEq must be\nextremely rigid.”\nAuthors’ reply:\nWe appreciate you raising this question. Indeed, the QEq model curve is completely flat,\nwith only negligible floating-point errors that are virtually undetectable. We would like to\nemphasize that this example represents a static calculation of a water molecule under different\nelectric field strengths, primarily designed to demonstrate the response of atomic polarization\nwithin the molecule to external electric fields. The definition of electric potential energy can\nbe expressed as,\n𝐸 = −∑𝑞 𝝐∙𝒓 .\nN 0 𝒊\n(R3)\nThe water molecule consistently remains in the yz-plane with the electric field direction\nconsistently along the x-axis. Given that the water dipole is perpendicular to the direction of\nthe electric field, the eletrostatic energy of QEq-based model is always zero.\nThis is a fundamental limitation of the traditional QEq-based model—it lacks the capability\nto respond to external electric fields when atomic positions remain fixed. In contrast, our\n24/43\npolarizable model captures this physical phenomenon correctly, as clearly demonstrated by\ntheir non-flat energy curves in Fig. 2d (now Fig. 2b). This response highlights an important\nadvantage of our model for simulations involving external electric fields or environments with\nvarying electrostatic characteristics.\nIn the Results section,\n… Due to the external electrostatic field direction consistently being orthogonal to the\nmolecule, the electrostatic energy in the QEq-based model (like 4G-HDNNP) is 0 and is\ncompletely independent of field strength (i.e., 𝝁⋅𝝐 = 0, where 𝝁, 𝝐 are dipole moment and\nelectric field, respectively). This is because the QEq model treats atoms in the molecule as\npoint charges or Gaussian charges, without distinguishing between core and shell\ncomponents…\nQ14. “Computational cost being linear: I find that hard to believe. Long range interactions in\nperiodic systems are typically treated using some version of the Ewald summation method or\nan equivalent approach. While there exist linear scaling implementations, there is no mention\nof these here.”\nAuthors’ reply:\nWe sincerely thank the reviewer for their correction. We regretfully acknowledge our\noversight regarding the computational scaling. Due to the matrix inversion operations in our\ncode, the scaling is indeed cubic rather than linear as we initially stated. We did not measure\nthe scaling, instead, we unintentionally took the apparent trend as the y-axis upper limit was\ntoo large but the x-axis upper limit was too small. To fix this misleading presentation, we extend\nthe system sizes being tested to reflect the cubic scaling behaviour. However, we would like to\nclarify that for the systems of tens of thousands of atoms that we have simulated, the\ncomputational cost remains acceptable, and we have still achieved high computational\nefficiency. This practical performance is sufficient for the material systems and properties\ninvestigated in this work.\nFor future work, we plan to improve our algorithm by implementing more efficient methods\nfor solving the linear equation systems, such as conjugate gradient methods with appropriate\npreconditioning, which could potentially reduce the scaling behaviour. We have corrected the\nmanuscript to accurately reflect the computational scaling of our current implementation and\nto properly contextualize its practical performance for the systems studied.\nAgain, we are extremely sorry for this misrepresentation.\nIn the Results section,\nDespite the incorporation of additional long-range interaction calculations, our model\nmaintains high computational efficiency. Performance benchmarks were conducted on a single\nNVIDIA H100 GPU to quantitatively assess the computational cost. As illustrated in Fig. S6,\nfor a periodic system comprising 2160 atoms, our model requires approximately 0.07 s per\n25/43\nmolecular dynamics time step, representing a significant improvement over conventional\nuniversal MLIPs which require approximately 0.91 s. This computational advantage extends\nto larger systems, as demonstrated by simulations of a 24000-atom system, where our model\nmaintains efficiency with only 1.49 s per time step. Such computational performance enables\nnanosecond-scale molecular dynamics simulations of systems containing tens of thousands of\natoms, making it practical for more realistic modelling of materials.\nAnd in SI,\nOur model introduces explicit polarizable long-range interactions while still maintaining\ngood computational efficiency. To highlight the performance of this model, benchmarking was\ndone on molecular dynamics simulations. Supercell models were generated from Li PO N\n2 2\n(mp-1020019) from Materials Project [8] and simulations were conducted. Additionally, the\nPyTorch-based universal machine learning interatomic potentials CHGNet [7] was used for\ncomparison. All molecular dynamics simulations were carried out under the canonical\nensemble simulations with a Nosé-Hoover thermostat [16, 17] at 300 K temperature. All the\ncomputations were performed with a single NVIDIA H100 GPU.\nCHGNet\n(w/o D3)\nOur model\n(w/ D3)\nFig. S6 Benchmarking on the computational efficiency involved in molecular dynamics\nsimulations and geometry optimizations. All the calculations are performed on a single\nNVIDIA H100 GPU. The time consumption per step during 100 steps of molecular dynamics\nsimulations for our model is presented. Results for CHGNet are shown in blue dots; however,\nthe calculation for systems with more than 10 thousand atoms was halted due to out of GPU\nmemory in the calculation environment.\nAs shown in Fig. S6, we present efficiency in terms of time consumption per simulation\nstep. Although our model demonstrates O(N³) scaling due to inversing the Eq. S3, it\nsuccessfully handles large-scale systems exceeding 20,000 atoms, highlighting its applicability\nto complex materials simulations. In our studied reaction systems containing tens of thousands\nof atoms, the inversion method provides acceptable computational efficiency. Nevertheless, in\nfuture work, we plan to further refine our algorithms, such as implementing iterative solution\n26/43\nmethods to reduce complexity, thereby achieving more computationally efficient\nimplementations and lowering resource consumption.\nQ15. “The sentence ‘as well as account for the free energy surfaces under the conditions of\nfinite-temperature thermodynamics’ does not make sense. The free energy is a state function\nand not a surface.”\nAuthors’ reply:\nWe thank the reviewer for this correction regarding our terminology. We realise that this\nsentence is scientifically imprecise. Free energy is indeed a state function, not a surface. In\nResults section, we have revised as,\n…Accurately simulating phase changes in ferroelectric substances requires precise potential\nenergy functions that can respond to the small atomic shifts and structural changes, as well as\naccount for the free energy landscape under the conditions of finite-temperature\nthermodynamics [56] …\nThe term “free energy landscape” is the standard terminology for describing how free\nenergy varies across configuration space. We appreciate the reviewer’s attention to this detail,\nwhich helps improve the scientific rigour of our manuscript.\nQ16. “‘The QEq-based machine learning models inherently struggle with atomic polarization,\nleading to non-physical results under external electric fields.’ there is no evidence for that. In\nthe results quoted by the authors on 4G-HDNNPs they may be less accurate but not unphysical.”\nAuthors’ reply:\nWe thank the reviewer for this important point regarding our characterization of QEq-based\nmachine learning models. We have revised this sentence in the Disscusion section,\n… QEq-based models inherently struggle with atomic polarization because they do not\ndistinguish between core and shell components. Due to their inability to handle atomic-level\nresponses to external fields, they may exhibit non-physical behaviours in some cases as\ndescribed in the reference [47].\nBecause the QEq model assumes the system is a continuous conductor, allowing charge\ntransfer to occur at any distance, it produces problematic results for insulators or\nsemiconductors, where long-range charge transfer should not occur in response to small\ndifferences in potential. As demonstrated in the reference 47 [J. Chem. Theory Comput. 18, 1,\n580-594 (2022) https://doi.org/10.1021/acs.jctc.1c00975], the QEq model indeed predicts non-\nphysical charge distributions. We tested our model under identical conditions and found that\nwater molecules do not exhibit non-physical charge distributions.\nIn the Results section,\n27/43\n…As illustrated in Fig. S5, we conducted additional investigations of two water molecules\nunder an external electric field of 0.25 V/Å. In charge distribution models, it is essential to\nallow systems to form an internal electric field opposing the external field. The QEq model\nassumes the system behaves as a perfect conductor without penalizing charge transfer as a\nfunction of distance. This limitation of the QEq method leads to non-physical charge transfer\nbetween molecules separated by large distances [47]. In contrast, within PQEq model, the\nindividual water molecules can polarize, generating internal electric fields that counteract\nthe external field, thereby capturing more realistic electrostatic responses in molecular systems.\nUnder QEq scheme, the water molecules accumulate non-zero net charges and migrate to the\ntop and bottom of the simulation box. In contrast, our model maintains charge neutrality of the\nwater molecules, which remain stationary despite the applied electric field, aligning with DFT\nresults…\nAnd in SI,\nWe also examine the limitations of the QEq method when applied to water molecules under\nan external electric field, particularly its tendency to predict physically unrealistic charge\ndistributions. To illustrate the distinctions between QEq and PQEq charge distribution methods\nin the presence of an electric field, we performed a simple demonstration following the protocol\noutlined in reference [15]. Fig. S5 a presents visualizations of water molecules using both QEq\nand PQEq methods. Fig. S5 b and c depict the temporal evolution of corresponding z-positions\nand molecular charges within the simulation box for both methods at an electric field of 0.25\nV/Å. Under QEq, the water molecules accumulate non-zero net charges and migrate to the top\nand bottom of the box [15]. In contrast, using our model, the water molecules maintain charge\nneutrality and remain stationary despite the applied electric field. This behaviour aligns with\nDFT calculations, which similarly predict no molecular displacement for the two water\nmolecules.\na\nM2\nM1\nb c\nM1\nM2\n28/43\nFig. S5 a, Visualization of two interacting water molecules (M1 and M2) under an applied\nelectric field of 0.25 V/Å oriented along the positive z-axis. The water molecules model\noriented in the yz-plane, with oxygen and hydrogen atoms shown in red and white, respectively.\nb, Temporal evolution of z-positions difference (|𝑅 −𝑅 |) of the two water molecules\nGO,M G1,M\nand c, charge dynamics of both water molecules calculated using QEq method (yellow curves),\nour universal model (blue dots and curves), and DFT reference (green curves). Solid and dash\nlines represent two molecules, respectively. The QEq and DFT results are extracted from\nreference [15].\nQ17. “‘outperforming conventional equivariant GNN-based MLIPs’ no evidence is showed.\n(Ref. 48. Is not an MLIP that fits this description).”\nAuthors’ reply:\nWe thank the reviewer for their important observation. Upon review, we realize that our\nstatement was imprecise and potentially misleading. What we intended to emphasize was the\ndifference in Na-Na/Na-Cl dimers resolution between our model and message-passing GNNs\nwithout long-range interactions that we trained ourselves, not the model in the reference 48.\nWe have modified our statement to be more precise and accurately reflect our comparative\nanalysis in the Disscussion section,\nOur pretrained universal model successfully reproduced DFT potential energy surfaces for\nNa-Cl/Na-Na dimers across extended interatomic distances, outperforming conventional\nmessage-passing MLIP beyond the cut-off regions.\nQ18. “Remarks on code availability:\nI haven’t had the chance to try it out, but it is available, and there are sufficient details to\ndownload/use.”\nAuthors’ reply:\nWe thank the reviewer for verifying the availability of our code and acknowledging that\nsufficient details are provided for downloading and usage. To further facilitate ease of use, we\nhave open-sourced our entire model framework in a new repository.\nCode availability\nThe JAX implementation code is available through https://github.com/reaxnet/reaxnet.\n29/43\nResponse to reviewer #3\nQ1. “### Summary\nThe manuscript presents a semi-heuristic approach to incorporating long-range interactions,\nsuch as polarizability, into machine learning potentials. The proposed framework integrates a\nphysically motivated polarizable charge equilibration scheme within an equivariant graph\nneural network-based short-range potential. While the idea is clever, it is not entirely new—\nsimilar approaches have been explored before, such as in *SpookyNet*, although the current\nwork extends these ideas further. The authors claim their model is universal and exhibits linear\nscaling in system size while accurately capturing long-range interactions.\nWhile the proposed approach is an interesting idea and appears to perform well on specific test\ndatasets, there are significant concerns regarding the validity of the claims made in the title and\nthroughout the manuscript. The manuscript lacks rigorous benchmarking against existing\nuniversal models, and the claim of linear scaling contradicts fundamental physical /\nmathematical limitations. As such, I cannot recommend publication in its current form as\noutlined below in more details.”\nAuthors’ reply:\nThank you for your thorough review of our manuscript. We deeply respect your expert\nassessment and have carefully considered all your points. We have made substantial revisions\nto address the concerns you raised, particularly by adding comprehensive benchmarking\nagainst existing universal models as you suggested.\n(1) Regarding SpookyNet, it directly incorporates the physically observable dipole moment\ninto the loss function rather than fitting partial charges. However, its fundamental\nrepresentation still relies on point charges, which significantly limits its ability to respond\naccurately to external fields.\nWe add discussions in the Introduction section,\n… Another way to solve the partition-dependent issue is to take dipole moment into the\nloss function, as done by SpookyNet [37]. However, it still uses point-like charges, making it\nbehaves similar to Qeq-based methods when dealing with polarization effects.\n(2) For the scaling assertion, we are deeply sorry that we unintentionally made incorrect\nclaim. We now provide a mathematical analysis that clarifies the actual cubic complexity\n(inversing the matrix), acknowledging the fundamental physical/mathematical limitations you\npointed out.\n(3) For the universality, we have included additional benchmarking against relevant\nuniversal models to substantiate our comparative performance claims. We believe these\nsubstantial revisions address the core concerns while maintaining the scientific value of our\ncontribution. We appreciate your critical feedback which has significantly improved the\naccuracy and rigour of our manuscript.\n30/43\nQ2. “### Major Concerns\n#### 1. **Lack of Rigorous Benchmarking for a “Universal” Model**\nThe authors make a very strong claim that their model is universal, but the provided evidence\ndoes not support this assertion. There are two critical issues with this claim:\n**Absence of Comparison with Established Models:** The manuscript does not compare the\nmodel’s performance against other well-known universal machine learning potentials (e.g.,\nMACE-MP-0, MATTERSim, SevenNet, ORB). Not only do the authors omit numerical\ncomparisons, but these models are not even cited. Instead, the error analysis is restricted to\ntheir own training dataset, which is insufficient to validate universality.”\nAuthors’ reply:\nWe thank the reviewer for this important observation regarding comparisons with\nestablished models. We acknowledge this limitation in our current manuscript and agree that a\nmore comprehensive comparison is necessary to substantiate our claims of universality. We\nhave cited the universal models mentioned by the reviewer in the Introduction section,\n… Moreover, universal models trained on the periodic table such as M3GNet [15], CHGNet\n[16], GNoME [17], MACE-MP-0 [18], MatterSim [19], SevenNet [20], and ORB [21] have\nemerged, showing remarkable prospects for materials discovery.\nWe now conduct benchmarks to compare our model with these advanced models. It is\nimportant to note that we only compare models trained on the MPtrj dataset (CHGNet and\nMACE-MP-0). Although SevenNet was also trained on the MPtrj dataset, it did not split the\ndata into training, validation, and test sets during training (see https://github.com/MDIL-\nSNU/SevenNet/tree/main/sevenn/pretrained_potentials/SevenNet_0__11Jul2024#sevennet-0-\n11july2024). We employ the challenging test sets proposed by MatterSim for benchmarking\nand found that our energy and stress errors are reasonable, and our force error is better than\nboth universal models on all sets. Since forces are the most relevant metrics for molecular\ndynamics simulations, this demonstrates the robustness of our model.\nIn the Results section,\nWe also perform benchmarking of our model with other universal models on the test sets\nprovided by the reference [19]. Table S5 demonstrates the exceptional performance of our\nmodel across diverse test sets in force prediction capabilities, consistently outperforming both\nCHGNet [16] and MACE-MP-0 [18], while maintaining reasonable energy and stress\npredictions. We excluded MatterSim [19] results from our comparative analysis as its training\ndataset incorporates high-temperature data.\nAnd in SI,\nWe also conducted benchmarking with other universal models on the test sets provided by\nthe reference [11]. For near-equilibrium evaluations, two randomly sampled collections are\nused: 1,000 structures from the MPtrj dataset [7] (MPtrj-random-1k) and 1,000 from\n31/43\nAlexandria dataset [12] (Alexandria-1k). These assess model accuracy on energetically\nfavourable configurations, which is crucial for predicting material stability. Three increasingly\nchallenging thermodynamic perturbation (TP) test sets evaluate performance under varied\ntemperature and pressure conditions with non-equilibrium atomic configurations … Table S5\ndemonstrates the performance of our model across diverse test sets, particularly highlighting\nits superior force prediction capabilities. Our model achieves remarkable force prediction\naccuracy, consistently outperforming both CHGNet [7] and MACE-MP-0 [14] (all three\nmodels were trained on MPtrj dataset) across all benchmark datasets with significantly lower\nMAEs – nearly half the error of competing models in many cases. Impressively, our model\ndemonstrated robust generalization capability across increasingly challenging thermodynamic\nperturbation scenarios. Although we observe slightly higher energy prediction MAEs\ncompared to the two models, its performance across the entire spectrum of materials conditions\nis improved particularly for force predictions critical for molecular dynamics simulations.\nTable S5. Comparison of mean absolute errors (MAE) in the test set generated by the reference\n[11].\nTest set MAE CHGNet MACE-MP-0 Our model\nEnergy (eV/atom) 0.027 0.015 0.019\nMPtrj-random-1k Force (eV/Å) 0.120 0.117 0.066\nStress (GPa) 0.290 0.468 0.271\nEnergy (eV/atom) 0.150 0.092 0.115\nAlexandria-1k Force (eV/Å) 0.108 0.095 0.064\nStress (GPa) 1.643 0.160 1.055\nEnergy (eV/atom) 0.250 1.351 0.280\nMPF-Alkali-TP Force (eV/Å) 1.636 15.819 0.717\nStress (GPa) 12.625 25.723 13.841\nEnergy (eV/atom) 0.254 256.340 0.521\nMPF-TP Force (eV/Å) 3.313 1506.854 0.768\nStress (GPa) 25.208 202.093 15.445\nEnergy (eV/atom) 0.506 9.184 3.338\nRandom-TP Force (eV/Å) 3.950 88.327 1.199\nStress (GPa) 7.230 19.224 33.282\nWe believe these additions significantly strengthen the manuscript and provide the\nnecessary context for evaluating our model’s universality claims.\nQ3. “2. **Inadequate Testing of Universality:** The authors evaluate their model on the\ncommon datasets published by Ko et al. But, the author performs their benchmarks bespoke\nmodels rather than testing their claimed universal model. If the model is indeed universal, it\nshould be evaluated directly on these datasets without additional training. A fair comparison\nwould involve evaluating both the universal and bespoke models to determine whether the\nuniversal model effectively learns long-range physics.”\n32/43\nAuthors’ reply:\nWe appreciate the reviewer’s important point regarding the evaluation of our model’s\nuniversality. We would like to clarify that due to differences in basis sets, functionals, and\ncomputational software, we cannot directly validate a universal model trained with\nVASP/PBE+U functional and pseudopotential basis sets on the datasets proposed by Ko et al.\nThe fundamental reason we adopted Ko et al.’s datasets was to verify the feasibility of our\nframework. Our universal model achieves a test error of 18 meV/atom on the MPtrj dataset,\nwhich is significantly larger than the ~ 1 meV/atom error that would be achieved by the model\ntrained on the specific datasets proposed by Ko et al. To the best of our knowledge, so far, no\nuniversal model can reach ~ 1 meV/atom of precision.\nA fairer comparison might be to progressively ablate our models, including benchmarking\nmodels containing only GNN components against our framework model, to demonstrate that\nthe introduction of long-range interactions further improves model accuracy. This approach\nallows us to isolate and quantify the specific contribution of the long-range interaction\ncomponents within our framework.\nIn the Results section,\nBenchmarking against different-charge-state dataset. To evaluate the capability of our\nframework in capturing different charge states and long-range interactions, we validated our\nframework against the dataset developed by Ko et al. [28] and Maruf et al. [38], which\nencompass various charge states and charge transfer systems … We demonstrate the advantage\nof integrating equivariant message-passing networks with explicitly polarizable long-range\ninteraction models by benchmarking against both 4G-HDNNP, a local MLIP with explicit\nlong-range interactions, and NequIP, which implicitly captures long-range effects through its\nmessage-passing mechanism. 4G-HDNNP differentiates charged systems by explicitly training\nDFT partial charges and incorporating them as descriptors into neural networks for short-range\ninteractions. We employed the NequIP models trained in the reference [38] (Maruf’s NequIP)\nas our baseline equivariant model. Our framework explicitly incorporates long-range\ninteractions while adopts Maruf’s NequIP as the E component without touching its settings,\n0\nthereby ensuring a fair comparative analysis. The results are presented in Table 1, with detailed\nneural network architectural specifications provided in the Table S3.\nTable 1. The comparison of root means square test error metrics on the different charge-\nstate datasets.\nDataset 4G-HDNNP Maruf’s NequIP This work\nEnergy (meV/atom) 1.19 1.33 0.44\nC H /C H +\n10 2 10 3\nForce (eV/Å) 0.078 0.071 0.023\nEnergy (meV/atom) 0.48 N. A. 0.23\nNa Cl +\n8/9 8\nForce (eV/Å) 0.033 N. A. 0.006\nEnergy (meV/atom) 1.32 498.39 4.87\nAg +/-\n33/43\nDataset 4G-HDNNP Maruf’s NequIP This work\nEnergy (meV/atom) 0.22 1.03 0.14\nAu -MgO\nEnergy (meV/atom) N. A. 0.48 0.30\nBTA-Cu\nForce (eV/Å) N. A. 0.008 0.004\nEnergy (meV/atom) N. A. 0.71 0.19\nBTA (H O)-Cu\n… In the linear C H molecule, … after explicitly incorporating physical long-range\n10 2\ninteractions, our framework achieved additional improvements in error metrics compared to\nNequIP. For positively charged NaCl clusters, the predicted energy and forces by our model\nare still more accurate compared to 4G-HDNNP, as well as the potential energy surface shown\nin Fig. S1 b. For Au adsorption on MgO surfaces and benzotriazole interactions with Cu-(111)\n2\nsurfaces in dry and aqueous environments, our model outperforms NequIP significantly.\nAdditionally, our model correctly predicted the preferential adsorption geometries for both\ndoped and undoped MgO, as listed in Table S1.\nAs for the charged Ag clusters, due to their small system size and lacking long-range effects,\nthe entire structure falls within the typical cutoff radius used by the MLIP. The energy\ndependence on the overall charge state of clusters leads to degeneracy between atomic\nstructures and potential energy surfaces, resulting in poor performance of charge-independent\nNequIP model. We acknowledge that for such small, highly charged systems, 4G-HDNNP\napproach to learning charge equilibration parameters is more convenient. Nonetheless, we also\nachieved reduced force error with the fixed physical parameters compared to 4G-HDNNP. We\nre-optimized the geometries of Ag +/− using our model, with results shown in Fig. S1 c. Our\n3\nmodel accurately predicts the structures, with root mean squared displacements of only 0.002\nÅ and 0.006 Å compared to DFT-optimized structures.\nOur comprehensive evaluation demonstrates that explicit incorporation of physical long-\nrange interactions significantly enhances the performance of MLIPs across diverse charge-state\nsystems. Notably, our proposed framework outperforms 4G-HDNNP in most cases for energies\nand forces without fitting partial charges. This advantage likely stems from a more physically\nmeaningful partition of the total energy.\nAfter verifying that our model could perform well on the challenging dataset, we then\nfocused on training a universal model. As mentioned above, we have supplemented our\nuniversal model with benchmarks to validate its generalizability across different datasets.\nThese additional evaluations demonstrate that our approach not only excels on specific test\ncases but also maintains good performance when applied broadly across diverse chemical\nsystems and properties.\nQ4. “#### 2. **Misleading Scaling Claims**\n34/43\nOne of the most problematic claims in the manuscript is that the model achieves linear scaling\nfor long-range interactions:\n- **Fundamental Computational Limits:** Long-range interactions require solving the Poisson\nequation, which fundamentally cannot be computed in linear time. The best possible scaling\nfor such a computation is **O(N log N)** using Ewald summation on a grid, while traditional\nsolvers scale as **O(N³/²)** or worse.\n- **Questionable Interpretation of Figure S4:** The claim of O(N) scaling is based on Figure\nS4, but the presented data does not conclusively support this claim. If the authors indeed\nachieve linear scaling, they must be using a cutoff solver, which is known to introduce severe\nerrors in energy and force computations.\nGiven these issues, the scaling behaviour of the proposed model needs significant clarification.”\nAuthors’ reply:\nWe sincerely appreciate the reviewer’s expert insights. We apologize for being misled by\nour original Fig. S4, where the large y-axis values and insufficient number of atoms incorrectly\nsuggested linear scaling. Given that our approach involves matrix inversion, it should exhibit\ncubic scaling. We did not measure the scaling, instead, we unintentionally took the apparent\nbut wrong trend. To fix this misleading presentation, we extend the system sizes to include up\nto 24,000 atoms being tested to reflect the cubic scaling behaviour.\nWe apologize again for this oversight. Nevertheless, for the reaction systems containing tens\nof thousands of atoms that we studied, the computational cost of cubic scaling remains\nacceptable, and the absolute computation time is still significantly better than CHGNet. In\nfuture work, we plan to update our algorithms by implementing iterative methods to solve the\nlinear equation systems, which will achieve lower computational complexity.\nIn the Results section,\nComputational cost. Despite the incorporation of additional long-range interaction\ncalculations, our model maintains high computational efficiency. Performance benchmarks\nwere conducted on a single NVIDIA H100 GPU to quantitatively assess the computational cost.\nAs illustrated in Fig. S6, for a periodic system comprising 2160 atoms, our model requires\napproximately 0.07 s per molecular dynamics time step, representing a significant\nimprovement over conventional universal MLIPs which require approximately 0.91 s. This\ncomputational advantage extends to larger systems, as demonstrated by simulations of a 24000-\natom system, where our model maintains efficiency with only 1.49 s per time step. Such\ncomputational performance enables nanosecond-scale molecular dynamics simulations of\nsystems containing tens of thousands of atoms, making it practical for more realistic modelling\nof materials.\nAnd in SI,\nOur model introduces explicit polarizable long-range interactions while still maintaining\ngood computational efficiency. To highlight the performance of this model, benchmarking was\n35/43\ndone on molecular dynamics simulations. Supercell models were generated from Li PO N\n2 2\n(mp-1020019) from Materials Project [8] and simulations were conducted. Additionally, the\nPyTorch-based universal machine learning interatomic potentials CHGNet [7] was used for\ncomparison. All molecular dynamics simulations were carried out under the canonical\nensemble simulations with a Nosé-Hoover thermostat [16, 17] at 300 K temperature. All the\ncomputations were performed with a single NVIDIA H100 GPU.\nCHGNet\n(w/o D3)\nOur model\n(w/ D3)\nFig. S6 Benchmarking on the computational efficiency involved in molecular dynamics\nsimulations and geometry optimizations. All the calculations are performed on a single\nNVIDIA H100 GPU. The time consumption per step during 100 steps of molecular dynamics\nsimulations for our model is presented. Results for CHGNet are shown in blue dots; however,\nthe calculation for systems with more than 10 thousand atoms was halted due to out of GPU\nmemory in the calculation environment.\nAs shown in Fig. S6, we present efficiency in terms of time consumption per simulation\nstep. Although our model demonstrates O(N³) scaling due to inversing the Eq. S3, it\nsuccessfully handles large-scale systems exceeding 20,000 atoms, highlighting its applicability\nto complex materials simulations. In our studied reaction systems containing tens of thousands\nof atoms, the inversion method provides acceptable computational efficiency. Nevertheless, in\nfuture work, we plan to further refine our algorithms, such as implementing iterative solution\nmethods to reduce complexity, thereby achieving more computationally efficient\nimplementations and lowering resource consumption.\nQ5. “- **Concerns Regarding Specific Applications**\nThe manuscript presents results on well-known test cases, but the methodology used for\ncomparison raises concerns:\n36/43\n- **Barium Titanate (BTO): ** The authors compare phase transition temperatures in Table 2,\nbut this comparison is misleading because the transition temperature is highly dependent on\nsystem size. The authors use a 10×10×10 cell, whereas the reference MLIP from Gigli et al.\nmodels use a 4×4×4 cell. The comparison should be made at equal system sizes or, at the very\nleast, the system size used should be explicitly stated in the table with a comprehensive\ndiscussion.”\nAuthors’ reply:\nWe appreciate the reviewer’s critical observation regarding the comparison of phase\ntransition temperatures for Barium Titanate in Table 2. The reviewer is correct that transition\ntemperatures can be significantly affected by system size, making direct comparisons\npotentially misleading without proper context.\nIn response to this valuable feedback, we have revised Table 2 to explicitly include the\nsystem size used for each model’s calculation, making the comparison transparent. Our\n10×10×10 supercell size is now clearly indicated alongside the 4×4×4 cell used in the\nGigli et al. reference.\nIn the Results section,\nTable 2. The phase transition temperatures of BaTiO3 obtained by different methods and\nexperiments. The size of supercell is also displayed.\nMethod 𝑇 (K) 𝑇 (K) 𝑇 (K)\n3,%PQ 3,QPR 3,RPS\nThis work (10×10×10) 145 205 285\nEffective Hamiltonian (12×12×12) [62] 150±10 195±5 265±5\nSecond Principles (16×16×16) [63] 140 180 224\nEffective Hamiltonian (16×16×16) [58] 119 158 257\nReaxFF (6×6×6) [64] N.A. N.A. 240\nMLIPs (4×4×4) [56] 18.6 91.4 182.4\nExperiments [65] 183 278 403\nWe have attempted to simulate phase transitions using our model with a 4×4×4 supercell to\nmatch the reference work by Gigli et al. However, we found that such small systems exhibit\nsignificant temperature fluctuations, making it challenging to accurately capture the free energy\nlandscape and precise transition temperatures. To address this issue, we have added the\nfollowing discussion to the manuscript.\nIn the Results section,\n… To detect lattice distortions and variations in free energy that distinguish different phases,\na 10×10×10 supercell of BaTiO was simulated. The use of larger simulation cells (10×10×10\n3\nin our work versus 4×4×4 in Gigli et al. [56]) may provide advantages for phase transition\nstudies. Larger supercells can effectively reduce the impact of temperature fluctuations,\nenabling better temperature control for sampling the free energy landscape, thereby yielding\nresults with greater statistical significance.\n37/43\nWe would like to clarify that our intention was not to directly compare our results with\nthose of Gigli et al. We do not consider results from different supercell sizes to be directly\ncomparable but rather aim to demonstrate that our model can reasonably predict phase\ntransition behaviour in complex material systems.\nIn the Results section,\n…It worth noting that the supercell size can affect phase transition temperatures, so that\ndirect comparison between different simulations is not expected…\nQ6. “- **Lithium Phosphorous Sulfide (LPS):** The claim that their model captures long-\nrange interactions in LPS is questionable. In this material, interactions are largely screened,\nand no significant long-range charge transfer is expected to occur. This system has already\nbeen simulated successfully with non-long-range potentials (e.g., *Phys. Rev. Materials 8,\n065403*), raising doubts about the necessity and impact of the proposed approach.”\nAuthors’ reply:\nWe thank the reviewer for the valuable insights. There may be a misunderstanding, as we\ndid not emphasize that our model captures long-range interactions in LPS, specifically in our\napplication examples of Li PS -Li and Li PS Cl-Li electrolyte-electrode interface reactive\n3 4 6 5\nmolecular dynamics. Regarding the reference [Phys. Rev. Materials 8, 065403 (2024)\nhttps://doi.org/10.1103/PhysRevMaterials.8.065403], it studied the bulk properties (thermal\nconductivity) of different LPS phases using MLIP, but not interfacial reactions.\nInstead, we used these examples to compare with experimentally observed phenomena and\nsome first-principles simulations merely to validate the robustness of our universal model –\nwhether a universal model trained solely on the MPtrj dataset (no interfacial data), without\nincorporating specific reaction data, can correctly describe reaction molecular dynamics and\neven reproduce certain experimental phenomena.\nAdditionally, our model can clearly reflect changes in partial charges during the molecular\ndynamics process, as shown in Fig. 5c, facilitating our analysis of SEI formation processes.\nThe distribution of Li partial charges gradually forms a stable plateau as the SEI forms.\nIn the Results section,\nWith the growth of the SEI layer, electrons were transferred from the Li anode to the Li PS\n3 4\nelectrolyte, which causes a gradual increase of the partial charges of Li as they migrate from\nthe anode, across the SEI, and into the electrolyte. This results in a transition from metal\nlithium to ions. The dynamic behaviour of partial charges at the anode-electrolyte interface is\na critical aspect showcased in Fig. 5 c. Initially, the partial charges of lithium changed almost\nlinearly from the anode to the electrolyte. Over time, this distribution evolved and showed a\ndistinct plateau. This phenomenon could be attributed to the ordered structuring of the SEI,\nsignifying the nucleation of crystalline structures. Ultimately, a stable lithium partial charges\nplateau was formed between the anode and the electrolyte, spanning a range of 15 to 21.5 nm.\n38/43\nThis is further supported by the visualized structure depicted in Fig. 5 b, where the plateau\nwithin the SEI aligns with the Li S crystal and amorphous regions.\n100, 2000, and 4000 ps states.\nThere have been previous examples [Energy Environ. Sci. 17, 2743-2752 (2024)\nhttps://doi.org/10.1039/D3EE03536K] using MLIP to study interfacial reactions, but these\nwere based on AIMD data of Li PS -Li interface reactions, whereas our model successfully\n3 4\nreproduces experimentally observed results through training on MPtrj dataset. Furthermore, we\nhave also used our universal model to study Li PS Cl-Li interfacial reactions, reproducing the\n6 5\npolycrystalline structure observed in experiments, as shown in Fig. S7. This represents an\nadvantage over specialized MLIP models that would typically require retraining to include\nthe Cl element.\nFig. S7 The formation of the Li PS Cl/Li SEI layer at different temperatures. The elements\n6 5\nare colour-coded for initial structures: Li in green, S in purple, P in pink, and Cl in light green.\nFor final structures at 300 K, 400 K, and 500 K, the partial charges on lithium ions are\nrepresented with colour coding to enhance the visibility of structural transformations during\nthe formation of SEI.\nQ7. “- **Lack of Fine-Tuning Capabilities:** Universal models are only widely useful if they\ncan be fine-tuned to specific systems. This is common practice in other domains (e.g., stable\ndiffusion for images or MATTERSim, MACE-MP0 for materials), yet the authors provide no\n39/43\nevidence that their model supports such fine-tuning. Examination of the code suggests that this\nfeature is not even planned.”\nAuthors’ reply:\nWe thank the reviewer for raising this important point about fine-tuning capabilities in\nuniversal models. We agree that the ability to fine-tune a universal model for specific systems\nis valuable for enhancing performance in targeted applications. We have included a\ndemonstration example in the supplementary materials showing how our model can be fine-\ntuned for high-temperature dataset, with quantitative improvements in property prediction\naccuracy.\nIn SI,\n… To evaluate our model’s capacity to learn high-temperature data representations, we fine-\ntuned our pretrained universal model on the MPF-TP dataset. During this fine-tuning process,\nonly the final two convolutional layers and the multilayer perceptron layers were set as\ntrainable parameters, while all remaining layers were frozen.\nOur fine-tuning results demonstrate remarkable improvements on high-temperature, non-\nequilibrium systems as shown in Table S6. The fine-tuned model achieves state-of-the-art\naccuracy across all prediction targets on the MPF-TP dataset. This suggests that while our\nuniversal model can maintain reasonable force accuracy under a wide range of conditions,\ntargeted fine-tuning is more likely to balance broad applicability with the superior accuracy\nrequired for specific applications.\nTable S6. Comparison of mean absolute errors in the MPF-TP test set of the universal\nmodels and finetuned model.\nEnergy (eV/atom) Force (eV/Å) Stress (GPa)\nCHGNet 0.254 3.313 25.208\nMACE-MP-0 256.340 1506.854 202.093\nMatterSim (M3GNet) 0.036 0.431 1.318\nMatterSim (Graphormer) 0.040 0.421 1.917\nOur universal model 0.521 0.768 15.445\nFine-tuned model 0.027 0.203 0.976\nWe also included a demonstration example in our code repository showing how our model\ncan be fine-tuned for specific material systems.\nWe believe this addition significantly enhances the practical utility of our universal model,\nand we appreciate the reviewer highlighting this important aspect of model development.\n40/43\nQ8. “### Minor Issues and Suggested Revisions\n1. **Dataset Clarifications:** Clearly state that MPtrj is the target training dataset in the main\ntext (Page 6, Line 169).”\nAuthors’ reply:\nWe thank the reviewer for this suggestion. We have revised it to clearly state that MPtrj is\nour target training dataset.\nIn the Results section,\n… We trained a universal model for all the periodic table elements up to Pu using MPtrj\ndataset [16] following our framework (our model), as described in the Methods section …\nAnd in the Methods section,\nTo train the universal equivariance neural network potential, we used the MPtrj dataset [16]\nsourced from Materials Project [48] as the training dataset …\nQ9. “2. **Consistency in Units:** Ensure uniform usage of units, likely V/Å, in all relevant\nsections (e.g., Page 7, Line 221).”\nAuthors’ reply:\nWe thank the reviewer for highlighting this important issue regarding consistency in units.\nWe have carefully reviewed the manuscript and standardized all electric field units to V/Å.\nIn the Results section,\n… As illustrated in Fig. S5, we conducted additional investigations of two water molecules\nunder an external electric field of 0.25 V/Å.\nIn SI,\n… Fig. S5 b and c depict the temporal evolution of corresponding z-positions and molecular\ncharges within the simulation box for both methods at an electric field of 0.25 V/Å.\nRegarding the units of PQEq parameters, to maintain consistency with original PQEq paper,\nwe have retained the conventional units where 𝜒 and 𝜂 are expressed in eV, while the spring\nconstant is expressed in kcal/mol/Å2. All other energy units throughout the work have been\nstandardized to eV (meV) for consistency. For forces and stresses units, eV/Å and GPa are\nused, respectively.\nQ10. “3. **Charge Values for Water Molecule:** Provide explicit charge values from Figure\n2d for easier comparison, similar to what is done in Figure S1a.”\n41/43\nAuthors’ reply:\nWe thank the reviewer for this helpful suggestion. We have added supplementary materials\nto illustrate the charge variations. We believe that due to differences in charge partitioning\nmethods, comparing absolute charge values may not be meaningful; however, the trends in\ncharge variations can effectively reflect our model’s response to external electric fields. The\nsupplementary information now includes a detailed analysis of how partial charges evolve\nunder different electric field strengths.\nIn SI,\n…. To illustrate this difference, we present in Fig. S4 the variations in O atom partial charge\ndifferences (𝑞 −𝑞 ) across different models as a function of the applied electric field\n!\"#\nstrengths …\nFig. S4 The variations in O atom partial charge differences (𝑞 −𝑞 ) between DFT\n!\"#\nHirshfeld calculations, our model (PQEq model), and the conventional QEq model.\nQ11. “4. **Thermostat Choice**: The authors use a Nosé-Hoover thermostat in their\nmolecular dynamics simulations, which does not properly sample the canonical ensemble for\nsmall systems and can artificially suppress fluctuations. A stochastic thermostat (e.g., Langevin\nor Bussi-Donadio-Parrinello) would be more appropriate to ensure correct statistical sampling.”\nAuthors’ reply:\nWe thank the reviewer for this insightful comment regarding thermostat choice in our\nmolecular dynamics simulations. We agree that the choice of thermostat can significantly\nimpact the quality of statistical sampling, especially for small systems. However, we would\nlike to clarify that the smallest system in this manuscript underwent MD simulations contains\n1,536 atoms. For systems of this size, we have conducted additional comparative simulations\nthat demonstrate minimal differences in results between the Nosé-Hoover thermostat and\nstochastic thermostats such as Langevin. Our supplementary tests confirm this equivalence, as\nshown in the Fig. R2 and Table R2.\n42/43\na b\nFig. R2 Evolution in temperature, potential energy, and pressure during molecular dynamics\nsimulations of 1536-atom LLZO using a, Nosé-Hoover, and b, Langevin thermostats at 1000K.\nTable R4 Confidence intervals of temperature, potential energy, and pressure during molecular\ndynamics simulations of 1536-atom LLZO with Nosé-Hoover and Langevin thermostats.\nThermostat Temperature (K) Potential energy (eV) Pressure (GPa)\nNosé-Hoover 1000.29 ± 21.40 −11627.45 ± 3.97 −0.01 ± 0.20\nLangevin 997.51 ± 20.60 −11628.02 ± 3.94 0.01 ± 0.20\nAs Fig. R2 and Table R4 demonstrate, the differences in temperature, potential energy, and\npressure between the two thermostats are negligible and well within statistical uncertainty for\nour systems. While we acknowledge that stochastic thermostats offer theoretical advantages\nfor canonical ensemble sampling in smaller systems, our specific applications involve\nsufficiently large simulation cells where the Nosé-Hoover thermostat performs adequately for\nthe properties of interest.\nQ12. “Remarks on code availability:\nThe provided repository includes a README file with installation instructions, and the code\nappears to be installable. However, the core model presented in the manuscript is not included\nin the repository, preventing direct verification of the results.”\nAuthors’ reply:\nWe thank the reviewer for their feedback regarding code availability. We have open-sourced\nthe complete model codebase in the new repository and supplemented it with additional\ntutorials and fine-tuning examples to facilitate implementation for users.\nCode availability\nThe JAX implementation code is available through https://github.com/reaxnet/reaxnet.\n43/43\nResponse to the reviewers\nWe sincerely thank the three reviewers for their continued time and effort in this second round\nof review. Their additional feedback has been invaluable, and we have revised the manuscript\nto address all remaining points. We have provided point-by-point responses to the questions\nraised by the reviewers in this document. The modifications to the manuscript are highlighted\nin dark red.\nTable of Contents\nResponse to reviewer #1 ............................................................................................................ 2\nResponse to reviewer #2 ............................................................................................................ 3\nResponse to reviewer #3 .......................................................................................................... 10\n1/18\nResponse to reviewer #1\n“The authors have addressed most of my concerns. One remaining point remains.\nQ6. “Figure 2 b-c: the authors should increase the x-axis range to 12 Å to actually see the long-\nrange scaling behaviour.”\nI do not see “other publications also do it like this” as a valid argument. Understandably, the\nsurrogate model cannot be more accurate than the data-generating DFT model. So, I’m not sure\nwhy the authors are trying to hide this. Thus, the authors should add the results. If not in the\nmain paper then at least in the supplementary information.”\nAuthors’ reply:\nWe sincerely thank the reviewer for acknowledging our previous revisions. We also\nappreciate your valuable suggestion, which has prompted us to further clarify an important\nscientific issue. We agree with your point that simply citing ‘standard practice’ of other\npublications is an insufficient justification, and we apologize for the lack of clarity in our\nprevious response. Your suggestion to extend the x-axis to 12 Å has indeed provided a much\nbetter perspective for a comprehensive evaluation of our model’s performance.\nIn the main text, we added the entire distance up to 12 Å for comparison,\n“Fig. 2 a, The interaction energies for Na-Na (green) and Na-Cl (purple) dimers predicted\nby our model, w/o-lr model, and reference calculations using DFT.\n…We analysed interaction energies for Na-Na and Na-Cl dimers across interatomic\ndistances ranging from 1.8 Å to 12.0 Å…”\nWe believe these additions fully address the reviewer’s point and have made the manuscript\nstronger. We thank the reviewer again for pushing us to include this important validation.\n2/18\nResponse to reviewer #2\n“Thanks to the authors for their reply and reviewing the manuscript. Factual errors have been\nmostly corrected, but there remains some open questions. These should be clarified and/or\ncorrected before considering this manuscript any further. Here are some comments:”\nAuthors’ reply:\nWe thank the reviewer for their feedback and for taking the time to review our revised\nmanuscript. We appreciate the acknowledgment that some factual errors have been corrected.\nWe have carefully considered the remaining questions and will address each point thoroughly\nin the updated manuscript to provide the necessary clarifications and corrections.\n“Table S5 does not represent any significant improvements over other methods. A selling point\ncould be increased transferability etc. but no evidence for that.”\nAuthors’ reply:\nWe thank the reviewer for this constructive and valuable feedback. We agree that the\noriginal manuscript did not sufficiently demonstrate the transferability of our method, which is\nindeed a key distinctive feature. Following the reviewer’s suggestion, we have performed\nadditional analyses and included substantial new evidence to explicitly showcase the superior\ntransferability of our approach. Specifically,\n1. Transferability of long-range physics. First, we demonstrate our model’s ability to\naccurately describe dimer interactions across multiple distance scales. As shown in the\nnew Fig. S7, our model correctly reproduces the potential energy surfaces for both Na-\nNa and Na-Cl dimers. In addition, our model outperforms CHGNet and MACE-MP-0\non both short-range and long-range distances. Second, we investigated the OH-OH\nsystem in both its neutral (weakly interacting) and ionic (charge-repulsive) states. We\nshowcase our model’s unique ability to distinguish between different charge states and\ncorrectly model their corresponding long-range physics—a fundamental challenge for\nother universal models.\n2. Transferability of polarization. To further provide evidence of our model’s\ntransferability in terms of polarization, we tested its ability to predict responses to\nexternal electric fields—a crucial test of physical fidelity that reveals limitations in other\nuniversal models. First, in a periodic system of 64 water molecules, our model accurately\nreproduced the quadratic energy dependence under a varying electric field, showing\noutstanding agreement with DFT reference calculations (Fig. S9). Notably, its\nperformance matched that of FIREANN, a state-of-the-art model specifically trained on\nthat system, whereas leading universal models like CHGNet and MACE-MP-0 were\ncompletely unresponsive to the field. In addition, we validated our model’s ability to\npredict molecular polarizability for 7,211 diverse molecules from the QM-7b dataset, a\ntask beyond the capabilities of other universal models. Our model achieved a low mean\nabsolute error of 4.57 a.u. against high-level theory reference values (Fig. S10).\n3/18\nThese additions provide the quantitative evidence for transferability of our foundation\nmodel. We are confident that these results now clearly establish a significant improvement over\nother methods, not just in the metrics of Table S5, but in the crucial dimension of physical\ntransferability.\nIn SI, we detail the transferability,\n“In comparison with existing universal models, our model demonstrates transferability in\ndescribing dimer interactions. As mentioned in the main text, regarding the potential energy\nsurface capturing of Na-Cl and Na-Na dimers, we have discovered that our model not only\naccurately reproduces the long-range behaviour predicted by DFT/PBE [18] calculations, but\nalso effectively captures the potential energy surfaces of these dimers in the short-range region\ncompared to CHGNet and MACE-MP-0, as illustrated in Fig. S7. The enhanced transferability\nof our model across different dimers suggests its potential utility for a broad spectrum of\napplications requiring accurate representation of interatomic forces across multiple distance\nscales.\nFig. S7 The interaction energies for Na-Na (green) and Na-Cl (purple) dimers predicted by\nour model, CHGNet [7], MACE-MP-0 (large version) [11], and reference DFT calculations.\nFor the neutral state (dashed lines), which represents the physical scenario of a weak, non-\nbonded interaction, the DFT reference shows a very shallow potential well with relative\nenergies. Our model demonstrates good agreement with the DFT reference, accurately\ncapturing the subtle attractive interaction. In contrast, MACE-MP-0 and CHGNet struggle to\ncapture this long-range interaction; they appear less transferable than our model as they fail to\nreproduce the precise shape and depth of the DFT curve. For the ionic state (solid lines), which\nrepresents a strong, long-range Coulombic repulsion between charged species, the DFT\nreference shows a significant positive (repulsive) energy that decays slowly with distance. Our\nmodel correctly identifies the ionic state and qualitatively reproduces the long-range repulsive\nenergy curve from DFT. Despite a slight overestimation of the energy, the model captures the\nphysical trend and magnitude of the interaction across the entire distance range. Conversely,\nMACE-MP-0 and CHGNet are completely unable to distinguish the ionic from the neutral state.\nThis represents a fundamental physical problem, indicating that their underlying architectures\ncannot capture charge-state-dependent long-range physics.\nFig. S8 Relative interaction energy versus distance for an OH-OH system in its neutral\n(dashed lines) and ionic (solid lines) states. Our model’s predictions are compared with DFT,\nMACE-MP-0, and CHGNet. The results demonstrate that our model can reproduce the DFT\nreference for both the weak interaction of the neutral state and, crucially, the long-range\nCoulombic repulsion of the ionic state. Universal models like MACE-MP-0 and CHGNet fail\nto capture the long-range physics of the ionic state.\n5/18\nA crucial test of a model’s physical fidelity is its response to external fields. We evaluated\nour model’s ability to predict the energy response of a system to a uniform electric field. We\nadopted a periodic system containing 64 water molecules with each cell length measured\n12.4185 Å, as described in reference [19]. All the water molecules are within yz-plane. The\nfirst and third layers are same while the second and fourth are same. An electric field was then\napplied along the x-axis with varying strength.\nFig. S9 Evaluation of polarizable interactions in the periodic water model. The plot\ncompares the performance of our model against DFT [19] (reference), FIREANN [19],\nCHGNet, and MACE-MP-0.\nThe results are presented in the Fig. S9, which plots the system energy as a function of field\nstrength. Our model demonstrates good agreement with the DFT reference, accurately\nreproducing the energy curve across the entire field range from −2 to 2 V/Å. Notably, the\nperformance of our foundation model is on par with that of FIREANN [19], a state-of-the-art\nmodel that was specifically trained on this water system. This highlights the transferability of\nour model, as it correctly captures complex polarization physics without explicit training for\nthis task. In contrast, the CHGNet and MACE-MP-0 models are entirely unresponsive to the\nfield, predicting a constant energy. This represents a fundamental failure to capture one of the\nmost essential physical interactions, underscoring the superior predictive power and physical\nrealism of our model.\nFig. S10 Comparison of the predicted isotropic polarizability (𝛼 ) against reference coupled-\n!\"#\ncluster theory including single and double excitations values for 7,211 molecules from the QM-\n7b dataset [20]. The model achieves a mean absolute error (MAE) of 4.57 a.u., demonstrating\nits strong capability to generalize and predict complex electronic response properties accurately.\nConventional universal machine learning potentials like MACE-MP-0 and CHGNet are\ntypically trained on field-free energies and forces, and their architectures lack the components\nto describe fundamental electrostatic properties like the dipole moment. Consequently, they are\ninherently incapable of predicting response properties such as polarizability, which is defined\nby the change in energy or dipole moment under an applied field. We further validated the\npowerful transferability of our model by assessing its ability to predict molecular polarizability,\na key electronic response property. To accomplish this, we curated a challenging test set\ncomprising 7,211 diverse small molecules from the QM-7b dataset [20-23] … This involved\ncalculating the system’s dipole response to a small, externally applied electric field—a physical\nprocess the model was not explicitly trained to reproduce. The results of this evaluation are\n6/18\npresented in the Fig. S10. The plot, which compares our model’s predictions against the\nreference values, demonstrates good agreement. The model achieves a MAE of just 4.57 a.u.\nacross this vast and chemically diverse dataset. This high level of accuracy, achieved without\ndirect training on polarizability data, provides compelling evidence for the model’s strong\ngeneralization capabilities and its ability to capture the fundamental physics governing\nmolecular electronic response, a capability that lies beyond the scope of other universal\nmachine learning potential models.”\n“Table S6. I have no stake in MACE-MP-0 but I can immediately tell that those results can’t\nbe right and I am sure the MACE-MP-0 authors will challenge those as soon as they see this\ntable. It seems extremely unlikely that in all the other metrics MACE-MP-0 does well and\nsuddenly the errors become huge.”\nAuthors’ reply:\nWe thank the reviewer for their careful scrutiny of the MACE-MP-0 results and for raising\nthis important point. We understand the concern, as the reported errors are indeed strikingly\nlarge. The large errors reported are indeed correct and stem from the specific nature of the\nbenchmark, which is designed to evaluate models under challenging thermodynamic\nperturbation datasets. The large version of MACE-MP-0 trained on the MPtrj dataset\nconsistently demonstrates limitations when evaluated on such extreme and challenging test\ncases.\nFurthermore, to ensure a rigorous and fair comparison, we conducted a thorough review of\nthe evaluation metrics. We discovered that the definition for MAE of forces and stresses used\nin the MatterSim paper (from which we cited their results) differs from the definition used\nduring their training. We have therefore re-evaluated all models using a single, unified\ndefinition to ensure consistency.\nTo guarantee full transparency and reproducibility, we have made our complete evaluation\nscript for MACE-MP-0 available on Figshare at the following link:\nhttps://figshare.com/s/7ecc56e3538d83c83010. We believe this clarification and the provided\nscript will fully address the reviewer’s concerns about the validity of our findings.\n“About the forces/stresses. Since the forces depend on the charges, which depend on the\npositions, one must take the chain rule and calculate the appropriate term. Table R2 is not\nconvincing at all - this compares to another implementation presumably using the same\nformulae. Did the authors compare their forces to forces obtained using finite differences?”\nAuthors’ reply:\nWe sincerely thank the reviewer for their insight on this critical point. After a careful re-\nderivation, we identified that it has a Hellmann-Feynman-like term which will vanish as you\nmentioned before. To provide a definitive and convincing validation, we have also followed\nyour suggestion and compared our auto-gradient forces against the forces obtained using finite\n7/18\ndifferences for a water molecule as described in Table R1. As depicted in Table R2, the auto-\ngradient forces and finite differences results show excellent agreement, confirming the\ncorrectness of our implementation. We are grateful for the reviewer’s suggestion, which has\nsignificantly improved the rigor of our work.\nThe complete and detailed derivation has been added to the SI,\n“The total forces of our model are,\n𝜕𝐸 𝜕𝐸 𝜕𝑞\n%#& %#& '\n𝑭\n$\n= −\n𝜕𝒓\n−(\n𝜕𝑞 𝜕𝒓\n, (S9)\n$ ' 𝒊\n'\nwhere 𝐸 , 𝒓 , and 𝑞 are total potential energy as described in Eq. 1, atomic coordinates, and\n%#& $ $\npartial charges, respectively. The Lagrange multiplier method (𝐸 = 𝐸 −𝜆∑𝑞 ) requires\n) %#& $\n𝜕𝐸\n)\n− = 0, (S10)\n𝜕𝑞\n'\nso that,\n𝜕𝐸\n%#&\n= 𝜆. (S11)\n𝜕𝑞\n'\nThus, the forces are simplified as,\n𝜕𝐸 𝜕𝑞\n%#& '\n𝑭\n$\n= −\n𝜕𝒓\n−𝜆(\n𝜕𝒓\n. (S12)\n$ 𝒊\n'\nAnd with the constraint,\n∑𝑞 = 𝑄, (S13)\n'\n*+\nwe have = 0 since the net charge 𝑄 is a constant. Finally, the atomic forces are,\n*𝒓𝒊\n𝜕𝐸\n%#&\n𝑭 = − . (S14)\n$ 𝜕𝒓\n$\n…”\nTable R1 Water molecule coordinates.\nx (Å) y (Å) z (Å)\nO −1.464 0.099 0.300\nH −1.956 0.624 −0.340\nH −1.797 −0.799 0.206\nTable R2 Comparison of the auto-gradient and finite differences forces on the water molecule.\nForces (eV/Å)\n𝐹./012!345 𝐹6/&# 𝐹./012!345 𝐹6/&# 𝐹./012!345 𝐹6/&#\n- - 7 7 8 8\nO −0.13566780 −0.13566779 −0.06141234 −0.06141233 −0.12068319 −0.12068319\nH 0.06315283 0.06315281 0.07240831 0.07240830 0.04431891 0.04431891\nH 0.07251500 0.07251498 −0.01099598 −0.01099597 0.07636429 0.07636429\n8/18\n“Regarding Fig. 2b. and R3 yes, it is clear now that due to the specific direction the electric\nfield was applied this is indeed the case. However, I disagree with the general statement “This\nis a fundamental limitation of the traditional QEq-based model—it lacks the capability to\nrespond to external electric fields when atomic positions remain fixed.” If the electric field had\ncomponents in the y or z directions, the energy would change.”\nAuthors’ reply:\nWe thank the reviewer for this insightful comment and for pointing out the imprecision in\nour original statement. The reviewer is absolutely correct. A traditional QEq-based model\nwould indeed show an energy change if the electric field applied has a component parallel to\nthe direction of molecular dipole moment. We are grateful to the reviewer for pointing this out.\nWe have revised it as in the main text,\n“…Conventional QEq-based models have limitations in their response to an external electric\nfield applied orthogonally to the molecular dipole…”\n““The term “free energy landscape” is the standard terminology for describing how free\nenergy varies across configuration space.” this is not quite right. This term describes the free\nenergy as a function of a collective variable or set of collective variables. But the manuscript\ntext is OK.”\nAuthors’ reply:\nWe thank the reviewer for the valuable clarification on the term “free energy landscape”\nand for the approval of our manuscript’s content. Reviewer’s rigorous approach and expertise\nin mathematical physics are, once again, truly impressive.\n9/18\nResponse to reviewer #3\n“The authors have addressed many of the concerns raised in my previous review, and the\nmanuscript has improved in several respects. For example, the addition of Table 1 is a valuable\ncontribution. However, I still find that the very strong claims made in the title and abstract—\nparticularly the notion of a universal long-range model—are not sufficiently supported by the\npresented evidence.”\nAuthors’ reply:\nWe sincerely thank the reviewer for acknowledging the improvements we made in the\nmanuscript, particularly for recognizing the value of the newly added Table 1. Our model can\ndemonstrate the capability to capture accurate long-range physics in certain scenarios. For\ninstance, it accurately models interactions across various distance scales, reproducing the\npotential energy surfaces for Na-Na and Na-Cl dimers (new Fig. S7) while outperforming\nbenchmarks like CHGNet and MACE-MP-0. Crucially, our model also addresses a key\nlimitation of other universal potentials by successfully distinguishing between different charge\nstates. This is highlighted by its ability to correctly model the distinct long-range physics of the\nOH-OH system in both its neutral (weakly interacting) and ionic (charge-repulsive) forms.\nHowever, we recognize that for more challenging and diverse systems, such as the dataset\nproposed by Ko et al., our model serves as a foundation model that requires further finetuning\nto achieve high accuracy. We will elaborate on this specific point in our response to the\nreviewer’s next question. We now position our work as a “foundation model”. This term more\naccurately reflects our contribution: a robust, physically-informed starting point that is\ndesigned to be efficiently finetuned to achieve ab initio accuracy for specific, challenging\napplications. This two-step approach—pretraining and finetuning—represents, in our view, a\nmore pragmatic and powerful paradigm for modern materials modelling. The revised title is\nnow,\n“A foundation machine learning potential with polarizable long-range interactions for\nmaterials modelling”\nAnd in the abstract, we further address the foundation model and finetuning,\n“…The foundation model, trained across the periodic table up to Pu, demonstrates strong\nperformance across key materials modelling challenges… Furthermore, we show that as a\nfoundation model, it can be efficiently finetuned to achieve high-level accuracy for specific\nchallenging systems.”\nDetailed physical test of long-range physics in SI,\n“In comparison with existing universal models, our model demonstrates transferability in\ndescribing dimer interactions. As mentioned in the main text, regarding the potential energy\nsurface capturing of Na-Cl and Na-Na dimers, we have discovered that our model not only\naccurately reproduces the long-range behaviour predicted by DFT/PBE [18] calculations, but\nalso effectively captures the potential energy surfaces of these dimers in the short-range region\ncompared to CHGNet and MACE-MP-0, as illustrated in Fig. S7. The enhanced transferability\nof our model across different dimers suggests its potential utility for a broad spectrum of\n10/18\napplications requiring accurate representation of interatomic forces across multiple distance\nscales.\nFig. S7 The interaction energies for Na-Na (green) and Na-Cl (purple) dimers predicted by\nour model, CHGNet [7], MACE-MP-0 (large version) [11], and reference DFT calculations.\nAnother example that showcases the transferability of our model is its ability to distinguish\nbetween different charge states and accurately simulate their corresponding long-range\ninteractions. We selected the OH-OH system as a case study, investigating the potential energy\ncurves when the system consists of two neutral hydroxyl radicals versus when it carries a single\nnegative charge. As illustrated in Fig. S8, the plot compares the performance of our model\nagainst DFT (the ground truth), MACE-MP-0, and CHGNet for the system in both “ionic” and\n“neutral” states as a function of distance.\nFig. S8 Relative interaction energy versus distance for an OH-OH system in its neutral\n(dashed lines) and ionic (solid lines) states. Our model’s predictions are compared with DFT,\nMACE-MP-0, and CHGNet. The results demonstrate that our model can reproduce the DFT\nreference for both the weak interaction of the neutral state and, crucially, the long-range\nCoulombic repulsion of the ionic state. Universal models like MACE-MP-0 and CHGNet fail\nto capture the long-range physics of the ionic state.\n11/18\nFor the neutral state (dashed lines), which represents the physical scenario of a weak, non-\nbonded interaction, the DFT reference shows a very shallow potential well with relative\nenergies. Our model demonstrates good agreement with the DFT reference, accurately\ncapturing the subtle attractive interaction. In contrast, MACE-MP-0 and CHGNet struggle to\ncapture this long-range interaction; they appear less transferable than our model as they fail to\nreproduce the precise shape and depth of the DFT curve. For the ionic state (solid lines), which\nrepresents a strong, long-range Coulombic repulsion between charged species, the DFT\nreference shows a significant positive (repulsive) energy that decays slowly with distance. Our\nmodel correctly identifies the ionic state and qualitatively reproduces the long-range repulsive\nenergy curve from DFT. Despite a slight and consistent overestimation of the energy, the model\ncaptures the physical trend and magnitude of the interaction across the entire distance range.\nConversely, MACE-MP-0 and CHGNet are completely unable to distinguish the ionic from\nthe neutral state. This represents a fundamental physical problem, indicating that their\nunderlying architectures cannot capture charge-state-dependent long-range physics.”\n“A central issue remains: while the model is trained to be universal, the authors do not\nconvincingly demonstrate that it has actually learned the correct long-range physics, as opposed\nto simply benefitting from increased expressivity due to the added long-range block. This point\nwas raised in my initial review, but the response did not adequately address it.\nThe authors argue that:\n> “we cannot directly validate a universal model trained with VASP/PBE+U functional and\npseudopotential basis sets on the datasets proposed by Ko et al.”\nWhile technically true, this is not a convincing justification for avoiding such validation. The\nKo et al. datasets contain only a few hundred small structures, and recomputing them using the\nauthors’ own DFT setup (VASP/PBE+U) is computationally feasible—especially in the\ncontext of training a large-scale universal model. Recomputing small datasets for consistent\nevaluation is standard practice in the field.\nThe further claim that: > “so far, no universal model can reach ~1 meV/atom of precision.”\nmisses the point. The Ko et al. datasets include physically motivated test cases, such as\ndisplacement scans of NaCl and Au₂ systems, which reveal whether a model correctly captures\nenergy and force trends with interatomic distance. The authors should at least reproduce such\nphysical validation tests with both their universal and bespoke models to support their claims.”\nAuthors’ reply:\nWe sincerely thank the reviewer for this critical and insightful feedback. We agree that our\nprevious response did not adequately address the central issue: whether our model has learned\nthe correct long-range physics or is merely benefiting from increased architectural expressivity.\nWe agree that performing validation on the potential energy surfaces of NaCl clusters and Au\n2\nsystems is indispensable.\nAs the reviewer previously pointed out, a pretrained foundation model should demonstrate\nits capacity for learning challenging physics through finetuning. We first evaluated various\nmodels on the potential energy surfaces (PES) of the Na Cl ⁺ and Na Cl ⁺ charged clusters\n8 8 9 8\nfrom the Ko et al. dataset. As illustrated in Fig. S11 (a, c, and e), all three foundation models\n12/18\n(ours, CHGNet, and MACE-MP-0) fail to reproduce the distinct DFT potential energy surfaces\nfor these two systems. This confirms the reviewer’s scepticism about the out-of-the-box\nperformance of universal models on such specific, challenging tasks. However, even at this\nstage, our model showed a nascent ability to distinguish between the two clusters, whereas the\nothers did not. We then finetuned all three models using a small subset (20%) of the\nconfigurations. The same algorithm was used for all models, and only the mean squared force\nloss was used to update the weights. After finetuning, PES of our finetuned model for both\nNa Cl ⁺ and Na Cl ⁺ achieve good agreements with the DFT reference (Fig. S11 b). In contrast,\n8 8 9 8\nwhile finetuning reduced the overall error for CHGNet and MACE-MP-0, their final predicted\nenergy curves settled into an intermediate position, averaging the two distinct DFT reference\ncurves (Fig. S11 d and f). They remain unable to differentiate between the two distinct clusters.\nWe may attribute this limitation to the absence of explicit descriptors for the system’s net\ncharge within their architectures, which prevents them from tweaking the unique energetic\nsignature of each cluster. We then extended our investigation to the complex periodic system\nof an Au dimer adsorbed on both Al-doped and undoped MgO surfaces. Following a similar\n2\nfinetuning protocol, we found that message-passing-based models are indeed capable of\nlearning these precise potential energy surfaces. As summarized in Table S6, both finetuned\nCHGNet and our model successfully reproduce the DFT equilibrium bond lengths on both\nsurfaces, with our model achieving slightly higher fidelity. The inferior performance of\nMACE-MP-0 in this case may be attributed to its limited effective receptive field (arising from\nonly two message-passing layers).\nOur new results not only validate this approach but also reveal the critical importance of\narchitectural design—specifically, the inclusion of net charge as a descriptor—in successfully\ncapturing complex electrostatic interactions.\nIn SI,\n“While our model provides a foundational understanding, achieving ab initio accuracy for\nspecific and challenging systems like charged clusters necessitates a refinement of the potential\nenergy surface. To rigorously assess the physical fidelity of the foundation models, we\nconducted an evaluation on the potential energy surfaces (PES) of charged sodium chloride\nclusters from the dataset proposed by Ko et al. [4]. The initial test involved scanning the PES\nof the Na Cl + cluster systems along a coordinate defined by varying the distance between\n8/9 8\ntwo sodium atoms. As illustrated in Fig. S11 (a, c, and e), a comparison with DFT calculations\nreveals that all three foundation models fail to accurately describe the PES for these systems.\nNotably, our model demonstrated a potential to distinguish between the two clusters due to the\ncharge state descriptor, whereas other models did not.\nWe acknowledge that capturing such a refined PES requires model specialization, so we\nperformed a finetuning procedure on all three models. A subset was created by randomly\nselecting 20% of the total 5000 configurations from the Na Cl + dataset. For our model and\n8/9 8\nCHGNet, a parameter-efficient approach was used where only the final two convolutional\nlayers, and the multilayer perceptron were updated. For MACE-MP-0, the entire model was\nsubjected to finetuning, as implemented in its own code. The finetuning protocol was as follows:\nthe subset was partitioned into training and validation sets with a 4:1 ratio. The optimization\nobjective was exclusively the atomic forces, using the Mean Squared Error (MSE) as the loss\n13/18\nfunction. All models were trained for 1000 epochs with a batch size of 10 and a learning rate\nof 5×10-4.\nAfter finetuning, PES of our model for Na Cl + clusters achieve a good agreement with the\n8/9 8\nDFT reference as shown in Fig. S11 b. This demonstrates that by accurately learning the forces,\nthe model implicitly reconstructs a highly precise and physically consistent potential energy\nsurface, capturing the subtle energetic details that were absent in the foundation model. On the\nother hand, following the finetuning procedure, both CHGNet and MACE-MP-0 exhibited a\nreduction in error as shown in Fig. S11 d and f. However, their final predicted potential energy\nsurfaces settled into an intermediate position, averaging the two distinct DFT reference surfaces\nfor two clusters. This indicates that while finetuning enables these foundation models to learn\na more refined potential energy surface to some extent, they remain unable to differentiate\nbetween the two distinct clusters. We may attribute this limitation to the absence of explicit\ndescriptors for the system’s net charge within their architectures, which prevents them from\ntweaking the unique energetic signature of each cluster.\n(a) (b)\n(c) (d)\n(e) (f)\nFig. S11 Potential energy surfaces of the Na Cl + cluster system. The surfaces were generated\n8/9 8\nby varying the position of the first Na atom relative to the second Na atom as described in Fig.\nS1. Left column: foundation models’ results of a, our model, c, CHGNet, and e, MACE-MP-0\n(large version). Right column: the finetuned models’ results of b, our model, d, CHGNet, and\nf, MACE-MP-0. Our model achieves good agreement with the DFT reference, demonstrating\nits superior capability to refine the potential energy surface for complex charged systems.\nHaving established the efficacy of finetuning for isolated clusters, we extended our\ninvestigation to more complex periodic systems to assess the models’ ability to capture subtle\n14/18\nsurface chemistry. We focused on the refined potential energy surface of an Au dimer\n2\nadsorbed on both Al-doped and undoped MgO surfaces, using the dataset provided by Ko et\nal. [4]. We found that message-passing-based foundation models, after finetuning, are indeed\ncapable of learning these precise potential energy surfaces. The finetuning protocol was\nconsistent with our previous tests: a subset comprising 20% of the configurations was randomly\nselected and partitioned into a 4:1 training-to-validation ratio. The models were then trained\nfor 1000 epochs with a 5×10⁻4 learning rate, exclusively minimizing the MSE loss on atomic\nforces. As summarized in Table S6, both finetuned CHGNet and our model successfully\nreproduce the DFT equilibrium bond lengths on both doped and undoped structures.\nNonetheless, our model achieves slightly higher fidelity with the DFT reference. In contrast,\nthe inferior performance of MACE-MP-0 may be attributed to its limited effective receptive\nfield arising from only two message-passing layers, which is likely insufficient to model the\ncomplex interactions modulated by the dopant.\nTable S6 The equilibrium bond lengths of Au dimer absorbed on both Al-doped and undoped\nAu on Al-doped MgO (Å) Au on MgO (Å)\n2 2\nTrue value of DFT 2.332 2.190\n4G-HDNNP +0.010 −0.003\nBespoke model −0.002 0.000\nOur foundation model −0.120 +0.038\nRelative\nFoundation CHGNet −0.098 +0.040\nvalue\ncompared\nFoundation MACE-MP-0 +0.028 +0.154\nto DFT\nFinetuned our model +0.002 −0.002\nFinetuned CHGNet −0.002 −0.004\nFinetuned MACE-MP-0 −0.068 +0.034\nIn summary, our work establishes a powerful two-step approach: first, pretraining on a\nuniversal dataset to build a robust foundation model, and second, performing targeted, data-\nefficient finetuning to achieve ab initio accuracy for specific applications. The success of this\nsecond step is, however, contingent on the inherent capabilities of the model’s architecture. A\nsuccessful foundation model must not only learn general physics but also possess the necessary\ndescriptors and structural flexibility—such as handling net charge and capturing long-range\neffects—to be effectively specialized for the diverse and challenging problems in modern\nmaterials science.”\nAnd the physical test of bespoke model of NaCl cluster is also provided in SI,\n15/18\nNa Cl +\n8 8\nNa Cl +\n9 8\n“Fig. S1 b, Potential energy surfaces of the Na Cl + cluster system, comparing DFT results\n8/9 8\nwith our model predictions. The surfaces were generated by varying the position of the first Na\natom relative to the second Na atom as shown in dashed lines, with Na and Cl atoms represented\nin yellow and green, respectively.”\nAs for the physical test of Au dimer absorption of the bespoke model, it has been already\n2\nprovided in Table S6. All the bespoke models achieve near-perfect agreements with DFT.\nWe demonstrated that our model is able to learn the correct underlying physics like in\ncharged NaCl clusters rather than simply fitting the data through increased expressivity. We\nare grateful to the reviewer for pushing us to perform these experiments, which have\nsignificantly strengthened the conclusions of our manuscript.\n“Additionally, while Table 1 shows the model’s performance on the 4G dataset, the\ncomparison is made only against a basic bespoke model from the original benchmark. Given\nthe field’s rapid progress, this is no longer a sufficient baseline. For instance, the recent model\nby Kim et al. (https://arxiv.org/abs/2412.15455) demonstrates significantly better performance\nfor some systems (e.g., an energy error of 0.073 meV/atom and force error of 0.008 eV/Å for\nthe gold dimer on MgO). A fair evaluation should include comparisons to such modern long-\nrange models.”\nAuthors’ reply:\nWe thank the reviewer for raising this important point. We agree that a fair evaluation\nrequires benchmarking against the state-of-the-art, and the comparison against a basic bespoke\nmodel is no longer sufficient. The work by Kim et al. is indeed a strong and relevant baseline.\nTo provide a more rigorous and contemporary evaluation, we have conducted benchmark\ncomparisons against the long-range capable model from Kim et al. (referred to as cace-lr) on\nKo et al. datasets. For the Ag +/− and Au -MgO systems, the polarizable charge equilibration\n3 2\n(PQEq) parameters for Ag and Au were reoptimized, allowing the model to learn the specific\nlong-range interactions more accurately. The results of this direct comparison are summarized\nin Table S10. In summary, the new benchmarks demonstrate that our methodology is highly\ncompetitive with the very recent reference although using the PQEq parameters invariant to\ngeometries. Future work could focus on two key avenues: first, developing dynamic models\nwhere PQEq parameters become a function of the local atomic environment, and second,\n16/18\nexploring non-charge-equilibration frameworks, such as those proposed by Kim et al., to\ncapture long-range physics with potentially greater efficiency and expressive power.\nIn SI,\n“Table S10 Comparison of our bespoke models and cace-lr models.\nDataset cace-lr This work\nEnergy (meV/atom) 0.73 0.44\nC H /C H +\n10 2 10 3\nForce (eV/Å) 0.037 0.023\nEnergy (meV/atom) 0.21 0.16\nNa Cl +\n8/9 8\nForce (eV/Å) 0.010 0.005\nEnergy (meV/atom) 0.162 4.87\nAg +/-\nEnergy (meV/atom) 0.073 0.13\nAu -MgO\nTable S11 Final PQEq parameters used for Ag and Au in the bespoke models for this section.\nAtom 𝜒9 (eV) 𝜂9 (eV) 𝑍 Radius (Å) 𝐾 (kcal/mol/Å2)\n:\nAg 3.50747 4.80698 1.00000 1.17602 221.26430\nAu 5.76708 6.91695 1.00000 0.55600 740.00000\n”\n“In summary, while the manuscript has improved, the main claims still require more rigorous\nand convincing validation. I would strongly recommend the authors perform standard\nphysicality tests and include comparisons to other state-of-the-art universal models to support\ntheir conclusions.”\nAuthors’ reply:\nWe sincerely thank the reviewer for their constructive feedback throughout this process. We\nfully agree that the main claims of a manuscript require rigorous and convincing validation. In\nresponse to the reviewer’s strong recommendation, we have performed an extensive series of\nexperiments to systematically address these points. We are confident that these additions now\nprovide the robust validation that was requested. Specifically, we have addressed the\nreviewer’s two primary recommendations including physical test and comparison with other\nuniversal models as follows.\n1. Fundamental interactions. We validated that our foundation model possesses the\ntransferability to describe correct long-range physics, successfully capturing potential\nenergy surfaces for both dimer interactions and systems with different charge states, as\ndescribed in the first question.\n2. Challenging validations. We then subjected our model to the challenging tests from the\nKo et al. dataset. Through finetuning, we demonstrated that our model, unlike others,\ncan learn the subtle differences in the potential energy surfaces of the Na Cl + and\n8 8\nNa Cl + clusters, proving it correctly interprets the system’s physics, as described in the\n9 8\nsecond question.\n17/18\n3. Unique polarization test. We introduced a novel physicality test that other universal\nmodels typically cannot address: response to an external electric field. Our model\nsuccessfully reproduces the DFT potential energy surface under varying field strengths,\ndemonstrating a genuine understanding of polarization.\nIn summary, we have systematically performed the requested physical tests and\ncomparisons with other universal models. These new results provide robust and convincing\nvalidation for our main claims. We are grateful to the reviewer for guidance, which has\nsignificantly strengthened the manuscript.\nFor the polarization test, details are provided in SI,\n“…The results are presented in the Fig. S9, which plots the system energy as a function of\nfield strength. Our model demonstrates good agreement with the DFT reference, accurately\nreproducing the energy curve across the entire field range from −2 to 2 V/Å. Notably, the\nperformance of our foundation model is on par with that of FIREANN [19], a state-of-the-art\nmodel that was specifically trained on this water system. This highlights the transferability of\nour model, as it correctly captures complex polarization physics without explicit training for\nthis task. In contrast, the CHGNet and MACE-MP-0 models are entirely unresponsive to the\nfield, predicting a constant energy. This represents a fundamental failure to capture one of the\nmost essential physical interactions, underscoring the superior predictive power and physical\nrealism of our approach.\nFig. S9 Evaluation of polarizable interactions in the periodic water model. The plot\ncompares the performance of our model against DFT [19] (reference), FIREANN [19],\nCHGNet, and MACE-MP-0.”\n18/18\nResponse to the reviewers\nWe are sincerely grateful to the two reviewers for their valuable time and insightful feedback\nduring this round of review. We have carefully considered all the comments and have provided\ndetailed, point-by-point responses in this document. The modifications to the manuscript are\nhighlighted in dark red.\nTable of Contents\nResponse to reviewer #2 ............................................................................................................ 2\nResponse to reviewer #3 ............................................................................................................ 3\n1/5\nResponse to reviewer #2\n“Remarks to the Author:\nThis revision of the manuscript has brought further improvements and the authors have now\nclarified most pertinent points raised by the reviewers. The manuscript is publishable.”\nAuthors’ reply:\nWe would like to express our gratitude for your positive and encouraging feedback on our\nrevised manuscript. We are very pleased to hear that you found our revisions effective and that\nthe manuscript is now considered publishable. Your expert comments and guidance have been\ninvaluable in improving the quality of our work, and we appreciate the time and effort you\ndedicated to the review process.\n“Remarks on code availability:\nI haven’t run the code but the repository exists with examples and instructions.”\nAuthors’ reply:\nWe appreciate you taking the time to check our code repository and for confirming that the\nexamples and instructions are available.\n2/5\nResponse to reviewer #3\n“The manuscript has improved substantially compared to the previous version. The results now\nlook mostly solid and trustworthy, and many of my earlier concerns have been addressed. I\nonly have a few remaining comments, including one point that still puzzles me.”\nAuthors’ reply:\nWe sincerely thank you for your thoughtful and constructive feedback on our revised\nmanuscript. We are very encouraged to hear that you find it has improved substantially and\nthat our results are now mostly solid and trustworthy. We have carefully noted your remaining\ncomments and will address them as follows.\n“Major comment:\nPhysical decay of long-range interactions (Fig. S8). I recommend plotting this figure on a log–\nlog scale in order to assess whether the correct power-law decay of interactions is reproduced.\nFrom visual inspection, it seems that the extracted exponent p is too small; for Coulombic\ninteractions, one expects p = 1. Could the authors comment on this discrepancy, especially\ngiven the claim that their model reproduces correct long-range physics?”\nAuthors’ reply:\nWe thank the reviewer for this insightful comment and valuable suggestion. The reviewer\nis correct that a log-log plot is the most appropriate way to analyse the power-law decay of\nlong-range interactions. Following this recommendation, we have generated a new log-log plot,\nnow included in the new Fig. S9.The reviewer correctly observed that the interaction does not\nfollow a simple 1/r decay across the entire distance range. This is expected in the distance range\nwe plotted in Fig. S8, as contributions from higher-order electrostatic terms cannot be ignored.\nTo show that our method captures the correct asymptotic behaviour, we extended the distance\nto more than 16 Å as shown in Fig. S9. Indeed, the curve from our model is tangent to the line\nwith slope = −1, confirming that it successfully reproduces the correct p = 1 power-law decay\nin the long-range limit.\nIn Supplementary Information, we have included the log-log plot of our model in\nextended distance range,\n3/5\n“Fig. S9 A log-log plot showing the interaction energy as a function of distance for an OH-\nOH dimer in an ionic state. The data from our model (solid green line) is compared with a\nreference line representing a pure 1/r decay (dashed black line).”\n“Minor comments:\n1. Fine-tuning strategy (S11). The additional comparison is very helpful and the results look\nconvincing. Could the authors comment on their fine-tuning strategy in more detail? For\nexample, is this full fine-tuning, multi-head fine-tuning, or a parameter-efficient method such\nas LoRA? In this context, I am also curious about the performance of the fine-tuned models on\nthe original training tasks—did the authors observe catastrophic forgetting?”\nAuthors’ reply:\nWe sincerely thank the reviewer for their positive feedback on our supplementary results\nand for asking these insightful questions regarding our finetuning strategy and the issue of\ncatastrophic forgetting.\nRegarding our finetuning strategy, we adopted “full finetuning” approach. For our model\nand CHGNet, we unfroze and updated the parameters of the final two convolutional layers and\nthe subsequent Multilayer Perceptron (MLP) layers. For the MACE-MP model, we finetuned\nall model parameters, following the officially recommended practice for this architecture. Our\ndecision to employ these methods, rather than using multi-head finetuning or parameter-\nefficient techniques such as LoRA, was driven by the high-precision nature of the Ko et al.\ndataset. We hypothesized that a more flexible approach—achieved by unfreezing some\nconvolutional and MLP layers—was necessary to provide the models with sufficient capacity\nto capture these subtle but critical features of the potential energy surfaces with high fidelity.\nIt is important to clarify that our primary goal in finetuning the foundation model on these\nchallenging systems was to demonstrate its inherent capability to capture these complex\ninteractions. That said, the high-precision nature of the potential energy surfaces in Ko et al.\ndatasets meant that we did observe a degree of catastrophic forgetting. However, the extent of\nthis effect was dependent on the similarity between the finetuning system and the domains\ncovered in the original training data. For the Au -MgO system, the task involves modelling a\n2\nsurface system, which does not include isolated, charged clusters. This represents a relatively\nmodest domain shift from our original MPtrj dataset. Consequently, we found that our model’s\nperformance on the original test set is acceptable after finetuning (MAE for forces and stresses\nwere 0.87 eV/Å and 0.64 GPa on the MPtrj test set, respectively). Conversely, the Na Cl +\n8/9 8\nsystems represent a significant departure from our original MPtrj dataset, as they involve\nisolated, charged clusters with unique electrostatic characteristics. Finetuning on this domain,\nwhich is substantially different, led to more pronounced catastrophic forgetting. Based on this\nanalysis, we conclude that for highly distinct systems like these charged clusters, the most\nreliable path to achieving a highly accurate description is to train a bespoke model.\nIn the Supplementary Information, we included details about the finetuning strategy,\n“…Regarding our finetuning strategy, we adopted a full finetuning approach. For our model\nand CHGNet, only the final two convolutional layers, and the multilayer perceptron were\nupdated. For MACE-MP-0, the entire model was subjected to finetuning, as implemented in\nits own code…”\n4/5\n“2. Comparison table (S10). Thank you for adding this comparison. The results look strong and\nconvincing. It also appears that the Ko dataset is becoming a widely used benchmark. I would\nlike to point out another recent paper (https://arxiv.org/abs/2507.19382) which further raises\nthe benchmark results; it may be worth discussing this in the context of related work.”\nAuthors’ reply:\nWe thank the reviewer for their positive feedback and for bringing the recent and highly\nrelevant LOREM paper [https://arxiv.org/abs/2507.19382] to our attention. We agree that its\nmethodology which treats charges as equivariant tensors provides an important contrast to\nother long-range approaches.Since the LOREM framework also does not require a charge\nequilibration process, we have integrated a description of this preprint into our manuscript’s\nDiscussion section, highlighting it among other methods that bypass this step.\nIn the main text,\n“…LOREM [84] models charges as equivariant tensors to capture orientation-dependent\ninteractions beyond the standard cutoff radius…\n84. Rumiantsev, E., et al., Learning Long-Range Representations with Equivariant Messages.\narXiv preprint arXiv:2507.19382, 2025.”\n5/5\nResponse to reviewer #3\n“I thank the authors for their revision and improvements. I am happy with the manuscript and\nsupport publication in its current form.”\nAuthors’ reply:\nWe sincerely thank the reviewer for their positive feedback and for endorsing the publication\nof our manuscript in its current form. We are grateful for their time and evaluation."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have addressed most of my concerns. One remaining point remains.\nQ6. “Figure 2 b-c: the authors should increase the x-axis range to 12 Å to actually see the long-range scaling behaviour.”\nI do not see \"other publications also do it like this\" as a valid argument. Understandably, the surrogate model cannot be more accurate than the data-generating DFT model. So, I'm not sure why the authors are trying to hide this.\nThus, the authors should add the results. If not in the main paper then at least in the supplementary information.", "reviewer_2": "Thanks to the authors for their reply and reviewing the manuscript. Factual errors have been mostly corrected, but there\nremains some open questions. These should be clarified and/or corrected before considering this manuscript any further.\nHere are some comments:\nTable S5 does not represent any significant improvements over other methods. A selling point could be increased\ntransferability etc. but no evidence for that.\nTable S6. I have no stake in MACE-MP-0 but I can immediately tell that those results can't be right and I am sure the MACE-\nMP-0 authors will challenge those as soon as they see this table. It seems extremely unlikely that in all the other metrics\nMACE-MP-0 does well and suddenly the errors become huge.\nAbout the forces/stresses. Since the forces depend on the charges, which depend on the positions, one must take the chain\nrule and calculate the appropriate term. Table R2 is not convincing at all - this compares to another implementation\npresumably using the same formulae. Did the authors compare their forces to forces obtained using finite differences?\nRegarding Fig. 2b. and R3 yes, it is clear now that due to the specific direction the electric field was applied this is indeed\nthe case. However, I disagree with the general statement \"This is a fundamental limitation of the traditional QEq-based\nmodel—it lacks the capability\nto respond to external electric fields when atomic positions remain fixed.\" If the electric field had components in the y or z\ndirections, the energy would change.\n\"The term “free energy landscape” is the standard terminology for describing how free\nenergy varies across configuration space.\" this is not quite right. This term describes the free energy as a function of a\ncollective variable or set of collective variables. But the manuscript text is OK.", "reviewer_3": "The authors have addressed many of the concerns raised in my previous review, and the manuscript has improved in\nseveral respects. For example, the addition of Table 1 is a valuable contribution. However, I still find that the very strong\nclaims made in the title and abstract—particularly the notion of a universal long-range model—are not sufficiently supported\nby the presented evidence.\nA central issue remains: while the model is trained to be universal, the authors do not convincingly demonstrate that it has\nactually learned the correct long-range physics, as opposed to simply benefitting from increased expressivity due to the\nadded long-range block. This point was raised in my initial review, but the response did not adequately address it.\nThe authors argue that:\n> “we cannot directly validate a universal model trained with VASP/PBE+U functional and pseudopotential basis sets on the\ndatasets proposed by Ko et al.”\nWhile technically true, this is not a convincing justification for avoiding such validation. The\nKo et al. datasets contain only a few hundred small structures, and recomputing them using the\nauthors’ own DFT setup (VASP/PBE+U) is computationally feasible—especially in the\ncontext of training a large-scale universal model. Recomputing small datasets for consistent\nevaluation is standard practice in the field.\nThe further claim that:\n> “so far, no universal model can reach ~1 meV/atom of precision.”\nmisses the point. The Ko et al. datasets include physically motivated test cases, such as displacement scans of NaCl and\nAu₂ systems, which reveal whether a model correctly captures energy and force trends with interatomic distance. The\nauthors should at least reproduce such physical validation tests with both their universal and bespoke models to support\ntheir claims.\nAdditionally, while Table 1 shows the model’s performance on the 4G dataset, the comparison is made only against a basic\nbespoke model from the original benchmark. Given the field’s rapid progress, this is no longer a sufficient baseline. For\ninstance, the recent model by Kim et al. (https://arxiv.org/abs/2412.15455) demonstrates significantly better performance\nfor some systems (e.g., an energy error of 0.073 meV/atom and force error of 0.008 eV/Å for the gold dimer on MgO). A fair\nevaluation should include comparisons to such modern long-range models.\nIn summary, while the manuscript has improved, the main claims still require more rigorous and convincing validation. I\nwould strongly recommend the authors perform standard physicality tests and include comparisons to other state-of-the-art\nuniversal models to support your conclusions."}, "author_response": "Response to the reviewers\nWe sincerely thank the three reviewers for their continued time and effort in this second round\nof review. Their additional feedback has been invaluable, and we have revised the manuscript\nto address all remaining points. We have provided point-by-point responses to the questions\nraised by the reviewers in this document. The modifications to the manuscript are highlighted\nin dark red.\nTable of Contents\nResponse to reviewer #1 ............................................................................................................ 2\nResponse to reviewer #2 ............................................................................................................ 3\nResponse to reviewer #3 .......................................................................................................... 10\n1/18\nResponse to reviewer #1\n“The authors have addressed most of my concerns. One remaining point remains.\nQ6. “Figure 2 b-c: the authors should increase the x-axis range to 12 Å to actually see the long-\nrange scaling behaviour.”\nI do not see “other publications also do it like this” as a valid argument. Understandably, the\nsurrogate model cannot be more accurate than the data-generating DFT model. So, I’m not sure\nwhy the authors are trying to hide this. Thus, the authors should add the results. If not in the\nmain paper then at least in the supplementary information.”\nAuthors’ reply:\nWe sincerely thank the reviewer for acknowledging our previous revisions. We also\nappreciate your valuable suggestion, which has prompted us to further clarify an important\nscientific issue. We agree with your point that simply citing ‘standard practice’ of other\npublications is an insufficient justification, and we apologize for the lack of clarity in our\nprevious response. Your suggestion to extend the x-axis to 12 Å has indeed provided a much\nbetter perspective for a comprehensive evaluation of our model’s performance.\nIn the main text, we added the entire distance up to 12 Å for comparison,\n“Fig. 2 a, The interaction energies for Na-Na (green) and Na-Cl (purple) dimers predicted\nby our model, w/o-lr model, and reference calculations using DFT.\n…We analysed interaction energies for Na-Na and Na-Cl dimers across interatomic\ndistances ranging from 1.8 Å to 12.0 Å…”\nWe believe these additions fully address the reviewer’s point and have made the manuscript\nstronger. We thank the reviewer again for pushing us to include this important validation.\n2/18\nResponse to reviewer #2\n“Thanks to the authors for their reply and reviewing the manuscript. Factual errors have been\nmostly corrected, but there remains some open questions. These should be clarified and/or\ncorrected before considering this manuscript any further. Here are some comments:”\nAuthors’ reply:\nWe thank the reviewer for their feedback and for taking the time to review our revised\nmanuscript. We appreciate the acknowledgment that some factual errors have been corrected.\nWe have carefully considered the remaining questions and will address each point thoroughly\nin the updated manuscript to provide the necessary clarifications and corrections.\n“Table S5 does not represent any significant improvements over other methods. A selling point\ncould be increased transferability etc. but no evidence for that.”\nAuthors’ reply:\nWe thank the reviewer for this constructive and valuable feedback. We agree that the\noriginal manuscript did not sufficiently demonstrate the transferability of our method, which is\nindeed a key distinctive feature. Following the reviewer’s suggestion, we have performed\nadditional analyses and included substantial new evidence to explicitly showcase the superior\ntransferability of our approach. Specifically,\n1. Transferability of long-range physics. First, we demonstrate our model’s ability to\naccurately describe dimer interactions across multiple distance scales. As shown in the\nnew Fig. S7, our model correctly reproduces the potential energy surfaces for both Na-\nNa and Na-Cl dimers. In addition, our model outperforms CHGNet and MACE-MP-0\non both short-range and long-range distances. Second, we investigated the OH-OH\nsystem in both its neutral (weakly interacting) and ionic (charge-repulsive) states. We\nshowcase our model’s unique ability to distinguish between different charge states and\ncorrectly model their corresponding long-range physics—a fundamental challenge for\nother universal models.\n2. Transferability of polarization. To further provide evidence of our model’s\ntransferability in terms of polarization, we tested its ability to predict responses to\nexternal electric fields—a crucial test of physical fidelity that reveals limitations in other\nuniversal models. First, in a periodic system of 64 water molecules, our model accurately\nreproduced the quadratic energy dependence under a varying electric field, showing\noutstanding agreement with DFT reference calculations (Fig. S9). Notably, its\nperformance matched that of FIREANN, a state-of-the-art model specifically trained on\nthat system, whereas leading universal models like CHGNet and MACE-MP-0 were\ncompletely unresponsive to the field. In addition, we validated our model’s ability to\npredict molecular polarizability for 7,211 diverse molecules from the QM-7b dataset, a\ntask beyond the capabilities of other universal models. Our model achieved a low mean\nabsolute error of 4.57 a.u. against high-level theory reference values (Fig. S10).\n3/18\nThese additions provide the quantitative evidence for transferability of our foundation\nmodel. We are confident that these results now clearly establish a significant improvement over\nother methods, not just in the metrics of Table S5, but in the crucial dimension of physical\ntransferability.\nIn SI, we detail the transferability,\n“In comparison with existing universal models, our model demonstrates transferability in\ndescribing dimer interactions. As mentioned in the main text, regarding the potential energy\nsurface capturing of Na-Cl and Na-Na dimers, we have discovered that our model not only\naccurately reproduces the long-range behaviour predicted by DFT/PBE [18] calculations, but\nalso effectively captures the potential energy surfaces of these dimers in the short-range region\ncompared to CHGNet and MACE-MP-0, as illustrated in Fig. S7. The enhanced transferability\nof our model across different dimers suggests its potential utility for a broad spectrum of\napplications requiring accurate representation of interatomic forces across multiple distance\nscales.\nFig. S7 The interaction energies for Na-Na (green) and Na-Cl (purple) dimers predicted by\nour model, CHGNet [7], MACE-MP-0 (large version) [11], and reference DFT calculations.\nFor the neutral state (dashed lines), which represents the physical scenario of a weak, non-\nbonded interaction, the DFT reference shows a very shallow potential well with relative\nenergies. Our model demonstrates good agreement with the DFT reference, accurately\ncapturing the subtle attractive interaction. In contrast, MACE-MP-0 and CHGNet struggle to\ncapture this long-range interaction; they appear less transferable than our model as they fail to\nreproduce the precise shape and depth of the DFT curve. For the ionic state (solid lines), which\nrepresents a strong, long-range Coulombic repulsion between charged species, the DFT\nreference shows a significant positive (repulsive) energy that decays slowly with distance. Our\nmodel correctly identifies the ionic state and qualitatively reproduces the long-range repulsive\nenergy curve from DFT. Despite a slight overestimation of the energy, the model captures the\nphysical trend and magnitude of the interaction across the entire distance range. Conversely,\nMACE-MP-0 and CHGNet are completely unable to distinguish the ionic from the neutral state.\nThis represents a fundamental physical problem, indicating that their underlying architectures\ncannot capture charge-state-dependent long-range physics.\nFig. S8 Relative interaction energy versus distance for an OH-OH system in its neutral\n(dashed lines) and ionic (solid lines) states. Our model’s predictions are compared with DFT,\nMACE-MP-0, and CHGNet. The results demonstrate that our model can reproduce the DFT\nreference for both the weak interaction of the neutral state and, crucially, the long-range\nCoulombic repulsion of the ionic state. Universal models like MACE-MP-0 and CHGNet fail\nto capture the long-range physics of the ionic state.\n5/18\nA crucial test of a model’s physical fidelity is its response to external fields. We evaluated\nour model’s ability to predict the energy response of a system to a uniform electric field. We\nadopted a periodic system containing 64 water molecules with each cell length measured\n12.4185 Å, as described in reference [19]. All the water molecules are within yz-plane. The\nfirst and third layers are same while the second and fourth are same. An electric field was then\napplied along the x-axis with varying strength.\nFig. S9 Evaluation of polarizable interactions in the periodic water model. The plot\ncompares the performance of our model against DFT [19] (reference), FIREANN [19],\nCHGNet, and MACE-MP-0.\nThe results are presented in the Fig. S9, which plots the system energy as a function of field\nstrength. Our model demonstrates good agreement with the DFT reference, accurately\nreproducing the energy curve across the entire field range from −2 to 2 V/Å. Notably, the\nperformance of our foundation model is on par with that of FIREANN [19], a state-of-the-art\nmodel that was specifically trained on this water system. This highlights the transferability of\nour model, as it correctly captures complex polarization physics without explicit training for\nthis task. In contrast, the CHGNet and MACE-MP-0 models are entirely unresponsive to the\nfield, predicting a constant energy. This represents a fundamental failure to capture one of the\nmost essential physical interactions, underscoring the superior predictive power and physical\nrealism of our model.\nFig. S10 Comparison of the predicted isotropic polarizability (𝛼 ) against reference coupled-\n!\"#\ncluster theory including single and double excitations values for 7,211 molecules from the QM-\n7b dataset [20]. The model achieves a mean absolute error (MAE) of 4.57 a.u., demonstrating\nits strong capability to generalize and predict complex electronic response properties accurately.\nConventional universal machine learning potentials like MACE-MP-0 and CHGNet are\ntypically trained on field-free energies and forces, and their architectures lack the components\nto describe fundamental electrostatic properties like the dipole moment. Consequently, they are\ninherently incapable of predicting response properties such as polarizability, which is defined\nby the change in energy or dipole moment under an applied field. We further validated the\npowerful transferability of our model by assessing its ability to predict molecular polarizability,\na key electronic response property. To accomplish this, we curated a challenging test set\ncomprising 7,211 diverse small molecules from the QM-7b dataset [20-23] … This involved\ncalculating the system’s dipole response to a small, externally applied electric field—a physical\nprocess the model was not explicitly trained to reproduce. The results of this evaluation are\n6/18\npresented in the Fig. S10. The plot, which compares our model’s predictions against the\nreference values, demonstrates good agreement. The model achieves a MAE of just 4.57 a.u.\nacross this vast and chemically diverse dataset. This high level of accuracy, achieved without\ndirect training on polarizability data, provides compelling evidence for the model’s strong\ngeneralization capabilities and its ability to capture the fundamental physics governing\nmolecular electronic response, a capability that lies beyond the scope of other universal\nmachine learning potential models.”\n“Table S6. I have no stake in MACE-MP-0 but I can immediately tell that those results can’t\nbe right and I am sure the MACE-MP-0 authors will challenge those as soon as they see this\ntable. It seems extremely unlikely that in all the other metrics MACE-MP-0 does well and\nsuddenly the errors become huge.”\nAuthors’ reply:\nWe thank the reviewer for their careful scrutiny of the MACE-MP-0 results and for raising\nthis important point. We understand the concern, as the reported errors are indeed strikingly\nlarge. The large errors reported are indeed correct and stem from the specific nature of the\nbenchmark, which is designed to evaluate models under challenging thermodynamic\nperturbation datasets. The large version of MACE-MP-0 trained on the MPtrj dataset\nconsistently demonstrates limitations when evaluated on such extreme and challenging test\ncases.\nFurthermore, to ensure a rigorous and fair comparison, we conducted a thorough review of\nthe evaluation metrics. We discovered that the definition for MAE of forces and stresses used\nin the MatterSim paper (from which we cited their results) differs from the definition used\nduring their training. We have therefore re-evaluated all models using a single, unified\ndefinition to ensure consistency.\nTo guarantee full transparency and reproducibility, we have made our complete evaluation\nscript for MACE-MP-0 available on Figshare at the following link:\nhttps://figshare.com/s/7ecc56e3538d83c83010. We believe this clarification and the provided\nscript will fully address the reviewer’s concerns about the validity of our findings.\n“About the forces/stresses. Since the forces depend on the charges, which depend on the\npositions, one must take the chain rule and calculate the appropriate term. Table R2 is not\nconvincing at all - this compares to another implementation presumably using the same\nformulae. Did the authors compare their forces to forces obtained using finite differences?”\nAuthors’ reply:\nWe sincerely thank the reviewer for their insight on this critical point. After a careful re-\nderivation, we identified that it has a Hellmann-Feynman-like term which will vanish as you\nmentioned before. To provide a definitive and convincing validation, we have also followed\nyour suggestion and compared our auto-gradient forces against the forces obtained using finite\n7/18\ndifferences for a water molecule as described in Table R1. As depicted in Table R2, the auto-\ngradient forces and finite differences results show excellent agreement, confirming the\ncorrectness of our implementation. We are grateful for the reviewer’s suggestion, which has\nsignificantly improved the rigor of our work.\nThe complete and detailed derivation has been added to the SI,\n“The total forces of our model are,\n𝜕𝐸 𝜕𝐸 𝜕𝑞\n%#& %#& '\n𝑭\n$\n= −\n𝜕𝒓\n−(\n𝜕𝑞 𝜕𝒓\n, (S9)\n$ ' 𝒊\n'\nwhere 𝐸 , 𝒓 , and 𝑞 are total potential energy as described in Eq. 1, atomic coordinates, and\n%#& $ $\npartial charges, respectively. The Lagrange multiplier method (𝐸 = 𝐸 −𝜆∑𝑞 ) requires\n) %#& $\n𝜕𝐸\n)\n− = 0, (S10)\n𝜕𝑞\n'\nso that,\n𝜕𝐸\n%#&\n= 𝜆. (S11)\n𝜕𝑞\n'\nThus, the forces are simplified as,\n𝜕𝐸 𝜕𝑞\n%#& '\n𝑭\n$\n= −\n𝜕𝒓\n−𝜆(\n𝜕𝒓\n. (S12)\n$ 𝒊\n'\nAnd with the constraint,\n∑𝑞 = 𝑄, (S13)\n'\n*+\nwe have = 0 since the net charge 𝑄 is a constant. Finally, the atomic forces are,\n*𝒓𝒊\n𝜕𝐸\n%#&\n𝑭 = − . (S14)\n$ 𝜕𝒓\n$\n…”\nTable R1 Water molecule coordinates.\nx (Å) y (Å) z (Å)\nO −1.464 0.099 0.300\nH −1.956 0.624 −0.340\nH −1.797 −0.799 0.206\nTable R2 Comparison of the auto-gradient and finite differences forces on the water molecule.\nForces (eV/Å)\n𝐹./012!345 𝐹6/&# 𝐹./012!345 𝐹6/&# 𝐹./012!345 𝐹6/&#\n- - 7 7 8 8\nO −0.13566780 −0.13566779 −0.06141234 −0.06141233 −0.12068319 −0.12068319\nH 0.06315283 0.06315281 0.07240831 0.07240830 0.04431891 0.04431891\nH 0.07251500 0.07251498 −0.01099598 −0.01099597 0.07636429 0.07636429\n8/18\n“Regarding Fig. 2b. and R3 yes, it is clear now that due to the specific direction the electric\nfield was applied this is indeed the case. However, I disagree with the general statement “This\nis a fundamental limitation of the traditional QEq-based model—it lacks the capability to\nrespond to external electric fields when atomic positions remain fixed.” If the electric field had\ncomponents in the y or z directions, the energy would change.”\nAuthors’ reply:\nWe thank the reviewer for this insightful comment and for pointing out the imprecision in\nour original statement. The reviewer is absolutely correct. A traditional QEq-based model\nwould indeed show an energy change if the electric field applied has a component parallel to\nthe direction of molecular dipole moment. We are grateful to the reviewer for pointing this out.\nWe have revised it as in the main text,\n“…Conventional QEq-based models have limitations in their response to an external electric\nfield applied orthogonally to the molecular dipole…”\n““The term “free energy landscape” is the standard terminology for describing how free\nenergy varies across configuration space.” this is not quite right. This term describes the free\nenergy as a function of a collective variable or set of collective variables. But the manuscript\ntext is OK.”\nAuthors’ reply:\nWe thank the reviewer for the valuable clarification on the term “free energy landscape”\nand for the approval of our manuscript’s content. Reviewer’s rigorous approach and expertise\nin mathematical physics are, once again, truly impressive.\n9/18\nResponse to reviewer #3\n“The authors have addressed many of the concerns raised in my previous review, and the\nmanuscript has improved in several respects. For example, the addition of Table 1 is a valuable\ncontribution. However, I still find that the very strong claims made in the title and abstract—\nparticularly the notion of a universal long-range model—are not sufficiently supported by the\npresented evidence.”\nAuthors’ reply:\nWe sincerely thank the reviewer for acknowledging the improvements we made in the\nmanuscript, particularly for recognizing the value of the newly added Table 1. Our model can\ndemonstrate the capability to capture accurate long-range physics in certain scenarios. For\ninstance, it accurately models interactions across various distance scales, reproducing the\npotential energy surfaces for Na-Na and Na-Cl dimers (new Fig. S7) while outperforming\nbenchmarks like CHGNet and MACE-MP-0. Crucially, our model also addresses a key\nlimitation of other universal potentials by successfully distinguishing between different charge\nstates. This is highlighted by its ability to correctly model the distinct long-range physics of the\nOH-OH system in both its neutral (weakly interacting) and ionic (charge-repulsive) forms.\nHowever, we recognize that for more challenging and diverse systems, such as the dataset\nproposed by Ko et al., our model serves as a foundation model that requires further finetuning\nto achieve high accuracy. We will elaborate on this specific point in our response to the\nreviewer’s next question. We now position our work as a “foundation model”. This term more\naccurately reflects our contribution: a robust, physically-informed starting point that is\ndesigned to be efficiently finetuned to achieve ab initio accuracy for specific, challenging\napplications. This two-step approach—pretraining and finetuning—represents, in our view, a\nmore pragmatic and powerful paradigm for modern materials modelling. The revised title is\nnow,\n“A foundation machine learning potential with polarizable long-range interactions for\nmaterials modelling”\nAnd in the abstract, we further address the foundation model and finetuning,\n“…The foundation model, trained across the periodic table up to Pu, demonstrates strong\nperformance across key materials modelling challenges… Furthermore, we show that as a\nfoundation model, it can be efficiently finetuned to achieve high-level accuracy for specific\nchallenging systems.”\nDetailed physical test of long-range physics in SI,\n“In comparison with existing universal models, our model demonstrates transferability in\ndescribing dimer interactions. As mentioned in the main text, regarding the potential energy\nsurface capturing of Na-Cl and Na-Na dimers, we have discovered that our model not only\naccurately reproduces the long-range behaviour predicted by DFT/PBE [18] calculations, but\nalso effectively captures the potential energy surfaces of these dimers in the short-range region\ncompared to CHGNet and MACE-MP-0, as illustrated in Fig. S7. The enhanced transferability\nof our model across different dimers suggests its potential utility for a broad spectrum of\n10/18\napplications requiring accurate representation of interatomic forces across multiple distance\nscales.\nFig. S7 The interaction energies for Na-Na (green) and Na-Cl (purple) dimers predicted by\nour model, CHGNet [7], MACE-MP-0 (large version) [11], and reference DFT calculations.\nAnother example that showcases the transferability of our model is its ability to distinguish\nbetween different charge states and accurately simulate their corresponding long-range\ninteractions. We selected the OH-OH system as a case study, investigating the potential energy\ncurves when the system consists of two neutral hydroxyl radicals versus when it carries a single\nnegative charge. As illustrated in Fig. S8, the plot compares the performance of our model\nagainst DFT (the ground truth), MACE-MP-0, and CHGNet for the system in both “ionic” and\n“neutral” states as a function of distance.\nFig. S8 Relative interaction energy versus distance for an OH-OH system in its neutral\n(dashed lines) and ionic (solid lines) states. Our model’s predictions are compared with DFT,\nMACE-MP-0, and CHGNet. The results demonstrate that our model can reproduce the DFT\nreference for both the weak interaction of the neutral state and, crucially, the long-range\nCoulombic repulsion of the ionic state. Universal models like MACE-MP-0 and CHGNet fail\nto capture the long-range physics of the ionic state.\n11/18\nFor the neutral state (dashed lines), which represents the physical scenario of a weak, non-\nbonded interaction, the DFT reference shows a very shallow potential well with relative\nenergies. Our model demonstrates good agreement with the DFT reference, accurately\ncapturing the subtle attractive interaction. In contrast, MACE-MP-0 and CHGNet struggle to\ncapture this long-range interaction; they appear less transferable than our model as they fail to\nreproduce the precise shape and depth of the DFT curve. For the ionic state (solid lines), which\nrepresents a strong, long-range Coulombic repulsion between charged species, the DFT\nreference shows a significant positive (repulsive) energy that decays slowly with distance. Our\nmodel correctly identifies the ionic state and qualitatively reproduces the long-range repulsive\nenergy curve from DFT. Despite a slight and consistent overestimation of the energy, the model\ncaptures the physical trend and magnitude of the interaction across the entire distance range.\nConversely, MACE-MP-0 and CHGNet are completely unable to distinguish the ionic from\nthe neutral state. This represents a fundamental physical problem, indicating that their\nunderlying architectures cannot capture charge-state-dependent long-range physics.”\n“A central issue remains: while the model is trained to be universal, the authors do not\nconvincingly demonstrate that it has actually learned the correct long-range physics, as opposed\nto simply benefitting from increased expressivity due to the added long-range block. This point\nwas raised in my initial review, but the response did not adequately address it.\nThe authors argue that:\n> “we cannot directly validate a universal model trained with VASP/PBE+U functional and\npseudopotential basis sets on the datasets proposed by Ko et al.”\nWhile technically true, this is not a convincing justification for avoiding such validation. The\nKo et al. datasets contain only a few hundred small structures, and recomputing them using the\nauthors’ own DFT setup (VASP/PBE+U) is computationally feasible—especially in the\ncontext of training a large-scale universal model. Recomputing small datasets for consistent\nevaluation is standard practice in the field.\nThe further claim that: > “so far, no universal model can reach ~1 meV/atom of precision.”\nmisses the point. The Ko et al. datasets include physically motivated test cases, such as\ndisplacement scans of NaCl and Au₂ systems, which reveal whether a model correctly captures\nenergy and force trends with interatomic distance. The authors should at least reproduce such\nphysical validation tests with both their universal and bespoke models to support their claims.”\nAuthors’ reply:\nWe sincerely thank the reviewer for this critical and insightful feedback. We agree that our\nprevious response did not adequately address the central issue: whether our model has learned\nthe correct long-range physics or is merely benefiting from increased architectural expressivity.\nWe agree that performing validation on the potential energy surfaces of NaCl clusters and Au\n2\nsystems is indispensable.\nAs the reviewer previously pointed out, a pretrained foundation model should demonstrate\nits capacity for learning challenging physics through finetuning. We first evaluated various\nmodels on the potential energy surfaces (PES) of the Na Cl ⁺ and Na Cl ⁺ charged clusters\n8 8 9 8\nfrom the Ko et al. dataset. As illustrated in Fig. S11 (a, c, and e), all three foundation models\n12/18\n(ours, CHGNet, and MACE-MP-0) fail to reproduce the distinct DFT potential energy surfaces\nfor these two systems. This confirms the reviewer’s scepticism about the out-of-the-box\nperformance of universal models on such specific, challenging tasks. However, even at this\nstage, our model showed a nascent ability to distinguish between the two clusters, whereas the\nothers did not. We then finetuned all three models using a small subset (20%) of the\nconfigurations. The same algorithm was used for all models, and only the mean squared force\nloss was used to update the weights. After finetuning, PES of our finetuned model for both\nNa Cl ⁺ and Na Cl ⁺ achieve good agreements with the DFT reference (Fig. S11 b). In contrast,\n8 8 9 8\nwhile finetuning reduced the overall error for CHGNet and MACE-MP-0, their final predicted\nenergy curves settled into an intermediate position, averaging the two distinct DFT reference\ncurves (Fig. S11 d and f). They remain unable to differentiate between the two distinct clusters.\nWe may attribute this limitation to the absence of explicit descriptors for the system’s net\ncharge within their architectures, which prevents them from tweaking the unique energetic\nsignature of each cluster. We then extended our investigation to the complex periodic system\nof an Au dimer adsorbed on both Al-doped and undoped MgO surfaces. Following a similar\n2\nfinetuning protocol, we found that message-passing-based models are indeed capable of\nlearning these precise potential energy surfaces. As summarized in Table S6, both finetuned\nCHGNet and our model successfully reproduce the DFT equilibrium bond lengths on both\nsurfaces, with our model achieving slightly higher fidelity. The inferior performance of\nMACE-MP-0 in this case may be attributed to its limited effective receptive field (arising from\nonly two message-passing layers).\nOur new results not only validate this approach but also reveal the critical importance of\narchitectural design—specifically, the inclusion of net charge as a descriptor—in successfully\ncapturing complex electrostatic interactions.\nIn SI,\n“While our model provides a foundational understanding, achieving ab initio accuracy for\nspecific and challenging systems like charged clusters necessitates a refinement of the potential\nenergy surface. To rigorously assess the physical fidelity of the foundation models, we\nconducted an evaluation on the potential energy surfaces (PES) of charged sodium chloride\nclusters from the dataset proposed by Ko et al. [4]. The initial test involved scanning the PES\nof the Na Cl + cluster systems along a coordinate defined by varying the distance between\n8/9 8\ntwo sodium atoms. As illustrated in Fig. S11 (a, c, and e), a comparison with DFT calculations\nreveals that all three foundation models fail to accurately describe the PES for these systems.\nNotably, our model demonstrated a potential to distinguish between the two clusters due to the\ncharge state descriptor, whereas other models did not.\nWe acknowledge that capturing such a refined PES requires model specialization, so we\nperformed a finetuning procedure on all three models. A subset was created by randomly\nselecting 20% of the total 5000 configurations from the Na Cl + dataset. For our model and\n8/9 8\nCHGNet, a parameter-efficient approach was used where only the final two convolutional\nlayers, and the multilayer perceptron were updated. For MACE-MP-0, the entire model was\nsubjected to finetuning, as implemented in its own code. The finetuning protocol was as follows:\nthe subset was partitioned into training and validation sets with a 4:1 ratio. The optimization\nobjective was exclusively the atomic forces, using the Mean Squared Error (MSE) as the loss\n13/18\nfunction. All models were trained for 1000 epochs with a batch size of 10 and a learning rate\nof 5×10-4.\nAfter finetuning, PES of our model for Na Cl + clusters achieve a good agreement with the\n8/9 8\nDFT reference as shown in Fig. S11 b. This demonstrates that by accurately learning the forces,\nthe model implicitly reconstructs a highly precise and physically consistent potential energy\nsurface, capturing the subtle energetic details that were absent in the foundation model. On the\nother hand, following the finetuning procedure, both CHGNet and MACE-MP-0 exhibited a\nreduction in error as shown in Fig. S11 d and f. However, their final predicted potential energy\nsurfaces settled into an intermediate position, averaging the two distinct DFT reference surfaces\nfor two clusters. This indicates that while finetuning enables these foundation models to learn\na more refined potential energy surface to some extent, they remain unable to differentiate\nbetween the two distinct clusters. We may attribute this limitation to the absence of explicit\ndescriptors for the system’s net charge within their architectures, which prevents them from\ntweaking the unique energetic signature of each cluster.\n(a) (b)\n(c) (d)\n(e) (f)\nFig. S11 Potential energy surfaces of the Na Cl + cluster system. The surfaces were generated\n8/9 8\nby varying the position of the first Na atom relative to the second Na atom as described in Fig.\nS1. Left column: foundation models’ results of a, our model, c, CHGNet, and e, MACE-MP-0\n(large version). Right column: the finetuned models’ results of b, our model, d, CHGNet, and\nf, MACE-MP-0. Our model achieves good agreement with the DFT reference, demonstrating\nits superior capability to refine the potential energy surface for complex charged systems.\nHaving established the efficacy of finetuning for isolated clusters, we extended our\ninvestigation to more complex periodic systems to assess the models’ ability to capture subtle\n14/18\nsurface chemistry. We focused on the refined potential energy surface of an Au dimer\n2\nadsorbed on both Al-doped and undoped MgO surfaces, using the dataset provided by Ko et\nal. [4]. We found that message-passing-based foundation models, after finetuning, are indeed\ncapable of learning these precise potential energy surfaces. The finetuning protocol was\nconsistent with our previous tests: a subset comprising 20% of the configurations was randomly\nselected and partitioned into a 4:1 training-to-validation ratio. The models were then trained\nfor 1000 epochs with a 5×10⁻4 learning rate, exclusively minimizing the MSE loss on atomic\nforces. As summarized in Table S6, both finetuned CHGNet and our model successfully\nreproduce the DFT equilibrium bond lengths on both doped and undoped structures.\nNonetheless, our model achieves slightly higher fidelity with the DFT reference. In contrast,\nthe inferior performance of MACE-MP-0 may be attributed to its limited effective receptive\nfield arising from only two message-passing layers, which is likely insufficient to model the\ncomplex interactions modulated by the dopant.\nTable S6 The equilibrium bond lengths of Au dimer absorbed on both Al-doped and undoped\nAu on Al-doped MgO (Å) Au on MgO (Å)\n2 2\nTrue value of DFT 2.332 2.190\n4G-HDNNP +0.010 −0.003\nBespoke model −0.002 0.000\nOur foundation model −0.120 +0.038\nRelative\nFoundation CHGNet −0.098 +0.040\nvalue\ncompared\nFoundation MACE-MP-0 +0.028 +0.154\nto DFT\nFinetuned our model +0.002 −0.002\nFinetuned CHGNet −0.002 −0.004\nFinetuned MACE-MP-0 −0.068 +0.034\nIn summary, our work establishes a powerful two-step approach: first, pretraining on a\nuniversal dataset to build a robust foundation model, and second, performing targeted, data-\nefficient finetuning to achieve ab initio accuracy for specific applications. The success of this\nsecond step is, however, contingent on the inherent capabilities of the model’s architecture. A\nsuccessful foundation model must not only learn general physics but also possess the necessary\ndescriptors and structural flexibility—such as handling net charge and capturing long-range\neffects—to be effectively specialized for the diverse and challenging problems in modern\nmaterials science.”\nAnd the physical test of bespoke model of NaCl cluster is also provided in SI,\n15/18\nNa Cl +\n8 8\nNa Cl +\n9 8\n“Fig. S1 b, Potential energy surfaces of the Na Cl + cluster system, comparing DFT results\n8/9 8\nwith our model predictions. The surfaces were generated by varying the position of the first Na\natom relative to the second Na atom as shown in dashed lines, with Na and Cl atoms represented\nin yellow and green, respectively.”\nAs for the physical test of Au dimer absorption of the bespoke model, it has been already\n2\nprovided in Table S6. All the bespoke models achieve near-perfect agreements with DFT.\nWe demonstrated that our model is able to learn the correct underlying physics like in\ncharged NaCl clusters rather than simply fitting the data through increased expressivity. We\nare grateful to the reviewer for pushing us to perform these experiments, which have\nsignificantly strengthened the conclusions of our manuscript.\n“Additionally, while Table 1 shows the model’s performance on the 4G dataset, the\ncomparison is made only against a basic bespoke model from the original benchmark. Given\nthe field’s rapid progress, this is no longer a sufficient baseline. For instance, the recent model\nby Kim et al. (https://arxiv.org/abs/2412.15455) demonstrates significantly better performance\nfor some systems (e.g., an energy error of 0.073 meV/atom and force error of 0.008 eV/Å for\nthe gold dimer on MgO). A fair evaluation should include comparisons to such modern long-\nrange models.”\nAuthors’ reply:\nWe thank the reviewer for raising this important point. We agree that a fair evaluation\nrequires benchmarking against the state-of-the-art, and the comparison against a basic bespoke\nmodel is no longer sufficient. The work by Kim et al. is indeed a strong and relevant baseline.\nTo provide a more rigorous and contemporary evaluation, we have conducted benchmark\ncomparisons against the long-range capable model from Kim et al. (referred to as cace-lr) on\nKo et al. datasets. For the Ag +/− and Au -MgO systems, the polarizable charge equilibration\n3 2\n(PQEq) parameters for Ag and Au were reoptimized, allowing the model to learn the specific\nlong-range interactions more accurately. The results of this direct comparison are summarized\nin Table S10. In summary, the new benchmarks demonstrate that our methodology is highly\ncompetitive with the very recent reference although using the PQEq parameters invariant to\ngeometries. Future work could focus on two key avenues: first, developing dynamic models\nwhere PQEq parameters become a function of the local atomic environment, and second,\n16/18\nexploring non-charge-equilibration frameworks, such as those proposed by Kim et al., to\ncapture long-range physics with potentially greater efficiency and expressive power.\nIn SI,\n“Table S10 Comparison of our bespoke models and cace-lr models.\nDataset cace-lr This work\nEnergy (meV/atom) 0.73 0.44\nC H /C H +\n10 2 10 3\nForce (eV/Å) 0.037 0.023\nEnergy (meV/atom) 0.21 0.16\nNa Cl +\n8/9 8\nForce (eV/Å) 0.010 0.005\nEnergy (meV/atom) 0.162 4.87\nAg +/-\nEnergy (meV/atom) 0.073 0.13\nAu -MgO\nTable S11 Final PQEq parameters used for Ag and Au in the bespoke models for this section.\nAtom 𝜒9 (eV) 𝜂9 (eV) 𝑍 Radius (Å) 𝐾 (kcal/mol/Å2)\n:\nAg 3.50747 4.80698 1.00000 1.17602 221.26430\nAu 5.76708 6.91695 1.00000 0.55600 740.00000\n”\n“In summary, while the manuscript has improved, the main claims still require more rigorous\nand convincing validation. I would strongly recommend the authors perform standard\nphysicality tests and include comparisons to other state-of-the-art universal models to support\ntheir conclusions.”\nAuthors’ reply:\nWe sincerely thank the reviewer for their constructive feedback throughout this process. We\nfully agree that the main claims of a manuscript require rigorous and convincing validation. In\nresponse to the reviewer’s strong recommendation, we have performed an extensive series of\nexperiments to systematically address these points. We are confident that these additions now\nprovide the robust validation that was requested. Specifically, we have addressed the\nreviewer’s two primary recommendations including physical test and comparison with other\nuniversal models as follows.\n1. Fundamental interactions. We validated that our foundation model possesses the\ntransferability to describe correct long-range physics, successfully capturing potential\nenergy surfaces for both dimer interactions and systems with different charge states, as\ndescribed in the first question.\n2. Challenging validations. We then subjected our model to the challenging tests from the\nKo et al. dataset. Through finetuning, we demonstrated that our model, unlike others,\ncan learn the subtle differences in the potential energy surfaces of the Na Cl + and\n8 8\nNa Cl + clusters, proving it correctly interprets the system’s physics, as described in the\n9 8\nsecond question.\n17/18\n3. Unique polarization test. We introduced a novel physicality test that other universal\nmodels typically cannot address: response to an external electric field. Our model\nsuccessfully reproduces the DFT potential energy surface under varying field strengths,\ndemonstrating a genuine understanding of polarization.\nIn summary, we have systematically performed the requested physical tests and\ncomparisons with other universal models. These new results provide robust and convincing\nvalidation for our main claims. We are grateful to the reviewer for guidance, which has\nsignificantly strengthened the manuscript.\nFor the polarization test, details are provided in SI,\n“…The results are presented in the Fig. S9, which plots the system energy as a function of\nfield strength. Our model demonstrates good agreement with the DFT reference, accurately\nreproducing the energy curve across the entire field range from −2 to 2 V/Å. Notably, the\nperformance of our foundation model is on par with that of FIREANN [19], a state-of-the-art\nmodel that was specifically trained on this water system. This highlights the transferability of\nour model, as it correctly captures complex polarization physics without explicit training for\nthis task. In contrast, the CHGNet and MACE-MP-0 models are entirely unresponsive to the\nfield, predicting a constant energy. This represents a fundamental failure to capture one of the\nmost essential physical interactions, underscoring the superior predictive power and physical\nrealism of our approach.\nFig. S9 Evaluation of polarizable interactions in the periodic water model. The plot\ncompares the performance of our model against DFT [19] (reference), FIREANN [19],\nCHGNet, and MACE-MP-0.”\n18/18\nResponse to the reviewers\nWe are sincerely grateful to the two reviewers for their valuable time and insightful feedback\nduring this round of review. We have carefully considered all the comments and have provided\ndetailed, point-by-point responses in this document. The modifications to the manuscript are\nhighlighted in dark red.\nTable of Contents\nResponse to reviewer #2 ............................................................................................................ 2\nResponse to reviewer #3 ............................................................................................................ 3\n1/5\nResponse to reviewer #2\n“Remarks to the Author:\nThis revision of the manuscript has brought further improvements and the authors have now\nclarified most pertinent points raised by the reviewers. The manuscript is publishable.”\nAuthors’ reply:\nWe would like to express our gratitude for your positive and encouraging feedback on our\nrevised manuscript. We are very pleased to hear that you found our revisions effective and that\nthe manuscript is now considered publishable. Your expert comments and guidance have been\ninvaluable in improving the quality of our work, and we appreciate the time and effort you\ndedicated to the review process.\n“Remarks on code availability:\nI haven’t run the code but the repository exists with examples and instructions.”\nAuthors’ reply:\nWe appreciate you taking the time to check our code repository and for confirming that the\nexamples and instructions are available.\n2/5\nResponse to reviewer #3\n“The manuscript has improved substantially compared to the previous version. The results now\nlook mostly solid and trustworthy, and many of my earlier concerns have been addressed. I\nonly have a few remaining comments, including one point that still puzzles me.”\nAuthors’ reply:\nWe sincerely thank you for your thoughtful and constructive feedback on our revised\nmanuscript. We are very encouraged to hear that you find it has improved substantially and\nthat our results are now mostly solid and trustworthy. We have carefully noted your remaining\ncomments and will address them as follows.\n“Major comment:\nPhysical decay of long-range interactions (Fig. S8). I recommend plotting this figure on a log–\nlog scale in order to assess whether the correct power-law decay of interactions is reproduced.\nFrom visual inspection, it seems that the extracted exponent p is too small; for Coulombic\ninteractions, one expects p = 1. Could the authors comment on this discrepancy, especially\ngiven the claim that their model reproduces correct long-range physics?”\nAuthors’ reply:\nWe thank the reviewer for this insightful comment and valuable suggestion. The reviewer\nis correct that a log-log plot is the most appropriate way to analyse the power-law decay of\nlong-range interactions. Following this recommendation, we have generated a new log-log plot,\nnow included in the new Fig. S9.The reviewer correctly observed that the interaction does not\nfollow a simple 1/r decay across the entire distance range. This is expected in the distance range\nwe plotted in Fig. S8, as contributions from higher-order electrostatic terms cannot be ignored.\nTo show that our method captures the correct asymptotic behaviour, we extended the distance\nto more than 16 Å as shown in Fig. S9. Indeed, the curve from our model is tangent to the line\nwith slope = −1, confirming that it successfully reproduces the correct p = 1 power-law decay\nin the long-range limit.\nIn Supplementary Information, we have included the log-log plot of our model in\nextended distance range,\n3/5\n“Fig. S9 A log-log plot showing the interaction energy as a function of distance for an OH-\nOH dimer in an ionic state. The data from our model (solid green line) is compared with a\nreference line representing a pure 1/r decay (dashed black line).”\n“Minor comments:\n1. Fine-tuning strategy (S11). The additional comparison is very helpful and the results look\nconvincing. Could the authors comment on their fine-tuning strategy in more detail? For\nexample, is this full fine-tuning, multi-head fine-tuning, or a parameter-efficient method such\nas LoRA? In this context, I am also curious about the performance of the fine-tuned models on\nthe original training tasks—did the authors observe catastrophic forgetting?”\nAuthors’ reply:\nWe sincerely thank the reviewer for their positive feedback on our supplementary results\nand for asking these insightful questions regarding our finetuning strategy and the issue of\ncatastrophic forgetting.\nRegarding our finetuning strategy, we adopted “full finetuning” approach. For our model\nand CHGNet, we unfroze and updated the parameters of the final two convolutional layers and\nthe subsequent Multilayer Perceptron (MLP) layers. For the MACE-MP model, we finetuned\nall model parameters, following the officially recommended practice for this architecture. Our\ndecision to employ these methods, rather than using multi-head finetuning or parameter-\nefficient techniques such as LoRA, was driven by the high-precision nature of the Ko et al.\ndataset. We hypothesized that a more flexible approach—achieved by unfreezing some\nconvolutional and MLP layers—was necessary to provide the models with sufficient capacity\nto capture these subtle but critical features of the potential energy surfaces with high fidelity.\nIt is important to clarify that our primary goal in finetuning the foundation model on these\nchallenging systems was to demonstrate its inherent capability to capture these complex\ninteractions. That said, the high-precision nature of the potential energy surfaces in Ko et al.\ndatasets meant that we did observe a degree of catastrophic forgetting. However, the extent of\nthis effect was dependent on the similarity between the finetuning system and the domains\ncovered in the original training data. For the Au -MgO system, the task involves modelling a\n2\nsurface system, which does not include isolated, charged clusters. This represents a relatively\nmodest domain shift from our original MPtrj dataset. Consequently, we found that our model’s\nperformance on the original test set is acceptable after finetuning (MAE for forces and stresses\nwere 0.87 eV/Å and 0.64 GPa on the MPtrj test set, respectively). Conversely, the Na Cl +\n8/9 8\nsystems represent a significant departure from our original MPtrj dataset, as they involve\nisolated, charged clusters with unique electrostatic characteristics. Finetuning on this domain,\nwhich is substantially different, led to more pronounced catastrophic forgetting. Based on this\nanalysis, we conclude that for highly distinct systems like these charged clusters, the most\nreliable path to achieving a highly accurate description is to train a bespoke model.\nIn the Supplementary Information, we included details about the finetuning strategy,\n“…Regarding our finetuning strategy, we adopted a full finetuning approach. For our model\nand CHGNet, only the final two convolutional layers, and the multilayer perceptron were\nupdated. For MACE-MP-0, the entire model was subjected to finetuning, as implemented in\nits own code…”\n4/5\n“2. Comparison table (S10). Thank you for adding this comparison. The results look strong and\nconvincing. It also appears that the Ko dataset is becoming a widely used benchmark. I would\nlike to point out another recent paper (https://arxiv.org/abs/2507.19382) which further raises\nthe benchmark results; it may be worth discussing this in the context of related work.”\nAuthors’ reply:\nWe thank the reviewer for their positive feedback and for bringing the recent and highly\nrelevant LOREM paper [https://arxiv.org/abs/2507.19382] to our attention. We agree that its\nmethodology which treats charges as equivariant tensors provides an important contrast to\nother long-range approaches.Since the LOREM framework also does not require a charge\nequilibration process, we have integrated a description of this preprint into our manuscript’s\nDiscussion section, highlighting it among other methods that bypass this step.\nIn the main text,\n“…LOREM [84] models charges as equivariant tensors to capture orientation-dependent\ninteractions beyond the standard cutoff radius…\n84. Rumiantsev, E., et al., Learning Long-Range Representations with Equivariant Messages.\narXiv preprint arXiv:2507.19382, 2025.”\n5/5\nResponse to reviewer #3\n“I thank the authors for their revision and improvements. I am happy with the manuscript and\nsupport publication in its current form.”\nAuthors’ reply:\nWe sincerely thank the reviewer for their positive feedback and for endorsing the publication\nof our manuscript in its current form. We are grateful for their time and evaluation."}, {"round": 3, "reviewer_comments": {"reviewer_2": "This revision of the manuscript has brought further improvements and the authors have now clarified most pertinent points raised by the reviewers. The manuscript is publishable.", "reviewer_3": "The manuscript has improved substantially compared to the previous version. The results now look mostly solid and trustworthy, and many of my earlier concerns have been addressed. I only have a few remaining comments, including one point that still puzzles me.\nMajor comment:\nPhysical decay of long-range interactions (Fig. S8). I recommend plotting this figure on a log–log scale in order to assess whether the correct power-law decay of interactions is reproduced. From visual inspection, it seems that the extracted exponent p is too small; for Coulombic interactions, one expects p = 1. Could the authors comment on this discrepancy, especially given the claim that their model reproduces correct long-range physics?\nMinor comments:\n1. Fine-tuning strategy (S11). The additional comparison is very helpful and the results look convincing. Could the authors comment on their fine-tuning strategy in more detail? For example, is this full fine-tuning, multi-head fine-tuning, or a parameter-efficient method such as LoRA? In this context, I am also curious about the performance of the fine-tuned models on the original training tasks—did the authors observe catastrophic forgetting?\n2. Comparison table (S10). Thank you for adding this comparison. The results look strong and convincing. It also appears that the Ko dataset is becoming a widely used benchmark. I would like to point out another recent paper (https://arxiv.org/abs/2507.19382) which further raises the benchmark results; it may be worth discussing this in the context of related work."}, "author_response": "Response to the reviewers\nWe are sincerely grateful to the two reviewers for their valuable time and insightful feedback\nduring this round of review. We have carefully considered all the comments and have provided\ndetailed, point-by-point responses in this document. The modifications to the manuscript are\nhighlighted in dark red.\nTable of Contents\nResponse to reviewer #2 ............................................................................................................ 2\nResponse to reviewer #3 ............................................................................................................ 3\n1/5\nResponse to reviewer #2\n“Remarks to the Author:\nThis revision of the manuscript has brought further improvements and the authors have now\nclarified most pertinent points raised by the reviewers. The manuscript is publishable.”\nAuthors’ reply:\nWe would like to express our gratitude for your positive and encouraging feedback on our\nrevised manuscript. We are very pleased to hear that you found our revisions effective and that\nthe manuscript is now considered publishable. Your expert comments and guidance have been\ninvaluable in improving the quality of our work, and we appreciate the time and effort you\ndedicated to the review process.\n“Remarks on code availability:\nI haven’t run the code but the repository exists with examples and instructions.”\nAuthors’ reply:\nWe appreciate you taking the time to check our code repository and for confirming that the\nexamples and instructions are available.\n2/5\nResponse to reviewer #3\n“The manuscript has improved substantially compared to the previous version. The results now\nlook mostly solid and trustworthy, and many of my earlier concerns have been addressed. I\nonly have a few remaining comments, including one point that still puzzles me.”\nAuthors’ reply:\nWe sincerely thank you for your thoughtful and constructive feedback on our revised\nmanuscript. We are very encouraged to hear that you find it has improved substantially and\nthat our results are now mostly solid and trustworthy. We have carefully noted your remaining\ncomments and will address them as follows.\n“Major comment:\nPhysical decay of long-range interactions (Fig. S8). I recommend plotting this figure on a log–\nlog scale in order to assess whether the correct power-law decay of interactions is reproduced.\nFrom visual inspection, it seems that the extracted exponent p is too small; for Coulombic\ninteractions, one expects p = 1. Could the authors comment on this discrepancy, especially\ngiven the claim that their model reproduces correct long-range physics?”\nAuthors’ reply:\nWe thank the reviewer for this insightful comment and valuable suggestion. The reviewer\nis correct that a log-log plot is the most appropriate way to analyse the power-law decay of\nlong-range interactions. Following this recommendation, we have generated a new log-log plot,\nnow included in the new Fig. S9.The reviewer correctly observed that the interaction does not\nfollow a simple 1/r decay across the entire distance range. This is expected in the distance range\nwe plotted in Fig. S8, as contributions from higher-order electrostatic terms cannot be ignored.\nTo show that our method captures the correct asymptotic behaviour, we extended the distance\nto more than 16 Å as shown in Fig. S9. Indeed, the curve from our model is tangent to the line\nwith slope = −1, confirming that it successfully reproduces the correct p = 1 power-law decay\nin the long-range limit.\nIn Supplementary Information, we have included the log-log plot of our model in\nextended distance range,\n3/5\n“Fig. S9 A log-log plot showing the interaction energy as a function of distance for an OH-\nOH dimer in an ionic state. The data from our model (solid green line) is compared with a\nreference line representing a pure 1/r decay (dashed black line).”\n“Minor comments:\n1. Fine-tuning strategy (S11). The additional comparison is very helpful and the results look\nconvincing. Could the authors comment on their fine-tuning strategy in more detail? For\nexample, is this full fine-tuning, multi-head fine-tuning, or a parameter-efficient method such\nas LoRA? In this context, I am also curious about the performance of the fine-tuned models on\nthe original training tasks—did the authors observe catastrophic forgetting?”\nAuthors’ reply:\nWe sincerely thank the reviewer for their positive feedback on our supplementary results\nand for asking these insightful questions regarding our finetuning strategy and the issue of\ncatastrophic forgetting.\nRegarding our finetuning strategy, we adopted “full finetuning” approach. For our model\nand CHGNet, we unfroze and updated the parameters of the final two convolutional layers and\nthe subsequent Multilayer Perceptron (MLP) layers. For the MACE-MP model, we finetuned\nall model parameters, following the officially recommended practice for this architecture. Our\ndecision to employ these methods, rather than using multi-head finetuning or parameter-\nefficient techniques such as LoRA, was driven by the high-precision nature of the Ko et al.\ndataset. We hypothesized that a more flexible approach—achieved by unfreezing some\nconvolutional and MLP layers—was necessary to provide the models with sufficient capacity\nto capture these subtle but critical features of the potential energy surfaces with high fidelity.\nIt is important to clarify that our primary goal in finetuning the foundation model on these\nchallenging systems was to demonstrate its inherent capability to capture these complex\ninteractions. That said, the high-precision nature of the potential energy surfaces in Ko et al.\ndatasets meant that we did observe a degree of catastrophic forgetting. However, the extent of\nthis effect was dependent on the similarity between the finetuning system and the domains\ncovered in the original training data. For the Au -MgO system, the task involves modelling a\n2\nsurface system, which does not include isolated, charged clusters. This represents a relatively\nmodest domain shift from our original MPtrj dataset. Consequently, we found that our model’s\nperformance on the original test set is acceptable after finetuning (MAE for forces and stresses\nwere 0.87 eV/Å and 0.64 GPa on the MPtrj test set, respectively). Conversely, the Na Cl +\n8/9 8\nsystems represent a significant departure from our original MPtrj dataset, as they involve\nisolated, charged clusters with unique electrostatic characteristics. Finetuning on this domain,\nwhich is substantially different, led to more pronounced catastrophic forgetting. Based on this\nanalysis, we conclude that for highly distinct systems like these charged clusters, the most\nreliable path to achieving a highly accurate description is to train a bespoke model.\nIn the Supplementary Information, we included details about the finetuning strategy,\n“…Regarding our finetuning strategy, we adopted a full finetuning approach. For our model\nand CHGNet, only the final two convolutional layers, and the multilayer perceptron were\nupdated. For MACE-MP-0, the entire model was subjected to finetuning, as implemented in\nits own code…”\n4/5\n“2. Comparison table (S10). Thank you for adding this comparison. The results look strong and\nconvincing. It also appears that the Ko dataset is becoming a widely used benchmark. I would\nlike to point out another recent paper (https://arxiv.org/abs/2507.19382) which further raises\nthe benchmark results; it may be worth discussing this in the context of related work.”\nAuthors’ reply:\nWe thank the reviewer for their positive feedback and for bringing the recent and highly\nrelevant LOREM paper [https://arxiv.org/abs/2507.19382] to our attention. We agree that its\nmethodology which treats charges as equivariant tensors provides an important contrast to\nother long-range approaches.Since the LOREM framework also does not require a charge\nequilibration process, we have integrated a description of this preprint into our manuscript’s\nDiscussion section, highlighting it among other methods that bypass this step.\nIn the main text,\n“…LOREM [84] models charges as equivariant tensors to capture orientation-dependent\ninteractions beyond the standard cutoff radius…\n84. Rumiantsev, E., et al., Learning Long-Range Representations with Equivariant Messages.\narXiv preprint arXiv:2507.19382, 2025.”\n5/5"}, {"round": 4, "reviewer_comments": {"reviewer_3": "I thank the authors for their revision and improvements. I am happy with the manuscript and support publication in its current form."}, "author_response": ""}]

Decision letter