The question that didn’t have an answer

Every charge-aware ML potential in the literature predicts a flat electrochemical potential across the electrode-electrolyte interface. Not flat as in “approximate” or “averaged” — flat as in mathematically guaranteed, regardless of how carefully the model is trained. The reason is the architecture: standard charge equilibration (QEq) enforces a single electrochemical potential across all atoms by construction. One μ, everywhere. The electric double layer — the potential drop from electrode surface into electrolyte that controls reaction rates, drives selectivity, and determines how a catalyst responds to applied voltage — cannot exist within that framework. Vondrák, Reuter, and Margraf confirmed this in 2025: the suppression is structural, not a training artifact, and it persists in every charge-aware MLIP that has been published.

This matters because electrochemical interfaces sit at the center of how we make green hydrogen, reduce CO₂ to fuels, and store energy in batteries and capacitors. The reaction rate at an electrode surface is not just a function of what atoms are there — it depends on the local potential, the charge of the surface, and the structure of the double layer. Simulate without those, and you are simulating something else.

When I joined LANL’s Theoretical Division (T-1) in early 2024, I came to work on DFT and catalyst mechanisms. Within the first few months it became clear that the charge-aware ML potentials the field was developing for electrochemical interfaces did not actually satisfy the minimum physics requirements for that task. So I built one that does.

What I built — HIPPIE-NN and Soft-FQEq

HIPPIE-NN is a charge-aware ML interatomic potential built on LANL’s existing HIPPYNN PyTorch framework. The design goal was to capture the electric double layer at reactive interfaces — meaning the full spatial potential drop from electrode into electrolyte, not just accurate per-atom charge values.

The core of what makes this different from prior frameworks is the charge solver: Soft-FQEq, short for Soft Fragmented Charge Equilibration. Every prior QEq-based MLIP inherits two architectural pathologies. First, a single electrochemical potential across all atoms suppresses the double layer. Second, with charge free to flow between any atoms, an electrolyte simulation allows water molecules to donate electrons to ions across vacuum — unphysical at the level of individual molecules. Soft-FQEq fixes both by enforcing charge conservation per molecular fragment instead of globally.

That fix creates a new problem: in a reactive simulation where proton transfers and O–O bonds are forming and breaking continuously, you cannot predefine which atoms belong to which molecule. The fragment membership itself is dynamic. Soft-FQEq solves this through a sigmoid-sharpened graph-Laplacian resolvent that identifies fragment membership from neural-network-predicted bond conductances. When a proton transfers, fragment assignments update smoothly — no discrete reassignment, no discontinuity in the potential energy surface.

The full architecture is a single shared HIPPYNN message-passing backbone feeding four MLP readout heads: short-range energy (E_short), atomic electronegativity (χ), source charges (q_c), and a bond-conductance correction (Δlog κ). These feed into an augmented-Lagrangian Uzawa solver that enforces per-fragment charge neutrality at each iteration. Forces come from direct autograd through every Uzawa iteration in float64 — I run K=3 iterations during training and K=5-10 in production. The envelope theorem does not apply here because constraint violation is nonzero at finite iteration count, so the charge–position coupling (dq/dR) is retained explicitly through the unroll — no shortcuts through the solver.

“Every charge-aware ML potential in the literature predicts a flat potential across the electrode-electrolyte interface. That’s not a bug people ignore — it’s a mathematical limitation baked into how they partition charge.”

The ablation that validates the architecture: when I train on DFT energies, forces, and DDEC6 reference charges on IrO₂/water/Na⁺/ClO₄⁻ slabs, the electric double layer emerges from charge-only training. With the per-fragment constraint in place, the potential drops from electrode surface into solution. Remove the constraint — same weights, same architecture, same training — and the potential collapses to uniform everywhere. The architecture is doing the electrochemistry.

The estimated speedup over DFT molecular dynamics is approximately 50,000× on GPU. That translates to nanosecond-scale reactive trajectories at controlled electrode potential on systems that DFT-MD cannot touch: too large, too slow, too expensive.

