The question: for hydrogen at a Ni(111) surface, what are the stable sites, the barriers between them, and the strain dependence — i.e. the thermodynamic and kinetic basis for superstitial/coherent-domain loading scenarios. We computed adsorption energies for all four surface sites and two subsurface interstitial sites, the strain derivative of each, and two full minimum-energy paths (fcc→octa, fcc→tetra) with climbing-image NEB. Method: GPAW (PBE, plane waves, 500 eV cutoff, spin-polarized), 4×4 Ni(111) four-layer slab + dipole correction, k-meshes 3×3×1 (relaxation/NEB) and 4×4×1 (final single points, all quoted numbers). Adsorption energies are referenced to ½H₂ in the gas phase. Validation: the site ladder reproduces four independent published studies (Bhatia–Sholl, Noordhoek et al., Shirazi et al., Henkelman-group anchors) within ≤0.02 eV. Every quoted number re-derives from raw output through an independent audit script (17/17 checks).
The six local minima for H, with adsorption energies Eads relative to ½H₂(g) at the 4×4×1 mesh:
| Site | Geometry | Eads | Interpretation |
|---|---|---|---|
| fcc | Threefold hollow, no atom directly beneath (ABC stacking) | −0.519 eV | Global minimum on clean Ni(111); anchor −0.50 ± 0.02 ✓ (Bhatia–Sholl, Noordhoek, Shirazi) |
| hcp | Threefold hollow above a second-layer atom | −0.503 eV | Quasi-degenerate with fcc (16 meV); standard for (111) surfaces |
| bridge | Twofold bridge | −0.377 eV | Metastable ledge; agrees with Shirazi −0.35 |
| top | On-top | +0.037 eV | Barely bound; agrees with Shirazi +0.05 |
| octa | Subsurface octahedral interstice, directly below the fcc hollow | +0.178 eV | Endothermic vs ½H₂ at zero coverage — the target state for subsurface loading |
| tetra | Subsurface tetrahedral interstice, laterally displaced (see §5) | +0.511 eV | Strongly endothermic; also laterally offset from the fcc column |
k-point bracket: k₄₄₁ vs k₆₆₁ differs by 17.6 meV on fcc → ±9 meV systematic on Eads. Coverage note: the 4×4 cell corresponds to 0.25 ML, which pushes subsurface states up relative to the dilute limit — direction and magnitude consistent with Shirazi's coverage dependence.
Main road, fcc→octa. Ea = 0.846 eV forward, 0.149 eV reverse (k₄₄₁). Literature anchor range for this barrier is 0.60–0.84 eV (Henkelman-group NEB; Bhatia–Sholl) — we sit at the top edge, expected from 0.25 ML coverage in the cell. Kinetic consequence, quantified: with attempt frequency ν ≈ 10¹²⁻¹³ s⁻¹, the residence time of an untrapped octa H is τ = ν⁻¹exp(Erev/kBT) ≈ 30 ps at 300 K and ≈ 1–2 ps at 650 K. The subsurface state is kinetically evaporating on clean Ni.
Side road, fcc→tetra. Ea ≈ 1.07 eV (3×3×1 converged climbing image; the k₄₄₁ refinement was skipped after an infrastructure failure at exit — bounded by the octa band's k₃₃₁→k₄₄₁ shift of ≲60 meV), reverse barrier 0.034 eV: the tetra interstice is a shelf, not a store. Additionally, a structural finding: the tetrahedral interstice below the fcc hollow is laterally displaced by 1.44 Å (the octa interstice is exactly in-column, displacement 0.000 Å). A naive vertical path therefore drives H through the close-packed Ni layer — our first band literally hit +10 eV Pauli walls at the layer plane — and the MEP must dive in-column, traverse the interlayer gap, and descend into the displaced cage. The tetra route is harder for two independent reasons: a higher barrier and a misaligned door. The purple stick in the 3D view shows the reconstructed path.
Each computed quantity maps onto a protocol parameter and an observable. The chain is: number → mechanism → protocol → what you would see.
| Protocol problem | Computed basis | Design consequence | Observable signature |
|---|---|---|---|
| Getting H in (absorption) | Ea = 0.85 eV (octa route): Arrhenius-attempt probability exp(−0.85/kBT) ≈ 10⁻⁸ at 300 K vs ≈ 10⁻⁴ at 700 K; pressure/electrochemical potential shifts the driving force, not the barrier | Load hot — 600–700 K — where thermal activation plus chemical-potential driving can move H over the hill; room-temperature loading of clean Ni cannot populate the subsurface at any practical feed rate | Uptake rate with ~Arrhenius T-dependence; activation energy ≈ 0.85 eV (or ≈ 1.07 eV if the tetra route dominated — it does not) |
| Keeping H in (retention) | Erev = 0.149 eV → τ ≈ 30 ps (300 K); unassisted retention is zero on experimental timescales | Continuous H₂ feed (needle-valve protocols) is mandatory and exactly compensates the leak; steady-state subsurface coverage is set by feed rate × residence time. For true accumulation you need the stabilization (strain, below) or traps (vacancies/voids — binding energies 0.4–0.6 eV in the literature, to be computed next) | Post-loading desorption transients on a ps-leak system look like instant loss at T; a trapped population would instead show staged, activated release (trapping energies set the release temperatures) |
| Making inside favorable (the switch) | dEads/dε ≈ −50 meV/% for subsurface vs −28 meV/% surface; exothermicity flip at ε ≈ +2%. Mechanism: interstitial volume enlargement. Practical strain sources: external tensile strain, or internal — the FM→PM magnetoelastic expansion at the Curie point, tunable via alloying (NiCu, NiFe), plus thermal expansion itself | Use tensile-strained foils/foams or a Curie-matched alloy so that the magnetic transition delivers the expansion at the operating temperature; expect the sign of uptake enthalpy to flip as the material crosses the threshold | A sharp uptake knee / exothermic burst when strain state crosses ≈ +2% (or the alloy equivalent); absence in the unstrained control |
| Coupled feedback | H lowers the Ni moment and (literature) the Curie temperature; the Curie transition changes the lattice; loading therefore feeds back on the switch | Treat loading fraction, strain state, and magnetic state as one coupled system during analysis — not as independent knobs | Hysteresis between loading and de-loading; uptake features tracking TC(loading) rather than fixed TC(0) |
(1) Paramagnetic/disordered-local-moment site ladder at experimental lattice constants vs T — quantifies whether the Curie transition alone delivers the ≈+2%-equivalent stabilization. (2) Vacancy and void trapping energies for H — the quantitative "unless trapped" term with release temperatures. (3) H-induced Curie shift and the loading↔magnetic feedback loop. (4) QE cross-validation of the two barrier numbers with the independent code.
| Tool | Type | Role here | Status |
|---|---|---|---|
| GPAW 26.7 | Plane-wave/real-space DFT | All campaign numbers (PBE-PW, spin-polarized, dipole-corrected slabs) | Validated — 4 literature anchors reproduced |
| ASE 3.29 | Atomistics framework | NEB/climbing-image machinery, FIRE optimization, structure I/O, image-parallel orchestration | Validated — campaign infrastructure |
| Quantum ESPRESSO + SSSP pseudos | Independent plane-wave DFT | Cross-validation of barriers/site energies with a second code — guards against implementation artifacts | Installed, staged for campaign 2 |
| LAMMPS | Classical MD | μs-scale diffusion and many-H loading configurations once a Ni–H potential is fitted to this campaign's DFT data | Installed, awaiting fitted potential |
| SymPy / SciPy | Symbolic & numeric oracles | Independent re-derivation of every quoted claim (mutation audit) | Active — 17/17 |