From ce5148aa34eba3d262ab3f3508f1cf192906ef09 Mon Sep 17 00:00:00 2001
From: logan-nc <nikolas.logan@columbia.edu>
Date: Mon, 19 Oct 2026 01:12:11 -0400
Subject: [PATCH 1/5] KineticForces - FEATURE! - Write the perturbation that
 drives the torque

KineticForces builds the Lagrangian delta|B|/B and div(xi_perp) spectra from the
PerturbedEquilibrium Clebsch displacements (`set_perturbation_data!`, the port of
Fortran pentrc set_peq) and then keeps them only as interpolants, so nothing of the
perturbation survives into gpec.h5. The torque is the only trace left of the field
that produced it.

That matters outside GPEC. The helical ripple delta_B/B on a flux surface, and its
poloidal spectrum, is the input another NTV code needs to be run on a GPEC field --
it is the quantity the ITPA TC-24 benchmarking activity asked for so NTVTOK could be
run on GPEC output. Fortran GPEC publishes it via `pmodb_flag`; Julia had no route to
it at all.

Retain the two mode tables alongside the interpolants and write them once per run to
`KineticForces/Perturbation/`:

    psi               PE psi_N grid
    m, n              mode numbers
    dBoB_lagrangian   Lagrangian delta|B|/B Fourier modes  (dimensionless)
    div_xi_perp       div(xi_perp) Fourier modes           [1/m]

Lagrangian, not Eulerian: Fortran fills the same spline from `lagbfun(j)/eqfun%f(1)`
(pentrc/inputs.f90:477), and the matrix route this code takes is its sibling branch,
`-(jbbdivx + jbbkapx)` at inputs.f90:855.

The perturbation is a property of the run, not of a species or a method, so it is
written by the call that has no `species_label` and guarded against a second write --
a multi-species run calls `write_to_hdf5!` once per species and once for the total.
Nothing in the torque calculation reads the new fields, so this is additive: no
existing output changes.
---
 src/GeneralizedPerturbedEquilibrium.jl    |  4 +-
 src/KineticForces/KineticForcesStructs.jl | 15 +++++++
 src/KineticForces/Output.jl               | 30 ++++++++++++-
 test/runtests_kinetic.jl                  | 52 +++++++++++++++++++++++
 4 files changed, 97 insertions(+), 4 deletions(-)

diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 6322d2c79..4dfcfb53d 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -1077,7 +1077,7 @@ function run_kinetic_forces(
             if kf_ctrl.write_outputs_to_HDF5
                 h5open(joinpath(result.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file
                     KineticForces.write_to_hdf5!(h5file, kf_state;
-                        dVdpsi_spline=result.equil.profiles.dVdpsi_spline)
+                        dVdpsi_spline=result.equil.profiles.dVdpsi_spline, perturbation=kf_intr)
                 end
             end
         else
@@ -1102,7 +1102,7 @@ function run_kinetic_forces(
                             dVdpsi_spline=result.equil.profiles.dVdpsi_spline, species_label=label)
                     end
                     KineticForces.write_to_hdf5!(h5file, kf_state;
-                        dVdpsi_spline=result.equil.profiles.dVdpsi_spline)
+                        dVdpsi_spline=result.equil.profiles.dVdpsi_spline, perturbation=kf_intr)
                 end
             end
         end
diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl
index ed9ef14a6..346074010 100644
--- a/src/KineticForces/KineticForcesStructs.jl
+++ b/src/KineticForces/KineticForcesStructs.jl
@@ -243,6 +243,14 @@ this struct.
     dbob_m::Any = nothing          # δB/B perturbation modes (CubicSeriesInterpolant)
     divx_m::Any = nothing          # ∇·ξ⊥ perturbation modes (CubicSeriesInterpolant)
 
+    # The same two quantities as the raw mode tables the interpolants were built from, kept
+    # so they can be written to the HDF5 file. δB/B is the Lagrangian perturbation of the
+    # field strength, which is what an NTV code needs as input; nothing in the torque
+    # calculation reads these fields back.
+    perturbation_psi::Vector{Float64} = Float64[]            # ψ_N of each row
+    dbob_m_data::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # (nψ, mpert)
+    divx_m_data::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # (nψ, mpert)
+
     # Sticky bracket-search hints for amortized lookup in `tpsi!`. When the outer
     # quadrature evaluates ψ non-monotonically, FastInterpolations falls back to a
     # broader search when the hint is stale — still a net win over a cold bracket
@@ -445,6 +453,13 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur
     kf_intr.dbob_m = cubic_interp(psi_grid, Series(dbob_m_data); itp_opts...)
     kf_intr.divx_m = cubic_interp(psi_grid, Series(divx_m_data); itp_opts...)
 
+    # Retain the mode tables for output. They are the physical δB/B and ∇·ξ⊥ spectra on the
+    # PE ψ grid, in the same poloidal angle and (m, n) convention as the PerturbedEquilibrium
+    # field spectra, and are of direct interest to other NTV codes run on GPEC fields.
+    kf_intr.perturbation_psi = copy(psi_grid)
+    kf_intr.dbob_m_data = dbob_m_data
+    kf_intr.divx_m_data = divx_m_data
+
     if kf_intr.verbose
         @info "set_perturbation_data!: built dbob_m/divx_m/xs_m interpolants " *
               "(npsi=$npsi, mpert=$mpert)"
diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl
index 75f527ae2..902cb259c 100644
--- a/src/KineticForces/Output.jl
+++ b/src/KineticForces/Output.jl
@@ -8,7 +8,7 @@ then write to gpec.h5 in a single pass.
 
 """
     write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spline=nothing,
-                   species_label=nothing)
+                   species_label=nothing, perturbation=nothing)
 
 Write KineticForces results to the "KineticForces" group in gpec.h5.
 
@@ -18,13 +18,16 @@ Write KineticForces results to the "KineticForces" group in gpec.h5.
 - `dVdpsi_spline`: Optional dV/dψ_N profile interpolant; when given, dV/dψ_N is
   written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ)
   is directly available
+- `perturbation`: optional `KineticForcesInternal`; its δB/B and ∇·ξ⊥ mode tables are written
+  once to `KineticForces/Perturbation/`. The perturbation is a property of the run, not of a
+  species, so it is written only on the call that has no `species_label`.
 - `species_label`: `nothing` writes the run total to `KineticForces/<method>/`; a label
   (e.g. `"ion_z1_m2"`, `"electron"`) writes one species' contribution to
   `KineticForces/PerSpecies/<label>/<method>/`. A multi-species run calls this once per
   species and once for the summed total, so the group is opened-or-created each time.
 """
 function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spline=nothing,
-    species_label::Union{Nothing,AbstractString}=nothing)
+    species_label::Union{Nothing,AbstractString}=nothing, perturbation=nothing)
     root = haskey(h5file, "KineticForces") ? h5file["KineticForces"] : create_group(h5file, "KineticForces")
     g = root
     if species_label !== nothing
@@ -94,8 +97,31 @@ function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spl
     for method_name in method_names
         Utilities.HDF5Annotations.annotate!(g[method_name], KF_METHOD_H5_ANNOTATIONS)
     end
+
+    # The perturbation driving the torque: one table per run, not per species or per method.
+    if species_label === nothing && perturbation !== nothing &&
+       !isempty(perturbation.perturbation_psi) && !haskey(root, "Perturbation")
+        pg = create_group(root, "Perturbation")
+        pg["psi"] = perturbation.perturbation_psi
+        pg["m"] = collect(perturbation.mlow:perturbation.mhigh)
+        pg["n"] = perturbation.nlow
+        pg["dBoB_lagrangian"] = perturbation.dbob_m_data
+        pg["div_xi_perp"] = perturbation.divx_m_data
+        Utilities.HDF5Annotations.annotate!(pg, KF_PERTURBATION_H5_ANNOTATIONS)
+    end
 end
 
