🔬 physics

samQEq: universal charge equilibration for molecules and conductors

This paper introduces samQEq, a unified LAMMPS package that enables simultaneous charge equilibration for both molecules and conductors by combining distinct neutrality constraints and self-energy corrections into a single, efficient, and stable matrix-free solver.

Original authors: Nicholas Solan, David Sun, Tod A. Pascal

Published 2026-11-04

📖 6 min read🧠 Deep dive

CC BY 4.0

Original authors: Nicholas Solan, David Sun, Tod A. Pascal

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

In the microscopic world of materials science, atoms are rarely static. They are constantly shifting, reacting to their neighbors, and reshaping their electrical personalities. When a molecule sits in a liquid, or when a metal surface meets water, the electrons within those atoms do not stay put; they flow to balance the electrical forces around them. This constant redistribution is what gives materials their ability to conduct electricity, dissolve salts, or catalyze chemical reactions. For decades, scientists have used computer simulations to watch these interactions unfold, but they faced a stubborn dilemma. To simulate a single molecule, they needed a method that kept the total charge of that molecule fixed, preventing electrons from leaking out into the void. To simulate a metal, they needed the opposite: a method that allowed charge to flow freely across the entire surface, just as it does in a real conductor. Until now, no single computer tool could handle both situations at once. Researchers had to choose one approach for molecules and a completely different one for metals, making it nearly impossible to study the complex interfaces where the two meet, such as a battery electrode soaking in an electrolyte.

A team of researchers at the University of California San Diego has built a new digital framework called samQEq that finally bridges this gap. Their work, published as a preprint, introduces a unified system capable of simulating both insulating molecules and conducting metals within the same calculation. The core achievement is a mathematical engine that treats the flow of electrical charge as a single, continuous problem, regardless of whether the atoms involved belong to a floating water molecule or a solid gold slab. In their tests, the team demonstrated that this system could keep a gold slab electrically stable while simultaneously simulating a film of water sitting on top of it. Crucially, the simulation allowed the water molecules to maintain their own electrical neutrality while the gold slab acted as a perfect conductor, screening electric fields just as a real metal would. This was achieved without forcing the water to behave like a metal or the gold to behave like an insulator; instead, the software applied the correct physical rules to each part of the system simultaneously.

The researchers verified their new tool by recreating known physical behaviors with high precision. They successfully reproduced the electrical properties of a specific water model that had been published previously, matching the results without needing to adjust any of the underlying numbers. They also showed that their system could mimic the charge calculations used in a popular reactive force field called ReaxFF, matching its results to within one part in a billion. These tests confirmed that the new framework is not just a theoretical idea but a practical tool that can be trusted to produce accurate data. The team further proved that the system remains stable even when the metal is pushed to its limits, simulating a scenario where the metal's ability to resist charge flow drops to near zero. In these extreme conditions, the system continued to function correctly, preventing the mathematical breakdowns that usually crash such simulations.

To make this work fast enough for real-world scientific discovery, the team optimized the code to run on modern graphics processing units, the powerful chips often used for video games and artificial intelligence. On a single computer node equipped with four of these graphics cards, the new system processed simulations roughly twice as fast as a standard computer node with 128 processor cores. This speedup is significant because it allows scientists to run larger simulations or observe longer periods of time, which is essential for understanding how materials behave over the long term. The researchers also built in a series of automatic checks that run before every simulation starts. These checks ensure that the mathematical rules being applied match the physical setup, catching potential errors that would otherwise go unnoticed and lead to incorrect results.

The implications of this work extend to the design of future energy technologies. By allowing scientists to simulate the exact interface between a metal electrode and a liquid electrolyte in one go, the tool opens the door to studying batteries and fuel cells with a level of detail that was previously out of reach. The researchers showed that their method could handle a system containing thousands of atoms, including copper electrodes and water molecules, while maintaining electrical stability over thousands of steps. They noted, however, that while the tool works exceptionally well on a single machine, scaling it up to run across many machines simultaneously remains a challenge due to the way the calculations communicate with each other. For now, the tool is best suited for running on a single, powerful computer node.

The success of samQEq lies in its ability to respect the distinct nature of different materials while solving them together. In the past, trying to simulate a metal and a molecule together often resulted in a chaotic mess where charge would flow uncontrollably between them, or the calculation would simply fail. The new approach uses a specific mathematical trick to keep the metal's part of the calculation stable, even when the metal is acting like a perfect conductor. This stability allows the simulation to capture the subtle dance of electrons at the boundary between a solid and a liquid. The researchers found that when they applied this method to a gold slab with a water film, the metal effectively screened a point charge placed above it, reducing the electrical influence by a constant amount at every height beyond a certain distance. This behavior matched the theoretical expectations for a real metal, confirming that the simulation was capturing the correct physics.

This work represents a significant step forward in computational materials science. It provides a single, consistent layer of physics that can be applied to any type of material, from simple gases to complex solid-liquid interfaces. By removing the need to switch between different simulation methods, it allows researchers to focus on the chemistry and physics of the system rather than the limitations of their tools. The team has made their code available to the scientific community, ensuring that others can build upon this foundation. While the current version has some limitations, such as not yet supporting constant-potential simulations or certain complex machine-learning potentials, it establishes a robust framework for the next generation of materials modeling. The ability to simulate the interface between a conductor and an insulator with such fidelity promises to accelerate the discovery of new materials for energy storage, catalysis, and electronics.

Technical Summary of samQEq: Universal Charge Equilibration for Molecules and Conductors

Problem Statement
Molecular simulations require accurate treatment of electrostatics, where atomic charges respond to their environment. Traditional fixed-charge models fail to capture polarization, while charge equilibration (QEq) methods allow charges to vary but face a fundamental dichotomy:

