qNEP: A highly efficient neuroevolution potential with dynamic charges for large-scale atomistic simulations

TL;DR

qNEP addresses the challenge of incorporating long-range electrostatics in machine-learned interatomic potentials by learning environment-dependent partial charges and integrating them with a short-range NEP core. It provides two evaluation modes (mode 1: real+reciprocal; mode 2: reciprocal only) and consistent force/virial expressions including dynamic-charge contributions, with a PPPM-based acceleration to achieve near-linear scaling. The method is implemented in GPUMD and validated on water, LLZO, BaTiO3, and Mg–water interfaces, delivering accurate energies, forces, stresses, dielectric properties, and polarization dynamics, at million-atom scales on consumer GPUs. The work enables direct computation of dielectric response, infrared spectra, and field-matter coupling in large-scale simulations and suggests future extensions to dispersion or other long-range interactions.

Abstract

Although electrostatics can be incorporated into machine-learned interatomic potentials, existing approaches are computationally very demanding, limiting large-scale, long-time simulations of electrostatics-driven phenomena such as dielectric response, infrared activity, and field-matter coupling. Here, we extend the neuroevolution potential (NEP), a highly efficient machine-learned interatomic potential, to a charge-aware framework (qNEP) by introducing explicit, environment-dependent partial charges. Each ionic partial charge is represented by a neural network as a function of the local descriptor vector, analogous to the NEP site-energy model. This formulation enables the direct prediction of the Born effective charge tensor for each ion and, consequently, the polarization. As a result, dielectric properties, infrared spectra, and coupling to external electric fields can be evaluated within a unified framework. We derive consistent expressions for the forces and virials that explicitly account for the position dependence of the partial charges. The qNEP method has been implemented in the free-and-open-source GPUMD package, with support for both Ewald summation and particle-particle particle-mesh treatments of electrostatics. We demonstrate the accuracy and efficiency of the qNEP approach through representative applications to water, Li7La3Zr2O12, BaTiO3, and a magnesium-water interface. These results show that qNEP enables accurate atomistic simulations with explicit long-range electrostatics, scalable to million-atom systems on nanosecond time scales using consumer-grade GPUs.

Figures (7)

• Figure 1: Schematic illustration of the qNEP framework. (a) Neural network architecture of a qNEP model with two outputs: a site energy $U_i$ and a partial charge $q_i$. The site energies are summed to yield the NEP contribution $U^\text{NEP}$ (\ref{['sect:nep']}) to the total energy $U^\text{tot}$, while interactions between the partial charges give rise to the electrostatic energy $U^\text{ES}$ (\ref{['sect:qnep']}). (b) Evaluation of the electrostatic energy and corresponding forces $\mathbf{F}^\text{ES}_i$ (\ref{['sect:energy-derivatives']}), including either both real- and reciprocal-space contributions (mode 1) or the reciprocal-space contribution only (mode 2). (c) Derived response properties obtained from the partial charges, such as the polarization $\mathbf{P}$, the Born effective charges $\mathbf{\mathcal{Z}}_i$, and forces induced by external electric fields (\ref{['sect:born-effective-charges']}). (d) Computational cost for evaluating the electrostatic contribution, comparing direct Ewald summation with the method (\ref{['sect:ewald-and-pppm']}), which offers superior computational performance, particularly for large systems (shown here for water; \ref{['sect:water']}).
• Figure 2: Performance of qNEP for water. (a) Relative validation of NEP and qNEP models. (b,c) Parity plots of Born effective charges for H and O, respectively, obtained from the qNEP model trained on the full reference data set using mode 2. (d,e) Performance comparison of qNEP models with , a CACE model with , and a deep potential model with long-range electrostatic interactions (DPLR), in terms of (d) computational speed and (e) the simulated time achievable within one day on a single Nvidia RTX4090 GPU using a time step of 0.5.
• Figure 3: Properties of water with NEP and qNEP models. (a,b) Partial radial distribution functions for (a) O--H and (b) H--H pairs in liquid water at 300. Solid and dashed lines correspond to quantum simulations performed using and classical simulations using standard , respectively, for the NEP (blue) and qNEP (red) models. Classical ab-initio (AIMD) reference data (dotted line) from Ref. XuHaoLia23 are included for comparison. (c,d) Infrared spectra obtained from classical simulations via the time of the ionic electric current $\dot{\mathbf{P}}$. (c) Spectra obtained using different combinations of the NEP and qNEP (mode 2) models for sampling the energy landscape ($E$) and the TNEP (from Ref. XuRosSch24) and qNEP models for computing the polarization ($P$) in comparison with experiment. (d) Temperature dependence of infrared spectra obtained using the qNEP model (mode 2) for both $E$ and $P$.
• Figure 4: Garnet-type llzo. (a) Relative training of and qNEP models, with values reported above the corresponding columns. (b,c) Temperature dependence of (b) lattice parameters and (c) heat capacity obtained from heating and cooling simulations. Vertical dashed lines indicate the corresponding phase transition temperatures. Experimental data for the lattice parameters from Ref. chen2015study. (d) Distribution of Li partial charges in the tetragonal (t-) and cubic (c-) phases after structural relaxation. The left-hand peaks in the distributions correspond to Li-ions occupying tetrahedral sites (see inset) while the right-hand peaks correspond to octahedral sites (see inset). (e) Arrhenius plot of the ionic conductivity $\sigma$ multiplied by temperature $T$. Activation energies extracted for the tetragonal and cubic phases are indicated. The vertical dashed line marks the average transition temperature of 900 obtained from heating and cooling runs.
• Figure 5: Ferroelectricity and dielectric response in barium titanate (BaTiO3). (a--c) Temperature dependence of (a) the lattice constant, (b) the polarization, and (c) the dielectric constant during heating and cooling, revealing the sequence of phase transitions from the rhombohedral (R) ground state through the orthorhombic (O) and tetragonal (T) phases to the high-temperature cubic (C) phase. Vertical dashed lines indicate the experimentally observed transition temperatures Merz1949. Experimental data points in (c) from Ref. Benedict1958. (d) Polarization--electric field ($P$--$E$) hysteresis loops and (e,f) the real and imaginary parts of the dielectric function at different temperatures, corresponding to all four phases (indicated by arrows in (b)). The gray diamond in (d) marks the experimental value for the spontaneous polarization at room temperature Wieder1955.