Bootstrap training data came from an automated interface-builder pipeline that generated and SOAP-FPS-selected roughly 7,000 structures across four chemistry categories — small molecules, water-ion boxes, bulk IrOₓ, and electrode-electrolyte interface slabs — from a much larger auto-generated pool of over 20,000 candidates. An active learning pipeline integrated with LANL’s ALF framework runs ongoing DFT labeling cycles on A100 GPUs via Parsl-orchestrated Slurm workflows, with a 5-model ensemble and four charge-aware uncertainty triggers (force dispersion, charge dispersion, fragment potential dispersion, EDL potential dispersion). A manuscript describing HIPPIE-NN is in preparation with mentors Blas Uberuaga and Travis Jones; an arXiv preprint is forthcoming pending LA-UR clearance.

Details on the software infrastructure — the interface builder, active learning pipeline, and deployment tooling — are on the /software page.

The plumbing — qe-tools and the ESM-RISM sweep engine

Before any of the DFT science described below was possible at the scale we ran it, I needed to solve a different problem: the Quantum ESPRESSO ESM-RISM workflow was powerful but not production-usable without substantial engineering.

ESM-RISM calculates in constant-charge mode: fix the surface charge, run the calculation, get an (energy, potential) pair. To build a continuous energy-vs-potential curve for a single surface configuration, you sweep across many charge states — each requiring its own calculation, its own convergence, its own recovery if the RISM solver fails or the HPC job times out. Then you do this for every configuration in a surface phase diagram, which for an OER study might be 60 configurations or more. The manual version of this workflow consumed 2–3 months of researcher attention per campaign.

I built a 12-tool Python automation framework that runs the whole pipeline unattended: bulk convergence, slab construction for rutile MO₂ and layered MOOH geometries, ESM-RISM constant-charge sweeps with an adaptive threshold ladder spanning 1×10⁻⁸ to 1×10⁻², NEB transition state searches with five-level automatic recovery, XAS simulation via XSpectra, zero-point energy corrections, and publication-ready electrochemical phase diagram generation. A custom Parsl execution provider — ChainedSlurmProvider — pre-queues successor jobs into the Slurm scheduler before the current job completes, cutting roughly 2.5 hours of queue wait per 10-restart convergence sequence. Seven error classes are auto-recovered without human intervention: walltime kills, SCF divergence, RISM convergence failures, stagnation, NEB image corruption, Hubbard card loss, and coordinate unit mismatch.

Result: a full OER campaign runs in 2–3 weeks of unattended compute instead of 2–3 months of researcher attention. Every LANL research project described in the next section ran on this infrastructure.

I also authored a 3-file, 14-line patch to Quantum ESPRESSO 7.5 enabling XSpectra calculations under ESM-RISM boundary conditions. Prior to this patch, operando XAS simulation with an implicit electrolyte was not possible in QE. The patch fixes all QE post-processing codes, not just XSpectra — it is the first implementation of its kind in the codebase.

What this has produced

β-NiOOH OER mechanism. NiOOH-based materials are the dominant anode catalysts in commercial alkaline and anion exchange membrane electrolyzers — a market currently valued around $2.2B and projected to reach $8–18B by 2032. For over a decade, the community has assumed OER on NiOOH proceeds through an oxyl radical (M–O·) coupling step. I constructed exhaustive DFT+U surface phase diagrams for the (101̄0) surface — approximately 64 configurations, PBE+U with U = 5.5 eV, 7-layer slabs — systematically varying hydroxyl, water, oxygen, and bridging coverage across the accessible potential window.

The data do not support the oxyl pathway. OER on β-NiOOH is driven by pseudocapacitive hydroxyl charge storage via water nucleophilic attack. A Brønsted–Frumkin kinetic model fits experimental Tafel slopes at R² = 0.91. On the spectroscopy side, I simulated O K-edge and Ni L-edge XAS using XSpectra and MULTIX; the computed spectra assign specific features in operando synchrotron data collected at BESSY II (Berlin) to specific surface terminations from the phase diagram. Three independent operando spectroscopy datasets from collaborators at FHI Berlin, MPI-CEC, HZB/BESSY II, and CNR Trieste corroborate the computational finding. I owned all computational deliverables in this 17-author, 3-country collaboration. The manuscript is under review at Nature Catalysis.

Pt₃O₄ metallic conductivity. PEM electrolyzers need conducting, corrosion-stable anode supports. Carbon corrodes above ~1.0 V vs. RHE — surface area can drop 70% within 24 hours. Most oxide alternatives are semiconducting or insulating. Pt₃O₄ is metallic, and understanding why matters for finding cheaper alternatives. The prior explanation invoked Pt(II)/Pt(IV) charge disproportionation (mixed valence), but systematic benchmarking tells a different story.