Molecular/Insulating Limit: Fluctuating-charge models typically enforce per-molecule neutrality to prevent unphysical charge transfer between molecules. However, standard global electronegativity equalization allows unlimited charge flow, behaving like a metal.
Conductor/Metallic Limit: Metals require a global neutrality constraint and near-zero hardness. However, as hardness approaches zero, the charge operator becomes indefinite (losing positive definiteness) unless a self-energy term is included. Without this, the solver fails to converge or returns a saddle point.
Integration Gap: Existing LAMMPS packages handle these limits separately. No single framework currently solves for molecules (with per-molecule constraints) and conductors (with global constraints and self-energy stabilization) simultaneously within one energy functional and one linear solve, particularly when long-range electrostatics are involved.

Methodology
The authors present samQEq (solids-and-molecules QEq), a unified charge-equilibration package for LAMMPS that resolves these limits within a single matrix-free solver.

Unified Energy Functional: The method minimizes a second-order energy functional based on electronegativity ( $\chi$ ), hardness ( $\eta$ ), and shielded Coulomb interactions ( $J_{ij}$ ). The functional is referenced to formal charges ( $q^0_i$ ) to handle non-neutral species.
Dual Constraints via Projector: An orthogonal projector is applied to the search direction in the Krylov solver. This allows the user to select run-time constraints:
- Per-fragment neutrality: Enforces charge conservation within each molecule (insulating limit), preventing inter-molecular charge transfer.
- Global neutrality: Enforces a single constraint for the entire system (conductor limit).
Conductor Stabilization (Self-Energy): To ensure the operator remains positive definite at near-zero hardness for conductors, samQEq adds a finite-width self-energy term to the on-site hardness for designated conductor types. This term is derived from a Gaussian charge distribution with a width ( $w$ ) narrower than the conductor's own Gaussian shielding width. This creates a rigorous lower bound on the smallest eigenvalue, preventing indefiniteness.
Matrix-Free Solver: The package employs a matrix-free approach where the Krylov solver (default: Projected Conjugate Gradient) receives only the action of the operator on a vector.
- Reciprocal-Space Integration: The reciprocal-space contribution (Particle-Mesh Ewald) is evaluated inside the matrix-vector product. This avoids an outer self-consistent field loop, allowing a single projected solve to yield fully long-range charges.
- Kernel Flexibility: Three shielding kernels are supported as run-time settings: Cube-root, Gaussian, and exact Slater-orbital integrals. The exact Slater integral allows for the reproduction of published models without refitting parameters.
Code Architecture: Implemented as a LAMMPS fix with a companion pair style. It supports fixed-charge, reactive (ReaxFF), and machine-learned potentials via hybrid/overlay. Device kernels (Kokkos) are available for GPU acceleration.

Key Contributions

Simultaneous Solvability: Demonstrates that a single linear solve can handle both insulating molecules (with per-molecule constraints) and metallic conductors (with global constraints and self-energy) in the same system without code path changes.
Mathematical Well-Posedness: Provides a rigorous sufficient condition (Eq. 9) for the conductor limit to remain positive definite at any hardness, provided the self-energy width is narrower than the shielding width.
Implementation of Exact Kernels: Enables the use of exact Slater-orbital integrals for shielding, facilitating the direct porting of published fluctuating-charge models (e.g., TIP4P-FQ) without parameter refitting.
Performance: Achieves approximately double the throughput of a 128-core CPU node when running on a single node with four NVIDIA V100 GPUs, utilizing GPU-aware MPI and device kernels for the matrix-vector product.

Results and Validation

Gold Slab with Water Film: On a system comprising a gold slab and a water film, the authors show that only the combination of per-molecule constraints (for water) and self-energy (for gold) keeps the operator positive definite down to the conductor limit ( $\eta \to 0$ ). Without self-energy, the operator becomes indefinite, leading to unbounded charge fluctuations and solver failure. The metal successfully screens a point charge by a constant amount beyond 4 Å.
Reproduction of Published Models:
- TIP4P-FQ Water: Using the exact Slater kernel and original parameters, samQEq reproduces the gas-phase dipole (1.8634 D vs. 1.85 D) and liquid dipole (2.630 D vs. 2.637 D) of the Rick et al. model without refitting.
- ReaxFF Charge Step: The package reproduces the native ReaxFF charge equilibration step to within $10^{-9} e$ for water and RDX systems, using ReaxFF's native parameters.
Numerical Stability: The solver maintains stability under dynamics. In a 5 ps simulation of water on a gold slab, the self-energy term prevented charge drift and kept site charges bounded, whereas the system without self-energy exhibited uncontrolled charge growth and instability.
Initialization Checks: The code includes robust checks to detect mismatches between the charge solver and the force field (e.g., kernel mismatches, cutoff inconsistencies) that would otherwise lead to silent errors.

Significance and Claims
The paper claims that samQEq provides a "universal charge-equilibration layer" that bridges the gap between molecular and metallic electrostatics. Its primary significance lies in:

Unification: It treats molecules, conductors, and their interfaces as a single electrostatic problem solvable within one framework.
Robustness: By embedding the reciprocal-space term in the matrix-vector product and using a projector, it ensures the solver is well-posed for both limits without requiring separate code paths or outer loops.
Extensibility: It serves as a foundation for future developments, such as constant-potential electrodes, redox reactions, and coupling with machine-learned potentials that predict electronegativities, without requiring derivatives of the charge solution.

The authors note current limitations, including the lack of constant-potential modes in this release, the requirement for manual parameter derivation for new electrode materials, and the fact that the device implementation currently scales well only within a single node due to latency-bound collectives. However, the framework establishes a physically consistent and computationally efficient method for simulating complex electrochemical interfaces.

More like this