Source-linked AI summary
qNEP: A highly efficient neuroevolution potential with dynamic charges for large-scale atomistic simulations
Zheyong Fan, Benrui Tang, Esmée Berger, Ethan Berger, Erik Fransson, Ke Xu, Zihan Yan, Zhoulin Liu, Zichen Song, Haikuan Dong, Shunda Chen, Lei Li, Ziliang Wang, Yizhou Zhu, Julia Wiktor, Paul Erhart
TL;DR
Existing electrostatic MLIPs are too computationally demanding for many large-scale, long-time simulations of polarization-driven phenomena. qNEP extends NEP with learned partial charges and consistent electrostatic derivatives, and achieves accurate applications with only 1.5 to 3 times the cost of comparable NEP models.
Problem
Existing electrostatic MLIPs constrain large-scale and long-time simulations of phenomena involving long-range electrostatics, polarization, dielectric response, spectroscopy, and field-driven dynamics.
Method
qNEP extends NEP with latent partial charges, charge-conservation regularization, polarization-derived BECs, and explicit electrostatics with Ewald or PPPM evaluation.
Results
qNEP is only 1.5 to 3 times slower than equivalently trained NEP models and supports accurate applications across water, ionic, ferroelectric, and reactive-interface systems.
Takeaways & Limitations
qNEP provides a practical route to simulations of charge-transfer- and polarization-driven phenomena across extended length and time scales.
Abstract
from arXiv · showhide
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.
I. INTRODUCTION
qNEP addresses the limitations of short-ranged and computationally demanding electrostatic MLIPs by extending NEP with charge-aware modeling for scalable simulations and charge-related observables.
- Short-ranged MLIPs become inadequate for sizable partial charges, weak screening, and explicit coupling to external electric fields.
- Existing electrostatic MLIPs remain computationally demanding for systems with hundreds of thousands to millions of atoms and nanosecond time scales.
- qNEP extends NEP with latent partial charges, charge-conservation regularization, and BECs derived from polarization.
- PPPM makes qNEP only 1.5 to 3 times slower than equivalently trained NEP models while supporting explicit long-range electrostatics and million-atom nanosecond simulations.
- A. The original NEP model architecture: NEP uses a feedforward neural network whose trainable parameters map local descriptors to site energies.
- A. The original NEP model architecture: The qNEP framework combines site-energy and partial-charge outputs, with electrostatic interactions contributing to total energy and derived response properties.
B. The qNEP model architecture
The qNEP architecture adds a partial-charge output to NEP and combines learned short-range energy with explicit electrostatic contributions, while supporting two electrostatic evaluation modes and response-property calculations.
- qNEP uses one neural network with an additional output node for environment-dependent partial charges rather than separate energy and charge networks.
- Periodic electrostatics are evaluated with Ewald decomposition into real-space, reciprocal-space, and self-energy terms.
- The total energy combines the NEP contribution with electrostatic energy computed from the predicted partial charges.
- Mode 1 includes real- and reciprocal-space electrostatics, whereas mode 2 includes only the reciprocal-space contribution because NEP may already capture short-range interactions.
- The framework derives forces, virials, BEC tensors, dielectric functions, infrared spectra, ionic conductivity, and external-field coupling from the electrostatic model.
- Training augments the NEP loss with charge-conservation penalties and optional BEC constraints, while implementation supports Ewald summation and PPPM in GPUMD.
1. Real-space electrostatic energy
The real-space electrostatic term is a cutoff-limited contribution of the Ewald formulation, with its convergence controlled by the inverse-length parameter α and coordinated with the NEP cutoff.
- In mode 1, the real-space electrostatic energy is evaluated using pairwise charges, interatomic distances, the complementary error function, and vacuum permittivity.
- The real-space contribution is truncated at cutoff radius r_c, chosen to match the associated NEP pairwise cutoff.
- α controls the relative convergence rates of real- and reciprocal-space contributions in the Ewald decomposition.
- The reciprocal-space contribution U_k is required in both electrostatic modes and is expressed through the Green’s function and structure factor.
- The reciprocal-space sum excludes k = 0, truncates at k_max, and achieves approximately 10^-5 accuracy with k_max = 2πα.
- Mode 1 includes a self-energy term that removes each charge’s unphysical interaction with its own screening cloud.
C. Energy derivatives
qNEP energy derivatives must include both ordinary Coulomb position dependence and additional contributions from configuration-dependent charges, yielding consistent forces and virials for the many-body model.
- Energy derivatives include static charge contributions from explicit 1/r Coulomb dependence and dynamic charge contributions from position-dependent charges.
- Because qNEP is a many-body potential, forces can be represented as pairwise sums of partial-force contributions while retaining weak Newton’s third law.
- The per-atom virial tensor is constructed from partial forces and the distance vector between atoms.
- Configuration-dependent charges add chain-rule force contributions to the real-space, reciprocal-space, and self-energy parts of the Ewald sum.
- The electrostatic contribution is treated componentwise alongside the NEP partial force when deriving forces and virials.
1. The real-space contribution
The real-space electrostatic force contains a pairwise contribution, while configuration-dependent charges add position-dependent terms obtained through charge derivatives. Reciprocal-space forces and virials likewise distinguish static-charge and dynamic-charge contributions.
- The real-space contribution: Static-charge real-space electrostatics is pairwise, but dynamic charges add a partial-force contribution through their position dependence.The charge derivative is evaluated using a chain rule with respect to relative position vectors.
- The reciprocal-space contribution: Reciprocal-space forces include separate static-charge and dynamic-charge contributions.The dynamic-charge term accounts for configuration-dependent charge variation.
- The reciprocal-space contribution: The reciprocal-space virial can be decomposed per atom using a k-space stress kernel that maps each k-mode contribution onto a second-rank tensor.The kernel is derived from the reciprocal-space energy's response to homogeneous cell strain.
- The reciprocal-space contribution: The stress tensor has six independent components because the reciprocal-space contribution derives from an underlying pairwise electrostatic interaction.The total virial can be evaluated more cheaply when a per-atom virial is unnecessary.
3. The self-energy contribution
The self-energy is force-neutral for fixed charges but contributes additional forces when charges depend on atomic positions. The predicted charges also provide polarization, Born effective charges, field-induced forces, currents, infrared spectra, and conductivity.
- The self-energy contribution: Configuration-dependent charges make the self-energy position-dependent, producing an additional force contribution.This contribution is evaluated using the charge-position derivative and the chain rule.
- Born effective charge and related properties: The output charges are scaled using the high-frequency relative permittivity before calculating polarization and Born effective charges.The permittivity accounts for electronic screening absent from the ionic degrees of freedom.
- Born effective charge and related properties: For non-periodic systems, polarization reduces to the dipole moment, whereas periodic systems require a translationally invariant formulation.Born effective charge tensors are obtained from derivatives of polarization.
- Born effective charge and related properties: Born effective charges combine with external electric fields to determine ionic forces and with atomic velocities to evaluate the time derivative of polarization.The polarization trajectory follows by integrating this derivative when its initial value is known.
E. Training of the models
qNEP trains shared potential-energy and charge outputs using a weighted loss over energies, forces, virials, Born effective charges, total charges, and regularization. Charge penalties encourage conservation, while a final correction enforces the target total charge.
- E. Training of the models: All descriptor, neural-network, and dielectric-permittivity parameters are optimized with the SNES method.The loss function collects the trainable parameters in an abstract vector.
- E. Training of the models: The loss combines RMSE terms for energies, forces, virials, Born effective charges, and total charges with L1 and L2 regularization.The total-charge term concerns each structure's total charge rather than individual partial charges.
- E. Training of the models: A final total-charge correction enforces strict charge conservation or neutrality before electrostatic energy and Born effective charges are evaluated.The correction is especially relevant to simulations with external electric fields.
- E. Training of the models: Reference Born effective charges are optional and can be supplied for only a subset of training structures.This limits the number of computationally demanding reference BEC calculations.
F. Accelerated calculation of the reciprocal-space contribution using PPPM
PPPM accelerates reciprocal-space electrostatics by replacing direct summation with FFT-based mesh calculations. Its near-linear practical scaling substantially lowers cost for large systems while retaining force and virial evaluation capabilities.
- F. Accelerated calculation of the reciprocal-space contribution using PPPM: PPPM uses FFT-based particle–mesh calculations to reduce reciprocal-space electrostatic cost relative to direct Ewald summation.Particle–mesh Ewald variants are mathematically related approaches for this acceleration.
- F. Accelerated calculation of the reciprocal-space contribution using PPPM: Charges are interpolated onto a regular mesh using an assignment function, and the structure factor is computed on that mesh.The mesh has dimensions Nx×Ny×Nz, with interpolation orders specified for the assignment functions.
- F. Accelerated calculation of the reciprocal-space contribution using PPPM: A mesh spacing below 1 Å and interpolation order P = 5 yield approximately 10^-4 accuracy in the implementation.The optimized Green's function is chosen consistently with the mesh interpolation scheme.
- F. Accelerated calculation of the reciprocal-space contribution using PPPM: Static-charge forces use ik differentiation, whereas dynamic-charge forces use analytical differentiation.Backward FFTs provide the Cartesian static-force components and the dynamic-charge force contribution.
- F. Accelerated calculation of the reciprocal-space contribution using PPPM: PPPM scales as O(N log N) formally and is near-linear in the studied systems, making its cost one to several orders of magnitude below Ewald summation.The method has a small prefactor over the system sizes considered.
- F. Accelerated calculation of the reciprocal-space contribution using PPPM: Water performance comparisons cover validation errors, Born effective-charge parity, computational speed, and one-day simulated time for qNEP and competing models.The one-day comparison uses a 0.5 fs time step on a single Nvidia RTX4090 GPU.
III. RESULTS
Across water, LLZO, BaTiO3, and a magnesium–water interface, qNEP improves accuracy while enabling electrostatic, polarization, and charge-dependent analyses at extended scales. The applications cover liquids, ionic crystals, ferroelectrics, and reactive interfaces.
- A. Liquid water: qNEP systematically improves water accuracy over regular NEP at modest additional computational cost.The models reproduce structural distributions and infrared spectra while incorporating long-range electrostatics.
- B. Lithium lanthanum zirconate crystal: LLZO simulations reproduce experimental thermal expansion and the tetragonal-to-cubic transition near 900 K.The model also resolves phase-dependent Li charge distributions and their relation to ionic transport.
- C. Barium titanate: qNEP reproduces BaTiO3 phase transitions, structural changes, polarization evolution, and temperature-dependent dielectric behavior.The model also supports polarization–electric-field hysteresis loops for coupling to external fields.
- D. Magnesium–water interface: At a magnesium–water interface, qNEP captures environment-dependent charge states and conversion of metallic Mg into hydroxylated and solvated species.Its efficiency enables simulations of highly reactive conditions over many nanoseconds.
A. Liquid water
The applications show that qNEP improves accuracy and supports electrostatic observables across water and LLZO. It reproduces water structure and infrared spectra, while LLZO simulations capture phase behavior, charge distributions, and transport changes.
- A. Liquid water: 1388 training and 500 validation structures of 384 atoms support the water models’ DFT-based training.The dataset provides energies, forces, and stresses for liquid-water configurations.
- A. Liquid water: qNEP systematically reduces energy, force, and stress errors relative to NEP for water.The reciprocal-space-only mode performs marginally better for energies and forces, while both qNEP variants reproduce BECs.
- A. Liquid water: √ϵ∞ = n = 1.77 and 1.53 for modes 1 and 2, respectively, compared with the experimental value 1.33 at ambient conditions.These fitted values are in reasonable agreement with experiment despite primarily serving as a training hyperparameter.
- A. Liquid water: 10^7 atom step/s and up to 40 ns of MD simulation per day are achieved on one RTX 4090, with qNEP and PPPM increasing cost by about twofold.The reported throughput applies to systems containing at least 10^4 atoms.
- A. Liquid water: Water radial distribution functions from NEP and qNEP are essentially indistinguishable and agree well with ab initio MD.This holds for both classical and path-integral molecular-dynamics simulations.
- A. Liquid water: Room-temperature infrared spectra compare well with experiment, while increasing temperature blueshifts O–H stretching and redshifts libration.The spectra are calculated from the time autocorrelation function of polarization or its time derivative.
- B. Lithium lanthanum zirconate crystal: qNEP reduces LLZO energy, force, and stress RMSEs by approximately 20 % to 30 % relative to NEP.The models are evaluated for an ionic crystal where electrostatic interactions influence lithium transport.
- B. Lithium lanthanum zirconate crystal: LLZO lattice parameters agree well with experiment, and qNEP captures the tetragonal-to-cubic transition at approximately 900 K.Heating and cooling simulations also show a hysteresis of about 55 K.
C. Barium titanate
For BaTiO3, qNEP reproduces the ferroelectric phase behavior and polarization while enabling field-dependent and dielectric observables. Its predicted Born effective charges also recover the long-range LO–TO phonon splitting that NEP misses at Γ.
- Phase transitions: 151 K, 235 K, and 390 K are the predicted transition temperatures for the four BaTiO3 phases, in good agreement with experiment.A hysteresis of up to 50 K occurs between heating and cooling runs.
- Polarization: qNEP resolves all four phases through their polarization, which vanishes above 407 K during heating and 390 K during cooling.Polarization increases from the rhombohedral through orthorhombic to tetragonal phases.
- External fields: qNEP produces switchable polarization–electric-field hysteresis loops and reproduces the room-temperature spontaneous polarization.The 500 MHz switching frequency is much higher than experimentally accessible frequencies, so coercive fields are not quantitatively comparable.
- Dielectric response: Approximately 3000 is the maximum dielectric constant near the high-temperature side of the tetragonal–cubic phase boundary, matching experiment in magnitude and temperature dependence.The dielectric response depends strongly on temperature and frequency, with resonances from 20 meV to 40 meV attributed to longitudinal optical modes.
- Phonons: qNEP preserves LO–TO separation at Γ in both harmonic and finite-temperature phonon dispersions, whereas NEP predicts coincident modes.The splitting arises from long-range Coulomb interactions and requires non-analytic corrections based on Born effective charges.
IV. SUMMARY AND CONCLUSIONS
qNEP extends NEP with latent environment-dependent charges, charge-conserving corrections, and explicit long-range electrostatics. Across diverse systems, it improves accuracy while providing charge- and polarization-related observables at practical computational cost, though the present work treats electrostatics as its only long-range interaction.
- qNEP learns partial charges as latent features, enforces charge conservation through regularization and total-charge correction, and derives polarization and Born effective charges consistently.
- Across liquid, ionic, ferroelectric, and reactive-interface systems, qNEP improves energies, forces, and stresses relative to NEP while exposing charge- and field-related observables.
- qNEP enables simulations of transport, dielectric response, spectroscopy, and electrochemical reactivity across extended length and time scales.
- The present framework includes electrostatic interactions as the long-range contribution but does not yet incorporate dispersion forces.Adding other long-range interactions is identified as future work.