Vantasner AG Hydride-Piezo Device Program Open-Source Simulation Toolchain

1 Why this toolchain exists

Dan's own density functional theory (DFT) findings frame the whole program: lattice strain, hydride phase transitions, and sharp mechanical transients are the levers that move deuteron behavior in a metal lattice. The device is designed for them, and the toolchain is built around them:

• Lattice strain. Strain shifts hydride site energies, diffusion barriers, and phase stability. Every model in the chain treats strain as an explicit input, never as an afterthought.
• Phase transitions. Palladium deuteride crosses an alpha-beta miscibility gap on the road from low loading to the D/Pd of about 0.7 that electrochemistry delivers; titanium and zirconium deuterides sit in their own fluorite-structure windows near D/Ti and D/Zr of about 2. The simulations resolve where those transitions sit under strain and how fast the fronts move.
• Sharp transients. Frictional force chains in the shear-thickening suspension deliver contact stress pulses, and the 20 kHz ultrasound drive delivers a periodic train of them. The simulations resolve what a single pulse does to a contact, a grain, and a gap.

Three questions the toolchain must answer, in order:

The concept rests on reported experimental anchors: the sub-keV D-D yield plateau with an effective screening potential rising to about 3.0 keV on palladium deuteride at 0.25 keV beam energy, reported by Karahadian, Colborne, Persaud, Schenkel and Munday in Nature Communications 17, 8845 (2026), doi:10.1038/s41467-026-74421-1; the metallic electron screening potentials of a few hundred electronvolts mapped by Czerski et al., Europhysics Letters 54, 449 (2001), doi:10.1209/epl/i2001-00265-7, and Raiola et al., European Physical Journal A 19, 283 (2004), doi:10.1140/epja/i2003-10125-0; and the crystal-field acceleration precedent of Naranjo, Gimzewski and Putterman, Nature 434, 1115 (2005), doi:10.1038/nature03575. Metzler, Hunt, Hagelstein and Galvanetto tally in New Journal of Physics 26, 101202 (2024), doi:10.1088/1367-2630/ad091c, that roughly 40 orders of magnitude of rate enhancement stand between bare-nucleus room-temperature fusion and observability: that arithmetic is exactly what our kernel codes carry around, point by point, for each candidate configuration. Chen et al., Nature 644, 640-645 (2025), doi:10.1038/s41586-025-09042-7, report that electrochemical loading of a metal target raises measured D-D fusion rates: the queue treats loading state as a first-class simulation variable for that reason. Kim, Esser-Kahn, Rowan and Jaeger, PNAS 120, e2310088120 (2023), doi:10.1073/pnas.2310088120, demonstrated piezoelectric nanoparticles generating charge exactly when frictional contact networks form in a shear-thickening suspension: our FEM and DEM chain models precisely that mechanism at device scale.

2 Density functional theory: where every parameter starts

2.1 Quantum ESPRESSO

• License and status. GNU General Public License; current stable release 7.4.1 (checked 2026-10-21). Plane-wave pseudopotential suite supporting norm-conserving, ultrasoft, and projector-augmented-wave (PAW) datasets; density-functional perturbation theory (DFPT) for phonons; pairs with phonopy (BSD-3-Clause) for post-processing and quasi-harmonic approximation (QHA) thermal expansion.
• What it computes for THIS device.
  • Cell relaxations and strain scans of PdDx across the alpha-beta transition region (x from 0.3 to 1.0), TiD2, ZrD2, and NiHx: energy-volume curves, bulk moduli, and formation energies as functions of hydrostatic and epitaxial strain over plus or minus 5 percent. This is Dan's strain finding turned into a production pipeline.
  • Tetrahedral versus octahedral site occupancy energies and site electron densities at each loading: the electron density at interstitial sites is the input the screening and tunneling kernels need.
  • Nudged elastic band (NEB) diffusion barriers for D hopping as a function of applied strain: loading protocol design.
  • DFPT phonons plus QHA thermal expansion and phonon softening around the hydride phase transitions: where the lattice goes mechanically soft, and what expansion mismatch the hydride-piezo pair accumulates on temperature excursions.
