From 6fba5886b4adb9c22a81a5decfa85024e73d1982 Mon Sep 17 00:00:00 2001
From: Matthew Pharr <m.pharr@protonmail.com>
Date: Fri, 23 Oct 2026 17:20:15 +0200
Subject: [PATCH 1/3]  Repo - API! - Delete the no-op KineticForces and
 PerturbedEquilibrium keys and stop runs that set deprecated keys

---
 .pre-commit-config.yaml                       |  2 +-
 benchmarks/benchmark_against_fortran_run.jl   |  2 -
 docs/development/hdf5-conventions.md          |  6 +-
 docs/development/toml-conventions.md          |  6 +-
 docs/src/kinetic_forces.md                    | 29 +-------
 .../DIIID-like_error_field_example/gpec.toml  |  3 -
 .../gpec.toml                                 |  2 -
 examples/DIIID-like_ideal_example/gpec.toml   |  3 -
 .../gpec.toml                                 |  2 -
 examples/Solovev_ideal_example/gpec.toml      |  2 -
 examples/Solovev_ideal_example_3D/gpec.toml   |  2 -
 .../Solovev_kinetic_NTV_example/gpec.toml     |  2 -
 .../gpec.toml                                 |  2 -
 .../gpec.toml                                 |  2 -
 src/Equilibrium/EquilibriumTypes.jl           |  8 +--
 src/GeneralizedPerturbedEquilibrium.jl        | 26 +++----
 src/KineticForces/BounceAveraging.jl          | 24 +++----
 .../CalculatedKineticMatrices.jl              |  8 +--
 src/KineticForces/Compute.jl                  |  3 -
 src/KineticForces/EnergyIntegration.jl        | 16 ++---
 src/KineticForces/KineticForces.jl            |  1 -
 src/KineticForces/KineticForcesStructs.jl     | 45 +-----------
 src/KineticForces/Output.jl                   | 71 -------------------
 src/KineticForces/PitchIntegration.jl         | 21 +++---
 src/KineticForces/Torque.jl                   | 15 ++--
 .../PerturbedEquilibrium.jl                   | 10 +--
 .../PerturbedEquilibriumStructs.jl            | 21 +-----
 src/PerturbedEquilibrium/Response.jl          |  2 +-
 src/PerturbedEquilibrium/ResponseMatrices.jl  |  4 +-
 src/PerturbedEquilibrium/SingularCoupling.jl  |  4 +-
 src/Rerun.jl                                  |  6 +-
 test/runtests_equil.jl                        | 23 +++---
 test/runtests_h5_schema.jl                    |  2 +-
 test/runtests_kinetic.jl                      | 16 ++---
 test/runtests_kinetic_validity.jl             |  4 +-
 test/runtests_rerun_from_h5.jl                | 17 ++---
 36 files changed, 106 insertions(+), 306 deletions(-)

diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 5efde01b0a..697ae70a2f 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -77,5 +77,5 @@ repos:
       - id: toml-no-deprecated-keys
         name: 'TOML conventions: no deprecated config keys'
         language: pygrep
-        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid)\s*='
+        entry: '^(mer_flag|force_wv_symmetry|ode_flag|cyl_flag|mat_flag|reform_eq_with_psilim|use_riccati|use_parallel|parallel_threads|populate_dense_xi|gal_flag|nstep|diagnose_ca|diagnose|ksing|power_bp|power_b|power_r|power_rc|use_galgrid|moment|eq_out|theta_out|xlmda_out|eqpsi_out|save_records|gpec_file|fnml_flag|ellip_flag|diag_flag|force_xialpha|fixed_boundary|output_eigenmodes|filter_modes|singular_point_method)\s*='
         files: ^(examples/.*\.toml|test/test_data/.*\.toml)$