I ran 12 levels of DFT theory — PBE and r2SCAN, with and without spin-orbit coupling, at three Hubbard U values — with COHP orbital decomposition, Bader and DDEC6 charge analysis, and ZSA framework classification. Across all 12 settings, the result is consistent: metallicity arises from persistent Pt–O covalency creating delocalized antibonding states at the Fermi level. Pt³O₄ is a negative-Δ charge-transfer metal in the ZSA framework. Covalency-driven metallicity, not charge disproportionation — a design principle for identifying which oxide compositions will stay conducting under electrochemical conditions. Published in Electronic Structure (2026) as first and co-corresponding author.

PtO₂ dopant screening. No systematic computational screening of transition metal dopants in PtO₂ for conducting oxide supports existed before this work. I screened 28 transition metals across alpha- and beta-PtO₂ at three compositions each — approximately 168 DFT systems total — using automated PBE-to-r2SCAN workflows combined with Pourbaix thermodynamics to assess dissolution stability under PEM electrolyzer conditions. The screen identifies three distinct conductivity mechanisms, providing transferable design rules beyond the specific PtO₂ host. First-author manuscript in preparation.

Cu nanoparticle CO₂RR electrochemical strain. Cu is the only heterogeneous catalyst that reduces CO₂ to multi-carbon products — ethylene, ethanol — with significant Faradaic efficiency. Operando XRD shows 50–100 femtometer lattice expansion in Cu nanoparticles going from cathodic to anodic potential. The naive expectation is contraction: removing electrons should shrink the lattice. The data show the opposite.

Using constant-charge ESM-RISM calculations with adsorbate-resolved surface phase diagrams, I decomposed the electrochemical strain into charge vs. adsorbate contributions for a collaborative study with a German experimental group. The finding: adsorbate-induced surface stress dominates. CO desorption and OH/carbonate adsorption at anodic potentials increases compressive surface stress, reducing Laplace pressure and driving expansion. DFT predicts 100–500 fm expansion — consistent with the 50–100 fm experimental observation. This redirects strain engineering strategy from potential tuning to adsorbate coverage control. Draft exists; co-first with the German experimental collaborators.

CoOOH neutral-pH OER. I extended the ESM-RISM methodology to CoOOH, a benchmark neutral-pH OER catalyst. Voltage-dependent surface phase diagrams show the surface fills with μ₁-OH species at OER-relevant potentials. Computed O K-edge XAS at 528 eV and 529 eV as a function of electrode potential — one peak rises then saturates, the other rises continuously — matching experimental operando trends and enabling direct spectral feature assignment. The ESM-RISM sweep engine made this feasible without extending the campaign timeline. Manuscript in preparation.

Trajectory

The infrastructure now exists. What comes next is deploying it.

The immediate target is production HIPPIE-NN runs at electrochemical interfaces: nanosecond-scale reactive trajectories on IrO₂/water/electrolyte systems at controlled electrode potential, with OER barrier calculations that include explicit solvent and the electric double layer. DFT cannot do this — at ~100 atoms and ~10 picoseconds, AIMD leaves every meaningful timescale unreached. HIPPIE-NN at ~50,000× speedup brings IrO₂ reaction kinetics into reach for the first time.

From there, the framework generalizes. The QE automation toolkit generates DFT training data for any catalyst system; the active learning pipeline builds the model; Soft-FQEq provides physically correct charge distributions and EDL representation. Extending to RuO₂, MnO₂, CoOₓHᵧ for OER screening, or to Cu systems for CO₂RR dynamics, requires new DFT data rather than new architecture.

The open questions this stack lets us address — how catalyst surfaces restructure under potential, how the double layer modulates reaction barriers, how adsorbate strain couples to selectivity — are not questions the field has answered with simulation. They remain open because the tools did not exist until now.

I am also working toward open-source release of the Soft-FQEq solver and QE automation toolkit. The solver architecture is general; every group building charge-aware ML potentials for electrochemical interfaces needs a fragment-constraint implementation that is differentiable, PyTorch-native, and compatible with reactive bond topology. That is what Soft-FQEq is.