• Typical job. 8 to 64 atom cells, PAW datasets, cutoffs near 50 to 80 Ry (estimate), k-meshes from 6x6x6 to 12x12x12 (estimate); 2x2x2 phonon supercells up to about 96 atoms with 100 to 200 displacements per geometry (estimate).
• Runtime (estimate). Relaxation of a 64-atom hydride cell: 2 to 8 hours on 32 cores. An 11-point strain scan for one hydride at one loading: 1 to 2 slot-days. Phonon + QHA sweep over 5 volumes for 4 compositions: 2 to 3 server-days amortized across the queue.
• Constrains. Equilibrium lattice constants and stiffness versus strain; hydride phase stability windows and transition pressures; thermal expansion mismatch strain between hydride grain, piezoelectric grain, and carrier; site electron densities feeding the screening kernels; strain-dependent D diffusivity for the loading protocol.

2.2 GPAW

• License and status. GNU General Public License v3 or later (checked 2026-10-21). Python-based DFT package built on the projector-augmented wave method with three wave-function representations: real-space grids, plane waves, and numerical atomic orbitals (LCAO mode for large cells at low cost). Hybrid functional support (HSE06 class) for selected single points.
• What it computes for THIS device. Metal hydride surfaces and the hydride-piezoelectric interface, which plane-wave bulk codes handle poorly:
  • Deuterated Pd(111), Ti(0001), and Zr(0001) slabs: deuterium adsorption and subsurface absorption energetics, work-function and surface-dipole changes under applied strain.
  • BaTiO3/PdD and PZT/TiD2 interface cells: charge transfer, band offsets, and the local field across a nanometer-scale contact gap, relaxed with the field switched on. This is the direct simulation of the concept's central object: a piezoelectrically charged surface looking across a gap at a deuteride surface.
  • Field-versus-gap-distance curves for capacitor geometries with finite screening: how much of the 10^8 to 10^9 V/m design field survives at 1 to 10 nanometer gaps.
• Typical job. Slabs of 80 to 250 atoms with 15 to 20 angstrom vacuum; interface cells of 150 to 400 atoms (estimate).
• Runtime (estimate). LCAO-mode relaxations: 6 to 24 hours on 32 cores. Hybrid functional single points: 24 to 72 hours on 48 cores.
• Constrains. The field-at-contact budget (design target 10^8 to 10^9 V/m per the concept); surface charge density per unit contact stress; whether deuteride surfaces hold their polarization state in the dry-bed carrier.

2.3 CP2K

• License and status. GNU General Public License v2 or later (checked 2026-10-21). Gaussian and plane-wave (GPW) mixed basis; mixed-precision options; efficient for large cells and for Born-Oppenheimer molecular dynamics (MD).
• What it computes for THIS device.
  • Grain-boundary and interface supercells of 500 to 2000 atoms: defect and boundary energetics where plane-wave codes waste their basis on vacuum.
  • Short ab initio MD trajectories (5 to 20 picoseconds, 1 femtosecond steps, estimate) of a hydride grain under a contact-stress pulse: the femtosecond-scale charge response at contacts, which is exactly the transient that survives ionic screening in an electrolyte carrier.
  • Reference data for training machine-learned interatomic potentials (see section 6).
• Runtime (estimate). 24 to 96 hours on 48 cores for a production trajectory.
• Constrains. Reference charges and energetics at grain boundaries and contacts where the reactive force field must be trustworthy; the sub-picosecond charge redistribution time constant that decides whether piezo gap fields survive in the electrolyte-carrier variant.

3 Molecular dynamics: transients are the device's native language

3.1 LAMMPS