diff --git a/benchmarks/benchmark_against_fortran_run.jl b/benchmarks/benchmark_against_fortran_run.jl
index 8419953015..a4de7bdcdb 100644
--- a/benchmarks/benchmark_against_fortran_run.jl
+++ b/benchmarks/benchmark_against_fortran_run.jl
@@ -439,7 +439,6 @@ function write_gpec_toml_coil(
         end
         println(io)
         println(io, "[PerturbedEquilibrium]")
-        println(io, "fixed_boundary            = false")
         println(io, "compute_response          = true")
         println(io, "compute_singular_coupling = true")
         println(io, "verbose                   = true")
@@ -491,7 +490,6 @@ function write_gpec_toml_file(
         println(io, "forcing_data_format = \"ascii\"")
         println(io)
         println(io, "[PerturbedEquilibrium]")
-        println(io, "fixed_boundary            = false")
         println(io, "compute_response          = true")
         println(io, "compute_singular_coupling = true")
         println(io, "verbose                   = true")
diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md
index 515186317d..01bfee8bb7 100644
--- a/docs/development/hdf5-conventions.md
+++ b/docs/development/hdf5-conventions.md
@@ -7,7 +7,7 @@ Conventions for the structure and naming of the `gpec.h5` output file. Every wri
 **The schema must be intuitive to a plasma physicist who is not a developer of this code.** The top level largely mirrors the TOML sections / major `src/` modules — which are themselves organized along physics lines — but **physics intuition wins whenever the two diverge**. Worked examples of that rule:
 
 - All per-rational-surface stability results consolidate under `SingularSurfaces/`, regardless of which algorithm produced them: the GGJ coefficients sit beside the Δ′/PEST-3 matching results, and the Riccati BVP and the Galerkin outer-region solve write the *same* `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_*` paths rather than each owning a producer-named subgroup.
-- `EulerLagrangeMatrices` over a bare `Matrices` — group names must say what the data *is*, not which array it came from. Same reasoning renamed `records/` → `EnergyIntegrals/` and `matrices_<method>/` → `KineticMatrices/`.
+- `EulerLagrangeMatrices` over a bare `Matrices` — group names must say what the data *is*, not which array it came from. Same reasoning renamed `matrices_<method>/` → `KineticMatrices/`.
 
 Physics-topic groups elevated to top level (rather than nested under their producer): `Info/`, `Input/`, `SingularSurfaces/`, `LocalStability/`, `SurfaceGeometries/`.
 
@@ -44,7 +44,7 @@ Top level (11 groups):
 | `LocalStability/` | Mercier `D_I`, resistive interchange `D_R`, `ballooning_Delta_prime` on `psi`; the ballooning α boundary on `ballooning_psi` |
 | `SingularSurfaces/` | Per-rational-surface data: `rational_psi`/`rational_q`/`rational_m`/`rational_n`, GGJ coefficients, `Delta_prime_matrix`/`Delta_prime_raw`/`Delta_coil`/`pest3_A`/`pest3_B`/`pest3_Gamma` (Riccati or Galerkin alike), the inner-layer `Match/` results, `Kinetic/` |
 | `PerturbedEquilibrium/` | `ForcingModes/`, `Response/`, `ResponseMatrices/`, `SingularCoupling/`, `Energies/`, control-surface spectra |
-| `KineticForces/` | `<method>/` (torque/energy profiles, `EnergyIntegrals/`, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` with the same per-method layout, summing to the top-level total |
+| `KineticForces/` | `<method>/` (torque/energy profiles, `KineticMatrices/`); multi-ion runs add `PerSpecies/<species>/<method>/` with the same per-method layout, summing to the top-level total |
 | `ErrorFields/` | `CoilSensitivities/` (per-coil-set control-surface spectra and their rigid shift/tilt derivatives, `DominantMode/` full-window projection); `MonteCarlo/` (intrinsic and corrected `\|δ\|` histograms over the sampled tolerances, per batch and averaged); `Risk/` (threshold density, `P(lock\|δ)`, locking probabilities, `ToleranceScan/`); `NTV/` (correction-coil overlap and NTV torque couplings per kAt) |
 | `Tearing/` | `PerSurface/` (+ `DpMatrix/`), `Roots/`, `LayerWidths/`, `Diagnostics/{ValidRoots,Poles,FilteredRoots}`, `Scan/Surface_<k>/` |
 | `SurfaceGeometries/` | `{Plasma,Wall}/{x,y,z}` point clouds |
@@ -69,7 +69,7 @@ Mechanism: writers stay table-driven — each writer keeps a `path => (; long_na
 - Complex quantities are stored as the native HDF5.jl compound type (readable by h5py as a compound dtype) — **never split into `*_real`/`*_imag` dataset pairs**. Sole sanctioned exception: `Input/RawInputs/ForcingTerms/amplitude_{real,imag}`, which mirrors the external forcing ingest-file format and keeps pre-existing snapshots replayable.
 - `NaN` is the not-computed sentinel in numeric datasets (e.g. auto-derived settings, rootless growth-rate entries).
 - A **zero-extent array** is the not-computed sentinel for whole datasets that a given run never produces (e.g. `SingularSurfaces/ca_left`/`ca_right` on kinetic or galerkin-matched runs, the free-boundary energies when `vac_flag=false`, the on-demand derivative stores). Never write unpopulated (`undef`) memory.
-- Ragged (variable-length) data uses the flat-plus-`offsets` companion pattern (`offsets[k+1] - offsets[k]` = length of row `k`) rather than HDF5 VLEN types, e.g. `KineticForces/<method>/EnergyIntegrals/` and `Tearing/Diagnostics/*`.
+- Ragged (variable-length) data uses the flat-plus-`offsets` companion pattern (`offsets[k+1] - offsets[k]` = length of row `k`) rather than HDF5 VLEN types, e.g. `Tearing/Diagnostics/*`.
 
 ## Back-compatibility policy
 
diff --git a/docs/development/toml-conventions.md b/docs/development/toml-conventions.md
index 5ee2c3d5d8..97139921ae 100644
--- a/docs/development/toml-conventions.md
+++ b/docs/development/toml-conventions.md
@@ -47,9 +47,9 @@ fails when its pattern matches a violation) lint-check the covered TOMLs on comm
 - `toml-no-decorative-dividers` — no `# ----`-style divider lines (all covered files).
 - `toml-inline-annotations` — every `key = value` line carries a trailing `# description`
   (`examples/*` and `test/test_data/*` only, per the rule-2 exception above).
-- `toml-no-deprecated-keys` — no deprecated config keys; its key list mirrors
-  `_DEPRECATED_FFS_KEYS`/`_DEPRECATED_EQUIL_KEYS` in `src/GeneralizedPerturbedEquilibrium.jl`,
-  so extend the hook's pattern whenever a key is deprecated there.
+- `toml-no-deprecated-keys` — no deprecated config keys; its key list mirrors the
+  `_DEPRECATED_*_KEYS` lists in `src/GeneralizedPerturbedEquilibrium.jl`, which stop any run that
+  sets one, so extend the hook's pattern whenever a key is added there.
 
 Run them manually with:
 
diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md
index ad866cf06b..3b2ecf7106 100644
--- a/docs/src/kinetic_forces.md
+++ b/docs/src/kinetic_forces.md
@@ -70,8 +70,7 @@ The summed total is written to `KineticForces/<method>/` exactly as in a single-
 run, and each species' own contribution to `KineticForces/PerSpecies/<label>/<method>/`
 — labels are `ion_z<z>_m<m>`, `impurity_z<z>_m<m>`, and `electron` (a numeric suffix
 is appended only if a run repeats a `(z, m)` pair). Every per-species group carries the
-same datasets and metadata as the total, so `dTdpsi`, `T`, and `EnergyIntegrals/` are
-available per species. The summed cumulative torque profile is a diagnostic (linear
+same datasets and metadata as the total, so `dTdpsi` and `T` are available per species. The summed cumulative torque profile is a diagnostic (linear
 interpolation onto the union grid); the summed `total_torque` scalar is the exact
 Gauss-Kronrod value.
 
@@ -171,35 +170,13 @@ dissipation ``-2n\,\mathrm{Im}\langle \xi, u_2\rangle/4\mu_0`` — two independe
 of the same quantity — agree to a fraction of a percent with the regularization off, and disagree
 by tens of percent with it on.
 
-## HDF5 outputs: complex torque convention and the EnergyIntegrals layout
+## HDF5 outputs: complex torque convention
 
 The method level of `KineticForces/<method>/` reports the two physical scalars a user
 wants first: `total_torque` = ``T_\phi`` (N·m) and `total_energy` = ``\delta W_k`` (J).
 Internally both are halves of one complex quantity ``T = T_\phi + 2in\,\delta W_k``
 (Logan 2013), and the ψ-profiles `dTdpsi` and `T` store that complex ``T`` directly —
-so `imag(T)` carries the ``2n`` factor while `total_energy` has it divided out. The
-per-record `EnergyIntegrals/torque` and `EnergyIntegrals/kinetic_energy` are separate
-complex diagnostics of the two integrand halves at each ``(\psi, \lambda, \ell)``
-evaluation, which is why they are not packed into one number there.
-
-`EnergyIntegrals/` stores the variable-length integration trajectories in the
-flat-plus-offsets ragged layout (chosen over HDF5 VLEN types for cross-language
-support; `Tearing/Diagnostics/*` uses the same pattern). Record `k` spans
-`offsets[k]+1 : offsets[k+1]` (Julia, 1-based) of each `*_all` array:
-
-```julia
-h5open("gpec.h5", "r") do f
-    g = f["KineticForces/fgar/EnergyIntegrals"]
-    off = read(g["trajectory_offsets"])
-    x_k = read(g["x_all"])[off[k]+1:off[k+1]]          # record k's abscissae
-    I_k = read(g["integrand_all"])[off[k]+1:off[k+1]]  # its complex integrand
-end
-```
-
-```python
-g = f["KineticForces/fgar/EnergyIntegrals"]            # h5py, 0-based
-x_k = g["x_all"][g["trajectory_offsets"][k]:g["trajectory_offsets"][k + 1]]
-```
+so `imag(T)` carries the ``2n`` factor while `total_energy` has it divided out.
 
 ```@autodocs
 Modules = [GeneralizedPerturbedEquilibrium.KineticForces]
diff --git a/examples/DIIID-like_error_field_example/gpec.toml b/examples/DIIID-like_error_field_example/gpec.toml
index a4f303562c..eee1c3e353 100644
--- a/examples/DIIID-like_error_field_example/gpec.toml
+++ b/examples/DIIID-like_error_field_example/gpec.toml
@@ -196,8 +196,6 @@ height = -1.5780                        # Hoop height: centroid Z of the pack cr
 currents = [55000.0]                    # 1 kA per turn × 55 turns; sensitivities scale linearly with current [A]
 
 [PerturbedEquilibrium]
-fixed_boundary = false                  # Use fixed boundary conditions
-output_eigenmodes = true                # Output eigenmode fields as b-fields
 compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = false                         # Enable verbose logging
@@ -250,7 +248,6 @@ impurity_mass = 12                      # Impurity mass [proton masses]
 electron = false                        # Ion-only NTV
 nutype = "harmonic"                     # Collision operator
 f0type = "maxwellian"                   # Equilibrium distribution
-moment = "pressure"                     # Pressure-moment NTV torque
 atol_pitch = 1e-9                       # Absolute tolerance for inner pitch + energy integrations
 rtol_pitch = 1e-5                       # Relative tolerance for inner pitch + energy integrations
 write_outputs_to_HDF5 = true            # Write outputs to the HDF5 file
diff --git a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
index 5398fd0d6b..f291a598cb 100644
--- a/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
+++ b/examples/DIIID-like_gal_resistive_pe_example/gpec.toml
@@ -93,8 +93,6 @@ name = "c"                              # DIII-D C-coil (6 coils); resolves to c
 currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0]  # n=1 cosine phasing, 1 kA peak [A]
 
 [PerturbedEquilibrium]
-fixed_boundary = false                  # Use fixed boundary conditions
-output_eigenmodes = true                # Output eigenmode fields as b-fields
 compute_response = false                # Compute plasma response to forcing
 compute_singular_coupling = false       # Compute singular layer coupling metrics
 verbose = true                          # Enable verbose logging
diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml
index 448ac5d078..cb9378287d 100644
--- a/examples/DIIID-like_ideal_example/gpec.toml
+++ b/examples/DIIID-like_ideal_example/gpec.toml
@@ -65,8 +65,6 @@ name = "c"                              # Coil set; resolves to coil_geometries/
 currents = [1000.0, 500.0, -500.0, -1000.0, -500.0, 500.0]  # n=1 cosine phasing, 1 kA peak [A]
 
 [PerturbedEquilibrium]
-fixed_boundary = false                  # Use fixed boundary conditions
-output_eigenmodes = true                # Output eigenmode fields as b-fields
 compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = true                          # Enable verbose logging
@@ -91,7 +89,6 @@ impurity_mass = 12                      # Impurity mass [proton masses]
 electron = false                        # Ion-only NTV
 nutype = "harmonic"                     # Collision operator
 f0type = "maxwellian"                   # Equilibrium distribution
-moment = "pressure"                     # Pressure-moment NTV torque
 atol_pitch = 1e-9                       # Absolute tolerance for inner pitch + energy integrations
 rtol_pitch = 1e-5                       # Relative tolerance for inner pitch + energy integrations
 axis_validity_suppression = true        # Suppress kinetic terms inside the near-axis orbit-width validity limit
diff --git a/examples/DIIID-like_ideal_multi_n_example/gpec.toml b/examples/DIIID-like_ideal_multi_n_example/gpec.toml
index 325a248545..898db615e9 100644
--- a/examples/DIIID-like_ideal_multi_n_example/gpec.toml
+++ b/examples/DIIID-like_ideal_multi_n_example/gpec.toml
@@ -67,8 +67,6 @@ name = "c"                              # Coil set; resolves to coil_geometries/
 currents = [1800.0, -50.0, -450.0, -800.0, -450.0, -50.0]  # Sum of n=1, 2, 3 cosine phasings with 1000, 500, 300 A peaks; n=3 lies outside nn_high [A]
 
 [PerturbedEquilibrium]
-fixed_boundary = false                  # Use fixed boundary conditions
-output_eigenmodes = true                # Output eigenmode fields as b-fields
 compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = false                         # Enable verbose logging
diff --git a/examples/Solovev_ideal_example/gpec.toml b/examples/Solovev_ideal_example/gpec.toml
index a42fe4cb02..87fb2ad36f 100644
--- a/examples/Solovev_ideal_example/gpec.toml
+++ b/examples/Solovev_ideal_example/gpec.toml
@@ -36,8 +36,6 @@ forcing_data_file = "forcing.dat"       # Manual mode table (n, m, amplitude_rea
 forcing_data_format = "ascii"           # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires)
 
 [PerturbedEquilibrium]
-fixed_boundary = false                  # Use fixed boundary conditions
-output_eigenmodes = true                # Output eigenmode fields as b-fields
 compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = true                          # Enable verbose logging
diff --git a/examples/Solovev_ideal_example_3D/gpec.toml b/examples/Solovev_ideal_example_3D/gpec.toml
index ab1909e15d..411b33d980 100644
--- a/examples/Solovev_ideal_example_3D/gpec.toml
+++ b/examples/Solovev_ideal_example_3D/gpec.toml
@@ -57,8 +57,6 @@ equal_arc_wall = false                  # Equal arc length distribution of nodes
 # [PerturbedEquilibrium]
 # forcing_data_file = "forcing.dat"      # Manual mode table (n, m, amplitude_real, amplitude_imag)
 # forcing_data_format = "ascii"          # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires)
-# fixed_boundary = false                 # Use fixed boundary conditions
-# output_eigenmodes = true               # Output eigenmode fields as b-fields
 # compute_response = true                # Compute plasma response to forcing
 # compute_singular_coupling = true       # Compute singular layer coupling metrics
 # verbose = true                         # Enable verbose logging
diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml
index cd02c559f5..afd921754c 100644
--- a/examples/Solovev_kinetic_NTV_example/gpec.toml
+++ b/examples/Solovev_kinetic_NTV_example/gpec.toml
@@ -38,8 +38,6 @@ forcing_data_file = "forcing.dat"       # Manual mode table (n, m, amplitude_rea
 forcing_data_format = "ascii"           # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires)
 
 [PerturbedEquilibrium]
-fixed_boundary = false                  # Use fixed boundary conditions
-output_eigenmodes = true                # Output eigenmode fields as b-fields
 compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = true                          # Enable verbose logging
diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml
index a6bce31b95..6e5bef8776 100644
--- a/examples/Solovev_kinetic_calculated_example/gpec.toml
+++ b/examples/Solovev_kinetic_calculated_example/gpec.toml
@@ -33,8 +33,6 @@ forcing_data_file = "forcing.dat"       # Manual mode table (n, m, amplitude_rea
 forcing_data_format = "ascii"           # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires)
 
 [PerturbedEquilibrium]
-fixed_boundary = false                  # Use fixed boundary conditions
-output_eigenmodes = true                # Output eigenmode fields as b-fields
 compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = false                         # Enable verbose logging
diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml
index bee59220ed..f66eeac158 100644
--- a/examples/Solovev_kinetic_multiion_example/gpec.toml
+++ b/examples/Solovev_kinetic_multiion_example/gpec.toml
@@ -38,8 +38,6 @@ forcing_data_file = "forcing.dat"       # Manual mode table (n, m, amplitude_rea
 forcing_data_format = "ascii"           # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires)
 
 [PerturbedEquilibrium]
-fixed_boundary = false                  # Use fixed boundary conditions
-output_eigenmodes = true                # Output eigenmode fields as b-fields
 compute_response = true                 # Compute plasma response to forcing
 compute_singular_coupling = true        # Compute singular layer coupling metrics
 verbose = true                          # Enable verbose logging
diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl
index c195b3b549..1038b33ccc 100644
--- a/src/Equilibrium/EquilibriumTypes.jl
+++ b/src/Equilibrium/EquilibriumTypes.jl
@@ -13,10 +13,10 @@ specified in the input.
   - `eq_type::String` - Type of equilibrium file ("efit", "solovev", "lar", etc.)
   - `eq_filename::String` - Path to equilibrium input file
   - `jac_type::String` - Jacobian coordinate type ("hamada", "pest", "equal_arc", "boozer", "park", "custom")
-  - `power_bp::Int` - Poloidal field power exponent for Jacobian (internal, derived from `jac_type`; deprecated as a TOML input key)
-  - `power_b::Int` - Total field power exponent for Jacobian (internal, derived from `jac_type`; deprecated as a TOML input key)
-  - `power_r::Int` - Major radius power exponent for Jacobian (internal, derived from `jac_type`; deprecated as a TOML input key)
-  - `power_rc::Int` - Minor radius (rfac = √((R-R₀)²+(Z-Z₀)²)) power exponent for Jacobian (internal, derived from `jac_type`; deprecated as a TOML input key)
+  - `power_bp::Int` - Poloidal field power exponent for Jacobian (internal, derived from `jac_type`; not a deck key)
+  - `power_b::Int` - Total field power exponent for Jacobian (internal, derived from `jac_type`; not a deck key)
+  - `power_r::Int` - Major radius power exponent for Jacobian (internal, derived from `jac_type`; not a deck key)
+  - `power_rc::Int` - Minor radius (rfac = √((R-R₀)²+(Z-Z₀)²)) power exponent for Jacobian (internal, derived from `jac_type`; not a deck key)
   - `jac_custom_power_bp::Int` - Poloidal-field exponent for a custom Jacobian `J = Bp^bp · B^b / (R^r · rfac^rc)`; read only when `jac_type = "custom"`, ignored for named types
   - `jac_custom_power_b::Int` - Total-field exponent for a custom Jacobian; read only when `jac_type = "custom"`, ignored for named types
   - `jac_custom_power_r::Int` - Major-radius exponent for a custom Jacobian; read only when `jac_type = "custom"`, ignored for named types
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 6af667a0a4..2367412004 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -96,17 +96,15 @@ const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_
     "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi",
     "gal_flag", "nstep", "diagnose_ca", "diagnose", "ksing")
 const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc", "use_galgrid")
-
-# Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing
-# instead of throwing an unknown-keyword error; warn so the removal is not silent.
-function _drop_deprecated_keys!(table, deprecated_keys, section::String)
-    for k in deprecated_keys
-        if haskey(table, k)
-            @warn "`$k` in [$section] is deprecated and ignored; please remove it from gpec.toml."
-            delete!(table, k)
-        end
-    end
-    return table
+const _DEPRECATED_KF_KEYS = ("moment", "eq_out", "theta_out", "xlmda_out", "eqpsi_out", "save_records", "gpec_file",
+    "fnml_flag", "ellip_flag", "diag_flag", "force_xialpha")
+const _DEPRECATED_PE_KEYS = ("fixed_boundary", "output_eigenmodes", "filter_modes", "singular_point_method")
+
+# Stop the run when a parsed gpec.toml section sets keys that no longer exist, so the deck is updated.
+function _stop_on_deprecated_keys(table, deprecated_keys, section::String)
+    found = [k for k in deprecated_keys if haskey(table, k)]
+    isempty(found) || error("[$section] sets $(join(("`$k`" for k in found), ", ")), which no longer exist; remove them from gpec.toml")
+    return nothing
 end
 
 function main(args::Vector{String}=String[]; dd::Union{IMASdd.dd,Nothing}=nothing)
@@ -146,7 +144,7 @@ function build_inputs_from_toml(path::String; dd::Union{IMASdd.dd,Nothing}=nothi
     inputs = TOML.parsefile(joinpath(path, "gpec.toml"))
 
     haskey(inputs, "Equilibrium") || error("No [Equilibrium] section in gpec.toml")
-    _drop_deprecated_keys!(inputs["Equilibrium"], _DEPRECATED_EQUIL_KEYS, "Equilibrium")
+    _stop_on_deprecated_keys(inputs["Equilibrium"], _DEPRECATED_EQUIL_KEYS, "Equilibrium")
     eq_config = Equilibrium.EquilibriumConfig(inputs["Equilibrium"], path)
 
     # An equilibrium is analytic (sol/lar/tj, parameters from its embedded section),
@@ -212,7 +210,9 @@ function main_from_inputs(
     # Build data structures from inputs
     intr = ForceFreeStatesInternal(; dir_path=path)
     ffs_table = inputs["ForceFreeStates"]
-    _drop_deprecated_keys!(ffs_table, _DEPRECATED_FFS_KEYS, "ForceFreeStates")
+    _stop_on_deprecated_keys(ffs_table, _DEPRECATED_FFS_KEYS, "ForceFreeStates")
+    haskey(inputs, "KineticForces") && _stop_on_deprecated_keys(inputs["KineticForces"], _DEPRECATED_KF_KEYS, "KineticForces")
+    haskey(inputs, "PerturbedEquilibrium") && _stop_on_deprecated_keys(inputs["PerturbedEquilibrium"], _DEPRECATED_PE_KEYS, "PerturbedEquilibrium")
     ctrl = ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in ffs_table)...)
 
     resolve_mode_space!(intr, ctrl)
diff --git a/src/KineticForces/BounceAveraging.jl b/src/KineticForces/BounceAveraging.jl
index 0f974a09c3..f824ee0024 100644
--- a/src/KineticForces/BounceAveraging.jl
+++ b/src/KineticForces/BounceAveraging.jl
@@ -34,7 +34,6 @@ Produced by `compute_bounce_data()`, consumed by pitch integration.
 struct BounceData
     nlmda::Int
     lambda::Vector{Float64}           # pitch angle grid points
-    dlambda::Vector{Float64}          # weights (dx/dnorm from powspace)
     sigma::Vector{Int}                # 0=trapped, 1=passing
     wb::Vector{Float64}               # bounce frequency ωb(λ) [rad/s]
     wd::Vector{Float64}               # precession drift ωd(λ) [rad/s]
@@ -243,7 +242,7 @@ function compute_bounce_data(
 
     # Build λ grid based on method variant (Fortran lines 569-585)
     method_char = method[1]
-    lambda, dlambda = _build_lambda_grid(method_char, lmdatpb, lmdamax, nlmda)
+    lambda = _build_lambda_grid(method_char, lmdatpb, lmdamax, nlmda)
 
     # Pre-allocate output arrays
     wb_arr = zeros(Float64, nlmda)
@@ -290,7 +289,7 @@ function compute_bounce_data(
         end
     end
 
-    return BounceData(nlmda, lambda, dlambda, sigma_arr,
+    return BounceData(nlmda, lambda, sigma_arr,
         wb_arr, wd_arr, dJdJ_arr, wmats_arr)
 end
 
@@ -416,21 +415,20 @@ end
 # ============================================================================
 
 """
-Build λ grid based on method character (f/t/p).
-Returns (lambda, dlambda) with endpoints excluded.
+Build λ grid based on method character (f/t/p), with endpoints excluded.
 """
 function _build_lambda_grid(method_char::Char, lmdatpb::Float64, lmdamax::Float64, nlmda::Int)
     lmdamin = 0.0
 
     if method_char == 't'
         # Trapped only: λ ∈ (lmdatpb, lmdamax), exclude endpoints
-        pts_inc, wts_inc = powspace(lmdatpb, lmdamax, 1, 2 + nlmda, "both")
-        return pts_inc[2:end-1], wts_inc[2:end-1]
+        pts_inc, _ = powspace(lmdatpb, lmdamax, 1, 2 + nlmda, "both")
+        return pts_inc[2:end-1]
 
     elseif method_char == 'p'
         # Passing only: λ ∈ (lmdamin, lmdatpb), exclude endpoints
-        pts_inc, wts_inc = powspace(lmdamin, lmdatpb, 1, 2 + nlmda, "both")
-        return pts_inc[2:end-1], wts_inc[2:end-1]
+        pts_inc, _ = powspace(lmdamin, lmdatpb, 1, 2 + nlmda, "both")
+        return pts_inc[2:end-1]
 
     else  # 'f' = full
         if lmdatpb ≈ lmdamax
@@ -439,12 +437,10 @@ function _build_lambda_grid(method_char::Char, lmdatpb::Float64, lmdamax::Float6
         # Passing half with refinement near upper boundary
         nhalf_p = nlmda ÷ 2
         nhalf_t = nlmda - nhalf_p
-        pts_p, wts_p = powspace(lmdamin, lmdatpb, 2, 2 + nhalf_p, "upper")
-        pts_t, wts_t = powspace(lmdatpb, lmdamax, 2, 2 + nhalf_t, "lower")
+        pts_p, _ = powspace(lmdamin, lmdatpb, 2, 2 + nhalf_p, "upper")
+        pts_t, _ = powspace(lmdatpb, lmdamax, 2, 2 + nhalf_t, "lower")
         # Exclude boundary points
-        lambda = vcat(pts_p[2:end-1], pts_t[2:end-1])
-        dlambda = vcat(wts_p[2:end-1], wts_t[2:end-1])
-        return lambda, dlambda
+        return vcat(pts_p[2:end-1], pts_t[2:end-1])
     end
 end
 
diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl
index eaa35bae7a..86007bea63 100644
--- a/src/KineticForces/CalculatedKineticMatrices.jl
+++ b/src/KineticForces/CalculatedKineticMatrices.jl
@@ -26,9 +26,7 @@ This routine reads equilibrium-derived profiles (q, dV/dψ, ⟨r⟩, ⟨R⟩) di
 from named splines on `equil.profiles` and `equil.geometry`, and kinetic
 profiles (n, T, ω_E, ν) from the `kinetic_profiles` argument, avoiding the
 former shadow-copy pattern in KineticForcesInternal. The perturbation-mode
-interpolants (`kf_intr.dbob_m`, `kf_intr.divx_m`) remain unwired and are
-tracked as follow-up work blocked on PR #196 — see the plan's "Out of scope"
-section.
+interpolants (`kf_intr.dbob_m`, `kf_intr.divx_m`) are not wired into this routine.
 
 # Arguments
 
@@ -42,7 +40,7 @@ section.
 # Keyword arguments
 
   - `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to
-    carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the
+    carry NTV-specific knobs (nl, zi, mi, wdfac, electron) that the
     KineticForces kernel needs but ForceFreeStatesControl does not expose.
   - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic-
     profile splines loaded via `Equilibrium.load_kinetic_profiles`.
@@ -153,7 +151,7 @@ function compute_calculated_kinetic_matrices(
                     fill!(block_t, 0)
                     compute_kinetic_matrices_at_psi!(
                         block_w, block_t, psi, n, ell,
-                        z_s, m_s, kf_ctrl.magnetic_drift_factor, kf_ctrl.divxi_factor,
+                        z_s, m_s, kf_ctrl.magnetic_drift_factor,
                         el_s, equil, intr_t, prof_s;
                         nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.collisionality_factor,
                         atol_xlmda=kf_ctrl.atol_pitch, rtol_xlmda=kf_ctrl.rtol_pitch,
diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl
index b110b67705..fdfd2325fc 100644
--- a/src/KineticForces/Compute.jl
+++ b/src/KineticForces/Compute.jl
@@ -283,7 +283,6 @@ function combine_species_states(states::AbstractVector{KineticForcesState})
             psi_grid=grid, dtdpsi=dtdpsi, t_cumulative=tcum,
             psi_nsteps=sum(r.psi_nsteps for r in results))
     end
-    combined.completed = true
     return combined
 end
 
@@ -407,6 +406,4 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor
             println("---------------------------------------------")
         end
     end
-
-    state.completed = true
 end
diff --git a/src/KineticForces/EnergyIntegration.jl b/src/KineticForces/EnergyIntegration.jl
index 191eb4bfb7..da95fc1cd2 100644
--- a/src/KineticForces/EnergyIntegration.jl
+++ b/src/KineticForces/EnergyIntegration.jl
@@ -31,7 +31,6 @@ struct EnergyParams
     nutype::String    # collision operator type
     f0type::String    # distribution function type
     nufac::Float64    # collisionality scaling factor
-    ximag::Float64    # deprecated: imaginary contour offset, no longer used
     qt::Bool          # heat flux calculation flag
 end
 
@@ -287,7 +286,7 @@ end
 """
     integrate_energy(wn, wt, we, wd, wb, nuk, ell, leff, n, psi, lambda, method;
                      nutype="harmonic", f0type="maxwellian", nufac=1.0,
-                     ximag=0.0, qt=false, atol=1e-7, rtol=1e-5) → ComplexF64
+                     qt=false, atol=1e-7, rtol=1e-5) → ComplexF64
 
 Integrate the kinetic resonance operator over normalized energy x = E/T.
 
@@ -304,9 +303,6 @@ its real-axis pole resolved analytically (see `_integrate_energy_resonant`).
 Collision operator types (`nutype`): `"zero"`, `"small"`, `"krook"`, `"harmonic"`.
 Distribution function types (`f0type`): `"maxwellian"`, `"park"`, `"cgl"`.
 
-`ximag` is accepted for backward compatibility but no longer used — resonance
-poles are now handled analytically rather than by contour deformation.
-
 # Returns
 - `ComplexF64`: energy integral value
 """
@@ -314,11 +310,11 @@ function integrate_energy(wn::Float64, wt::Float64, we::Float64, wd::Float64,
                               wb::Float64, nuk::Float64, ell::Int, leff::Float64,
                               n::Int, psi::Float64, lambda::Float64, method::String;
                               nutype::String="harmonic", f0type::String="maxwellian",
-                              nufac::Float64=1.0, ximag::Float64=0.0, qt::Bool=false,
+                              nufac::Float64=1.0, qt::Bool=false,
                               atol::Float64=1e-7, rtol::Float64=1e-5, segbuf=nothing)::ComplexF64
 
     p = EnergyParams(wn, wt, we, wd, wb, nuk, leff, n,
-                     nutype, f0type, nufac, ximag, qt)
+                     nutype, f0type, nufac, qt)
 
     # CGL has no resonance denominator and no pole — integrate the physical
     # x-space integrand directly over the half line (QuadGK maps [0,∞) itself).
@@ -335,7 +331,7 @@ end
 """
     evaluate_energy_integrand(x_grid; wn, wt, we, wd, wb, nuk, leff, n,
                                nutype="harmonic", f0type="maxwellian",
-                               nufac=1.0, ximag=0.0, qt=false) → Vector{ComplexF64}
+                               nufac=1.0, qt=false) → Vector{ComplexF64}
 
 Diagnostic convenience: evaluate the physical x-space energy integrand
 N(x)·exp(-x)/denom(x) at specified x = E/T values. Returns the integrand value
@@ -355,7 +351,7 @@ function evaluate_energy_integrand(x_grid::AbstractVector{Float64};
                                     wd::Float64, wb::Float64, nuk::Float64,
                                     leff::Float64, n::Int,
                                     nutype::String="harmonic", f0type::String="maxwellian",
-                                    nufac::Float64=1.0, ximag::Float64=0.0, qt::Bool=false)
-    p = EnergyParams(wn, wt, we, wd, wb, nuk, leff, n, nutype, f0type, nufac, ximag, qt)
+                                    nufac::Float64=1.0, qt::Bool=false)
+    p = EnergyParams(wn, wt, we, wd, wb, nuk, leff, n, nutype, f0type, nufac, qt)
     return [energy_integrand_scalar(x, p) for x in x_grid]
 end
diff --git a/src/KineticForces/KineticForces.jl b/src/KineticForces/KineticForces.jl
index cd27590edb..183297ada8 100644
--- a/src/KineticForces/KineticForces.jl
+++ b/src/KineticForces/KineticForces.jl
@@ -66,7 +66,6 @@ const twopi = 2π
 const μ₀    = 4e-7 * π
 const rad2deg = 180 / π
 const deg2rad = π / 180
-const iunit = 1im               # equivalent to Fortran's (0,1)
 
 # Treat-as-zero threshold for singular denominators and degenerate-coefficient
 # checks. Used by the resonance solver, the Sokhotski-Plemelj pole shift, and
diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl
index 9014dd9a4a..068fb3ee1b 100644
--- a/src/KineticForces/KineticForcesStructs.jl
+++ b/src/KineticForces/KineticForcesStructs.jl
@@ -95,9 +95,6 @@ builds a second control for its differing tolerance).
     equilibrium and kinetic profiles; there are no tuning parameters. Default `true`.
 """
 @kwdef struct KineticForcesControl
-    # Moment type
-    moment::String = "pressure"     # "heat" or "pressure"
-
     # Method flags (which methods to calculate)
     fgar_flag::Bool = true          # Full general-aspect-ratio
     tgar_flag::Bool = false         # Trapped particle GAR
@@ -163,26 +160,12 @@ builds a second control for its differing tolerance).
     # Diagnostic parameters
     psilims::Vector{Float64} = [0.0, 1.0]  # Integration limits in psi
 
-    # Diagnostic output flags
-    eq_out::Bool = false            # Output equilibrium profiles
-    theta_out::Bool = false         # Output theta-dependent quantities
-    xlmda_out::Bool = false         # Output pitch angle profiles
-    eqpsi_out::Bool = false         # Output equilibrium psi profiles
-
     # Output configuration
     write_outputs_to_HDF5::Bool = true
     HDF5_filename::String = "gpec.h5"
-    save_records::Bool = false      # Save detailed integration trajectories
 
     # Input files (relative to dir_path)
     kinetic_file::String = "kinetic.dat"
-    gpec_file::String = "gpec.h5"
-
-    # Diagnostic flags
-    fnml_flag::Bool = false         # Fourier mode coupling diagnostics
-    ellip_flag::Bool = false        # Elliptic integral diagnostics
-    diag_flag::Bool = false         # General diagnostics
-    force_xialpha::Bool = false     # Force xi_alpha formulation
 
     # Debugging
     verbose::Bool = false
@@ -198,7 +181,7 @@ end
 # `cubic_interp` call so the 7-parameter `CubicSeriesInterpolant{...}` need not be
 # spelled out. Used to concretely type the smats..zmats fields below so the per-ψ
 # `intr.smats(psi)` reads in `Torque.jl` are statically dispatched, not boxed
-# `Any` calls (issue #247).
+# `Any` calls.
 const GeomInterp = typeof(cubic_interp(collect(0.0:0.25:1.0), Series(zeros(ComplexF64, 5, 1)); extrap=ExtendExtrap()))
 
 """
@@ -280,14 +263,8 @@ this struct.
     # Clebsch displacement vectors for mode-coupled dW contraction
     xs_m::Any = nothing            # Vector of 3 CubicSeriesInterpolants: [ξ_ψ, ξ_+, ξ_-]
 
-    # Integration results
-    tphi::ComplexF64 = 0.0 + 0.0im    # Total torque/energy
-    tsurf::ComplexF64 = 0.0 + 0.0im   # Surface torque/energy
-    teq::ComplexF64 = 0.0 + 0.0im     # Equilibrium grid result
-
     # Method tracking — available methods/docs come from `METHOD_REGISTRY`.
     method::String = ""
-    wtw::Array{ComplexF64,3} = Array{ComplexF64}(undef, 0, 0, 0)
 
     verbose::Bool = false
 end
@@ -485,24 +462,6 @@ end
 # Result Structs
 # ============================================================================
 
-"""
-    EnergyIntegrationResult
-
-Results from a single energy-space integration at one (ψ, λ, ℓ) point.
-Trajectory fields are only populated when `ctrl.save_records=true`.
-"""
-@kwdef struct EnergyIntegrationResult
-    psi::Float64 = 0.0
-    lambda::Float64 = 0.0
-    ell::Int = 0
-    leff::Float64 = 0.0
-    torque::ComplexF64 = 0.0 + 0.0im
-    kinetic_energy::ComplexF64 = 0.0 + 0.0im
-    x_trajectory::Vector{Float64} = Float64[]
-    integrand_trajectory::Vector{ComplexF64} = ComplexF64[]
-    integral_trajectory::Vector{ComplexF64} = ComplexF64[]
-end
-
 """
     MethodResult
 
@@ -514,7 +473,6 @@ Results for one NTV computation method across all flux surfaces.
     torque_profile::Any = nothing     # Interpolant of dT/dψ(ψ) from ψ integration
     total_torque::ComplexF64 = 0.0 + 0.0im  # complex T: Re = torque T_φ [N·m], Im = 2n·δW_k (Logan 2013 Eq. 19)
     total_energy::Float64 = 0.0              # real kinetic energy δW_k = Im(T)/(2n) [J]
-    records::Vector{EnergyIntegrationResult} = EnergyIntegrationResult[]
     # Per-step ψ profile from outer quadrature
     psi_grid::Vector{Float64} = Float64[]
     dtdpsi::Vector{ComplexF64} = ComplexF64[]
@@ -535,5 +493,4 @@ Written to gpec.h5 under the "KineticForces" group.
     method_results::Dict{String,MethodResult} = Dict{String,MethodResult}()
     # Block-diagonal kinetic matrices: key=method, value=(numpert_total, numpert_total, 6)
     kinetic_matrices::Dict{String,Array{ComplexF64,3}} = Dict{String,Array{ComplexF64,3}}()
-    completed::Bool = false
 end
diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl
index ffd6d67c45..d19a3cf668 100644
--- a/src/KineticForces/Output.jl
+++ b/src/KineticForces/Output.jl
@@ -71,10 +71,6 @@ function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spl
                 mg["dVdpsi"] = [dVdpsi_spline(p) for p in result.psi_grid]
             end
         end
-
-        if !isempty(result.records)
-            write_integration_records!(mg, result.records)
-        end
     end
 
     # Write the six drift-kinetic coefficient matrices if present
@@ -116,16 +112,6 @@ const KF_METHOD_H5_ANNOTATIONS = [
     "dTdpsi" => (; long_name="complex torque density dT/dψ_N = dT_φ/dψ_N + 2i·n·dδW_k/dψ_N; divide by dVdpsi for dT/dV", units="N*m", dims=("psi",), attach=(1 => "psi",)),
     "T" => (; long_name="cumulative complex torque T(ψ_N) = T_φ + 2i·n·δW_k (trapezoidal)", units="N*m", dims=("psi",), attach=(1 => "psi",)),
     "dVdpsi" => (; long_name="flux-surface volume derivative dV/dψ_N at quadrature points", units="m^3", dims=("psi",), attach=(1 => "psi",)),
-    "EnergyIntegrals/psi" => (; long_name="ψ_N of each energy-integration record"),
-    "EnergyIntegrals/lambda" => (; long_name="pitch λ = μB0/E of each record"),
-    "EnergyIntegrals/ell" => (; long_name="bounce harmonic ℓ of each record"),
-    "EnergyIntegrals/leff" => (; long_name="effective bounce harmonic ℓ_eff = ℓ + n·q (circulating) or ℓ (trapped)"),
-    "EnergyIntegrals/torque" => (; long_name="complex torque contribution of the record", units="N*m"),
-    "EnergyIntegrals/kinetic_energy" => (; long_name="complex kinetic energy contribution of the record", units="J"),
-    "EnergyIntegrals/trajectory_offsets" => (; long_name="ragged-array offsets: record k spans offsets[k]+1:offsets[k+1] of the *_all arrays"),
-    "EnergyIntegrals/x_all" => (; long_name="normalized energy x = E/T abscissae of all integration trajectories (concatenated)"),
-    "EnergyIntegrals/integrand_all" => (; long_name="complex energy-space torque integrand along all trajectories (concatenated)"),
-    "EnergyIntegrals/integral_all" => (; long_name="complex cumulative energy-space integral along all trajectories (concatenated)"),
     "KineticMatrices/A" => (; long_name="Logan 2015 drift-kinetic coefficient matrix A = W_Z†W_Z (energy normalization)", dims=("mode_row", "mode_col")),
     "KineticMatrices/B" => (; long_name="Logan 2015 drift-kinetic coefficient matrix B = W_Z†W_X (energy normalization)", dims=("mode_row", "mode_col")),
     "KineticMatrices/C" => (; long_name="Logan 2015 drift-kinetic coefficient matrix C = W_Z†W_Y (energy normalization)", dims=("mode_row", "mode_col")),
@@ -134,63 +120,6 @@ const KF_METHOD_H5_ANNOTATIONS = [
     "KineticMatrices/H" => (; long_name="Logan 2015 drift-kinetic coefficient matrix H = W_Y†W_Y (energy normalization)", dims=("mode_row", "mode_col"))
 ]
 
-"""
-    write_integration_records!(mg::HDF5.Group, records::Vector{EnergyIntegrationResult})
-
-Write variable-length integration trajectory records using offset-indexed concatenated arrays.
-This is the standard HDF5 ragged array pattern for storing variable-length data.
-
-# Arguments
-
-  - `mg::HDF5.Group`: HDF5 group for this method
-  - `records::Vector{EnergyIntegrationResult}`: Integration records to write
-"""
-function write_integration_records!(mg::HDF5.Group, records::Vector{EnergyIntegrationResult})
-    rg = create_group(mg, "EnergyIntegrals")
-
-    # Scalar fields per record
-    rg["psi"] = [r.psi for r in records]
-    rg["lambda"] = [r.lambda for r in records]
-    rg["ell"] = [r.ell for r in records]
-    rg["leff"] = [r.leff for r in records]
-    rg["torque"] = [r.torque for r in records]
-    rg["kinetic_energy"] = [r.kinetic_energy for r in records]
-
-    # Variable-length trajectories: concatenate all, store offsets for indexing
-    lengths = [length(r.x_trajectory) for r in records]
-    offsets = cumsum([0; lengths])
-    rg["trajectory_offsets"] = offsets
-
-    if sum(lengths) > 0
-        rg["x_all"] = vcat([r.x_trajectory for r in records]...)
-        rg["integrand_all"] = vcat([r.integrand_trajectory for r in records]...)
-        rg["integral_all"] = vcat([r.integral_trajectory for r in records]...)
-    end
-end
-
-
-"""
-    print_summary(state::KineticForcesState; verbose::Bool=false)
-
-Print a summary of KineticForces results to stdout.
-
-# Arguments
-
-  - `state::KineticForcesState`: Accumulated computation results
-  - `verbose::Bool`: Print detailed per-surface results
-"""
-function print_summary(state::KineticForcesState; verbose::Bool=false)
-    for (method_name, result) in state.method_results
-        @printf("%-8s  T_phi = %11.3e   2n*dW_k = %11.3e\n",
-            method_name, real(result.total_torque), imag(result.total_torque))
-    end
-    if verbose
-        for (method_name, _) in state.kinetic_matrices
-            println("  Kinetic matrices stored for: $method_name")
-        end
-    end
-end
-
 """
     write_validity!(h5file, kf_ctrl, species, kinetic_profiles, equil)
 
diff --git a/src/KineticForces/PitchIntegration.jl b/src/KineticForces/PitchIntegration.jl
index 7b80954d39..f0ae9a4052 100644
--- a/src/KineticForces/PitchIntegration.jl
+++ b/src/KineticForces/PitchIntegration.jl
@@ -37,7 +37,6 @@ struct PitchGARParams{F}
     nutype::String     # collision operator type
     f0type::String     # distribution function type
     nufac::Float64     # collisionality scaling factor
-    ximag::Float64     # imaginary contour offset
     qt::Bool           # heat flux flag
     energy_atol::Float64
     energy_rtol::Float64
@@ -78,13 +77,13 @@ function integrate_pitch_gar_quadgk(
     nqty::Int, ell::Int, n::Int,
     rex::Float64, imx::Float64, psi::Float64, method::String;
     nutype::String="harmonic", f0type::String="maxwellian",
-    nufac::Float64=1.0, ximag::Float64=0.0, qt::Bool=false,
+    nufac::Float64=1.0, qt::Bool=false,
     energy_atol::Float64=1e-9, energy_rtol::Float64=1e-6,
     pitch_atol::Float64=1e-9, pitch_rtol::Float64=1e-6
 )
     params = PitchGARParams(
         wn, wt, we, nuk, bobmax, epsr, q, ell, n, psi, method,
-        nutype, f0type, nufac, ximag, qt,
+        nutype, f0type, nufac, qt,
         energy_atol, energy_rtol,
         rex, imx, nqty, fbnce, fbnce_norm, Ref(1),
         Vector{ComplexF64}(undef, 2 + nqty),
@@ -127,19 +126,19 @@ function _pitch_gar_kernel_quadgk!(out::Vector{ComplexF64}, lambda, p::PitchGARP
         xint_co = integrate_energy(p.wn, p.wt, p.we, wd, wb, nueff,
                                         p.ell, leff, p.n, p.psi, lambda, p.method;
                                         nutype=p.nutype, f0type=p.f0type,
-                                        nufac=p.nufac, ximag=p.ximag, qt=p.qt,
+                                        nufac=p.nufac, qt=p.qt,
                                         atol=p.energy_atol, rtol=p.energy_rtol, segbuf=p.esegbuf)
         xint_counter = integrate_energy(p.wn, p.wt, p.we, wd, -wb, nueff,
                                              p.ell, leff, p.n, p.psi, lambda, p.method;
                                              nutype=p.nutype, f0type=p.f0type,
-                                             nufac=p.nufac, ximag=p.ximag, qt=p.qt,
+                                             nufac=p.nufac, qt=p.qt,
                                              atol=p.energy_atol, rtol=p.energy_rtol, segbuf=p.esegbuf)
         xint = xint_co + xint_counter
     else
         xint = integrate_energy(p.wn, p.wt, p.we, wd, wb, nueff,
                                      p.ell, leff, p.n, p.psi, lambda, p.method;
                                      nutype=p.nutype, f0type=p.f0type,
-                                     nufac=p.nufac, ximag=p.ximag, qt=p.qt,
+                                     nufac=p.nufac, qt=p.qt,
                                      atol=p.energy_atol, rtol=p.energy_rtol)
     end
 
@@ -174,14 +173,14 @@ function integrate_pitch_gar_quadgk_wt(
     nqty::Int, ell::Int, n::Int,
     psi::Float64, method::String;
     nutype::String="harmonic", f0type::String="maxwellian",
-    nufac::Float64=1.0, ximag::Float64=0.0, qt::Bool=false,
+    nufac::Float64=1.0, qt::Bool=false,
     energy_atol::Float64=1e-9, energy_rtol::Float64=1e-6,
     pitch_atol::Float64=1e-9, pitch_rtol::Float64=1e-6
 )
     # rex/imx unused in the dual kernel; carry 1.0 placeholders for the struct.
     params = PitchGARParams(
         wn, wt, we, nuk, bobmax, epsr, q, ell, n, psi, method,
-        nutype, f0type, nufac, ximag, qt,
+        nutype, f0type, nufac, qt,
         energy_atol, energy_rtol,
         1.0, 1.0, nqty, fbnce, fbnce_norm, Ref(1),
         Vector{ComplexF64}(undef, 2 + nqty),
@@ -225,19 +224,19 @@ function _pitch_gar_kernel_quadgk_wt!(out::Vector{ComplexF64}, lambda, p::PitchG
         xint_co = integrate_energy(p.wn, p.wt, p.we, wd, wb, nueff,
                                         p.ell, leff, p.n, p.psi, lambda, p.method;
                                         nutype=p.nutype, f0type=p.f0type,
-                                        nufac=p.nufac, ximag=p.ximag, qt=p.qt,
+                                        nufac=p.nufac, qt=p.qt,
                                         atol=p.energy_atol, rtol=p.energy_rtol, segbuf=p.esegbuf)
         xint_counter = integrate_energy(p.wn, p.wt, p.we, wd, -wb, nueff,
                                              p.ell, leff, p.n, p.psi, lambda, p.method;
                                              nutype=p.nutype, f0type=p.f0type,
-                                             nufac=p.nufac, ximag=p.ximag, qt=p.qt,
+                                             nufac=p.nufac, qt=p.qt,
                                              atol=p.energy_atol, rtol=p.energy_rtol, segbuf=p.esegbuf)
         xint = xint_co + xint_counter
     else
         xint = integrate_energy(p.wn, p.wt, p.we, wd, wb, nueff,
                                      p.ell, leff, p.n, p.psi, lambda, p.method;
                                      nutype=p.nutype, f0type=p.f0type,
-                                     nufac=p.nufac, ximag=p.ximag, qt=p.qt,
+                                     nufac=p.nufac, qt=p.qt,
                                      atol=p.energy_atol, rtol=p.energy_rtol)
     end
 
diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl
index 6c8c24eefa..f1478fe454 100644
--- a/src/KineticForces/Torque.jl
+++ b/src/KineticForces/Torque.jl
@@ -417,7 +417,7 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo,
                        ymat=nothing, zmat=nothing,
                        nlmda::Int=128, ntheta::Int=128,
                        nutype::String="harmonic", f0type::String="maxwellian",
-                       nufac::Float64=1.0, ximag::Float64=0.0, qt::Bool=false,
+                       nufac::Float64=1.0, qt::Bool=false,
                        energy_atol::Float64=1e-7, energy_rtol::Float64=1e-5,
                        pitch_atol::Float64=1e-9, pitch_rtol::Float64=1e-6,
                        rex_override::Union{Nothing,Float64}=nothing,
@@ -490,7 +490,7 @@ function calculate_gar(psi, n, l, q, epsr, wdian, wdiat, welec, nuk, bo,
     lxint = integrate_pitch_gar_quadgk(
         wdian, wdiat, welec, nuk, bo / bmax, epsr, q,
         fbnce, fbnce_norm, nqty, l, n, rex, imx, psi, method;
-        nutype, f0type, nufac, ximag, qt,
+        nutype, f0type, nufac, qt,
         energy_atol, energy_rtol, pitch_atol, pitch_rtol)
 
     # Scalar torque, Eq. (19) of Logan et al., Phys. Plasmas 20, 122507 (2013).
@@ -738,7 +738,7 @@ function kinetic_energy_matrices_for_euler_lagrange!(
     intr::KineticForcesInternal;
     nlmda::Int=128, ntheta::Int=128,
     nutype::String="harmonic", f0type::String="maxwellian",
-    nufac::Float64=1.0, ximag::Float64=0.0,
+    nufac::Float64=1.0,
     energy_atol::Float64=1e-9, energy_rtol::Float64=1e-6,
     pitch_atol::Float64=1e-9, pitch_rtol::Float64=1e-6
 )
@@ -804,7 +804,7 @@ function kinetic_energy_matrices_for_euler_lagrange!(
         state.wdian, state.wdiat, state.welec, state.nuk,
         bo / state.bmax, state.epsr, state.q,
         fbnce, fbnce_norm, nqty, l, n, psi, "fwmm";
-        nutype, f0type, nufac, ximag, qt=false,
+        nutype, f0type, nufac, qt=false,
         energy_atol, energy_rtol, pitch_atol, pitch_rtol)
 
     # Torque → energy normalization, Eq. (19) of Logan et al., Phys. Plasmas 20, 122507 (2013).
@@ -860,7 +860,7 @@ end
 
 """
     compute_kinetic_matrices_at_psi!(kwmat, ktmat, psi, n, l, zi, mi,
-        wdfac, divxfac, electron, equil, intr, kinetic_profiles)
+        wdfac, electron, equil, intr, kinetic_profiles)
 
 Compute the six kinetic Euler-Lagrange coefficient matrices at a single
 flux surface and split them into `kwmat` and `ktmat` in the convention
@@ -887,8 +887,7 @@ Fortran `dcon/sing.f:967-1075` for non-Hermitian B_k, C_k, E_k.
 # Arguments
 - `kwmat::Array{ComplexF64,3}`: Output (mpert×mpert×6), fwmm half, zeroed on entry
 - `ktmat::Array{ComplexF64,3}`: Output (mpert×mpert×6), ftmm half, zeroed on entry
-- `psi, n, l, zi, mi, wdfac, divxfac, electron`: Same as `tpsi!` (divxfac unused
-  on the matrix path — retained for call-site compatibility)
+- `psi, n, l, zi, mi, wdfac, electron`: Same as `tpsi!`
 - `equil`: PlasmaEquilibrium
 - `intr::KineticForcesInternal`: Internal state with mode indexing, geometric
   matrices (smats/tmats/xmats/ymats/zmats), and per-surface θ-grid buffers
@@ -900,7 +899,7 @@ function compute_kinetic_matrices_at_psi!(
     kwmat::Array{ComplexF64,3},
     ktmat::Array{ComplexF64,3},
     psi::Float64, n::Int, l::Int,
-    zi::Int, mi::Int, wdfac::Float64, _divxfac::Float64,
+    zi::Int, mi::Int, wdfac::Float64,
     electron::Bool, equil::Equilibrium.PlasmaEquilibrium, intr::KineticForcesInternal,
     kinetic_profiles::Equilibrium.KineticProfileSplines;
     nutype::String="harmonic", f0type::String="maxwellian", nufac::Float64=1.0,
diff --git a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
index 50b23548c2..9df339eac2 100644
--- a/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
+++ b/src/PerturbedEquilibrium/PerturbedEquilibrium.jl
@@ -17,7 +17,6 @@ import ..ForcingTerms
 import ..ForcingTerms: ForcingMode, CoilSet, load_forcing_data!, convert_forcing_normalization!
 import ..Utilities
 import ..Utilities.FourierTransforms
-import DelimitedFiles: readdlm
 
 # Include module files
 include("PerturbedEquilibriumStructs.jl")
@@ -84,7 +83,7 @@ function compute_perturbed_equilibrium(
               "does not publish yet; set compute_response = false and compute_singular_coupling = false, or use a forward solve")
     end
 
-    # Step 0: Initialize mode arrays for convenient indexing
+    # Mode arrays for convenient indexing
     initialize_mode_arrays!(intr, ffs)
 
     # A run has at most one ξ solution, already closed at the rationals and carrying populated
@@ -103,14 +102,14 @@ function compute_perturbed_equilibrium(
         isempty(coil_sets) || (intr.coil_sets = coil_sets)
     end
 
-    # Step 2: Compute plasma response
+    # Plasma response
     if ctrl.compute_response &&
        ForceFreeStates.require(ffs, :free_boundary, "plasma response calculation") &&
        ForceFreeStates.require_solution(ffs, "plasma response calculation")
         compute_plasma_response!(state, equil, solution, ffs.free_boundary.wt0, mtheta_vacuum, ffs, intr, ctrl, ffs.metric, mats)
     end
 
-    # Step 3: Compute singular coupling metrics
+    # Singular coupling metrics
     if ctrl.compute_singular_coupling &&
        ForceFreeStates.require(ffs, :free_boundary, "singular coupling calculation") &&
        ForceFreeStates.require_solution(ffs, "singular coupling calculation")
@@ -118,9 +117,6 @@ function compute_perturbed_equilibrium(
         compute_dominant_coupling!(state, ctrl)
     end
 
-    # Step 4: Output eigenmode fields (integrated into HDF5 output)
-    # This is handled by write_outputs_to_HDF5 in main()
-
     return state
 end
 
diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
index 18cd783f31..3b08154394 100644
--- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
+++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
@@ -7,10 +7,8 @@ User-facing control parameters from TOML [PerturbedEquilibrium] section.
 
 Note: Forcing data file settings are now in [ForcingTerms] section.
 
-High Priority (MWE):
+Calculation:
 
-  - `fixed_boundary::Bool` - Fixed boundary flag (default: false)
-  - `output_eigenmodes::Bool` - Output mode fields as b-fields (default: true)
   - `compute_response::Bool` - Compute plasma response (default: true)
   - `compute_singular_coupling::Bool` - Compute singular coupling metrics (default: true)
   - `verbose::Bool` - Enable verbose logging (default: true)
@@ -20,11 +18,6 @@ Output Settings:
   - `output_filename::String` - Combined output file with ForceFreeStates results (default: uses ForceFreeStates HDF5_filename)
   - `write_outputs_to_HDF5::Bool` - Write perturbed equilibrium outputs to HDF5 (default: true)
 
-Medium Priority (defer for MWE):
-
-  - `filter_modes::Bool` - Enable mode filtering (default: false)
-  - `singular_point_method::String` - Method for singular point treatment (default: "standard")
-
 Regularization:
 
   - `regularization_width::Float64` - Regularization width for singular surface smoothing (default: 0.05). Set to 0
@@ -32,9 +25,7 @@ Regularization:
     resonant singularity to smooth.
 """
 @kwdef struct PerturbedEquilibriumControl
-    # High Priority (MWE)
-    fixed_boundary::Bool = false
-    output_eigenmodes::Bool = true
+    # Calculation
     compute_response::Bool = true
     compute_singular_coupling::Bool = true
     verbose::Bool = true
@@ -43,10 +34,6 @@ Regularization:
     output_filename::String = ""  # Empty means use ForceFreeStates filename
     write_outputs_to_HDF5::Bool = true
 
-    # Medium Priority (include but simple for MWE)
-    filter_modes::Bool = false
-    singular_point_method::String = "standard"
-
     # Regularization width for singular surface smoothing (matches Fortran gpec.f reg_spot).
     # Set to 0 to disable regularization. Must be non-negative.
     regularization_width::Float64 = 5e-2
@@ -64,7 +51,6 @@ Internal state variables for perturbed equilibrium calculations.
   - `coil_sets::Vector{CoilSet}` - Coil geometry used (when `forcing_data_format = "coil"`);
     captured for the gpec.h5 rerun snapshot
   - `plasma_response::Matrix{ComplexF64}` - Plasma response matrix
-  - `singular_coupling_metrics::Dict{String,Float64}` - Coupling metrics at singular surfaces
   - `m_modes::Vector{Int}` - Poloidal mode numbers for each index i in 1:numpert_total
   - `n_modes::Vector{Int}` - Toroidal mode numbers for each index i in 1:numpert_total
 """
@@ -73,7 +59,6 @@ Internal state variables for perturbed equilibrium calculations.
     forcing_modes::Vector{ForcingMode} = ForcingMode[]
     coil_sets::Vector{CoilSet} = CoilSet[]
     plasma_response::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)
-    singular_coupling_metrics::Dict{String,Float64} = Dict{String,Float64}()
     m_modes::Vector{Int} = Int[]
     n_modes::Vector{Int} = Int[]
     # ForceFreeStates-provided B_pen per (match surface × coil-drive column) from inner layer.
@@ -139,7 +124,7 @@ representations (all tesla; no flux/weber is stored):
   - `forcing_b_area`/`response_b_area` - area-weighted field b̄ (= S·b̃; flux is Φ = A·b̄)
 
 Control surface matrices [numpert_total × numpert_total], stored in the coordinate-invariant
-root-area-weighted field (b̃) space (issue #233 / Pharr 2026). Writing `S ≡ rootarea_to_area_weight`
+root-area-weighted field (b̃) space (Pharr 2026). Writing `S ≡ rootarea_to_area_weight`
 (b̃→b̄) and `A ≡ surface_area`, the brief internal flux-conform operator is `R = S·A` (Φ = R·b̃).
 Recover the area-weighted (b̄) forms by conforming with `S` (e.g. `L_b̄ = S·L̃·S†`); recover flux
 with `Φ = A·b̄`.
diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl
index 2100e52c1d..120847e73a 100644
--- a/src/PerturbedEquilibrium/Response.jl
+++ b/src/PerturbedEquilibrium/Response.jl
@@ -51,7 +51,7 @@ function compute_plasma_response!(
     intr.plasma_response = permeability
 
     # Conform the control-surface matrices to the coordinate-invariant root-area-weighted
-    # field (b̃) space for output (issue #233 / Pharr 2026). Store the b̃→b̄ operator S = Σ/√A
+    # field (b̃) space for output (Pharr 2026). Store the b̃→b̄ operator S = Σ/√A
     # and the scalar surface area A so users can recover the area-weighted field (b̄ = S·b̃) or
     # flux (Φ = A·b̄) — see Utils.jl output docs.
     rootarea_to_area_weight, surface_area = build_control_surface_rootarea_to_area_weight(equil, ffs)
diff --git a/src/PerturbedEquilibrium/ResponseMatrices.jl b/src/PerturbedEquilibrium/ResponseMatrices.jl
index 405b3bebfa..a9224c2b0d 100644
--- a/src/PerturbedEquilibrium/ResponseMatrices.jl
+++ b/src/PerturbedEquilibrium/ResponseMatrices.jl
@@ -195,10 +195,10 @@ function build_flux_matrix(
     ffs::ForceFreeStatesResult
 )::Matrix{ComplexF64}
 
-    # Step 1: Extract boundary displacements and equilibrium quantities
+    # Boundary displacements and equilibrium quantities
     boundary_data = extract_boundary_displacements(equil, solution, ffs)
 
-    # Step 2: Compute normal magnetic field at plasma boundary
+    # Normal magnetic field at the plasma boundary
     # This is the actual implementation of GPEC's bwp_mn calculation
     # bwp_mn[i,j] = i * (dΨ/dρ) * (m[i] - n*q_boundary) * ξ_ψ[i,j]
     flxmats = compute_normal_magnetic_field(boundary_data, ffs)
diff --git a/src/PerturbedEquilibrium/SingularCoupling.jl b/src/PerturbedEquilibrium/SingularCoupling.jl
index e8c5c9e09c..05a97049d3 100644
--- a/src/PerturbedEquilibrium/SingularCoupling.jl
+++ b/src/PerturbedEquilibrium/SingularCoupling.jl
@@ -425,7 +425,7 @@ function compute_singular_coupling_metrics!(
                 state.C_penetrated_area_weighted_field[row, :] = pen_row
             end
 
-            # LHS normalization audit (#233) — output scalar coordinate-invariance per row:
+            # LHS normalization audit — output scalar coordinate-invariance per row:
             #  - Δ' (1/length): the resonant-surface jump in ∂b^ψ/∂ψ over 2π·χ₁; the tearing index is
             #    coordinate-invariant (its sign/zero-crossing set the stability boundary) [Glasser 2016].
             #  - resonant (shielding) current: j_c already integrates jac·|∇ψ| over the surface, so the
@@ -474,7 +474,7 @@ function compute_singular_coupling_metrics!(
     state.forcing_solution_weights = C_coeffs * forcing_flux
 
     # Conform the stored coupling-matrix input basis to the coordinate-invariant root-area-weighted
-    # field (b̃) space (#233 / Pharr 2026): C̃ = C·R, so each stored row acts on the applied field
+    # field (b̃) space (Pharr 2026): C̃ = C·R, so each stored row acts on the applied field
     # b̃_x (Φ_x = R·b̃_x) and its singular values become coordinate-invariant. The conform operator
     # R = S·A (Σ·√A) is the only place flux briefly appears; it is built from the b̃→b̄ operator S and
     # the scalar surface area A. Done after the applied-vector evaluation above so those physical
diff --git a/src/Rerun.jl b/src/Rerun.jl
index 0723921603..29bca18887 100644
--- a/src/Rerun.jl
+++ b/src/Rerun.jl
@@ -290,11 +290,11 @@ end
 The equilibrium configuration and the input `setup_equilibrium` consumes, rebuilt from a stored
 TOML `inputs` dict and the equilibrium `ingest` read from a `gpec.h5`. Analytic kinds regenerate
 from their TOML section; file-based kinds rebuild splines from the stored ingest. A file-based
-run with no ingest can only come from a pre-ingest `gpec.h5` and is an error. Mutates `inputs`
-only to drop deprecated equilibrium keys.
+run with no ingest can only come from a pre-ingest `gpec.h5` and is an error. Stops on
+deprecated equilibrium keys.
 """
 function rebuild_equilibrium_inputs(inputs::Dict{String,Any}, ingest, output_dir::AbstractString)
-    _drop_deprecated_keys!(inputs["Equilibrium"], _DEPRECATED_EQUIL_KEYS, "Equilibrium")
+    _stop_on_deprecated_keys(inputs["Equilibrium"], _DEPRECATED_EQUIL_KEYS, "Equilibrium")
     # Clear eq_filename on a copy: unused on replay, a stale absolute path could mislead
     # downstream code, and `inputs` itself is re-serialized into the rerun's gpec_toml_raw.
     equil_dict = merge(inputs["Equilibrium"], Dict{String,Any}("eq_filename" => ""))
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index 4bfa677f93..7f776e1f14 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -164,15 +164,20 @@ using Statistics
         @test_throws r"stopped early at t = .* with retcode DtNaN" Equil.tj_analytic_find_nu(Equil.TJAnalyticShapeParams(tj), tj.qa; abstol=NaN)
     end
 
-    @testset "Deprecated TOML keys are dropped, not fatal" begin
-        # Removed control knobs must keep old gpec.toml decks (and older gpec.h5 replays,
-        # whose stored TOML blob goes through the same path) parsing with a warning.
-        ffs_table = Dict{String,Any}("reform_eq_with_psilim" => false, "nn_low" => 1)
-        @test_logs (:warn, r"reform_eq_with_psilim") GeneralizedPerturbedEquilibrium._drop_deprecated_keys!(
-            ffs_table, GeneralizedPerturbedEquilibrium._DEPRECATED_FFS_KEYS, "ForceFreeStates")
-        @test !haskey(ffs_table, "reform_eq_with_psilim")
-        @test GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(;
-            (Symbol(k) => v for (k, v) in ffs_table)...) isa GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl
+    @testset "Deprecated TOML keys stop the run" begin
+        # A listed key no longer exists, so the run stops and names it instead of ignoring it.
+        GPEC = GeneralizedPerturbedEquilibrium
+        @test_throws r"\[ForceFreeStates\] sets `reform_eq_with_psilim`" GPEC._stop_on_deprecated_keys(
+            Dict{String,Any}("reform_eq_with_psilim" => false, "nn_low" => 1), GPEC._DEPRECATED_FFS_KEYS, "ForceFreeStates")
+        @test_throws r"`moment`, `save_records`" GPEC._stop_on_deprecated_keys(
+            Dict{String,Any}("moment" => "pressure", "save_records" => false, "fgar_flag" => true), GPEC._DEPRECATED_KF_KEYS, "KineticForces")
+        @test_throws r"\[PerturbedEquilibrium\] sets `fixed_boundary`" GPEC._stop_on_deprecated_keys(
+            Dict{String,Any}("fixed_boundary" => false, "compute_response" => false), GPEC._DEPRECATED_PE_KEYS, "PerturbedEquilibrium")
+        @test GPEC._stop_on_deprecated_keys(Dict{String,Any}("nn_low" => 1), GPEC._DEPRECATED_FFS_KEYS, "ForceFreeStates") === nothing
+        # A listed deck key must not still be a control field, or listing it would stop runs that set a live key.
+        @test !any(k -> hasfield(GPEC.ForceFreeStates.ForceFreeStatesControl, Symbol(k)), GPEC._DEPRECATED_FFS_KEYS)
+        @test !any(k -> hasfield(GPEC.KineticForces.KineticForcesControl, Symbol(k)), GPEC._DEPRECATED_KF_KEYS)
+        @test !any(k -> hasfield(GPEC.PerturbedEquilibrium.PerturbedEquilibriumControl, Symbol(k)), GPEC._DEPRECATED_PE_KEYS)
     end
 
     @testset "Inverse grid types (original, rho)" begin
diff --git a/test/runtests_h5_schema.jl b/test/runtests_h5_schema.jl
index c743d6a7d6..edf256dae7 100644
--- a/test/runtests_h5_schema.jl
+++ b/test/runtests_h5_schema.jl
@@ -15,7 +15,7 @@ function _group_name_ok(parent_path, name)
     parent_path == "KineticForces" && return true
     # Multi-ion runs nest KineticForces/PerSpecies/<species>/<method>/: the species label
     # (ion_z1_m2, impurity_z6_m12, electron) and the method token are both data-driven.
-    # Exactly those two levels — anything deeper (EnergyIntegrals/, KineticMatrices/) is
+    # Exactly those two levels — anything deeper (KineticMatrices/) is
     # a normal group and still has to be CamelCase.
     occursin(r"^KineticForces/PerSpecies(/[^/]+)?$", parent_path) && return true
     occursin(r"^Surface_\d+$", name) && return true
diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl
index 4d141f6269..ecb58bd0aa 100644
--- a/test/runtests_kinetic.jl
+++ b/test/runtests_kinetic.jl
@@ -94,7 +94,7 @@
             # In CGL mode, fx = cx^2.5 * exp(-cx) / (i*n) — no resonance denominator
             p = KF.EnergyParams(
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1,
-                "zero", "cgl", 1.0, 0.0, false
+                "zero", "cgl", 1.0, false
             )
             x = 1.0
             val = KF.energy_integrand_scalar(x, p)
@@ -107,7 +107,7 @@
             # At large x, the integrand should decay exponentially
             p = KF.EnergyParams(
                 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1,
-                "zero", "maxwellian", 1.0, 0.0, false
+                "zero", "maxwellian", 1.0, false
             )
             val_large = KF.energy_integrand_scalar(50.0, p)
             # At x=50, exp(-50) ≈ 1.9e-22, so integrand should be negligible
@@ -118,11 +118,11 @@
             # With qt=true, integrand is multiplied by (x - 2.5)
             p_no_qt = KF.EnergyParams(
                 1.0, 0.0, 1.0, 0.5, 0.5, 0.0, 1.0, 1,
-                "zero", "maxwellian", 1.0, 0.0, false
+                "zero", "maxwellian", 1.0, false
             )
             p_qt = KF.EnergyParams(
                 1.0, 0.0, 1.0, 0.5, 0.5, 0.0, 1.0, 1,
-                "zero", "maxwellian", 1.0, 0.0, true
+                "zero", "maxwellian", 1.0, true
             )
             x = 3.0
             val_no = KF.energy_integrand_scalar(x, p_no_qt)
@@ -395,7 +395,7 @@
             # CGL mode: ∫₀^∞ x^2.5·exp(-x)/(i·n) dx = Γ(3.5)/(i·1) = -i·(15/8)√π.
             result = KF.integrate_energy(
                 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0, 0.0, 1, 0.5, 0.5, "fcgl";
-                nutype="zero", f0type="cgl", nufac=1.0, ximag=0.0, qt=false,
+                nutype="zero", f0type="cgl", nufac=1.0, qt=false,
                 atol=1e-12, rtol=1e-10
             )
             gamma_3_5 = 15/8 * sqrt(π)
@@ -418,7 +418,7 @@
             # independent reference for the u-substitution + pole-extraction result.
             wn, wt, we, wd, wb, nuk, leff, n = 0.5, 0.8, -2.0, 0.5, 1.0, 0.3, 1.0, 1
             p = KF.EnergyParams(wn, wt, we, wd, wb, nuk, leff, n,
-                                "harmonic", "maxwellian", 1.0, 0.0, false)
+                                "harmonic", "maxwellian", 1.0, false)
             x_res = KF.find_resonance_energies(leff, wb, n, we, wd)
             @test length(x_res) == 1   # this case has exactly one resonance
             reference, _ = KF.quadgk(x -> KF.energy_integrand_scalar(x, p),
@@ -541,7 +541,7 @@
         @test all(isfinite, vals)
         # Element-wise match against the scalar form documented in the docstring.
         p = KF.EnergyParams(0.5, 0.8, -2.0, 0.5, 1.0, 0.3, 1.0, 1,
-                            "harmonic", "maxwellian", 1.0, 0.0, false)
+                            "harmonic", "maxwellian", 1.0, false)
         @test vals ≈ [KF.energy_integrand_scalar(x, p) for x in x_grid]
     end
 
@@ -654,7 +654,6 @@
 
     @testset "KineticForcesState" begin
         state = KF.KineticForcesState()
-        @test !state.completed
         @test isempty(state.method_results)
         @test isempty(state.kinetic_matrices)
     end
@@ -663,7 +662,6 @@
         mr = KF.MethodResult(method="fgar", nn=1)
         @test mr.total_torque == 0.0 + 0.0im
         @test mr.total_energy == 0.0
-        @test isempty(mr.records)
     end
 
     # =========================================================================
diff --git a/test/runtests_kinetic_validity.jl b/test/runtests_kinetic_validity.jl
index a9e41c17ee..cca53fa453 100644
--- a/test/runtests_kinetic_validity.jl
+++ b/test/runtests_kinetic_validity.jl
@@ -107,7 +107,7 @@
         @test default_reg != 0.0   # otherwise the override above proves nothing
 
         # Unrelated kwargs pass through untouched.
-        out = PE.kinetic_regularization_kwargs(kinetic, pairs((; regularization_width=0.05, filter_modes=true)))
-        @test out.filter_modes == true
+        out = PE.kinetic_regularization_kwargs(kinetic, pairs((; regularization_width=0.05, verbose=false)))
+        @test out.verbose == false
     end
 end
diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl
index 7aea77bab6..6a67801cab 100644
--- a/test/runtests_rerun_from_h5.jl
+++ b/test/runtests_rerun_from_h5.jl
@@ -109,20 +109,11 @@ end
             end
         end
 
-        # A gpec.h5 written before the integrator refactor embeds retired keys such as
-        # `use_parallel` in its gpec_toml_raw, so a replay must warn and drop them rather than
-        # failing on an unknown keyword. The override flag stands in for that stale blob.
-        @testset "retired ForceFreeStates keys are warned and ignored on replay" begin
+        # A stored deck that sets a retired key stops the replay before any solve; the override
+        # flag stands in for a stale gpec_toml_raw.
+        @testset "retired ForceFreeStates keys stop a replay" begin
             mktempdir() do out_dir
-                inputs, = GeneralizedPerturbedEquilibrium.build_inputs_from_h5([source_h5, "--output-dir", out_dir, "--override", "ForceFreeStates.use_parallel=true"])
-                @test inputs["ForceFreeStates"]["use_parallel"] == true
-
-                @test_logs (:warn,) GeneralizedPerturbedEquilibrium._drop_deprecated_keys!(inputs["ForceFreeStates"], GeneralizedPerturbedEquilibrium._DEPRECATED_FFS_KEYS, "ForceFreeStates")
-                @test !haskey(inputs["ForceFreeStates"], "use_parallel")
-
-                # With the key gone the control struct builds again, on the default integrator.
-                ctrl = GeneralizedPerturbedEquilibrium.ForceFreeStates.ForceFreeStatesControl(; (Symbol(k) => v for (k, v) in inputs["ForceFreeStates"])...)
-                @test ctrl.integrator in ("forward", "riccati")
+                @test_throws r"`use_parallel`" GeneralizedPerturbedEquilibrium.main([source_h5, "--output-dir", out_dir, "--override", "ForceFreeStates.use_parallel=true"])
             end
         end
 

From 05d1faaba5dda8e9b7e048d39faa174bf890101e Mon Sep 17 00:00:00 2001
From: Matthew Pharr <m.pharr@protonmail.com>
Date: Mon, 26 Oct 2026 15:09:38 +0100
Subject: [PATCH 2/3] Repo - REFACTOR! - Stop a gpec.h5 replay with retired
 keys and point it to the original deck

A replay checks its stored deck before overrides and says the file was written by an older GPEC version, since overrides cannot delete keys. The stop message pluralizes on the number of keys found. Tests cover the stop through main on a staged deck and on a replayed gpec.h5.
---
 src/GeneralizedPerturbedEquilibrium.jl |  9 ++++++---
 src/Rerun.jl                           |  5 +++++
 test/runtests_equil.jl                 | 14 ++++++++++++--
 test/runtests_rerun_from_h5.jl         | 18 ++++++++++++++++--
 4 files changed, 39 insertions(+), 7 deletions(-)

diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index bac8319cc4..4eb4ec2c22 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -99,12 +99,15 @@ const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc", "u
 const _DEPRECATED_KF_KEYS = ("moment", "eq_out", "theta_out", "xlmda_out", "eqpsi_out", "save_records", "gpec_file",
     "fnml_flag", "ellip_flag", "diag_flag", "force_xialpha")
 const _DEPRECATED_PE_KEYS = ("fixed_boundary", "output_eigenmodes", "filter_modes", "singular_point_method")
+const _DEPRECATED_KEYS = ("Equilibrium" => _DEPRECATED_EQUIL_KEYS, "ForceFreeStates" => _DEPRECATED_FFS_KEYS,
+    "KineticForces" => _DEPRECATED_KF_KEYS, "PerturbedEquilibrium" => _DEPRECATED_PE_KEYS)
 
 # Stop the run when a parsed gpec.toml section sets keys that no longer exist, so the deck is updated.
-function _stop_on_deprecated_keys(table, deprecated_keys, section::String)
+function _stop_on_deprecated_keys(table, deprecated_keys, section::String; remedy::String="remove from gpec.toml")
     found = [k for k in deprecated_keys if haskey(table, k)]
-    isempty(found) || error("[$section] sets $(join(("`$k`" for k in found), ", ")), which no longer exist; remove them from gpec.toml")
-    return nothing
+    isempty(found) && return nothing
+    listed = join(("`$k`" for k in found), ", ")
+    error(length(found) == 1 ? "[$section] sets $listed, which is no longer a key; $remedy" : "[$section] sets $listed, which are no longer keys; $remedy")
 end
 
 function main(args::Vector{String}=String[]; dd::Union{IMASdd.dd,Nothing}=nothing)
diff --git a/src/Rerun.jl b/src/Rerun.jl
index e6cd618799..f280fe1e09 100644
--- a/src/Rerun.jl
+++ b/src/Rerun.jl
@@ -247,6 +247,11 @@ function build_inputs_from_h5(args::Vector{String})
     # Parse TOML blob, replay the frozen forcing file, then fold in overrides: file overrides
     # first, then CLI flags, so an overridden [ForcingTerms] still wins.
     inputs = TOML.parse(toml_raw)
+    # Overrides can set keys but not delete them, so a stored deck with retired keys cannot be replayed.
+    for (section, keys) in _DEPRECATED_KEYS
+        haskey(inputs, section) && _stop_on_deprecated_keys(inputs[section], keys, section;
+            remedy="$source_h5 was written by an older GPEC version; re-run from the original deck instead of replaying it")
+    end
     if forcing_snapshot !== nothing
         forcing_file = splitext(output_name)[1] * "_forcing.h5"
         modes, normalization = forcing_snapshot
diff --git a/test/runtests_equil.jl b/test/runtests_equil.jl
index 7ed9db7614..70d75d845f 100644
--- a/test/runtests_equil.jl
+++ b/test/runtests_equil.jl
@@ -1,4 +1,5 @@
 using Statistics
+using TOML
 
 @testset "Equilibrium Unit Tests" begin
 
@@ -169,11 +170,20 @@ using Statistics
         GPEC = GeneralizedPerturbedEquilibrium
         @test_throws r"\[ForceFreeStates\] sets `reform_eq_with_psilim`" GPEC._stop_on_deprecated_keys(
             Dict{String,Any}("reform_eq_with_psilim" => false, "nn_low" => 1), GPEC._DEPRECATED_FFS_KEYS, "ForceFreeStates")
-        @test_throws r"`moment`, `save_records`" GPEC._stop_on_deprecated_keys(
+        @test_throws r"`moment`, `save_records`, which are no longer keys" GPEC._stop_on_deprecated_keys(
             Dict{String,Any}("moment" => "pressure", "save_records" => false, "fgar_flag" => true), GPEC._DEPRECATED_KF_KEYS, "KineticForces")
-        @test_throws r"\[PerturbedEquilibrium\] sets `fixed_boundary`" GPEC._stop_on_deprecated_keys(
+        @test_throws r"\[PerturbedEquilibrium\] sets `fixed_boundary`, which is no longer a key" GPEC._stop_on_deprecated_keys(
             Dict{String,Any}("fixed_boundary" => false, "compute_response" => false), GPEC._DEPRECATED_PE_KEYS, "PerturbedEquilibrium")
         @test GPEC._stop_on_deprecated_keys(Dict{String,Any}("nn_low" => 1), GPEC._DEPRECATED_FFS_KEYS, "ForceFreeStates") === nothing
+        # The stop runs through main, before any equilibrium work.
+        mktempdir() do deck_dir
+            template_dir = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example")
+            foreach(name -> cp(joinpath(template_dir, name), joinpath(deck_dir, name)), readdir(template_dir))
+            inputs = TOML.parsefile(joinpath(deck_dir, "gpec.toml"))
+            inputs["KineticForces"] = Dict{String,Any}("moment" => 1)
+            open(io -> TOML.print(io, inputs), joinpath(deck_dir, "gpec.toml"), "w")
+            @test_throws r"\[KineticForces\] sets `moment`, which is no longer a key; remove from gpec.toml" GPEC.main([deck_dir])
+        end
         # A listed deck key must not still be a control field, or listing it would stop runs that set a live key.
         @test !any(k -> hasfield(GPEC.ForceFreeStates.ForceFreeStatesControl, Symbol(k)), GPEC._DEPRECATED_FFS_KEYS)
         @test !any(k -> hasfield(GPEC.KineticForces.KineticForcesControl, Symbol(k)), GPEC._DEPRECATED_KF_KEYS)
diff --git a/test/runtests_rerun_from_h5.jl b/test/runtests_rerun_from_h5.jl
index 6a67801cab..f97079c78e 100644
--- a/test/runtests_rerun_from_h5.jl
+++ b/test/runtests_rerun_from_h5.jl
@@ -109,14 +109,28 @@ end
             end
         end
 
-        # A stored deck that sets a retired key stops the replay before any solve; the override
-        # flag stands in for a stale gpec_toml_raw.
+        # An override that sets a retired key stops the replay before any solve.
         @testset "retired ForceFreeStates keys stop a replay" begin
             mktempdir() do out_dir
                 @test_throws r"`use_parallel`" GeneralizedPerturbedEquilibrium.main([source_h5, "--output-dir", out_dir, "--override", "ForceFreeStates.use_parallel=true"])
             end
         end
 
+        # A gpec.h5 from an older GPEC stores a deck with retired keys, which the replay cannot delete.
+        @testset "a stored deck with retired keys stops the replay" begin
+            mktempdir() do out_dir
+                old_h5 = joinpath(out_dir, "old.h5")
+                cp(source_h5, old_h5)
+                h5open(old_h5, "r+") do h5
+                    inputs = TOML.parse(read(h5, "Input/gpec_toml_raw"))
+                    inputs["PerturbedEquilibrium"]["fixed_boundary"] = false
+                    delete_object(h5, "Input/gpec_toml_raw")
+                    h5["Input/gpec_toml_raw"] = sprint(TOML.print, inputs)
+                end
+                @test_throws r"`fixed_boundary`, which is no longer a key; .*older GPEC version" GeneralizedPerturbedEquilibrium.main([old_h5, "--output-dir", out_dir])
+            end
+        end
+
         @testset "rerun refuses to overwrite its own source" begin
             # Default output is `<basename>_rerun.h5`, so clobbering only happens if the user
             # routes the rerun back onto the source file; exercise that guard.

From 5786d2d84e13a7bf00731e721f01aa05ff93c033 Mon Sep 17 00:00:00 2001
From: Matthew Pharr <m.pharr@protonmail.com>
Date: Tue, 27 Oct 2026 09:16:23 +0100
Subject: [PATCH 3/3] Test - TEST - Build equilibria with PlasmaEquilibrium and
 promote the re-formed analytic grid to a passing test

setup_equilibrium is gone from develop; three tests still called it. The equilibrium now stores its analytic config, so the deck/API parity test's known drift no longer occurs.
---
 test/runtests_coils.jl          | 2 +-
 test/runtests_localstability.jl | 2 +-
 test/runtests_solve_api.jl      | 6 +++---
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/runtests_coils.jl b/test/runtests_coils.jl
index 255c1d65a8..f5fde2f383 100644
--- a/test/runtests_coils.jl
+++ b/test/runtests_coils.jl
@@ -347,7 +347,7 @@ end
 # ---------------------------------------------------------------------------
 @testset "CoilFourier: B·n̂ file conversion onto the run's harmonics" begin
     eq_config = Equilibrium.EquilibriumConfig(; eq_type="sol", psilow=0.01, psihigh=0.995)
-    equil = Equilibrium.setup_equilibrium(eq_config, Equilibrium.SolovevConfig())
+    equil = Equilibrium.PlasmaEquilibrium(eq_config, Equilibrium.SolovevConfig())
     n, m_out, mtheta = 1, -10:14, 256
     convert(modes) = ForcingTerms.convert_forcing_normalization(modes, "normal_field_T", equil, n, m_out)
     spectrum(ms, amp) = [ForcingMode(; n, m, amplitude=amp(m)) for m in ms]
diff --git a/test/runtests_localstability.jl b/test/runtests_localstability.jl
index 2d2866bbb8..7ee46667e4 100644
--- a/test/runtests_localstability.jl
+++ b/test/runtests_localstability.jl
@@ -8,7 +8,7 @@ using TOML
     LS = GPEC.LocalStability
     template = joinpath(@__DIR__, "test_data", "regression_solovev_ideal_example")
     deck = TOML.parsefile(joinpath(template, "gpec.toml"))
-    equil = GPEC.Equilibrium.setup_equilibrium(
+    equil = GPEC.Equilibrium.PlasmaEquilibrium(
         GPEC.Equilibrium.EquilibriumConfig(deck["Equilibrium"], template),
         GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"])
     )
diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl
index 7237c017fc..56e9f34fef 100644
--- a/test/runtests_solve_api.jl
+++ b/test/runtests_solve_api.jl
@@ -377,7 +377,7 @@ using HDF5
             ffs_overrides = Dict{String,Any}("vac_flag" => false)
             mktempdir() do dir
                 reference = GPEC.main([_stage_deck(dir; eq_overrides=eq_overrides, ffs_overrides=ffs_overrides)]).ffs
-                pass1 = GPEC.Equilibrium.setup_equilibrium(
+                pass1 = GPEC.Equilibrium.PlasmaEquilibrium(
                     GPEC.Equilibrium.EquilibriumConfig(merge(deck["Equilibrium"], eq_overrides), template),
                     GPEC.Equilibrium.SolovevConfig(deck["SOL_INPUT"])
                 )
@@ -385,8 +385,8 @@ using HDF5
 
                 kwargs = merge(ffs_kwargs, Dict(Symbol(k) => v for (k, v) in ffs_overrides))
                 api = solve(pass1, Forward(); nn=1, wall=wall, dir_path=dir, kwargs...)
-                # Known drift: the equilibrium does not store its analytic config, so `solve` cannot re-form it.
-                @test_broken api.equil.rzphi_xs == reference.equil.rzphi_xs
+                # The equilibrium stores its analytic config, so `solve` re-forms it on the deck's grid.
+                @test api.equil.rzphi_xs == reference.equil.rzphi_xs
             end
         end