+# Metadata for KineticForces/Perturbation/. δB/B and ∇·ξ⊥ are the JBB-deweighted physical
+# spectra built in `set_perturbation_data!` (Fortran pentrc set_peq), in the same poloidal
+# angle and sign convention as the PerturbedEquilibrium field spectra.
+const KF_PERTURBATION_H5_ANNOTATIONS = [
+    "psi" => (; long_name="normalized poloidal flux ψ_N of the perturbation grid", scale="psi"),
+    "m" => (; long_name="poloidal mode number m", scale="m"),
+    "n" => (; long_name="toroidal mode number n of the perturbation"),
+    "dBoB_lagrangian" => (; long_name="Lagrangian δ|B|/B Fourier modes, δB/B(θ) = Σ_m value·exp(2πi·m·θ) with θ ∈ [0,1); dimensionless", dims=("psi", "m"), attach=(1 => "psi", 2 => "m")),
+    "div_xi_perp" => (; long_name="∇·ξ_⊥ Fourier modes, same angle and sign convention as dBoB_lagrangian", units="1/m", dims=("psi", "m"), attach=(1 => "psi", 2 => "m"))
+]
+
 # Storage slice → Logan 2015 matrix letter (Eqs 7.30-7.35); see BounceAveraging.jl
 # for the Hermitian/full packing details.
 const _KINETIC_MATRIX_LETTERS = (("A", 1), ("B", 2), ("C", 3), ("D", 4), ("E", 5), ("H", 6))
diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl
index 4f47db9fd..f2545f7bf 100644
--- a/test/runtests_kinetic.jl
+++ b/test/runtests_kinetic.jl
@@ -1,3 +1,4 @@
+using HDF5
 
 @testset "KineticForces Unit Tests" begin
 
@@ -666,6 +667,57 @@
         @test isempty(mr.records)
     end
 
+    @testset "perturbation output" begin
+        # KineticForces/Perturbation/ carries the δB/B and ∇·ξ⊥ spectra that drive the torque,
+        # so another NTV code can be run on the same perturbation. It is a property of the run,
+        # not of a species or a method, so it is written once by the call with no species_label.
+        @test isempty(KF.KineticForcesInternal().perturbation_psi)
+        @test isempty(KF.KineticForcesInternal().dbob_m_data)
+
+        npsi, mpert = 4, 3
+        intr = KF.KineticForcesInternal(; mlow=-1, mhigh=1, mpert=mpert, nlow=3, nhigh=3,
+            perturbation_psi=collect(range(0.1, 0.9; length=npsi)),
+            dbob_m_data=ComplexF64[0.1 * i + 0.01im * j for i in 1:npsi, j in 1:mpert],
+            divx_m_data=ComplexF64[0.2 * i - 0.02im * j for i in 1:npsi, j in 1:mpert])
+        st = KF.KineticForcesState()
+        st.method_results["fgar"] = KF.MethodResult(; method="fgar", nn=3,
+            psi_grid=collect(0.0:0.25:1.0), dtdpsi=fill(ComplexF64(1.0), 5),
+            t_cumulative=zeros(ComplexF64, 5))
+
+        h5 = tempname() * ".h5"
+        HDF5.h5open(h5, "cw") do f
+            KF.write_to_hdf5!(f, st; species_label="electron", perturbation=intr)  # must NOT write it
+            KF.write_to_hdf5!(f, st; perturbation=intr)
+        end
+        HDF5.h5open(h5, "r") do f
+            @test !haskey(f["KineticForces"]["PerSpecies"]["electron"], "Perturbation")
+            pg = f["KineticForces"]["Perturbation"]
+            @test read(pg["m"]) == [-1, 0, 1]
+            @test read(pg["n"]) == 3
+            @test read(pg["psi"]) ≈ intr.perturbation_psi
+            # HDF5.jl round-trips the Julia shape, so (npsi, mpert) here; h5py and MATLAB see
+            # the transpose, matching the field spectra written by PerturbedEquilibrium.
+            @test size(read(pg["dBoB_lagrangian"])) == (npsi, mpert)
+            @test read(pg["dBoB_lagrangian"]) ≈ intr.dbob_m_data
+            @test read(pg["div_xi_perp"]) ≈ intr.divx_m_data
+            for d in ("dBoB_lagrangian", "div_xi_perp")
+                a = HDF5.attributes(pg[d])
+                @test haskey(a, "long_name") && haskey(a, "dims")
+            end
+        end
+        rm(h5; force=true)
+
+        # Nothing is written when the perturbation was never built (PE did not run).
+        h5b = tempname() * ".h5"
+        HDF5.h5open(h5b, "cw") do f
+            KF.write_to_hdf5!(f, st; perturbation=KF.KineticForcesInternal())
+        end
+        HDF5.h5open(h5b, "r") do f
+            @test !haskey(f["KineticForces"], "Perturbation")
+        end
+        rm(h5b; force=true)
+    end
+
     # =========================================================================
     # Kinetic profile I/O: HDF5 schema, ASCII back-compat, writer round-trip
     # =========================================================================

From a69bed782f301f13457fa866d43fd537c6d5e4a6 Mon Sep 17 00:00:00 2001
From: logan-nc <nikolas.logan@columbia.edu>
Date: Mon, 19 Oct 2026 20:17:48 -0400
Subject: [PATCH 2/5] Docs - DOCS - Record the KineticForces perturbation group
 in the schema table

The schema table in hdf5-conventions.md enumerates the top-level groups and what each
holds, so a new group that is not listed there is invisible to anyone reading the schema
rather than the writer. `KineticForces/Perturbation/` is written once per run rather than
per species or per method, which is the part of its placement that is not guessable from
the group name.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FabuzB6rkkv638ZUCzaAA2
---
 docs/development/hdf5-conventions.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md