• License and status. GNU General Public License v2 or later; stable release cadence on the project site (checked 2026-10-21). Embedded atom method (EAM) and ReaxFF (reactive force field) parameter families ship with the code; published parameterizations exist across the literature for Pd-H, Ni-H, Ti-H, and Zr-H, with coverage and transferability that vary, so the queue runs a benchmark phase before production use and fits its own parameter sets where coverage is missing.
• What it computes for THIS device.
  • EAM runs for Pd-H and Ni-H: elastic response under strain, H site hopping rates versus strain and temperature, alpha-beta coexistence sampled directly, and contact-impact pulses reproducing force-chain contact stresses in the 0.1 to 1 GPa design range reported for dense shear-thickening suspensions.
  • Shock and transient response by velocity-ramp loading: Hugoniot states, spall thresholds, and pulse duration spectra for stress waves crossing a 1 to 10 micrometer hydride grain.
  • Driven-oscillation studies of hydride slabs at MD-resolvable frequencies (MHz to GHz) to extract strain-rate-dependent strength and damping; those constitutive parameters are handed up to the continuum codes, which carry the actual 20 kHz drive. The honest time-scale boundary: a 20 kHz period is 50 microseconds, far beyond direct MD; MD computes the per-pulse physics, the finite element method (FEM) codes compute the periodic drive.
  • ReaxFF runs for Ti-H and Zr-H and for hydride-carrier interfaces where charge transfer and surface chemistry matter: fracture and friction at grain contacts, surface chemistry relevant to fractoemission backgrounds (light and charged particles from cracking, reported in the fracture literature, and a diagnosed artifact class for the lab, not a nuclear channel).
  • Boundary statement: MD tracks atoms at eV kinetic energies. Deuteron energetics in the 0.1 to 1 keV class across piezo gaps are computed by the screening and tunneling kernels (section 5), not by MD.
• Typical job. 10^5 to 10^7 atoms, 10^6 to 10^9 timesteps at 0.1 to 1 femtoseconds (estimate).
• Runtime (estimate). EAM transient of 10^6 atoms: about 96 core-hours. ReaxFF mid-size cell: about 768 core-hours.
• Constrains. Contact-pulse amplitude and duration distributions (feed the piezo charge kernel); maximum safe drive amplitude; particle size effects on strain transfer; H diffusivity versus strain for the loading protocol.

3.2 ASE (Atomic Simulation Environment)

• License and status. GNU Lesser General Public License (checked 2026-10-21).
• Role. The workflow glue and input generator for everything above: one Python object model (Atoms, calculators, filters) generates Quantum ESPRESSO, GPAW, CP2K, and LAMMPS inputs; runs NEB barrier scans; post-processes trajectories; and emits the queue submission records. The agents on Marbell speak to the simulation codes exclusively through ASE factories plus template files, so every input in the database is reproducible from a code commit hash plus a configuration hash.
• Typical job. Milliseconds to minutes of single-core Python per generated input (estimate).
• Constrains. Nothing physical itself; it constrains error rates. Uniform input generation is what makes the results database trustworthy at thousands of jobs per month.

4 Phase-field: hydride phase fronts under strain

• PRISMS-PF. GNU Lesser General Public License v2.1 (checked 2026-10-21). Massively parallel finite element phase-field framework from the University of Michigan PRISMS center, matrix-free solvers on the deal.II library.
• MOOSE. GNU Lesser General Public License v2.1 (checked 2026-10-21). Multiphysics Object-Oriented Simulation Environment from Idaho National Laboratory: phase field, solid mechanics, and heat transfer modules; the ecosystem carries a nuclear-materials heritage, including zirconium hydride precipitation models in the Bison fuel-performance code, which is exactly the mathematical neighborhood our hydride precipitation problem lives in.
• Custom FEniCSx kernel, if needed. See section 6 for the license. If the coupled electro-chemo-mechanical term (piezoelectric field entering the diffusion chemical potential) proves to matter beyond what the two applications expose, a FEniCSx kernel implements it directly.
• What they compute for THIS device.
  • Hydride precipitation fronts: Cahn-Hilliard free-energy landscapes fitted to the DFT formation energies of section 2; alpha-beta interface mobility; front motion under oscillatory mechanical stress from the 20 kHz drive (drive-on versus drive-off).
  • Hydride sleeve growth on powder particles during ex-situ gas loading: loading ramp protocol design (temperature and pressure schedules).
  • Strain as a thermodynamic variable in the diffusion potential: Dan's finding operationalized. Precipitates nucleate where tensile strain concentrates; the simulations answer whether the oscillatory strain of the drive keeps phase fronts mobile or locks them in place.
• Typical job. 2D grids near 1024 by 1024 with coupled mechanics; 3D grids of 192 to 256 cells per side (estimate).
• Runtime (estimate). 2D mechanics-coupled run: about 384 core-hours. 3D run: about 2300 core-hours.
• Constrains. Loading ramp protocol; duty cycle (drive-on fraction) choice; phase-front timescales relative to the drive period; expected precipitate sizes for comparison with microscopy once the lab exists.

5 Screening and tunneling kernels (small custom Python)

• What they compute.
  • The D-D Gamow penetrability exp(-31.29 / sqrt(E_keV)) and the standard adiabatic screening enhancement in the form f(Ue) = exp(31.29 x Ue / (sqrt(E) x (E+Ue)^1.5)), with E and Ue in keV, as consolidated for metallic environments by Huke, Czerski, Heide, Ruprecht, Targosz and Zebrowski, Physical Review C 78, 015803 (2008), doi:10.1103/PhysRevC.78.015803.
  • Astrophysical S-factor integration over the measured sub-keV plateau; branching between the proton and neutron channels per measured ratio.
  • Scan grids: screening potential Ue from 0 to 3 keV crossed with effective deuteron energy from 10 eV to 10 keV; each point yields rates per cubic meter per second for the local deuteron density of a loaded hydride (order 5e28 per cubic meter at PdD0.7, order-of-magnitude context from the concept brief).
  • Fitting chains: DFT site electron densities (section 2) to configuration-specific screening estimates, benchmarked against the measured metallic screening potentials of a few hundred electronvolts (Czerski et al. 2001; Raiola et al. 2004) and against the reported sub-keV plateau enhancement above 1e18 relative to bare-nucleus rates below about 2 keV on palladium deuteride (Karahadian et al. 2026). Piezo-gap field maps (section 6) to effective deuteron energy spectra.
  • Signal-class kernels: discrete field-accelerated deuteron acceleration across a nanometer gap predicts monoenergetic proton lines near 3.02 MeV; heat-driven distributions predict continua. The kernels emit predicted spectra so the lab can plan detector exposures before the first run.
• Job shape. Single core, seconds to minutes per scan point (estimate); millions of points per campaign (estimate). Cost per point is about CHF 0.01 (estimate, see queue document).
• Constrains. Expected deuteron energy spectra; minimum detectable yield targets and detector exposure planning; the sensitivity of observable rates to Ue, gap field, and gap width: the three numbers the whole device concept trades on.

6 FEM acoustics and piezoelectric coupling

6.1 Elmer

• License and status. GNU General Public License (checked 2026-10-21). Open-source multiphysics FEM suite developed mainly by CSC, the IT Center for Science in Finland; solver library covers elasticity, electrostatics, acoustics (Helmholtz/wave equation), and heat transfer, coupled through a multi-solver framework.
• What it computes for THIS device.
  • Ultrasonic horn and suspension cell acoustics: eigenfrequency scans of the horn-stack plus cell assembly, mode shapes, quality factors, and stress-amplitude distributions versus drive frequency, with the suspension represented as an effective medium whose modulus and viscosity come from the rheology fits (section 7).
  • Resonance tracking versus temperature and loading state: how far the cell drifts as the suspension shear-thickens, and where the drive control loop must hold it.
  • Power ceiling versus cavitation onset margin: the acoustic amplitude at which bubble activity begins, an artifact class the lab diagnoses (sonoluminescence, neutron-counting bursts from cavitation in the prior record, for example Jorne, Fusion Technology 29, 1996), and that the design simply stays below.
  • Electrostatics of piezoelectric grains in representative contact cells: stress-to-boundary-charge conversion with the measured d33 of each candidate piezoelectric (PZT, barium titanate, zinc oxide, lithium niobate), and field distribution across gaps with finite screening set by carrier chemistry.
  • Piezoelectric coupling: Elmer couples its elasticity and electrostatics solvers to represent piezoelectric constitutive behavior; where a monolithic piezo formulation with time-domain transients is needed, the FEniCSx kernel below takes over.
• Typical job. 10^5 to 2x10^6 elements (estimate).
• Runtime (estimate). 1 to 48 hours on 16 to 32 cores per model point.
• Constrains. Drive frequency selection (20 kHz versus 40 kHz class); horn and cell geometry; input power ceiling; expected contact stress distribution in the suspension.

6.2 FEniCSx / DOLFINx