index c4631c984..82a32c12a 100644
--- a/docs/development/hdf5-conventions.md
+++ b/docs/development/hdf5-conventions.md
@@ -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), `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, `EnergyIntegrals/`, `KineticMatrices/`); `Perturbation/` (the Lagrangian δ\|B\|/B and ∇·ξ_⊥ mode tables the torque is evaluated on, written once per run rather than per species or per method); 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 |

From d98788ec7116dd0c9cb0d22e67de8fadb75f19eb Mon Sep 17 00:00:00 2001
From: logan-nc <nikolas.logan@columbia.edu>
Date: Tue, 20 Oct 2026 19:48:16 -0400
Subject: [PATCH 3/5] PerturbedEquilibrium - FEATURE - Write the perturbation
 that drives neoclassical transport

The Lagrangian delta|B|/B and div.xi_perp spectra are what an NTV calculation is
driven by, and GPEC already formed them -- but only inside KineticForces, on the
way to the torque, so getting them out meant running our own NTV model. Another
code wanting to run on GPEC's perturbation had no way to read it.

Nothing kinetic enters them. They are the PE Clebsch displacement put through the
S,T,X,Y,Z geometric matrices and deweighted by J.B^2, so they belong to the
perturbed equilibrium -- which is where the Fortran parent puts them: gpout_pmodb
runs in gpec under pmodb_flag, and pentrc reads that file back (read_pmodb). The
Julia port fused the two executables and the output landed on the wrong side of a
seam that is no longer visible.

So move the construction to PerturbedEquilibrium.compute_modb_modes, write it to
PerturbedEquilibrium/Response/{dBoB_lagrangian,div_xi_perp} beside the
displacement it comes from, and let KineticForces consume it. An ideal run with no
[KineticForces] block now carries the spectra; KineticForces gets smaller.

Harmonics only, as everywhere else in Response/: Analysis.modes_to_theta
reconstructs theta-space from any (psi, numpert_total) dataset, so the real-space
forms are a post-processing call rather than stored arrays. Only the two
independent spectra are written -- the curvature term is
-(dBoB_lagrangian + div_xi_perp), and the Eulerian delta|B|/B differs from the
Lagrangian by xi.grad|B|/|B|, which xi_clebsch_* and the equilibrium geometry
already determine.

Single toroidal mode only, as before: the geometric matrices are assembled at
n = nlow. Multi-n runs write nothing rather than something wrong.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FabuzB6rkkv638ZUCzaAA2
---
 docs/development/hdf5-conventions.md          |   4 +-
 src/GeneralizedPerturbedEquilibrium.jl        |   8 +-
 src/KineticForces/KineticForcesStructs.jl     | 140 +++---------------
 src/KineticForces/Output.jl                   |  30 +---
 .../FieldReconstruction.jl                    | 129 ++++++++++++++++
 .../PerturbedEquilibriumStructs.jl            |   4 +
 src/PerturbedEquilibrium/Response.jl          |  16 +-
 src/PerturbedEquilibrium/Utils.jl             |  20 ++-
 test/runtests_h5_schema.jl                    |  39 +++++
 test/runtests_kinetic.jl                      |  52 -------
 10 files changed, 234 insertions(+), 208 deletions(-)

diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md
index b66a3d3d8..93df52d43 100644
--- a/docs/development/hdf5-conventions.md
+++ b/docs/development/hdf5-conventions.md
@@ -43,8 +43,8 @@ Top level (11 groups):
 | `ForceFreeStates/` | `Solutions/ForwardIntegration/` (u-solutions), `Solutions/GalerkinIntegration/` (closed ξ profiles in the shared layout, the gal surface list, debug-gated `Basis/`), `EulerLagrangeMatrices/{Ideal,Kinetic}`, `FreeBoundaryStability/`, `EdgeScan/` |
 | `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/`); `Perturbation/` (the Lagrangian δ\|B\|/B and ∇·ξ_⊥ mode tables the torque is evaluated on, written once per run rather than per species or per method); multi-ion runs add `PerSpecies/<species>/<method>/` with the same per-method layout, summing to the top-level total |
+| `PerturbedEquilibrium/` | `ForcingModes/`, `Response/` (displacement and field spectra, plus the `dBoB_lagrangian` and `div_xi_perp` spectra that drive neoclassical transport), `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 |
 | `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 |
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 061a5afd6..add66e67a 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -1094,7 +1094,7 @@ function run_kinetic_forces(
     else
         # kf_ctrl and kinetic_profiles were loaded once before the equilibrium was re-formed.
         kf_intr = KineticForces.KineticForcesInternal(result.equil; verbose=kf_ctrl.verbose)
-        KineticForces.set_perturbation_data!(kf_intr, pe_state, result, result.equil, result.metric)
+        KineticForces.set_perturbation_data!(kf_intr, pe_state, result)
 
         if species === nothing
             kf_state = KineticForces.KineticForcesState()
@@ -1102,7 +1102,7 @@ function run_kinetic_forces(
             if kf_ctrl.write_outputs_to_HDF5
                 h5open(joinpath(result.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file
                     KineticForces.write_to_hdf5!(h5file, kf_state;
-                        dVdpsi_spline=result.equil.profiles.dVdpsi_spline, perturbation=kf_intr)
+                        dVdpsi_spline=result.equil.profiles.dVdpsi_spline)
                 end
             end
         else
@@ -1127,7 +1127,7 @@ function run_kinetic_forces(
                             dVdpsi_spline=result.equil.profiles.dVdpsi_spline, species_label=label)
                     end
                     KineticForces.write_to_hdf5!(h5file, kf_state;
-                        dVdpsi_spline=result.equil.profiles.dVdpsi_spline, perturbation=kf_intr)
+                        dVdpsi_spline=result.equil.profiles.dVdpsi_spline)
                 end
             end
         end
@@ -1345,7 +1345,7 @@ function efc_couplings(
         pe_intr.forcing_modes = modes
         pe_state = PerturbedEquilibrium.compute_perturbed_equilibrium(ffs, ForcingTerms.RMPField(ForcingTerms.ForcingTermsControl()), quiet, pe_intr)
         kf_intr = KineticForces.KineticForcesInternal(equil; verbose=verbose)
-        KineticForces.set_perturbation_data!(kf_intr, pe_state, ffs, equil, ffs.metric)
+        KineticForces.set_perturbation_data!(kf_intr, pe_state, ffs)
         return kf_intr
     end
     # Torque and cumulative torque profile (resampled onto the kinetic ψ grid) of one perturbation at one rotation.
diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl
index 7a9959983..327ac9d2f 100644
--- a/src/KineticForces/KineticForcesStructs.jl
+++ b/src/KineticForces/KineticForcesStructs.jl
@@ -243,14 +243,6 @@ this struct.
     dbob_m::Any = nothing          # δB/B perturbation modes (CubicSeriesInterpolant)
     divx_m::Any = nothing          # ∇·ξ⊥ perturbation modes (CubicSeriesInterpolant)
 
-    # The same two quantities as the raw mode tables the interpolants were built from, kept
-    # so they can be written to the HDF5 file. δB/B is the Lagrangian perturbation of the
-    # field strength, which is what an NTV code needs as input; nothing in the torque
-    # calculation reads these fields back.
-    perturbation_psi::Vector{Float64} = Float64[]            # ψ_N of each row
-    dbob_m_data::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # (nψ, mpert)
-    divx_m_data::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)  # (nψ, mpert)
-
     # Sticky bracket-search hints for amortized lookup in `tpsi!`. When the outer
     # quadrature evaluates ψ non-monotonically, FastInterpolations falls back to a
     # broader search when the hint is stale — still a net win over a cold bracket
@@ -329,25 +321,21 @@ function KineticForcesInternal(equil; verbose::Bool=false)
 end
 
 """
-    set_perturbation_data!(kf_intr, pe_state, ffs, equil, metric)
+    set_perturbation_data!(kf_intr, pe_state, ffs)
 
 Populate perturbation data from PerturbedEquilibriumState into KineticForcesInternal.
 
-Builds three interpolant sets from PE Clebsch displacements:
-1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ
-2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq)
-3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting
+Builds three interpolant sets over ψ:
+1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] from the PE Clebsch displacements
+2. `dbob_m` — Lagrangian δ|B|/B, taken from `pe_state.modb_modes`
+3. `divx_m` — ∇·ξ⊥, taken from `pe_state.modb_modes`
 
-The JBB deweighting algorithm (Fortran pentrc/inputs.f90:828-868):
-1. Apply geometric matrices S,T,X,Y,Z in m-space
-2. Inverse DFT to θ-space
-3. Divide by J·B² at each θ
-4. Forward DFT back to m-space
+The perturbation spectra themselves are a property of the perturbed equilibrium and are
+derived and written there (`PerturbedEquilibrium.compute_modb_modes`); this function only
+interpolates them onto the torque quadrature.
 """
 function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::PerturbedEquilibrium.PerturbedEquilibriumState,
-                                ffs::ForceFreeStates.ForceFreeStatesResult,
-                                equil::Equilibrium.PlasmaEquilibrium,
-                                metric::ForceFreeStates.MetricData)
+                                ffs::ForceFreeStates.ForceFreeStatesResult)
     ffs.npert == 1 || error("KineticForces supports a single toroidal mode; this run has n = $(ffs.nlow):$(ffs.nhigh).")
 
     # Copy mode numbers from FFS
@@ -386,79 +374,19 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur
     xs_m_3 = cubic_interp(psi_grid, Series(clebsch_alpha_mat); itp_opts...)
     kf_intr.xs_m = [xs_m_1, xs_m_2, xs_m_3]
 
-    # Build geometric matrices (S,T,X,Y,Z) for JBB deweighting
-    geom_mats = ForceFreeStates.build_kinetic_metric_matrices(equil, ffs, metric)
-
-    # Build FourierTransform for the JBB deweighting DFT round-trip
-    mthsurf = kf_intr.mthsurf
-    ft = Utilities.FourierTransforms.FourierTransform(mthsurf, mpert, ffs.mlow)
-
-    # JBB deweighting: convert Clebsch modes → physical δB/B and ∇·ξ⊥ modes
-    dbob_m_data = zeros(ComplexF64, npsi, mpert)
-    divx_m_data = zeros(ComplexF64, npsi, mpert)
-
-    # Pre-allocate buffers
-    smat_flat = Vector{ComplexF64}(undef, mpert^2)
-    tmat_flat = Vector{ComplexF64}(undef, mpert^2)
-    xmat_flat = Vector{ComplexF64}(undef, mpert^2)
-    ymat_flat = Vector{ComplexF64}(undef, mpert^2)
-    zmat_flat = Vector{ComplexF64}(undef, mpert^2)
-    jbb_kapx = Vector{ComplexF64}(undef, mpert)
-    jbb_divx = Vector{ComplexF64}(undef, mpert)
-    jbb_dbob = Vector{ComplexF64}(undef, mpert)
-    theta_buf = Vector{ComplexF64}(undef, mthsurf)
-
-    hint_s = Ref(1)
-    hint_t = Ref(1)
-    hint_x = Ref(1)
-    hint_y = Ref(1)
-    hint_z = Ref(1)
-
-    for ipsi in 1:npsi
-        psi = psi_grid[ipsi]
-
-        # Get Clebsch displacement vectors at this ψ
-        xsp  = view(xi_modes.clebsch_psi,  ipsi, :)       # ξ^ψ [mpert]
-        xmp1 = view(xi_modes.clebsch_psi1, ipsi, :)       # ∂ξ^ψ/∂ψ [mpert]
-        xms  = view(clebsch_alpha_mat, ipsi, :)            # ξ^α [mpert]
-
-        # Evaluate geometric matrices at ψ → mpert² flat vectors, reshape to mpert×mpert
-        geom_mats.smats(smat_flat, psi; hint=hint_s)
-        geom_mats.tmats(tmat_flat, psi; hint=hint_t)
-        geom_mats.xmats(xmat_flat, psi; hint=hint_x)
-        geom_mats.ymats(ymat_flat, psi; hint=hint_y)
-        geom_mats.zmats(zmat_flat, psi; hint=hint_z)
-
-        smat = reshape(smat_flat, mpert, mpert)
-        tmat = reshape(tmat_flat, mpert, mpert)
-        xmat = reshape(xmat_flat, mpert, mpert)
-        ymat = reshape(ymat_flat, mpert, mpert)
-        zmat = reshape(zmat_flat, mpert, mpert)
-
-        # Apply geometric matrices in m-space. xs_m ordering: (1)=∂ξ^ψ/∂ψ, (2)=ξ^ψ, (3)=ξ^α.
-        #   jbb_kapx = smat·xsp + tmat·xms;  jbb_divx = xmat·xmp1 + ymat·xsp + zmat·xms
-        mul!(jbb_kapx, smat, xsp)
-        mul!(jbb_kapx, tmat, xms, 1.0 + 0.0im, 1.0 + 0.0im)   # += tmat * xms
-        mul!(jbb_divx, xmat, xmp1)
-        mul!(jbb_divx, ymat, xsp,  1.0 + 0.0im, 1.0 + 0.0im)  # += ymat * xsp
-        mul!(jbb_divx, zmat, xms,  1.0 + 0.0im, 1.0 + 0.0im)  # += zmat * xms
-        @. jbb_dbob = -(jbb_divx + jbb_kapx)
-
-        # Inverse DFT to θ-space, divide by J·B², forward DFT back
-        _jbb_deweight!(view(dbob_m_data, ipsi, :), jbb_dbob, ft, psi, equil, mthsurf, theta_buf)
-        _jbb_deweight!(view(divx_m_data, ipsi, :), jbb_divx, ft, psi, equil, mthsurf, theta_buf)
-    end
+    # δB/B and ∇·ξ⊥ come from PerturbedEquilibrium, which derives them from the same
+    # displacement (`compute_modb_modes`, Fortran gpout_pmodb / set_peq) and writes them to
+    # `PerturbedEquilibrium/Response/`. Nothing kinetic enters them, so they are not rebuilt
+    # here; a run whose npert > 1 has already failed the single-mode check above.
+    modb = pe_state.modb_modes
+    modb === nothing && error("set_perturbation_data!: the perturbed equilibrium carries no " *
+                              "δ|B|/B spectra; KineticForces cannot evaluate the torque without them.")
+    size(modb.dbob) == (npsi, mpert) || error("set_perturbation_data!: δ|B|/B spectra are " *
+        "$(size(modb.dbob)), expected ($npsi, $mpert).")
 
     # Build CubicSeriesInterpolants over ψ for dbob_m and divx_m
-    kf_intr.dbob_m = cubic_interp(psi_grid, Series(dbob_m_data); itp_opts...)
-    kf_intr.divx_m = cubic_interp(psi_grid, Series(divx_m_data); itp_opts...)
-
-    # Retain the mode tables for output. They are the physical δB/B and ∇·ξ⊥ spectra on the
-    # PE ψ grid, in the same poloidal angle and (m, n) convention as the PerturbedEquilibrium
-    # field spectra, and are of direct interest to other NTV codes run on GPEC fields.
-    kf_intr.perturbation_psi = copy(psi_grid)
-    kf_intr.dbob_m_data = dbob_m_data
-    kf_intr.divx_m_data = divx_m_data
+    kf_intr.dbob_m = cubic_interp(psi_grid, Series(modb.dbob); itp_opts...)
+    kf_intr.divx_m = cubic_interp(psi_grid, Series(modb.divx); itp_opts...)
 
     if kf_intr.verbose
         @info "set_perturbation_data!: built dbob_m/divx_m/xs_m interpolants " *
@@ -466,34 +394,6 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur
     end
 end
 
-"""
-    _jbb_deweight!(out, jbb_modes, ft, psi, equil, mthsurf, theta_buf)
-
-JBB deweighting step: inverse DFT → divide by J(ψ,θ)·B(ψ,θ)² → forward DFT.
-
-Matches Fortran set_peq lines 859-868: transforms JBB-weighted m-space data
-to θ-space, removes the J·B² weighting at each poloidal angle, and transforms back.
-"""
-function _jbb_deweight!(out::AbstractVector{ComplexF64}, jbb_modes::Vector{ComplexF64},
-                        ft::Utilities.FourierTransforms.FourierTransform,
-                        psi::Float64, equil::Equilibrium.PlasmaEquilibrium,
-                        mthsurf::Int, theta_buf::Vector{ComplexF64})
-    # Inverse DFT: m-space → θ-space
-    theta_buf .= Utilities.FourierTransforms.inverse(ft, jbb_modes)
-
-    # Divide by J·B² at each θ point
-    for j in 1:mthsurf
-        theta_norm = (j - 1) / mthsurf   # θ ∈ [0, 1)
-        pt = (psi, theta_norm)
-        jac = equil.rzphi_jac(pt)
-        B = equil.eqfun_B(pt)
-        theta_buf[j] /= jac * B^2
-    end
-
-    # Forward DFT: θ-space → m-space
-    out .= ft(theta_buf)
-    return nothing
-end
 
 
 # ============================================================================
diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl
index 902cb259c..75f527ae2 100644
--- a/src/KineticForces/Output.jl
+++ b/src/KineticForces/Output.jl
@@ -8,7 +8,7 @@ then write to gpec.h5 in a single pass.
 
 """
     write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spline=nothing,
-                   species_label=nothing, perturbation=nothing)
+                   species_label=nothing)
 
 Write KineticForces results to the "KineticForces" group in gpec.h5.
 
@@ -18,16 +18,13 @@ Write KineticForces results to the "KineticForces" group in gpec.h5.
 - `dVdpsi_spline`: Optional dV/dψ_N profile interpolant; when given, dV/dψ_N is
   written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ)
   is directly available
-- `perturbation`: optional `KineticForcesInternal`; its δB/B and ∇·ξ⊥ mode tables are written
-  once to `KineticForces/Perturbation/`. The perturbation is a property of the run, not of a
-  species, so it is written only on the call that has no `species_label`.
 - `species_label`: `nothing` writes the run total to `KineticForces/<method>/`; a label
   (e.g. `"ion_z1_m2"`, `"electron"`) writes one species' contribution to
   `KineticForces/PerSpecies/<label>/<method>/`. A multi-species run calls this once per
   species and once for the summed total, so the group is opened-or-created each time.
 """
 function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spline=nothing,
-    species_label::Union{Nothing,AbstractString}=nothing, perturbation=nothing)
+    species_label::Union{Nothing,AbstractString}=nothing)
     root = haskey(h5file, "KineticForces") ? h5file["KineticForces"] : create_group(h5file, "KineticForces")
     g = root
     if species_label !== nothing
@@ -97,31 +94,8 @@ function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; dVdpsi_spl
     for method_name in method_names
         Utilities.HDF5Annotations.annotate!(g[method_name], KF_METHOD_H5_ANNOTATIONS)
     end
-
-    # The perturbation driving the torque: one table per run, not per species or per method.
-    if species_label === nothing && perturbation !== nothing &&
-       !isempty(perturbation.perturbation_psi) && !haskey(root, "Perturbation")
-        pg = create_group(root, "Perturbation")
-        pg["psi"] = perturbation.perturbation_psi
-        pg["m"] = collect(perturbation.mlow:perturbation.mhigh)
-        pg["n"] = perturbation.nlow
-        pg["dBoB_lagrangian"] = perturbation.dbob_m_data
-        pg["div_xi_perp"] = perturbation.divx_m_data
-        Utilities.HDF5Annotations.annotate!(pg, KF_PERTURBATION_H5_ANNOTATIONS)
-    end
 end
 
-# Metadata for KineticForces/Perturbation/. δB/B and ∇·ξ⊥ are the JBB-deweighted physical
-# spectra built in `set_perturbation_data!` (Fortran pentrc set_peq), in the same poloidal
-# angle and sign convention as the PerturbedEquilibrium field spectra.
-const KF_PERTURBATION_H5_ANNOTATIONS = [
-    "psi" => (; long_name="normalized poloidal flux ψ_N of the perturbation grid", scale="psi"),
-    "m" => (; long_name="poloidal mode number m", scale="m"),
-    "n" => (; long_name="toroidal mode number n of the perturbation"),
-    "dBoB_lagrangian" => (; long_name="Lagrangian δ|B|/B Fourier modes, δB/B(θ) = Σ_m value·exp(2πi·m·θ) with θ ∈ [0,1); dimensionless", dims=("psi", "m"), attach=(1 => "psi", 2 => "m")),
-    "div_xi_perp" => (; long_name="∇·ξ_⊥ Fourier modes, same angle and sign convention as dBoB_lagrangian", units="1/m", dims=("psi", "m"), attach=(1 => "psi", 2 => "m"))
-]
-
 # Storage slice → Logan 2015 matrix letter (Eqs 7.30-7.35); see BounceAveraging.jl
 # for the Hermitian/full packing details.
 const _KINETIC_MATRIX_LETTERS = (("A", 1), ("B", 2), ("C", 3), ("D", 4), ("E", 5), ("H", 6))
diff --git a/src/PerturbedEquilibrium/FieldReconstruction.jl b/src/PerturbedEquilibrium/FieldReconstruction.jl
index c768e3473..310470582 100644
--- a/src/PerturbedEquilibrium/FieldReconstruction.jl
+++ b/src/PerturbedEquilibrium/FieldReconstruction.jl
@@ -1055,3 +1055,132 @@ function _apply_rzphi_transform(
 
     return R_modes, Z_modes, phi_modes
 end
+
+"""
+    compute_modb_modes(xi_modes, psi_grid, equil, ffs, metric) -> (; dbob, divx)
+
+Lagrangian δ|B|/B and ∇·ξ_⊥ Fourier modes, `[npsi, mpert]`, on `psi_grid`.
+
+These are the two spectra a neoclassical-transport calculation is driven by: NTV is a
+functional of the Lagrangian perturbation of the field strength, not of the normal field.
+They are a property of the perturbed equilibrium alone — no kinetic quantity enters — so they
+are built here and merely consumed by `KineticForces`, which is what lets another NTV code be
+run on GPEC's perturbation without running ours. The Fortran parent writes them from the same
+stage (`gpec` `gpout_pmodb` under `pmodb_flag`, read back by `pentrc` `read_pmodb`); the
+algorithm below is the one `pentrc` `set_peq` uses on the Clebsch displacements.
+
+The JBB deweighting (Fortran `pentrc/inputs.f90:828-868`):
+
+ 1. apply the geometric matrices S,T,X,Y,Z in m-space, giving the J·B²-weighted κ·ξ and ∇·ξ_⊥
+ 2. inverse DFT to θ-space
+ 3. divide by J·B² at each θ
+ 4. forward DFT back to m-space
+
+Only the two independent spectra are returned. The curvature term follows from
+δB/B = −(∇·ξ_⊥ + κ·ξ), and the Eulerian δ|B|/B differs from the Lagrangian by ξ·∇|B|/|B|,
+which `Response/xi_clebsch_psi`, `Response/xi_clebsch_alpha` and the equilibrium geometry
+already determine — so neither is written.
+
+Single toroidal mode only: `build_kinetic_metric_matrices` assembles mpert×mpert matrices at
+`n = ffs.nlow`.
+"""
+function compute_modb_modes(xi_modes::NamedTuple, psi_grid::Vector{Float64},
+    equil::Equilibrium.PlasmaEquilibrium,
+    ffs::ForceFreeStatesResult,
+    metric::MetricData
+)
+    npsi = length(psi_grid)
+    mpert = ffs.mpert
+
+    # Geometric matrices (S,T,X,Y,Z) of the JBB-weighted Clebsch-to-physical transformation.
+    geom_mats = ForceFreeStates.build_kinetic_metric_matrices(equil, ffs, metric)
+
+    # The deweighting round-trip runs on the equilibrium's own poloidal grid.
+    mthsurf = length(equil.rzphi_ys) - 1
+    ft = Utilities.FourierTransforms.FourierTransform(mthsurf, mpert, ffs.mlow)
+
+    dbob = zeros(ComplexF64, npsi, mpert)
+    divx = zeros(ComplexF64, npsi, mpert)
+
+    smat_flat = Vector{ComplexF64}(undef, mpert^2)
+    tmat_flat = Vector{ComplexF64}(undef, mpert^2)
+    xmat_flat = Vector{ComplexF64}(undef, mpert^2)
+    ymat_flat = Vector{ComplexF64}(undef, mpert^2)
+    zmat_flat = Vector{ComplexF64}(undef, mpert^2)
+    jbb_kapx = Vector{ComplexF64}(undef, mpert)
+    jbb_divx = Vector{ComplexF64}(undef, mpert)
+    jbb_dbob = Vector{ComplexF64}(undef, mpert)
+    theta_buf = Vector{ComplexF64}(undef, mthsurf)
+
+    # Serial over ψ: the rows are independent, but the sticky bracket hints below are shared
+    # mutable state, and a threaded reassociation here would make the spectra depend on the
+    # thread count.
+    hint_s = Ref(1)
+    hint_t = Ref(1)
+    hint_x = Ref(1)
+    hint_y = Ref(1)
+    hint_z = Ref(1)
+
+    for ipsi in 1:npsi
+        psi = psi_grid[ipsi]
+
+        xsp = view(xi_modes.clebsch_psi, ipsi, :)        # ξ^ψ [mpert]
+        xmp1 = view(xi_modes.clebsch_psi1, ipsi, :)      # ∂ξ^ψ/∂ψ [mpert]
+        xms = view(xi_modes.clebsch_alpha, ipsi, :)      # ξ^α/χ₁ [mpert]
+
+        geom_mats.smats(smat_flat, psi; hint=hint_s)
+        geom_mats.tmats(tmat_flat, psi; hint=hint_t)
+        geom_mats.xmats(xmat_flat, psi; hint=hint_x)
+        geom_mats.ymats(ymat_flat, psi; hint=hint_y)
+        geom_mats.zmats(zmat_flat, psi; hint=hint_z)
+
+        smat = reshape(smat_flat, mpert, mpert)
+        tmat = reshape(tmat_flat, mpert, mpert)
+        xmat = reshape(xmat_flat, mpert, mpert)
+        ymat = reshape(ymat_flat, mpert, mpert)
+        zmat = reshape(zmat_flat, mpert, mpert)
+
+        # jbb_kapx = smat·xsp + tmat·xms;  jbb_divx = xmat·xmp1 + ymat·xsp + zmat·xms
+        mul!(jbb_kapx, smat, xsp)
+        mul!(jbb_kapx, tmat, xms, 1.0 + 0.0im, 1.0 + 0.0im)
+        mul!(jbb_divx, xmat, xmp1)
+        mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im)
+        mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im)
+        @. jbb_dbob = -(jbb_divx + jbb_kapx)
+
+        _jbb_deweight!(view(dbob, ipsi, :), jbb_dbob, ft, psi, equil, mthsurf, theta_buf)
+        _jbb_deweight!(view(divx, ipsi, :), jbb_divx, ft, psi, equil, mthsurf, theta_buf)
+    end
+
+    return (; dbob, divx)
+end
+
+"""
+    _jbb_deweight!(out, jbb_modes, ft, psi, equil, mthsurf, theta_buf)
+
+JBB deweighting step: inverse DFT → divide by J(ψ,θ)·B(ψ,θ)² → forward DFT.
+
+Matches Fortran `set_peq` lines 859-868: transforms JBB-weighted m-space data to θ-space,
+removes the J·B² weighting at each poloidal angle, and transforms back.
+"""
+function _jbb_deweight!(out::AbstractVector{ComplexF64}, jbb_modes::Vector{ComplexF64},
+    ft::Utilities.FourierTransforms.FourierTransform,
+    psi::Float64, equil::Equilibrium.PlasmaEquilibrium,
+    mthsurf::Int, theta_buf::Vector{ComplexF64}
+)
+    # Inverse DFT: m-space → θ-space
+    theta_buf .= Utilities.FourierTransforms.inverse(ft, jbb_modes)
+
+    # Divide by J·B² at each θ point
+    for j in 1:mthsurf
+        theta_norm = (j - 1) / mthsurf   # θ ∈ [0, 1)
+        pt = (psi, theta_norm)
+        jac = equil.rzphi_jac(pt)
+        B = equil.eqfun_B(pt)
+        theta_buf[j] /= jac * B^2
+    end
+
+    # Forward DFT: θ-space → m-space
+    out .= ft(theta_buf)
+    return nothing
+end
diff --git a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
index c433733db..9e7cc9b10 100644
--- a/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
+++ b/src/PerturbedEquilibrium/PerturbedEquilibriumStructs.jl
@@ -96,6 +96,8 @@ Response fields (mode space):
   - `b_modes::Union{Nothing, NamedTuple}` - Magnetic field; psi=b^ψ, b_psi_area_weighted=b^ψ/⟨J·|∇ψ|⟩_θ, theta/zeta=unregularized, theta_reg/zeta_reg=regularized [npsi, numpert_total]
   - `b_n_modes::Union{Nothing, Matrix{ComplexF64}}` - Physical normal field b_n [npsi, numpert_total]
   - `xi_n_modes::Union{Nothing, Matrix{ComplexF64}}` - Physical normal displacement xi_n [npsi, numpert_total]
+  - `modb_modes::Union{Nothing, NamedTuple}` - Perturbation spectra driving neoclassical transport;
+    dbob=Lagrangian δ|B|/B, divx=∇·ξ_⊥ [npsi, mpert]. Single-n runs only
 
 Coupling matrices [n_rational × numpert_total] — one row per resonant (surface, n) pair.
 Each row maps the full applied field to the resonant response at that surface.
@@ -173,6 +175,8 @@ well-conditioned flux-space inductances L, Λ:
     b_modes::Union{Nothing,NamedTuple} = nothing
     b_n_modes::Union{Nothing,Matrix{ComplexF64}} = nothing  # physical normal field b_n [npsi, numpert_total]
     xi_n_modes::Union{Nothing,Matrix{ComplexF64}} = nothing  # physical normal displacement xi_n [npsi, numpert_total]
+    # Lagrangian δ|B|/B and ∇·ξ_⊥ spectra (see compute_modb_modes); the input an NTV code needs
+    modb_modes::Union{Nothing,NamedTuple} = nothing
 
     # Coupling matrices [n_rational × numpert_total]
     C_resonant_area_weighted_field::Matrix{ComplexF64} = zeros(ComplexF64, 0, 0)
diff --git a/src/PerturbedEquilibrium/Response.jl b/src/PerturbedEquilibrium/Response.jl
index 2100e52c1..7a944e603 100644
--- a/src/PerturbedEquilibrium/Response.jl
+++ b/src/PerturbedEquilibrium/Response.jl
@@ -120,7 +120,21 @@ function compute_plasma_response!(
     state.b_n_modes = b_n_modes
     state.xi_n_modes = xi_n_modes
 
+    # Lagrangian δ|B|/B and ∇·ξ_⊥: the perturbation a neoclassical-transport calculation is
+    # driven by. Derived here, where the displacement is, rather than inside KineticForces --
+    # it contains nothing kinetic, and another NTV code needs it without running ours.
+    # Restricted to a single toroidal mode by the geometric matrices it is built from.
+    modb_dt = 0.0
+    if ffs.npert == 1
+        modb_dt = @elapsed state.modb_modes = compute_modb_modes(
+            xi_modes, state.psi_grid, equil, ffs, metric
+        )
+    end
+
     if ctrl.verbose
-        @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = $(@sprintf("%.3e", maximum(abs.(response_flux))))"
+        @info "Response complete: $(length(intr.forcing_modes)) forcing modes, max amplitude = " *
+              "$(@sprintf("%.3e", maximum(abs.(response_flux))))" *
+              (state.modb_modes === nothing ? "; no δ|B|/B spectra (npert > 1)" :
+               "; δ|B|/B and ∇·ξ_⊥ spectra in $(@sprintf("%.2f", modb_dt)) s")
     end
 end
diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl
index d2064aee4..5d78bd997 100644
--- a/src/PerturbedEquilibrium/Utils.jl
+++ b/src/PerturbedEquilibrium/Utils.jl
@@ -70,7 +70,9 @@ PerturbedEquilibrium/
 │   ├── b_n            # Physical normal field b_n (ComplexF64 [npsi, numpert_total])
 │   ├── xi_n           # Physical normal displacement xi_n (ComplexF64 [npsi, numpert_total])
 │   ├── Jb_theta
-│   └── Jb_zeta
+│   ├── Jb_zeta
+│   ├── dBoB_lagrangian  # Lagrangian δ|B|/B (ComplexF64 [npsi, numpert_total]) -- NTV input
+│   └── div_xi_perp      # ∇·ξ_⊥ (ComplexF64 [npsi, numpert_total]) -- NTV input
 ├── ResponseMatrices/         # [numpert_total × numpert_total], root-area-weighted field (b̃) space; R = S·A
 │   ├── plasma_inductance     # Λ̃ = R⁻¹·Λ·R⁻†
 │   ├── surface_inductance    # L̃ = R⁻¹·L·R⁻†
@@ -166,6 +168,14 @@ function write_outputs_to_HDF5(
             response_group["xi_clebsch_alpha"] = state.xi_modes.clebsch_alpha
         end
 
+        # Lagrangian δ|B|/B and ∇·ξ_⊥ (gpout_pmodb): the perturbation a neoclassical-transport
+        # calculation is driven by, written whether or not this run also computes the torque.
+        # Empty for npert > 1, which the geometric matrices behind them do not support.
+        if !isnothing(state.modb_modes)
+            response_group["dBoB_lagrangian"] = state.modb_modes.dbob
+            response_group["div_xi_perp"] = state.modb_modes.divx
+        end
+
         # Contravariant displacement (from gpeq_contra, all J-weighted)
         if have_xi
             response_group["Jxi_psi"] = state.xi_modes.psi_J
@@ -291,6 +301,14 @@ const PE_H5_ANNOTATIONS = [
         (; long_name="regularized ψ_N derivative of ξ^ψ (× singfac²/(singfac²+regularization_width²))", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
     "Response/xi_clebsch_alpha" =>
         (; long_name="Clebsch displacement component ξ^α/χ₁ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
+    "Response/dBoB_lagrangian" =>
+        (; long_name="Lagrangian perturbation of the field strength δ|B|/B (gpout_pmodb lagbpar/|B|), the NTV drive; " *
+            "δB/B(θ) = Σ_m value·exp(2πi·m·θ), θ ∈ [0,1). Eulerian form = this − ξ·∇|B|/|B| from xi_clebsch_*",
+            dims=("psi", "mode"), attach=(1 => "Response/psi",)),
+    "Response/div_xi_perp" =>
+        (; long_name="∇·ξ_⊥ (gpout_pmodb divxprp/|B|), the second NTV input; same angle and sign convention as " *
+            "dBoB_lagrangian, with κ·ξ = −(dBoB_lagrangian + div_xi_perp)", units="1/m",
+            dims=("psi", "mode"), attach=(1 => "Response/psi",)),
     "Response/xi_n" => (; long_name="physical normal displacement ξ_n", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
     "Response/xi_R" => (; long_name="cylindrical displacement component ξ_R (mode space)", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
     "Response/xi_Z" => (; long_name="cylindrical displacement component ξ_Z (mode space)", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
diff --git a/test/runtests_h5_schema.jl b/test/runtests_h5_schema.jl
index c743d6a7d..11a6788c2 100644
--- a/test/runtests_h5_schema.jl
+++ b/test/runtests_h5_schema.jl
@@ -105,3 +105,42 @@ end
         end
     end
 end
+
+@testset "PerturbedEquilibrium δ|B|/B output contract" begin
+    PE = GeneralizedPerturbedEquilibrium.PerturbedEquilibrium
+    npsi, mpert = 4, 3
+
+    # Written when the spectra exist. The NTV codes this serves read them as (m, ψ) — the
+    # transpose of the Julia shape — exactly like every other Response spectrum.
+    state = PE.PerturbedEquilibriumState(;
+        psi_grid=collect(range(0.1, 0.9; length=npsi)),
+        modb_modes=(dbob=ComplexF64[0.1i + 0.01im * j for i in 1:npsi, j in 1:mpert],
+                    divx=ComplexF64[0.2i - 0.02im * j for i in 1:npsi, j in 1:mpert]))
+    f1 = tempname() * ".h5"
+    PE.write_outputs_to_HDF5(state, PE.PerturbedEquilibriumInternal(), f1)
+    HDF5.h5open(f1, "r") do f
+        g = f["PerturbedEquilibrium/Response"]
+        @test size(read(g["dBoB_lagrangian"])) == (npsi, mpert)
+        @test read(g["dBoB_lagrangian"]) ≈ state.modb_modes.dbob
+        @test read(g["div_xi_perp"]) ≈ state.modb_modes.divx
+        for d in ("dBoB_lagrangian", "div_xi_perp")
+            a = HDF5.attributes(g[d])
+            @test haskey(a, "long_name") && haskey(a, "dims")
+            # Same ψ dimension scale as every other Response spectrum, which is what lets
+            # Analysis.PerturbedEquilibriumModes.modes_to_theta read these datasets unchanged.
+            @test HDF5.API.h5ds_is_attached(g[d], g["psi"], 0)
+        end
+        @test size(read(g["dBoB_lagrangian"])) == size(read(g["xi_clebsch_psi"])) ||
+              isempty(read(g["xi_clebsch_psi"]))
+    end
+    rm(f1; force=true)
+
+    # Absent, not empty, when the run could not produce them (npert > 1).
+    f2 = tempname() * ".h5"
+    PE.write_outputs_to_HDF5(PE.PerturbedEquilibriumState(), PE.PerturbedEquilibriumInternal(), f2)
+    HDF5.h5open(f2, "r") do f
+        @test !haskey(f["PerturbedEquilibrium/Response"], "dBoB_lagrangian")
+        @test !haskey(f["PerturbedEquilibrium/Response"], "div_xi_perp")
+    end
+    rm(f2; force=true)
+end
diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl
index 59db7b58e..4d141f626 100644
--- a/test/runtests_kinetic.jl
+++ b/test/runtests_kinetic.jl
@@ -1,4 +1,3 @@
-using HDF5
 
 @testset "KineticForces Unit Tests" begin
 
@@ -667,57 +666,6 @@ using HDF5
         @test isempty(mr.records)
     end
 
-    @testset "perturbation output" begin
-        # KineticForces/Perturbation/ carries the δB/B and ∇·ξ⊥ spectra that drive the torque,
-        # so another NTV code can be run on the same perturbation. It is a property of the run,
-        # not of a species or a method, so it is written once by the call with no species_label.
-        @test isempty(KF.KineticForcesInternal().perturbation_psi)
-        @test isempty(KF.KineticForcesInternal().dbob_m_data)
-
-        npsi, mpert = 4, 3
-        intr = KF.KineticForcesInternal(; mlow=-1, mhigh=1, mpert=mpert, nlow=3, nhigh=3,
-            perturbation_psi=collect(range(0.1, 0.9; length=npsi)),
-            dbob_m_data=ComplexF64[0.1 * i + 0.01im * j for i in 1:npsi, j in 1:mpert],
-            divx_m_data=ComplexF64[0.2 * i - 0.02im * j for i in 1:npsi, j in 1:mpert])
-        st = KF.KineticForcesState()
-        st.method_results["fgar"] = KF.MethodResult(; method="fgar", nn=3,
-            psi_grid=collect(0.0:0.25:1.0), dtdpsi=fill(ComplexF64(1.0), 5),
-            t_cumulative=zeros(ComplexF64, 5))
-
-        h5 = tempname() * ".h5"
-        HDF5.h5open(h5, "cw") do f
-            KF.write_to_hdf5!(f, st; species_label="electron", perturbation=intr)  # must NOT write it
-            KF.write_to_hdf5!(f, st; perturbation=intr)
-        end
-        HDF5.h5open(h5, "r") do f
-            @test !haskey(f["KineticForces"]["PerSpecies"]["electron"], "Perturbation")
-            pg = f["KineticForces"]["Perturbation"]
-            @test read(pg["m"]) == [-1, 0, 1]
-            @test read(pg["n"]) == 3
-            @test read(pg["psi"]) ≈ intr.perturbation_psi
-            # HDF5.jl round-trips the Julia shape, so (npsi, mpert) here; h5py and MATLAB see
-            # the transpose, matching the field spectra written by PerturbedEquilibrium.
-            @test size(read(pg["dBoB_lagrangian"])) == (npsi, mpert)
-            @test read(pg["dBoB_lagrangian"]) ≈ intr.dbob_m_data
-            @test read(pg["div_xi_perp"]) ≈ intr.divx_m_data
-            for d in ("dBoB_lagrangian", "div_xi_perp")
-                a = HDF5.attributes(pg[d])
-                @test haskey(a, "long_name") && haskey(a, "dims")
-            end
-        end
-        rm(h5; force=true)
-
-        # Nothing is written when the perturbation was never built (PE did not run).
-        h5b = tempname() * ".h5"
-        HDF5.h5open(h5b, "cw") do f
-            KF.write_to_hdf5!(f, st; perturbation=KF.KineticForcesInternal())
-        end
-        HDF5.h5open(h5b, "r") do f
-            @test !haskey(f["KineticForces"], "Perturbation")
-        end
-        rm(h5b; force=true)
-    end
-
     # =========================================================================
     # Kinetic profile I/O: HDF5 schema, ASCII back-compat, writer round-trip
     # =========================================================================

From 9016e64317ca6fd98c24a35e0ac584af9ba55489 Mon Sep 17 00:00:00 2001
From: logan-nc <nikolas.logan@columbia.edu>
Date: Tue, 20 Oct 2026 21:15:35 -0400
Subject: [PATCH 4/5] Test - TEST - Fix the perturbation-output assertions

Two of my own mistakes, caught by running them. The full-run block asserted the new
Response datasets appear in regression_solovev_ideal_example, which has no
[PerturbedEquilibrium] section at all -- no test deck does, so that assertion was
unsatisfiable by construction and the computation is covered by the harness case
solovev_kinetic_ntv instead. And the dimension-scale check used C index 0: H5DS
indexes file (row-major) dimensions, so Julia axis 1 of a 2-D dataset is C index 1,
per Utilities.HDF5Annotations.attach_scale!.

Kept as its own commit so the harness stamp from the parent still applies: src/ is
untouched here.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FabuzB6rkkv638ZUCzaAA2
---
 test/runtests_h5_schema.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/test/runtests_h5_schema.jl b/test/runtests_h5_schema.jl
index 11a6788c2..35c0478fd 100644
--- a/test/runtests_h5_schema.jl
+++ b/test/runtests_h5_schema.jl
@@ -128,10 +128,10 @@ end
             @test haskey(a, "long_name") && haskey(a, "dims")
             # Same ψ dimension scale as every other Response spectrum, which is what lets
             # Analysis.PerturbedEquilibriumModes.modes_to_theta read these datasets unchanged.
-            @test HDF5.API.h5ds_is_attached(g[d], g["psi"], 0)
+            # H5DS indexes file (row-major) dimensions, so Julia axis 1 of a 2-D dataset is
+            # C index 1 — see Utilities.HDF5Annotations.attach_scale!.
+            @test HDF5.API.h5ds_is_attached(g[d], g["psi"], 1)
         end
-        @test size(read(g["dBoB_lagrangian"])) == size(read(g["xi_clebsch_psi"])) ||
-              isempty(read(g["xi_clebsch_psi"]))
     end
     rm(f1; force=true)
 

From 136b0fb9f7d05ec7aa24463b99e0b3716ac67943 Mon Sep 17 00:00:00 2001
From: logan-nc <nikolas.logan@columbia.edu>
Date: Tue, 20 Oct 2026 22:39:47 -0400
Subject: [PATCH 5/5] PerturbedEquilibrium/KineticForces - MINOR - Format this
 PR's own new code

Runs the repository's pinned JuliaFormatter (1.0.62, the pre-commit rev) over
the lines this PR introduces: the relocated docstring's list, the two annotation
entries and one test continuation. The format workflow is not a PR gate -- it is
disabled pending a repo-wide cleanup -- so the three pre-existing hunks the
formatter also wanted in FieldReconstruction.jl and Response.jl are left alone
rather than carried as churn in a feature PR (it would have rewritten
`((in - 1) * mpert + 1)` as `((in-1)*mpert+1)`, which is not an improvement).

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01FabuzB6rkkv638ZUCzaAA2
---
 src/KineticForces/KineticForcesStructs.jl | 11 ++++++-----
 src/PerturbedEquilibrium/Utils.jl         | 10 ++++++----
 test/runtests_h5_schema.jl                |  2 +-
 3 files changed, 13 insertions(+), 10 deletions(-)

diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl
index 2dc594397..2711ab9ef 100644
--- a/src/KineticForces/KineticForcesStructs.jl
+++ b/src/KineticForces/KineticForcesStructs.jl
@@ -324,16 +324,17 @@ end
 Populate perturbation data from PerturbedEquilibriumState into KineticForcesInternal.
 
 Builds three interpolant sets over ψ:
-1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] from the PE Clebsch displacements
-2. `dbob_m` — Lagrangian δ|B|/B, taken from `pe_state.modb_modes`
-3. `divx_m` — ∇·ξ⊥, taken from `pe_state.modb_modes`
+
+ 1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] from the PE Clebsch displacements
+ 2. `dbob_m` — Lagrangian δ|B|/B, taken from `pe_state.modb_modes`
+ 3. `divx_m` — ∇·ξ⊥, taken from `pe_state.modb_modes`
 
 The perturbation spectra themselves are a property of the perturbed equilibrium and are
 derived and written there (`PerturbedEquilibrium.compute_modb_modes`); this function only
 interpolates them onto the torque quadrature.
 """
 function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::PerturbedEquilibrium.PerturbedEquilibriumState,
-                                ffs::ForceFreeStates.ForceFreeStatesResult)
+    ffs::ForceFreeStates.ForceFreeStatesResult)
     ffs.npert == 1 || error("KineticForces supports a single toroidal mode; this run has n = $(ffs.nlow):$(ffs.nhigh).")
 
     # Copy mode numbers from FFS
@@ -380,7 +381,7 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur
     modb === nothing && error("set_perturbation_data!: the perturbed equilibrium carries no " *
                               "δ|B|/B spectra; KineticForces cannot evaluate the torque without them.")
     size(modb.dbob) == (npsi, mpert) || error("set_perturbation_data!: δ|B|/B spectra are " *
-        "$(size(modb.dbob)), expected ($npsi, $mpert).")
+                                              "$(size(modb.dbob)), expected ($npsi, $mpert).")
 
     # Build CubicSeriesInterpolants over ψ for dbob_m and divx_m
     kf_intr.dbob_m = cubic_interp(psi_grid, Series(modb.dbob); itp_opts...)
diff --git a/src/PerturbedEquilibrium/Utils.jl b/src/PerturbedEquilibrium/Utils.jl
index 5d78bd997..7549517a8 100644
--- a/src/PerturbedEquilibrium/Utils.jl
+++ b/src/PerturbedEquilibrium/Utils.jl
@@ -302,12 +302,14 @@ const PE_H5_ANNOTATIONS = [
     "Response/xi_clebsch_alpha" =>
         (; long_name="Clebsch displacement component ξ^α/χ₁ (PENTRC input, gpout_xclebsch convention)", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
     "Response/dBoB_lagrangian" =>
-        (; long_name="Lagrangian perturbation of the field strength δ|B|/B (gpout_pmodb lagbpar/|B|), the NTV drive; " *
-            "δB/B(θ) = Σ_m value·exp(2πi·m·θ), θ ∈ [0,1). Eulerian form = this − ξ·∇|B|/|B| from xi_clebsch_*",
+        (;
+            long_name="Lagrangian perturbation of the field strength δ|B|/B (gpout_pmodb lagbpar/|B|), the NTV drive; " *
+                      "δB/B(θ) = Σ_m value·exp(2πi·m·θ), θ ∈ [0,1). Eulerian form = this − ξ·∇|B|/|B| from xi_clebsch_*",
             dims=("psi", "mode"), attach=(1 => "Response/psi",)),
     "Response/div_xi_perp" =>
-        (; long_name="∇·ξ_⊥ (gpout_pmodb divxprp/|B|), the second NTV input; same angle and sign convention as " *
-            "dBoB_lagrangian, with κ·ξ = −(dBoB_lagrangian + div_xi_perp)", units="1/m",
+        (;
+            long_name="∇·ξ_⊥ (gpout_pmodb divxprp/|B|), the second NTV input; same angle and sign convention as " *
+                      "dBoB_lagrangian, with κ·ξ = −(dBoB_lagrangian + div_xi_perp)", units="1/m",
             dims=("psi", "mode"), attach=(1 => "Response/psi",)),
     "Response/xi_n" => (; long_name="physical normal displacement ξ_n", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
     "Response/xi_R" => (; long_name="cylindrical displacement component ξ_R (mode space)", units="m", dims=("psi", "mode"), attach=(1 => "Response/psi",)),
diff --git a/test/runtests_h5_schema.jl b/test/runtests_h5_schema.jl
index 35c0478fd..cecbc464a 100644
--- a/test/runtests_h5_schema.jl
+++ b/test/runtests_h5_schema.jl
@@ -115,7 +115,7 @@ end
     state = PE.PerturbedEquilibriumState(;
         psi_grid=collect(range(0.1, 0.9; length=npsi)),
         modb_modes=(dbob=ComplexF64[0.1i + 0.01im * j for i in 1:npsi, j in 1:mpert],
-                    divx=ComplexF64[0.2i - 0.02im * j for i in 1:npsi, j in 1:mpert]))
+            divx=ComplexF64[0.2i - 0.02im * j for i in 1:npsi, j in 1:mpert]))
     f1 = tempname() * ".h5"
     PE.write_outputs_to_HDF5(state, PE.PerturbedEquilibriumInternal(), f1)
     HDF5.h5open(f1, "r") do f