• License and status. GNU Lesser General Public License v3 or later (checked 2026-10-21). Next-generation FEniCS problem-solving environment.
• What it computes for THIS device. Custom coupled elastodynamics plus electrostatics on representative force-chain cells: a piezoelectric grain between two hydride grains, driven at 20 kHz, with the strain-charge constitutive law solved monolithically in the time domain. Charge at every contact versus time; field-per-gap histograms; wave transit checks at the millimeter cell scale (microsecond regime, well resolved).
• Typical job. 10^4 to 10^6 elements (estimate).
• Runtime (estimate). Minutes to hours on 4 to 32 cores.
• Constrains. The field-at-gap budget versus carrier choice (dry bed versus deuterated oil versus electrolyte, with Debye-length arithmetic deciding what survives in each); duty cycle; the contact-by-contact charge series that feeds the screening kernels.

7 Optional adjuncts

• LIGGGHTS-PUBLIC (discrete element method, DEM). GPL (checked 2026-10-21), released by DCS Computing within the CFDEM project; the vendor's commercial successor Aspherix exists alongside the open-source code. DEM for the dry granular bed and for dense suspension contact networks: force-chain statistics, contact-stress histograms versus packing fraction and piezoelectric volume fraction (the concept baseline is about 10 volume percent PZT with about 50 volume percent hydride), with lubrication-adjusted contact laws for the fluid-loaded cases.
• Suspension rheology models. Python and FEniCSx effective-medium fits to measured storage and loss moduli versus applied stress; outputs feed Elmer's acoustic effective medium and the phase-field stress inputs. First-principles DEM with lubrication laws brackets the dry and wet limits.
• Machine-learned interatomic potentials (MLIP). MACE (Apache License 2.0, checked 2026-10-21; Batatia, Kovacs, Simm, Ortner and Csanyi, Advances in Neural Information Processing Systems 35, 2022). Trained on CP2K and Quantum ESPRESSO reference trajectories, it runs DFT-accuracy transient dynamics at force-field cost and closes the accuracy gap between ReaxFF and ab initio MD for shock-like contact pulses. Dataset generation is CP2K class K6 in the queue taxonomy; training is class K12.
• Visualization. ParaView (BSD 3-Clause, checked 2026-10-21) and OVITO (OVITO Basic binaries under an MIT license, source under GPLv3, OVITO Pro is a commercial product; checked 2026-10-21) for trajectories, phase-field isosurfaces, and FEM fields.
• Uncertainty quantification. SALib (BSD, checked 2026-10-21) for Sobol sensitivity analysis over the coupled model chain before the design-lock month.

8 Orchestration glue

ASE job factories generate inputs; Dask and Dask distributed (BSD 3-Clause, checked 2026-10-21) execute them across the Marbell compute pools; Valkey (BSD 3-Clause, Linux Foundation project, protocol-compatible with Redis, checked 2026-10-21) holds the queues; PostgreSQL holds the shared results database; Prometheus and Grafana watch the fleet. The full architecture, job taxonomy, throughput model, and campaign plan live in simulation/MARBELL-SIMULATION-QUEUE.md.

9 Validation ladder

Every tool enters production only after it clears its rung:

1. DFT anchors. Reproduce face-centered-cubic palladium lattice constant near 3.89 angstrom, the beta-PdD lattice expansion on the order of 3 to 4 percent on full loading (reported experimental range), TiH2 and ZrH2 fluorite stability, and bulk moduli against literature values.
2. MD anchors. Elastic constants and H solution energies versus the DFT results of rung 1; Hugoniot states against tabulated data.
3. Phase-field anchors. Dilute-limit chemical potential against DFT; interface width and mobility against MD.
4. Kernel anchors. Bare Gamow and S-factor arithmetic against standard compilations; enhancement curves against the Huke et al. 2008 consolidation.
5. FEM anchors. Horn impedance spectrum against network-analyzer measurements of the physical horn once hardware exists (lab link, months 1 to 3 of the build plan).

Every results row in the shared database carries provenance: code version, container digest, input factory commit, configuration hash, and checkpoint path.

10 License and version table

11 Fit to the Marbell queue

Each tool maps to a job class with a fixed core-hour budget in the queue taxonomy (K0 to K12 in simulation/MARBELL-SIMULATION-QUEUE.md). The campaign plan there schedules which tool runs in which month and states, for each campaign, which lab parameter it narrows: hydride choice, carrier choice, piezoelectric volume fraction, particle size window, drive frequency, power ceiling, duty cycle, loading ramp, and detector exposure.