diff --git a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md
index 1a3bc1fc1..9d17a0bdd 100644
--- a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md
+++ b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md
@@ -5,3 +5,4 @@
 - [KineticForces (NTV) Audit Checklist](kinetic_ntv_map.md) — pentrc->KineticForces map, Logan 2015 matrices, what to verify
 - [InnerLayer (Resistive) Audit Checklist](resistive_layer_map.md) — rmatch->InnerLayer map, GGJ Wasow basis / Δ′, what to verify
 - [Galerkin Δ′ Assembly Map](galerkin_assembly_map.md) — gal.f<->GalerkinAssembly.jl; resonant sign chain verified; PASS
+- [reg_spot vs singfac_min audit](reg_spot_regularization.md) — reg_spot field-recon smoothing CORRECT; distinct from singfac_min ODE gate
diff --git a/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md b/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md
index 5f5dfb564..c1724bd49 100644
--- a/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md
+++ b/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md
@@ -31,3 +31,60 @@ the Fortran `pentrc/` sources (file mapping in fortran_correspondence_map.md).
 - **Normalization**: torque Im(T)=2n·δW_k (Eq. 19) and diamagnetic-frequency sign/factor conventions match Fortran `torque.F90`.
 - **Mode indexing**: m, n ranges and block packing over n stay consistent with ForceFreeStates.
 - **Method variants**: FGAR/TGAR/PGAR/RLAR/CLAR/FCGL/TMM/WMM each present, not stubbed to a single fallback.
+
+## Multi-ion (D-T) NTV — composition/collisionality (2026-07)
+Both Fortran PENTRC (`inputs.f90:236-243`, `read_kin`) and Julia (`KineticProfiles.jl:261-278`)
+support ONE main ion (zi,mi) + ONE impurity (zimp,mimp) per run. Correct multi-main-ion NTV is
+an EXTENSION beyond Fortran, but consistent with Logan-Park 2013's pitch-angle (Lorentz) model.
+- Zeff = Σ_s Z_s² n_s / n_e; quasineutrality n_e = Σ_s Z_s n_s (all ions incl. impurity).
+- Bug in current D-T split (run D, then T, each with ni=Ni/2): the `ni` column sets BOTH species
+  density AND Zeff via z = zimp-(n_i/n_e)zi(zimp-zi). With ni=Ni/2 it treats the missing half as
+  high-Z impurity → Zeff≈3.7 instead of true ≈1–1.4. Corrupts zpitch → nue,nui.
+- Correct collisionality: ν_a ∝ Z_a² lnΛ · Σ_b n_b Z_b² /(√m_a T_a^{3/2}) = ∝ Z_a² lnΛ·(n_e·Zeff_true).
+  The zpitch(Zeff) polynomial is a main-ion+impurity closure (momentum-restoring correction); the
+  MINIMAL fix is to feed the TRUE (full-composition) Zeff into the existing zpitch/ν formulas.
+- Additivity: τ = Σ_s τ_s. Lorentz operator is additive over field species; species couple only
+  through shared δB, shared ω_E, shared Zeff/ν. Additive at this theory's order.
+- What changes single→multi: Zeff, zpitch, nue, nui (→nueff). UNCHANGED for equal-shape D-T:
+  wdian/wdiat (log-derivative, density factor cancels), wtran/wbhat/wdhat/wgyro (already per-species),
+  and the resonant-density prefactor (Ni/2 per species is correct).
+- ASIDE (separate fidelity bug, same code block): Fortran `inputs.f90:238` uses natural log for lnΛ;
+  Julia `KineticProfiles.jl:270` uses log10 — diverges away from the n=1e20,T=1keV reference point.
+
+## Multi-ion NTV — full-composition species set (2026-07 audit, PASS-with-caveats)
+Reviewed `resolve_ntv_species` (KineticProfiles.jl ~214-266) + `compute_calculated_kinetic_matrices`
+(CalculatedKineticMatrices.jl ~96-151). Verdict: physics is sound.
+- ν_s field-density RECONCILIATION (supersedes the earlier "should be n_e·Zeff" note above):
+  code uses `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ/(√m_s·T_i^1.5)`, i.e. field density = zpitch·n_main
+  (n_main = Σ MAIN-ion densities, no impurity, unweighted). This is CORRECT and MORE faithful to
+  single-ion PENTRC than n_e·Zeff: PENTRC's design is zpitch·n_i, NOT Σ_b n_b Z_b². The two are NOT
+  numerically equal (Zeff=1.5, C6, n_i/n_e=0.9 → zpitch·n_i = 1.225·n_e vs n_e·Zeff = 1.5·n_e, ~18%
+  apart); that gap IS the intended momentum-restoring design difference that zpitch(Zeff) carries.
+  n_main vs n_e·Zeff is a deliberate PENTRC-fidelity choice, exact in the z=1 main-ion domain.
+- z_s² test-particle factor: CORRECT and correctly placed (deflection freq ∝ test charge²). Single-ion
+  had no z² only because zi=1. Reduces EXACTLY to single-ion nui for one z=1,fraction=1 ion (verified).
+- Impurity as its own test species: field density zpitch·n_main is NOT undercounting — the impurity's
+  z_imp² contribution is already folded into Zeff inside zpitch. CAVEAT: zpitch is strictly a main-ion
+  momentum-restoring closure; reusing it for the impurity/electron ν is an approximation beyond the
+  single-ion theory. Acceptable (impurity δW ∝ n_imp is small); worth a one-line annotation.
+- Electron descriptor passes `ns[1]` (first ion's density) as its `ni_spline` — HARMLESS: the kernel
+  `_setup_surface_state` (Torque.jl:659-664) reads `ne_spline` (full shared n_e) when electron=true,
+  never ni_spline. Electron n_s = full n_e (correct). Cosmetic smell only.
+- Self-consistent δW summation (CalculatedKineticMatrices.jl:112-148): kw_flat/kt_flat are PURELY the
+  kinetic matrices (Logan 7.30-7.35), every term ∝ species phase-space density n_s·f0_s — NO species-
+  independent baseline. Fluid F,K,G added ONCE downstream in _compute_fkg_matrices!, outside the
+  species loop. So `+=` over species is clean additivity; NOTHING in kw is wrongly ×species-count.
+  → TC-24 n=3 δW +0.066(D) → −0.10(D+T+C+e) sign flip is PLAUSIBLE physics near marginal stability,
+  NOT a double-count. Dominant driver: the newly-ON electron channel (full n_e, opposite precession).
+  DECISIVE cheap diagnostic to confirm: run D(½)+T(½) with electron=OFF, impurity absent — should
+  return ≈ +0.066 (single-ion D). If D+T alone ≈ +0.066, the whole shift is electron+impurity = physical.
+
+## Single-ion nui vs multi-species z_s² (#339, 2026-07 decision)
+Fortran PENTRC `inputs.f90:240-241` single-ion nui = (zpitch/3.5e17)·n_i·lnΛ/(√mi·T_i^1.5) has
+NO explicit zi² (implicitly assumes zi=1, main ion hydrogenic). Julia `load_kinetic_profiles`
+(~L416) is a faithful exact port. Multi-species `_nu` (~L253) adds explicit test-particle z_s²
+(the physically correct pitch-angle/Lorentz form; zpitch is the field-side momentum-restoring
+factor of Zeff, independent of test charge — no double count). RECOMMENDATION: option (a) ADD
+zi² to single-ion nui. It is numerically identical for the default zi=1 (regression byte-identical),
+makes single-ion the true 1-species special case of `_nu`, and only "deviates" from Fortran in the
+exotic zi≠1 case where Fortran is physically wrong anyway. Annotate as a documented improvement.
diff --git a/.claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md b/.claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md
new file mode 100644
index 000000000..137301e28
--- /dev/null
+++ b/.claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md
@@ -0,0 +1,24 @@
+---
+name: reg_spot vs singfac_min regularization audit
+description: Verdict on GPEC-Julia reg_spot (field-reconstruction smoothing) port and its distinction from singfac_min (ODE crossing gate)
+metadata:
+  type: project
+---
+
+## Verdict (audited FieldReconstruction.jl, no local Fortran repo available)
+`reg_spot` port judged CORRECT-WITH-CAVEATS (caveat = Fortran verified from GPEC `xm*` convention + physics + known default, not from on-disk source).
+
+## reg_spot — field-reconstruction smoothing (PerturbedEquilibrium)
+- Factor form: `reg_factor = singfac²/(singfac²+reg_spot²)`, singfac = m - n·q. Correct GPEC form. → 0 at rational surface, → 1 away. Applied in singfac-space, NOT ψ-space.
+- Default 5e-2 = 0.05 (PerturbedEquilibriumStructs.jl:50). Matches GPEC gpec_input namelist default.
+- Quantities regularized (matches GPEC `xm*` modified-quantity convention):
+  - xmp1 = ξ^ψ' (FieldReconstruction.jl:393-395)
+  - xms (clebsch_alpha) computed FROM regularized xmp1 via xms=-A⁻¹(B·xmp1_reg+C·xsp) (line 398-415, gpeq_sol)
+  - xmt/xmz = regularized ξ^θ,ξ^ζ (line 584-594, gpeq_contra); xwt/xwz kept unregularized
+  - b^θ_reg,b^ζ_reg inherit reg via xmp1/xms (compute_modified_field_modes)
+- LEFT UNREGULARIZED: ξ^ψ primitive (xsp / clebsch_psi = copy, line 356; xwp uses raw xsp). CORRECT — ξ^ψ (normal disp) is finite at rational surface; only 1/singfac-divergent tangential/derivative quantities are softened.
+
+## singfac_min — DISTINCT mechanism (ForceFreeStates ODE crossing gate)
+- Default 1e-4 (ForceFreeStatesStructs.jl:267), "Matches Fortran STRIDE". Gates singular-surface crossings during EL ODE integration (EulerLagrange.jl:466-519, GeneralizedPerturbedEquilibrium.jl:426).
+- SingularCoupling.jl:204 uses `spot = 5e-4` as an offset (5e-4/|n·q1|) to evaluate bwp1 at lpsi/rpsi — labeled "matches Fortran default singfac_min". This is the ψ-offset for one-sided derivative eval, a THIRD distinct use, not reg_spot.
+- reg_spot (0.05, singfac-space smoothing of reconstructed fields) and singfac_min (1e-4, ODE crossing gate) are separate; both ported.
diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md
index 9787d9698..7f7e6f490 100644
--- a/docs/src/kinetic_forces.md
+++ b/docs/src/kinetic_forces.md
@@ -16,7 +16,10 @@ extension:
     `chi_phi` (toroidal momentum diffusivity ``\chi_\phi``). Each dataset
     carries a `units` attribute; the root carries `schema_version` and
     `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s,
-    diffusivities m²/s. Write files with `write_kinetic_h5`.
+    diffusivities m²/s. Additional datasets named `n_*` (e.g. `n_D`, `n_T`)
+    are named per-species densities for multi-ion runs; names outside `n_*`
+    are reserved for future schema fields and ignored. Write files with
+    `write_kinetic_h5` (which round-trips the per-species densities).
   - **ASCII (`.gpeckf`/`.kin`/`.dat`)** — legacy six-column whitespace table
     `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward
     compatibility. Header rows are skipped.
@@ -25,6 +28,49 @@ The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_phi`,
 when present, are carried for the resistive-layer (SLAYER) analysis and ignored
 here.
 
+## Multi-ion runs
+
+A plasma may declare an arbitrary list of main-ion species; the NTV is computed
+per species under one shared full-composition ``Z_\mathrm{eff}`` and summed
+(``\tau = \sum_s \tau_s``) over the main ions, the quasineutrality-closing
+impurity, and (with `electron = true`) the electrons. This applies to both NTV
+paths: the post-PE ψ-quadrature diagnostic and the self-consistent
+`kinetic_source = "calculated"` matrices.
+
+```toml
+[KineticForces]
+kinetic_file = "kinetic.h5"    # n_i column/dataset = TOTAL main-ion density
+electron = true                # add electron NTV in addition to the ion species
+zimp = 6                       # impurity charge (closes quasineutrality)
+mimp = 12                      # impurity mass
+
+[[KineticForces.ion_species]]
+z = 1
+m = 2
+fraction = 0.5                 # this species' share of the total n_i
+[[KineticForces.ion_species]]
+z = 1
+m = 3
+density = "n_T"                # or: explicit n_* dataset from the HDF5 kinetic file
+```
+
+Each species sets exactly one of `fraction` (share of the file's total `n_i`)
+or `density` (a named `n_*` dataset). An all-fraction list must sum to 1; in a
+mixed list the impurity content is set by the file's `n_i`/`n_e` deficit, not
+by a fraction shortfall, and fractions may sum below (never above) 1. Every
+main-ion charge must satisfy `z < zimp`. An empty `ion_species` list runs the
+single main ion from `zi`/`mi` (with `electron = true` still adding the
+electron species — the electron flag always means *in addition to* the ions).
+
+Per-species results are written to HDF5 groups `kinetic_forces_<label>`
+(e.g. `kinetic_forces_ion1_z1_m2`, `kinetic_forces_impurity_z6_m12`,
+`kinetic_forces_electron`) alongside the summed total in `kinetic_forces`. 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.
+
+The profile-scaling knobs below are not supported together with a multi-ion
+`ion_species` list (they error).
+
 ## Profile Scaling Knobs
 
 Seven scaling factors are available on `KineticForcesControl` to modify kinetic
diff --git a/examples/Solovev_kinetic_calculated_example/gpec.toml b/examples/Solovev_kinetic_calculated_example/gpec.toml
index 149f40af7..ddebd1400 100644
--- a/examples/Solovev_kinetic_calculated_example/gpec.toml
+++ b/examples/Solovev_kinetic_calculated_example/gpec.toml
@@ -1,7 +1,9 @@
 # Solovev analytical equilibrium — n=1 ideal stability with calculated kinetic matrices.
-# Exercises the real KineticForces (NTV) physics path at full strength; kinetic regression
-# fixture (et[1]≈1.894-1.525i). The equilibrium is generated analytically from the embedded
-# [SOL_INPUT] section, with n·T tied to the Solovev P₀≈4.27e4 Pa via the local kinetic.dat.
+# Exercises the real KineticForces (NTV) physics path at full strength, with the main ion
+# declared as a 50/50 D-T mix so the self-consistent kinetic matrices are accumulated over
+# the multi-ion species loop (pure plasma, n_i = n_e, so no impurity species arises). The
+# equilibrium is generated analytically from the embedded [SOL_INPUT] section, with n·T tied
+# to the Solovev P₀≈4.27e4 Pa via the local kinetic.dat.
 
 [Equilibrium]
 eq_type = "sol"                         # Type of the input 2D equilibrium file
@@ -48,6 +50,19 @@ ucrit = 1e3                   # Maximum fraction of solutions allowed before re-
 write_outputs_to_HDF5 = true  # Write outputs to the HDF5 file
 verbose = false               # Enable verbose logging
 
+[KineticForces]
+kinetic_file = "kinetic.dat"  # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns
+
+[[KineticForces.ion_species]]
+z = 1                         # Species charge (fundamental units)
+m = 2                         # Species mass (proton masses) — Deuterium
+fraction = 0.5                # This species' share of the total main-ion density n_i
+
+[[KineticForces.ion_species]]
+z = 1                         # Species charge (fundamental units)
+m = 3                         # Species mass (proton masses) — Tritium
+fraction = 0.5                # This species' share of the total main-ion density n_i
+
 # Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium.
 [SOL_INPUT]
 mr = 128       # Number of radial grid zones
diff --git a/examples/Solovev_kinetic_multiion_example/forcing.dat b/examples/Solovev_kinetic_multiion_example/forcing.dat
new file mode 100644
index 000000000..42eb94703
--- /dev/null
+++ b/examples/Solovev_kinetic_multiion_example/forcing.dat
@@ -0,0 +1,5 @@
+# Forcing data for perturbed equilibrium calculations
+# normalization: normal_field_T
+# Format: n  m  amplitude_real  amplitude_imag
+# Single mode test case: n=1, m=2, amplitude=1e-4 T
+1  2  1e-4  0.0
diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml
new file mode 100644
index 000000000..8f140127a
--- /dev/null
+++ b/examples/Solovev_kinetic_multiion_example/gpec.toml
@@ -0,0 +1,109 @@
+# Solovev analytical equilibrium — n=1 ideal stability + perturbed equilibrium + NTV torque.
+# Extends the Solovev ideal example with a [KineticForces] section, so after the plasma
+# response is computed the neoclassical toroidal viscosity torque is integrated over ψ as a
+# post-processing diagnostic (adaptive quadrature paneled at the rational surfaces).
+# The kinetic profiles (kinetic.dat) are constructed to be consistent with this equilibrium's
+# pressure; the equilibrium is generated analytically from the embedded [SOL_INPUT] section.
+
+[Equilibrium]
+eq_type = "sol"                         # Type of the input 2D equilibrium file
+jac_type = "pest"                       # Coordinate system (hamada, pest, boozer, equal_arc, park, other)
+grid_type = "ldp"                       # Radial grid packing type
+psilow = 1e-4                           # Lower limit of normalized poloidal flux
+psihigh = 0.9995                        # Upper limit of normalized poloidal flux
+mpsi = 128                              # Number of radial grid points (0 = auto-compute from psi_accuracy)
+mtheta = 256                            # Number of poloidal grid points
+newq0 = 0                               # Override for on-axis safety factor (0 = use input value)
+etol = 1e-7                             # Error tolerance for equilibrium solver
+force_termination = false               # Terminate after equilibrium setup (skip stability calculations)
+
+[Wall]
+# Close conformal wall is required to stabilize this Solovev fixture's n=1 external kink:
+# with nowall, et[1] = -6.8 (strongly unstable); with this wall, et[1] = +0.24 (barely stable).
+# The plasma is near marginal stability, so the BVP Δ' matrix values are pathological
+# (dpm magnitudes ~ 10¹¹, |Im/Re| ≫ 1). This fixture's role is integration-pipeline
+# smoke testing + et[1] regression, NOT BVP Δ' regression — DIIID-like is the canonical
+# Δ'-matrix fixture (stable et[1] = +1.6, clean BVP Δ').
+shape = "conformal"                     # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath)
+a = 0.2415                              # Distance from plasma (conformal) or shape parameter
+aw = 0.05                               # Half-thickness parameter for Dee-shaped walls
+bw = 1.5                                # Elongation parameter for wall shapes
+cw = 0                                  # Offset of wall center from major radius
+dw = 0.5                                # Triangularity parameter for wall shapes
+tw = 0.05                               # Sharpness of wall corners (try 0.05 as initial value)
+equal_arc_wall = true                   # Equal arc length distribution of nodes on wall
+
+[ForcingTerms]
+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)
+
+[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
+write_outputs_to_HDF5 = true            # Write perturbed equilibrium outputs to HDF5
+reg_spot = 0.05                         # Regularization width for singular surfaces (0 = disabled)
+
+[ForceFreeStates]
+local_stability_flag = true   # Perform local stability analysis (Mercier and ballooning) across the ψ profile
+mat_flag = true               # Construct coefficient matrices for diagnostic purposes
+ode_flag = true               # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC)
+vac_flag = true               # Compute plasma, vacuum, and total energies for free-boundary modes
+
+psiedge = 0.99                # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable
+qlow = 1.02                   # Integration initiated at q determined by min(q0, qlow)
+qhigh = 1e3                   # Integration terminated at q limit determined by min(qa, qhigh)
+sing_start = 0                # Start integration at the sing_start'th rational from the axis (psilow)
+
+nn_low = 1                    # Smallest toroidal mode number to include
+nn_high = 1                   # Largest toroidal mode number to include
+delta_mlow = 8                # Expands lower bound of Fourier harmonics
+delta_mhigh = 8               # Expands upper bound of Fourier harmonics
+mthvac = 960                  # Number of points used in splines over poloidal angle at the plasma-vacuum interface
+
+kinetic_source = "fixed"      # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model
+kinetic_factor = 0.0          # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode)
+eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations
+singfac_min = 1e-4            # Fractional distance from rational q at which ideal jump enforced
+ucrit = 1e3                   # Column-norm threshold that triggers solution renormalization
+force_wv_symmetry = true      # Enforce symmetry of the vacuum energy matrix
+save_interval = 3             # Save every Nth ODE step (1=all). Always saves near rational surfaces.
+
+# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details)
+use_parallel          = true   # Run parallel FM-propagator BVP path (unlocks singular/delta_prime_matrix)
+parallel_threads      = 2      # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates)
+populate_dense_xi     = true   # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section
+set_psilim_via_dmlim  = false  # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge
+dmlim                 = 0.2    # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true)
+# Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium.
+[SOL_INPUT]
+mr = 128       # Number of radial grid zones
+mz = 128       # Number of axial grid zones
+ma = 128       # Number of flux grid zones
+e = 1.6        # Elongation
+a = 0.33       # Minor radius
+r0 = 1.0       # Major radius
+q0 = 1.9       # Safety factor at the magnetic axis (O-point)
+p0fac = 1      # Scale on-axis pressure (β changes; Φ, q fixed)
+b0fac = 1      # Scale toroidal field at constant β (Bt changes; shape, β fixed)
+f0fac = 1      # Scale toroidal field at constant pressure (β, q change; Φ, p, Bp fixed)
+
+[KineticForces]
+# Multi-ion NTV: the analytic Solovev main ion is modelled as a 50/50 D-T mix (n_i = total),
+# so the summed NTV runs over D + T + the neutrality-closing impurity + electrons.
+kinetic_file = "kinetic.dat"  # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns
+electron = true               # Add electron NTV in addition to the ion species
+zimp = 6                      # Impurity charge (fundamental units); closes quasineutrality
+mimp = 12                     # Impurity mass (proton masses)
+
+[[KineticForces.ion_species]]
+z = 1                         # Species charge (fundamental units)
+m = 2                         # Species mass (proton masses) — Deuterium
+fraction = 0.5                # This species' share of the total main-ion density n_i
+
+[[KineticForces.ion_species]]
+z = 1                         # Species charge (fundamental units)
+m = 3                         # Species mass (proton masses) — Tritium
+fraction = 0.5                # This species' share of the total main-ion density n_i
diff --git a/examples/Solovev_kinetic_multiion_example/kinetic.dat b/examples/Solovev_kinetic_multiion_example/kinetic.dat
new file mode 100644
index 000000000..96fd9fa3e
--- /dev/null
+++ b/examples/Solovev_kinetic_multiion_example/kinetic.dat
@@ -0,0 +1,23 @@
+# Kinetic profiles consistent with this Solovev equilibrium's pressure.
+# Constraint: P(psi_n) = n_e*T_e + n_i*T_i must equal the equilibrium pressure
+#   P(psi_n) = P0*(1 - psi_n),  P0 = 4.273e4 Pa  (linear, zero at the edge).
+# Construction (quasineutral, n_e=n_i, T_e=T_i):
+#   T(psi_n) = 100 + 2600*(1 - psi_n) eV      (core 2700 eV, edge 100 eV)
+#   n(psi_n) = P0*(1 - psi_n) / (2*e*T)        (core ~4.9e19; product n0*T0 = P0/2e fixed)
+# Density carries the edge falloff to a small floor at psi_n=1 so T stays finite
+# (avoids T->0 in the collisionality / drift-frequency formulas); the residual edge
+# pressure from the floor is ~1e-6 of P0. omega_E is a free (non-pressure) profile.
+psi_n   n_i_m3       n_e_m3       T_i_eV    T_e_eV    omega_E_rads
+0.00    4.939e19     4.939e19     2700.0    2700.0    1.0e4
+0.10    4.919e19     4.919e19     2440.0    2440.0    9.0e3
+0.20    4.894e19     4.894e19     2180.0    2180.0    8.0e3
+0.30    4.862e19     4.862e19     1920.0    1920.0    7.0e3
+0.40    4.820e19     4.820e19     1660.0    1660.0    6.0e3
+0.50    4.763e19     4.763e19     1400.0    1400.0    5.0e3
+0.60    4.679e19     4.679e19     1140.0    1140.0    4.0e3
+0.70    4.546e19     4.546e19     880.0     880.0     3.0e3
+0.80    4.302e19     4.302e19     620.0     620.0     2.0e3
+0.90    3.705e19     3.705e19     360.0     360.0     1.0e3
+0.95    2.899e19     2.899e19     230.0     230.0     5.0e2
+0.99    1.058e19     1.058e19     126.0     126.0     1.0e2
+1.00    2.000e18     2.000e18     100.0     100.0     0.0
diff --git a/examples/Solovev_kinetic_multiion_example/run_example.jl b/examples/Solovev_kinetic_multiion_example/run_example.jl
new file mode 100644
index 000000000..fa0a9e53d
--- /dev/null
+++ b/examples/Solovev_kinetic_multiion_example/run_example.jl
@@ -0,0 +1,4 @@
+using Pkg;
+Pkg.activate(joinpath(@__DIR__, "../.."))
+using GeneralizedPerturbedEquilibrium
+GeneralizedPerturbedEquilibrium.main([dirname(@__FILE__)])
diff --git a/regression-harness/cases/solovev_kinetic_calculated.toml b/regression-harness/cases/solovev_kinetic_calculated.toml
index e5ebb29ef..022312cc1 100644
--- a/regression-harness/cases/solovev_kinetic_calculated.toml
+++ b/regression-harness/cases/solovev_kinetic_calculated.toml
@@ -1,10 +1,12 @@
 # Regression case: Solovev n=1 with the "calculated" kinetic path at full strength.
-# Guards the real KineticForces (PENTRC-style) physics path against drift. Each
-# [quantities.*] block names an HDF5 path in the run output, how to extract it, and the
-# noise floor below which a difference is treated as zero.
+# Guards the real KineticForces (PENTRC-style) physics path against drift, with the main
+# ion declared as a 50/50 D-T mix so the self-consistent kinetic-matrix species loop
+# (multi-ion accumulation) is pinned end-to-end. Each [quantities.*] block names an HDF5
+# path in the run output, how to extract it, and the noise floor below which a difference
+# is treated as zero.
 [case]
 name = "solovev_kinetic_calculated"
-description = "Solovev, n=1, kinetic_source=calculated, kinetic_factor=1.0, mode-converged band (delta_m=8, mpert=32). et[1]≈1.894-1.525i with the pressure-consistent kinetic.dat (n·T tied to the Solovev P₀≈4.27e4 Pa)."
+description = "Solovev, n=1, kinetic_source=calculated, kinetic_factor=1.0, 50/50 D-T multi-ion kinetic matrices, mode-converged band (delta_m=8, mpert=32); pressure-consistent kinetic.dat (n·T tied to the Solovev P₀≈4.27e4 Pa)."
 example_dir = "examples/Solovev_kinetic_calculated_example"
 
 # Energies — leading eigenvalues
diff --git a/regression-harness/cases/solovev_kinetic_multiion.toml b/regression-harness/cases/solovev_kinetic_multiion.toml
new file mode 100644
index 000000000..df97e6740
--- /dev/null
+++ b/regression-harness/cases/solovev_kinetic_multiion.toml
@@ -0,0 +1,67 @@
+# Regression case: Solovev n=1 + PE + MULTI-ION NTV torque. The Solovev main ion is modelled as a
+# 50/50 D-T mix, so the KineticForces quadrature sums over D + T + electrons (this fixture is a pure
+# plasma, n_e = n_i, so no impurity species is present). Tracks the summed total AND each per-species
+# contribution, pinning the multi-species loop (resolve_ntv_species + combine_species_states) and the
+# full-composition Zeff. The impurity-species path is covered by the unit tests + the TC-24 case.
+[case]
+name = "solovev_kinetic_multiion"
+description = "Solovev n=1, ideal + PE + multi-ion NTV (D + T + electron summed)"
+example_dir = "examples/Solovev_kinetic_multiion_example"
+
+# Summed total NTV over all species — the headline multi-ion quantity.
+[quantities.ntv_total]
+h5path = "kinetic_forces/fgar/total_torque"
+type = "real_scalar"
+extract = "value"
+label = "NTV total (D+T+imp+e) [N·m]"
+noise_threshold = 1e-8
+order = 10
+
+[quantities.ntv_D]
+h5path = "kinetic_forces_ion1_z1_m2/fgar/total_torque"
+type = "real_scalar"
+extract = "value"
+label = "NTV Deuterium [N·m]"
+noise_threshold = 1e-8
+order = 11
+
+[quantities.ntv_T]
+h5path = "kinetic_forces_ion2_z1_m3/fgar/total_torque"
+type = "real_scalar"
+extract = "value"
+label = "NTV Tritium [N·m]"
+noise_threshold = 1e-8
+order = 12
+
+[quantities.ntv_electron]
+h5path = "kinetic_forces_electron/fgar/total_torque"
+type = "real_scalar"
+extract = "value"
+label = "NTV electron [N·m]"
+noise_threshold = 1e-8
+order = 14
+
+# Stability anchors — confirm the upstream FFS/PE stages feeding the NTV diagnostic.
+[quantities.et_real]
+h5path = "FreeBoundaryStability/eigenmode_energies"
+type = "complex_vector"
+extract = "real_first"
+label = "total energy Re(et[1])"
+noise_threshold = 1e-10
+order = 20
+
+[quantities.msing]
+h5path = "singular/msing"
+type = "int_scalar"
+extract = "value"
+label = "# singular surfaces"
+noise_threshold = 0
+order = 21
+
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0
+order = 90
diff --git a/regression-harness/cases/solovev_kinetic_nuzero.toml b/regression-harness/cases/solovev_kinetic_nuzero.toml
index 078b9c363..65a8f4544 100644
--- a/regression-harness/cases/solovev_kinetic_nuzero.toml
+++ b/regression-harness/cases/solovev_kinetic_nuzero.toml
@@ -1,12 +1,12 @@
 # Regression case: Solovev n=1 collisionless kinetic run (nutype = "zero").
-# Reuses the Solovev_kinetic_calculated_example deck via [overrides] at a
-# mode-converged band; locks the collisionless real-x-space energy-integral
-# eigenvalue plus equilibrium and singular-surface anchors.
+# Reuses the Solovev_kinetic_calculated_example deck (50/50 D-T multi-ion) via
+# [overrides] at a mode-converged band; locks the collisionless real-x-space
+# energy-integral eigenvalue plus equilibrium and singular-surface anchors.
 # Each [quantities.*] block names an HDF5 path in the run output, how to extract it,
 # and the noise floor below which a difference is treated as zero.
 [case]
 name = "solovev_kinetic_nuzero"
-description = "Solovev, n=1, kinetic_source=calculated, nutype=zero (collisionless energy integral, issue #281), mode-converged band (delta_m=8, mpert=32). Reuses the calculated example deck via override; et[1]≈2.333-2.174i with the pressure-consistent kinetic.dat."
+description = "Solovev, n=1, kinetic_source=calculated, nutype=zero (collisionless energy integral), 50/50 D-T multi-ion kinetic matrices, mode-converged band (delta_m=8, mpert=32). Reuses the calculated example deck via override; pressure-consistent kinetic.dat."
 example_dir = "examples/Solovev_kinetic_calculated_example"
 
 # Run the shared calculated deck collisionless (no separate nuzero example directory).
diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl
index 3c7170767..196a9498e 100644
--- a/src/Equilibrium/KineticProfiles.jl
+++ b/src/Equilibrium/KineticProfiles.jl
@@ -12,6 +12,44 @@ shimming.
 using DelimitedFiles
 using HDF5
 
+# Physical constants and collisionality normalizations shared by `load_kinetic_profiles` and
+# `resolve_ntv_species` (Krook collision frequency, Logan & Park 2013 Eq. 6).
+const _EV_J = 1.602e-19          # eV → J
+const _E_CHARGE = 1.602e-19      # elementary charge [C] (PENTRC-parity value, not CODATA)
+const _NKIN = 100                # kinetic-profile resample: uniform ψ_n grid of _NKIN+1 points
+const _KEV_J = 1.602e-16         # 1 keV in J (temperature normalization in the collision frequency)
+const _NU_PREFAC = 3.5e17        # Krook collision-frequency prefactor
+const _MP = 1.672_614e-27        # proton mass [kg]
+const _ME = 9.109_1e-31          # electron mass [kg]
+
+# NRL-formulary Coulomb logarithm (natural log); n_e in m⁻³, T_e in J.
+_coulomb_log(ne, Te) = 17.3 - 0.5 * log(ne / 1.0e20) + 1.5 * log(Te / _KEV_J)
+
+# Hard stop on an unphysical negative density (bad input or a cubic-resample overshoot).
+function _assert_nonneg_density(kinetic_file, label, arr, psi_reg)
+    i = findfirst(<(0), arr)
+    i === nothing ||
+        error("kinetic file '$kinetic_file': negative $label = $(arr[i]) at ψ_n = $(psi_reg[i]) — unphysical")
+end
+
+"""
+    ResolvedNTVSpecies{P}
+
+One fully-resolved species in the NTV sum — a main ion, the neutrality-closing impurity, or the
+electrons. Fields: charge `z`, mass `m` (proton masses), `electron` flag, `label`, and `profiles`
+(a `KineticProfileSplines` view carrying this species' density as `ni_spline` and its
+full-composition collision frequency as `nui_spline`). Produced by `resolve_ntv_species`; both NTV
+paths — the KineticForces ψ-quadrature and the self-consistent kinetic-matrix build — iterate this
+and sum τ = Σ_s τ_s.
+"""
+struct ResolvedNTVSpecies{P}
+    z::Int
+    m::Int
+    electron::Bool
+    label::String
+    profiles::P
+end
+
 """
     KineticProfileData
 
@@ -21,15 +59,17 @@ present; every profile field is `nothing` when the source omits it, so each
 consumer validates only the fields it needs. Units: densities m⁻³,
 temperatures eV, frequencies rad/s, diffusivities m²/s.
 
-| field       | meaning                              |
-|-------------|--------------------------------------|
-| `n_i`       | main-ion density                     |
-| `n_e`       | electron density                     |
-| `T_i`/`T_e` | ion / electron temperature           |
-| `omega_E`   | ExB rotation                         |
-| `omega_tor` | toroidal rotation (optional)         |
-| `chi_e`     | perpendicular heat diffusivity χ⊥    |
-| `chi_phi`   | toroidal momentum diffusivity χ_φ    |
+| field               | meaning                                                                                                           |
+|:------------------- |:----------------------------------------------------------------------------------------------------------------- |
+| `n_i`               | main-ion density                                                                                                  |
+| `n_e`               | electron density                                                                                                  |
+| `T_i`/`T_e`         | ion / electron temperature                                                                                        |
+| `omega_E`           | ExB rotation                                                                                                      |
+| `omega_tor`         | toroidal rotation (optional)                                                                                      |
+| `chi_e`             | perpendicular heat diffusivity χ⊥                                                                                 |
+| `chi_phi`           | toroidal momentum diffusivity χ_φ                                                                                 |
+| `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras |
+| `provenance`        | short string recording the source file/format                                                                     |
 """
 struct KineticProfileData
     psi::Vector{Float64}
@@ -41,14 +81,16 @@ struct KineticProfileData
     omega_tor::Union{Nothing,Vector{Float64}}
     chi_e::Union{Nothing,Vector{Float64}}
     chi_phi::Union{Nothing,Vector{Float64}}
+    species_densities::Union{Nothing,Dict{String,Vector{Float64}}}
     provenance::String
 end
 
 function KineticProfileData(; psi, n_i=nothing, n_e=nothing, T_i=nothing, T_e=nothing,
-    omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing, provenance="")
+    omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing,
+    species_densities=nothing, provenance="")
     _v(x) = x === nothing ? nothing : Float64.(collect(x))
     return KineticProfileData(Float64.(collect(psi)), _v(n_i), _v(n_e), _v(T_i), _v(T_e),
-        _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), String(provenance))
+        _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), species_densities, String(provenance))
 end
 
 const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5")
@@ -72,15 +114,19 @@ function read_kinetic_file(path::AbstractString; group::AbstractString="/")
     return ext in _KINETIC_H5_EXTS ? _read_kinetic_h5(path; group=group) : _read_kinetic_ascii(path)
 end
 
-"""Read a legacy 6-column ASCII kinetic table into a `KineticProfileData`."""
+"""
+Read a legacy 6-column ASCII kinetic table into a `KineticProfileData`.
+"""
 function _read_kinetic_ascii(path::AbstractString)
     psi, ni, ne, Ti, Te, omegaE = _read_kinetic_table(path)
     return KineticProfileData(; psi=psi, n_i=ni, n_e=ne, T_i=Ti, T_e=Te, omega_E=omegaE,
         provenance="ASCII 6-column: $(basename(path))")
 end
 
-"""Read the GPEC HDF5 kinetic schema. Only `psi` is required; other datasets are
-optional and surfaced as-is."""
+"""
+Read the GPEC HDF5 kinetic schema. Only `psi` is required; other standard datasets are
+optional and surfaced as-is. Non-standard `n_*` datasets become named per-species densities.
+"""
 function _read_kinetic_h5(path::AbstractString; group::AbstractString="/")
     h5open(path, "r") do f
         g = group == "/" ? f : f[group]
@@ -89,10 +135,20 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/")
         prov = ""
         ats = attributes(g)
         haskey(ats, "provenance") && (prov = string(read(ats["provenance"])))
+        # Non-standard datasets named `n_*` (e.g. "n_D", "n_T") are named per-species density
+        # profiles for explicit multi-ion input; other names are reserved for future schema fields.
+        standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_phi")
+        extras = Dict{String,Vector{Float64}}()
+        for k in keys(g)
+            (k in standard || !startswith(k, "n_")) && continue
+            v = read(g[k])
+            v isa AbstractArray && (extras[k] = Float64.(vec(v)))
+        end
         return KineticProfileData(; psi=Float64.(vec(read(g["psi"]))),
             n_i=rd("n_i"), n_e=rd("n_e"), T_i=rd("T_i"), T_e=rd("T_e"),
             omega_E=rd("omega_E"), omega_tor=rd("omega_tor"),
-            chi_e=rd("chi_e"), chi_phi=rd("chi_phi"), provenance=prov)
+            chi_e=rd("chi_e"), chi_phi=rd("chi_phi"),
+            species_densities=(isempty(extras) ? nothing : extras), provenance=prov)
     end
 end
 
@@ -101,7 +157,8 @@ end
 
 Write a `KineticProfileData` to the GPEC HDF5 kinetic schema. Each present field
 is written as a dataset with a `units` attribute; absent (`nothing`) fields are
-skipped. The group carries `schema_version` and `provenance` attributes.
+skipped. Named per-species densities (`species_densities`, `n_*` datasets) round-trip
+with `read_kinetic_file`. The group carries `schema_version` and `provenance` attributes.
 """
 function write_kinetic_h5(path::AbstractString, data::KineticProfileData;
     group::AbstractString="/", schema_version::AbstractString="1.0",
@@ -122,12 +179,190 @@ function write_kinetic_h5(path::AbstractString, data::KineticProfileData;
         put("omega_tor", data.omega_tor)
         put("chi_e", data.chi_e)
         put("chi_phi", data.chi_phi)
+        # Named per-species densities (multi-ion input); dataset names must be `n_*` to round-trip.
+        if data.species_densities !== nothing
+            for (name, v) in data.species_densities
+                startswith(name, "n_") ||
+                    error("species density dataset `$name` must be named `n_*` to be readable as a per-species profile")
+                g[name] = collect(Float64, v)
+                attributes(g[name])["units"] = "m^-3"
+            end
+        end
         attributes(g)["schema_version"] = schema_version
         attributes(g)["provenance"] = provenance
     end
     return path
 end
 
+"""
+    multi_ion_composition(zs, ns, ne; zimp, mimp) -> (zeff, zpitch, n_main, n_imp)
+
+Full-composition effective charge and pitch-angle enhancement for a plasma with an
+arbitrary number of main-ion species. Inputs are per-point scalars (`zs`, `ns` are per-species
+values at one ψ; the caller loops grid points):
+
+  - `zs`, `ns`: main-ion charges and number densities (one entry per species)
+  - `ne`: electron density
+  - `zimp`, `mimp`: single trailing impurity charge / mass; its density closes quasineutrality
+
+Composition (per point):
+
+    n_main = Σ_s n_s                       total main-ion number density
+    n_imp  = (n_e − Σ_s z_s n_s) / z_imp   impurity density from quasineutrality
+    Zeff   = (Σ_s z_s² n_s + z_imp² n_imp) / n_e
+    zpitch = 1 + (1+m_imp)/(2 m_imp)·z_imp·(Zeff−1)/(z_imp−Zeff)   momentum-restoring closure
+
+Reduces **exactly** to the single-ion form `Zeff = z_imp − (n_i/n_e)·z_i·(z_imp−z_i)` for one
+z-species (verified algebraically). The per-species pitch-angle collision frequency built from
+this is `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s · (T_s)^{3/2})` — shared `zpitch`, `n_main`,
+`Zeff`, `lnΛ`; the test species contributes only its own `z_s²`, `m_s`, `T_s`.
+"""
+function multi_ion_composition(zs::AbstractVector, ns::AbstractVector, ne::Real;
+    zimp::Real, mimp::Real)
+    n_main = sum(ns)
+    # Vacuum point (ne ≤ 0): neutral composition, no impurity — avoids the zpitch pole at Zeff = zimp.
+    ne > 0 || return (1.0, 1.0, n_main, 0.0)
+    charge_density = sum(z * n for (z, n) in zip(zs, ns))   # Σ z_s n_s
+    n_imp = (ne - charge_density) / zimp
+    zeff = (sum(z^2 * n for (z, n) in zip(zs, ns)) + zimp^2 * n_imp) / ne
+    zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (zeff - 1.0) / (zimp - zeff)
+    return (zeff, zpitch, n_main, n_imp)
+end
+
+"""
+    resolve_ntv_species(kinetic_file, ion_species; electron, zimp, mimp, ...) -> Vector{ResolvedNTVSpecies}
+
+Resolve the **full NTV species set** for a multi-ion run: the main ions, the neutrality-closing
+impurity, and (if `electron`) the electrons. Returns `ResolvedNTVSpecies` descriptors — the
+single set that BOTH NTV paths (the KineticForces quadrature and the self-consistent kinetic-matrix
+build) loop over and sum. Errors on any negative density (unphysical input / resample overshoot).
+
+Each descriptor's `profiles` view carries **that species' resonant density** `n_s` (as `ni_spline`)
+and its **full-composition collision frequency** `ν_s` (as `nui_spline`); `ne/Te/ωE/Zeff/lnΛ` (and
+the electron `nue`) are shared. Main-ion density is resolved from `IonSpecies`: `fraction` ⇒
+`n_s = fraction · n_i` (the kinetic file's `n_i` = total main-ion density); `density` ⇒ an explicit
+`n_*`-named dataset in the HDF5 kinetic file (e.g. `"n_D"`), cubic-resampled to the working grid.
+In a mixed fraction/density list the impurity content is set by the file's n_i/n_e deficit, not by
+a fraction shortfall. Composition (Zeff, zpitch, `n_imp`) comes from `multi_ion_composition`;
+per-species `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s·(T_i)^{3/2})`. The `zpitch·n_main` field
+term is exact for z=1 main ions (D-T, D-H) and approximate otherwise (a z>1 main ion's charge
+enters only through Zeff in the impurity closure). The impurity (`zimp`,`mimp`, density = the
+quasineutrality `n_imp`) is included as its own resonant species. Reduces to
+`load_kinetic_profiles` for one z=1 main ion with `fraction=1`.
+"""
+function resolve_ntv_species(kinetic_file::AbstractString, ion_species::AbstractVector;
+    electron::Bool=false, zimp::Integer=6, mimp::Integer=12,
+    density_factor::Float64=1.0, temperature_factor::Float64=1.0,
+    ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0)
+
+    (density_factor == 1.0 && temperature_factor == 1.0 && ExB_rotation_factor == 1.0 &&
+     toroidal_rotation_factor == 1.0) ||
+        error("profile scaling factors are not yet supported together with a multi-ion ion_species list")
+
+    data = read_kinetic_file(kinetic_file)
+    _need(f, n) = f === nothing ? error("kinetic file '$kinetic_file' missing required '$n'") : f
+    psi_in = data.psi
+    ne_in = _need(data.n_e, "n_e")
+    Ti_in = _need(data.T_i, "T_i")
+    Te_in = _need(data.T_e, "T_e")
+    omE_in = _need(data.omega_E, "omega_E")
+    ni_in = data.n_i === nothing ? copy(ne_in) : data.n_i   # TOTAL main-ion density
+
+    nkin = _NKIN
+    psi_reg = collect(0:nkin) ./ nkin
+    ne = _cubic_resample(psi_in, ne_in, psi_reg)
+    ni_total = _cubic_resample(psi_in, ni_in, psi_reg)
+    Ti = _cubic_resample(psi_in, Ti_in, psi_reg) .* _EV_J
+    Te = _cubic_resample(psi_in, Te_in, psi_reg) .* _EV_J
+    omegaE = _cubic_resample(psi_in, omE_in, psi_reg)
+    _assert_nonneg_density(kinetic_file, "n_e", ne, psi_reg)
+    _assert_nonneg_density(kinetic_file, "n_i (total)", ni_total, psi_reg)
+
+    # Fraction bookkeeping: fractions are shares of the total n_i, so they can never sum above 1,
+    # and when every species is fraction-specified they must account for all of n_i. (In a mixed
+    # fraction/density list no full-sum check is possible; the n_imp ≥ 0 quasineutrality check below
+    # still catches over-allocation in physical units.)
+    fracs = [s.fraction for s in ion_species if !isnan(s.fraction)]
+    if !isempty(fracs)
+        fsum = sum(fracs)
+        fsum <= 1.0 + 1e-6 ||
+            error("ion_species fractions sum to $fsum > 1 — fractions are shares of the total main-ion density n_i")
+        length(fracs) == length(ion_species) && abs(fsum - 1.0) > 1e-6 &&
+            error("ion_species fractions sum to $fsum ≠ 1 — with all species fraction-specified the shares must account for the full n_i " *
+                  "(the impurity content is set by the file's n_i/n_e deficit, not by a fraction shortfall)")
+    end
+
+    # The momentum-restoring zpitch closure assumes the declared main ions are lighter/lower-z than
+    # the trailing impurity; z_s ≥ z_imp can drive Zeff to the zpitch pole at Zeff = z_imp.
+    all(s.z < zimp for s in ion_species) ||
+        error("every ion_species charge must satisfy z < zimp = $zimp (declare heavier species as the impurity)")
+
+    # Resolve each species' resonant density: `fraction` ⇒ share of the total n_i; `density` ⇒
+    # an explicit named profile from the (HDF5) kinetic file (resampled to the working grid).
+    zs = [Int(s.z) for s in ion_species]
+    ns = Vector{Vector{Float64}}(undef, length(ion_species))
+    for (si, s) in enumerate(ion_species)
+        has_frac = !isnan(s.fraction)
+        has_dens = !isempty(s.density)
+        (has_frac ⊻ has_dens) || error("ion_species[$si]: specify exactly one of `fraction` or `density`")
+        if has_dens
+            (data.species_densities !== nothing && haskey(data.species_densities, s.density)) ||
+                error(
+                    "ion_species[$si]: density profile `$(s.density)` not found in kinetic file " *
+                    "(available: $(data.species_densities === nothing ? "none — ASCII files carry only n_i" : join(keys(data.species_densities), ", ")))"
+                )
+            ns[si] = _cubic_resample(psi_in, data.species_densities[s.density], psi_reg)
+        else
+            ns[si] = s.fraction .* ni_total
+        end
+        _assert_nonneg_density(kinetic_file, "ion_species[$si] density", ns[si], psi_reg)
+    end
+
+    # Shared composition + collisionality (natural-log Coulomb log), per grid point.
+    npts = length(psi_reg)
+    zeff = zeros(npts)
+    zpitch = zeros(npts)
+    loglam = zeros(npts)
+    n_main = zeros(npts)
+    n_imp = zeros(npts)
+    nue = zeros(npts)
+    for i in 1:npts
+        z, zp, nm, nimp = multi_ion_composition(zs, [ns[si][i] for si in eachindex(ion_species)], ne[i]; zimp=zimp, mimp=mimp)
+        zeff[i] = z
+        zpitch[i] = zp
+        n_main[i] = nm
+        n_imp[i] = nimp
+        loglam[i] = (ne[i] > 0 && Te[i] > 0) ? _coulomb_log(ne[i], Te[i]) : 0.0
+        nue[i] = Te[i] > 0 ? (zp / _NU_PREFAC) * ne[i] * loglam[i] / (sqrt(_ME / _MP) * (Te[i] / _KEV_J)^1.5) : 0.0
+    end
+    # A negative impurity density means the main ions over-neutralize (Σ z_s n_s > n_e) — unphysical.
+    _assert_nonneg_density(kinetic_file, "impurity density (n_e < Σ z_s n_s)", n_imp, psi_reg)
+    # Approaching the zpitch pole at Zeff = zimp signals a dilute main-ion mix taken up by impurity.
+    maximum(zeff) > 0.9 * zimp &&
+        @warn "max(Zeff) = $(maximum(zeff)) is within 10% of zimp = $zimp — zpitch closure near its pole; check ion_species densities/fractions"
+
+    # Per-species Krook ν of Logan & Park (2013) Eq. (6); numerical form ported from Fortran
+    # PENTRC inputs.f90:238-244 with an added test-particle z_s². Shared zpitch/n_main/Zeff/lnΛ,
+    # species-specific z²/m/T_i; `zpitch` is a main-ion closure, applied approximately to the
+    # impurity/electron test species.
+    _nu(zsp, msp) = [Ti[i] > 0 ? (zpitch[i] / _NU_PREFAC) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * (Ti[i] / _KEV_J)^1.5) : 0.0 for i in 1:npts]
+    _view(dens, nu) = KineticProfileSplines(psi_reg, dens, ne, Ti, Te, omegaE, loglam, nu, nue, zeff)
+
+    # Full NTV species set: main ions, the neutrality-closing impurity, and (optionally) electrons.
+    species = ResolvedNTVSpecies[]
+    for (si, s) in enumerate(ion_species)
+        push!(species, ResolvedNTVSpecies(Int(s.z), Int(s.m), false, "ion$(si)_z$(s.z)_m$(s.m)", _view(ns[si], _nu(s.z, s.m))))
+    end
+    if any(>(0), n_imp)
+        push!(species, ResolvedNTVSpecies(Int(zimp), Int(mimp), false, "impurity_z$(zimp)_m$(mimp)", _view(n_imp, _nu(zimp, mimp))))
+    end
+    if electron
+        # The electron view's `ni_spline` is unused (the electron path reads `ne_spline`); carry ne.
+        push!(species, ResolvedNTVSpecies(-1, 1, true, "electron", _view(ne, nue)))
+    end
+    return identity.(species)   # narrow to a concrete-eltype Vector{ResolvedNTVSpecies{...}}
+end
+
 """
     load_kinetic_profiles(kinetic_file::AbstractString;
                           zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12,
@@ -162,11 +397,12 @@ six-column whitespace-separated ASCII table (header rows are filtered out):
 # Scaling sequence
 
 When any of `density_factor`, `temperature_factor`, `toroidal_rotation_factor` differ from 1.0:
-1. Build first-pass cubic splines from unscaled profiles (for derivatives)
-2. Compute diamagnetic frequencies `wdian`, `wdiat` and total toroidal rotation `wphi`
-3. Scale: `wdian_new = temperature_factor * wdian`, `wdiat_new = temperature_factor * wdiat` (`density_factor` cancels in `T*(dn/dψ)/n`)
-4. Reform: `omegaE = toroidal_rotation_factor * wphi - wdian_new - wdiat_new`
-5. Scale density/temperature arrays: `ni *= density_factor`, `Ti *= temperature_factor`, etc.
+
+ 1. Build first-pass cubic splines from unscaled profiles (for derivatives)
+ 2. Compute diamagnetic frequencies `wdian`, `wdiat` and total toroidal rotation `wphi`
+ 3. Scale: `wdian_new = temperature_factor * wdian`, `wdiat_new = temperature_factor * wdiat` (`density_factor` cancels in `T*(dn/dψ)/n`)
+ 4. Reform: `omegaE = toroidal_rotation_factor * wphi - wdian_new - wdiat_new`
+ 5. Scale density/temperature arrays: `ni *= density_factor`, `Ti *= temperature_factor`, etc.
 
 Then `ExB_rotation_factor` is applied independently: `omegaE *= ExB_rotation_factor`.
 
@@ -192,11 +428,7 @@ function load_kinetic_profiles(kinetic_file::AbstractString;
     omegaE_input = _need(data.omega_E, "omega_E")
     ni_input = data.n_i === nothing ? copy(ne_input) : data.n_i
 
-    eV_to_J = 1.602e-19
-    mp = 1.672_614e-27
-    me = 9.109_1e-31
-
-    nkin = 100
+    nkin = _NKIN
     psi_reg = collect(0:nkin) ./ nkin
 
     # Match Fortran pentrc/inputs.f90:215-232 — cubic-spline-then-resample, NOT
@@ -207,8 +439,8 @@ function load_kinetic_profiles(kinetic_file::AbstractString;
     # denominator). See feedback_kf_kin_profile_linear_interp.md.
     ni = _cubic_resample(psi_input, ni_input, psi_reg)
     ne = _cubic_resample(psi_input, ne_input, psi_reg)
-    Ti = _cubic_resample(psi_input, Ti_input_eV, psi_reg) .* eV_to_J
-    Te = _cubic_resample(psi_input, Te_input_eV, psi_reg) .* eV_to_J
+    Ti = _cubic_resample(psi_input, Ti_input_eV, psi_reg) .* _EV_J
+    Te = _cubic_resample(psi_input, Te_input_eV, psi_reg) .* _EV_J
     omegaE = _cubic_resample(psi_input, omegaE_input, psi_reg)
 
     needs_rotation_reform = density_factor != 1.0 || temperature_factor != 1.0 || toroidal_rotation_factor != 1.0
@@ -228,8 +460,7 @@ function load_kinetic_profiles(kinetic_file::AbstractString;
         dTi_dpsi = deriv1(Ti_spl)
 
         # Compute original wdian, wdiat, wphi and reform omegaE at each grid point
-        echarge = 1.602e-19
-        chrg_ion = zi * echarge
+        chrg_ion = zi * _E_CHARGE
         for i in eachindex(omegaE)
             ψ = psi_reg[i]
             wdian_i = ni[i] > 0 ? -2π * Ti[i] * dni_dpsi(ψ) / (chrg_ion * chi1 * ni[i]) : 0.0
@@ -258,7 +489,7 @@ function load_kinetic_profiles(kinetic_file::AbstractString;
     nue = zeros(Float64, nkin + 1)
     zeff = zeros(Float64, nkin + 1)
 
-    for i in 1:(nkin + 1)
+    for i in 1:(nkin+1)
         n_i = ni[i]
         n_e = ne[i]
         T_i = Ti[i]
@@ -267,14 +498,16 @@ function load_kinetic_profiles(kinetic_file::AbstractString;
         z = n_e > 0 ? zimp - (n_i / n_e) * zi * (zimp - zi) : Float64(zimp)
         zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (z - 1.0) / (zimp - z)
 
-        # Coulomb logarithm (NRL formulary form), natural log; n_e in 1e20 m^-3, T_e in keV.
-        ll = 17.3 - 0.5 * log(n_e / 1.0e20) + 1.5 * log(T_e / 1.602e-16)
+        ll = _coulomb_log(n_e, T_e)
         loglam[i] = ll
 
+        # Krook ν of Logan & Park (2013) Eq. (6); numerical form ported from Fortran PENTRC
+        # inputs.f90:238-244 with an added test-particle z² (generalizes the implicit zi=1 form);
+        # identical to the single-species limit of resolve_ntv_species.
         nui[i] = T_i > 0 ?
-            (zpitch / 3.5e17) * n_i * ll / (sqrt(1.0 * mi) * (T_i / 1.602e-16)^1.5) : 0.0
+                 (zpitch / _NU_PREFAC) * zi^2 * n_i * ll / (sqrt(1.0 * mi) * (T_i / _KEV_J)^1.5) : 0.0
         nue[i] = T_e > 0 ?
-            (zpitch / 3.5e17) * n_e * ll / (sqrt(me / mp) * (T_e / 1.602e-16)^1.5) : 0.0
+                 (zpitch / _NU_PREFAC) * n_e * ll / (sqrt(_ME / _MP) * (T_e / _KEV_J)^1.5) : 0.0
         zeff[i] = z
     end
 
diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl
index 27279f51c..1e46591ac 100755
--- a/src/GeneralizedPerturbedEquilibrium.jl
+++ b/src/GeneralizedPerturbedEquilibrium.jl
@@ -222,6 +222,10 @@ function main_from_inputs(
         KineticForces.KineticForcesControl()
 
     kinetic_profiles = nothing
+    # Multi-species NTV set (main ions + neutrality impurity + electrons), resolved once and
+    # summed by BOTH the self-consistent kinetic-matrix build and the KineticForces quadrature.
+    # `nothing` ⇒ single-species from kf_ctrl.zi/mi (unchanged behaviour).
+    kf_species = nothing
     needs_kinetic_profiles = haskey(inputs, "KineticForces") ||
                              (ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated")
     if needs_kinetic_profiles
@@ -233,6 +237,17 @@ function main_from_inputs(
             density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor,
             ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor,
             chi1=2π * equil.psio)
+        # `electron=true` always means electrons IN ADDITION to the ions (the Fortran PENTRC
+        # either/or was an IO limitation, not physics), so it routes through the species
+        # resolver even with an empty ion_species list (single main ion from zi/mi).
+        if !isempty(kf_ctrl.ion_species) || kf_ctrl.electron
+            speclist = isempty(kf_ctrl.ion_species) ?
+                       [KineticForces.IonSpecies(; z=kf_ctrl.zi, m=kf_ctrl.mi, fraction=1.0)] : kf_ctrl.ion_species
+            kf_species = Equilibrium.resolve_ntv_species(kinetic_file, speclist;
+                electron=kf_ctrl.electron, zimp=kf_ctrl.zimp, mimp=kf_ctrl.mimp,
+                density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor,
+                ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor)
+        end
     end
 
     # Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry,
@@ -425,7 +440,7 @@ function main_from_inputs(
         calculated_cb = (c, e, i, m, f) ->
             KineticForces.compute_calculated_kinetic_matrices(
                 c, e, i, m, f;
-                kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles)
+                kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, species=kf_species)
         make_kinetic_matrix(ctrl, equil, ffit, intr, metric;
             calculated_source=calculated_cb)
 
@@ -646,12 +661,37 @@ function main_from_inputs(
             kf_intr = KineticForces.KineticForcesInternal(equil; verbose=kf_ctrl.verbose)
             KineticForces.set_perturbation_data!(kf_intr, pe_state, intr, equil, metric)
 
-            kf_state = KineticForces.KineticForcesState()
-            KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, kinetic_profiles)
-
-            if kf_ctrl.write_outputs_to_HDF5
-                h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file
-                    KineticForces.write_to_hdf5!(h5file, kf_state)
+            if kf_species === nothing
+                # Single main ion (from zi/mi), no electrons — unchanged path.
+                kf_state = KineticForces.KineticForcesState()
+                KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, kinetic_profiles)
+                if kf_ctrl.write_outputs_to_HDF5
+                    h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file
+                        KineticForces.write_to_hdf5!(h5file, kf_state)
+                    end
+                end
+            else
+                # Multi-species: summed over main ions + neutrality impurity + electrons
+                # (kf_species, resolved once above and shared with the kinetic-matrix build).
+                states = KineticForces.KineticForcesState[]
+                labels = String[]
+                for sp in kf_species
+                    sctrl = deepcopy(kf_ctrl)
+                    sctrl.zi = sp.z; sctrl.mi = sp.m; sctrl.electron = sp.electron
+                    sctrl.ion_species = KineticForces.IonSpecies[]
+                    st = KineticForces.KineticForcesState()
+                    KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, equil, sp.profiles)
+                    push!(states, st); push!(labels, sp.label)
+                end
+                @info "Multi-species NTV: summed $(length(kf_species)) species ($(join((s.label for s in kf_species), ", ")))"
+                kf_state = KineticForces.combine_species_states(states)
+                if kf_ctrl.write_outputs_to_HDF5
+                    h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file
+                        for (st, lbl) in zip(states, labels)
+                            KineticForces.write_to_hdf5!(h5file, st; group_name="kinetic_forces_$lbl")
+                        end
+                        KineticForces.write_to_hdf5!(h5file, kf_state)  # "kinetic_forces" = total
+                    end
                 end
             end
         end
diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl
index d36fead99..ec0c0997a 100644
--- a/src/KineticForces/CalculatedKineticMatrices.jl
+++ b/src/KineticForces/CalculatedKineticMatrices.jl
@@ -31,22 +31,25 @@ tracked as follow-up work blocked on PR #196 — see the plan's "Out of scope"
 section.
 
 # Arguments
-- `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`)
-- `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines
-- `ffs_intr`: ForceFreeStatesInternal (mode indexing)
-- `metric`: MetricData (provides ψ grid via `metric.xs`)
-- `ffit`: FourFitVars (used only for `numpert_total` cross-check)
+
+  - `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`)
+  - `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines
+  - `ffs_intr`: ForceFreeStatesInternal (mode indexing)
+  - `metric`: MetricData (provides ψ grid via `metric.xs`)
+  - `ffit`: FourFitVars (used only for `numpert_total` cross-check)
 
 # Keyword arguments
-- `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to
-  carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, 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`.
+
+  - `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to
+    carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, 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`.
 
 # Returns
-- `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)`
-- `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)`
+
+  - `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)`
+  - `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)`
 """
 function compute_calculated_kinetic_matrices(
     _ffs_ctrl,
@@ -54,8 +57,9 @@ function compute_calculated_kinetic_matrices(
     ffs_intr,
     metric,
     ffit;
-    kf_ctrl::KineticForcesControl = KineticForcesControl(),
+    kf_ctrl::KineticForcesControl=KineticForcesControl(),
     kinetic_profiles::Equilibrium.KineticProfileSplines,
+    species::Union{Nothing,AbstractVector{<:Equilibrium.ResolvedNTVSpecies}}=nothing
 )
     xs = metric.xs
     mpsi = length(xs)
@@ -95,44 +99,58 @@ function compute_calculated_kinetic_matrices(
     # are read-only and safely shared through deepcopy semantics.
     nl = kf_ctrl.nl
     nthreads = Threads.maxthreadid()
-    thread_intrs   = [deepcopy(kf_intr) for _ in 1:nthreads]
-    thread_full_w  = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
-    thread_full_t  = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
+    thread_intrs = [deepcopy(kf_intr) for _ in 1:nthreads]
+    thread_full_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
+    thread_full_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
     thread_block_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
     thread_block_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads]
 
-    Threads.@threads for ipsi in 1:mpsi
-        tid     = Threads.threadid()
-        intr_t  = thread_intrs[tid]
-        full_w  = thread_full_w[tid]
-        full_t  = thread_full_t[tid]
-        block_w = thread_block_w[tid]
-        block_t = thread_block_t[tid]
-        psi     = xs[ipsi]
-        for in_idx in 1:npert
-            n = ffs_intr.nlow + in_idx - 1
-            fill!(full_w, 0)
-            fill!(full_t, 0)
-            for ell in -nl:nl
-                fill!(block_w, 0)
-                fill!(block_t, 0)
-                compute_kinetic_matrices_at_psi!(
-                    block_w, block_t, psi, n, ell,
-                    kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.wdfac, kf_ctrl.divxfac,
-                    kf_ctrl.electron, equil, intr_t, kinetic_profiles;
-                    nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.nufac,
-                    atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda
-                )
-                full_w .+= block_w
-                full_t .+= block_t
-            end
-
-            # Place the n-block on the diagonal of the full np×np matrix and flatten.
-            row_offset = (in_idx - 1) * mpert
-            for k in 1:6, j in 1:mpert, i in 1:mpert
-                idx = (row_offset + j - 1) * np + (row_offset + i)
-                kw_flat[ipsi, idx, k] = full_w[i, j, k]
-                kt_flat[ipsi, idx, k] = full_t[i, j, k]
+    # Species set summed into the kinetic matrices: the resolved multi-species set
+    # (main ions + neutrality impurity + electrons) when provided, else a single species
+    # from kf_ctrl. The kinetic W/torque matrices are additive over species, so we
+    # accumulate (+=) each species' block-diagonal contribution.
+    splist = species === nothing ?
+             [Equilibrium.ResolvedNTVSpecies(kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.electron, "single", kinetic_profiles)] : species
+
+    for sp in splist
+        # Hoisted per-species scalars: the closure then captures concrete values, not a
+        # union-typed `sp` (the ternary above mixes the resolved vector with the fallback).
+        z_s, m_s, el_s, prof_s = sp.z, sp.m, sp.electron, sp.profiles
+        # :static pins one task per thread so the threadid() buffer indexing below is sound
+        # (dynamic scheduling may migrate tasks at yield points, silently sharing buffers).
+        Threads.@threads :static for ipsi in 1:mpsi
+            tid = Threads.threadid()
+            intr_t = thread_intrs[tid]
+            full_w = thread_full_w[tid]
+            full_t = thread_full_t[tid]
+            block_w = thread_block_w[tid]
+            block_t = thread_block_t[tid]
+            psi = xs[ipsi]
+            for in_idx in 1:npert
+                n = ffs_intr.nlow + in_idx - 1
+                fill!(full_w, 0)
+                fill!(full_t, 0)
+                for ell in -nl:nl
+                    fill!(block_w, 0)
+                    fill!(block_t, 0)
+                    compute_kinetic_matrices_at_psi!(
+                        block_w, block_t, psi, n, ell,
+                        z_s, m_s, kf_ctrl.wdfac, kf_ctrl.divxfac,
+                        el_s, equil, intr_t, prof_s;
+                        nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.nufac,
+                        atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda
+                    )
+                    full_w .+= block_w
+                    full_t .+= block_t
+                end
+
+                # Place the n-block on the diagonal of the full np×np matrix and accumulate.
+                row_offset = (in_idx - 1) * mpert
+                for k in 1:6, j in 1:mpert, i in 1:mpert
+                    idx = (row_offset + j - 1) * np + (row_offset + i)
+                    kw_flat[ipsi, idx, k] += full_w[i, j, k]
+                    kt_flat[ipsi, idx, k] += full_t[i, j, k]
+                end
             end
         end
     end
diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl
index 4175a62b4..d0fdd591b 100644
--- a/src/KineticForces/Compute.jl
+++ b/src/KineticForces/Compute.jl
@@ -41,7 +41,7 @@ Warn when the ψ torque quadrature terminated without satisfying the requested t
 silent-garbage scenario for weak applied fields, since NTV scales as δB².
 """
 function check_psi_quadrature_convergence(total::ComplexF64, quad_err::Float64,
-                                          ctrl::KineticForcesControl, method::String)
+    ctrl::KineticForcesControl, method::String)
     if quad_err > max(ctrl.atol_psi, ctrl.rtol_psi * abs(total))
         @warn "ψ torque quadrature ($method) did not converge within maxevals_psi=$(ctrl.maxevals_psi): " *
               "error estimate $quad_err vs |T|=$(abs(total)) N·m. Raise maxevals_psi or loosen rtol_psi."
@@ -63,14 +63,16 @@ diagnostic T(ψ) profile at no extra cost (the values are computed anyway
 — we just keep them).
 
 # Returns
+
 NamedTuple with:
-- `total::ComplexF64`: Total integrated torque
-- `torque_profile`: NamedTuple of (psi, dtdpsi, t_cumulative) from evaluation points
-- `matrix_integrated`: Trapezoidal-integrated mpert×mpert×6 matrix (if matrix method)
-- `psi_nsteps::Int`: Number of integrand evaluations
-- `psi_quad_error::Float64`: Quadrature error estimate for the total torque
-- `panel_psis::Vector{Float64}`: Quadrature panel boundaries actually used (bounds + interior resonant surfaces)
-- `resonance_psis::Vector{Float64}`: Located kinetic-resonance ψ surfaces (Ω_ℓ(x=1)=0), for diagnostics/plotting
+
+  - `total::ComplexF64`: Total integrated torque
+  - `torque_profile`: NamedTuple of (psi, dtdpsi, t_cumulative) from evaluation points
+  - `matrix_integrated`: Trapezoidal-integrated mpert×mpert×6 matrix (if matrix method)
+  - `psi_nsteps::Int`: Number of integrand evaluations
+  - `psi_quad_error::Float64`: Quadrature error estimate for the total torque
+  - `panel_psis::Vector{Float64}`: Quadrature panel boundaries actually used (bounds + interior resonant surfaces)
+  - `resonance_psis::Vector{Float64}`: Located kinetic-resonance ψ surfaces (Ω_ℓ(x=1)=0), for diagnostics/plotting
 
 The integral is paneled at the rational-surface ψ locations (`intr.sing_psis`) and capped
 at `ctrl.maxevals_psi` evaluations; a warning is emitted if the quadrature fails to reach
@@ -95,7 +97,7 @@ function integrate_psi_quadgk(
 
     if x0 >= xout
         return (total=ComplexF64(0.0), torque_profile=nothing, matrix_integrated=nothing, psi_nsteps=0, psi_quad_error=0.0,
-                panel_psis=Float64[], resonance_psis=Float64[])
+            panel_psis=Float64[], resonance_psis=Float64[])
     end
 
     # The outer ψ-integral (the QuadGK batch / ψ-node loop) stays serial: QuadGK's refine
@@ -133,9 +135,9 @@ function integrate_psi_quadgk(
                 w = is_matrix_method ? thread_wtw[tid] : nothing
                 is_matrix_method && fill!(w, 0)
                 tpsi!(thread_tpsi[tid], psi, n, l, zi, mi, wdfac, divxfac,
-                      electron, method, equil, thread_intrs[tid], kinetic_profiles;
-                      op_wmats=w,
-                      atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda)
+                    electron, method, equil, thread_intrs[tid], kinetic_profiles;
+                    op_wmats=w,
+                    atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda)
                 harm_vals[ell_idx] = thread_tpsi[tid][]
                 is_matrix_method && (harm_elems[ell_idx] .= w)
             end
@@ -205,7 +207,7 @@ function integrate_psi_quadgk(
           "$(length(pts) - 1) panels ($n_rational rational + $n_resonance kinetic resonance surfaces)"
 
     return (total=total, torque_profile=torque_profile, matrix_integrated=matrix_integrated, psi_nsteps=npts, psi_quad_error=quad_err,
-            panel_psis=pts, resonance_psis=sort(resonance_psis))
+        panel_psis=pts, resonance_psis=sort(resonance_psis))
 end
 
 
@@ -213,6 +215,63 @@ end
 # High-level orchestration
 # ============================================================================
 
+"""
+    combine_species_states(states) -> KineticForcesState
+
+Sum per-species [`KineticForcesState`](@ref) results into a single total (τ = Σ_s τ_s). Per
+method: the scalar `total_torque`/`total_energy` are summed exactly; the dT/dψ profile is summed
+by linearly interpolating each species' own `(psi_grid, dtdpsi)` arrays onto the sorted union of
+the species ψ grids (zero outside a species' range), and the cumulative T(ψ) is re-integrated
+(trapezoid) from that summed profile. All species share the same ψ-integration range
+(`ctrl.psilims`), so the grids differ only in adaptive nodes. The interpolated/trapezoid
+`t_cumulative` is a diagnostic profile — its endpoint need not equal the exactly-summed
+Gauss-Kronrod `total_torque`, especially near sharp resonances. The combined `MethodResult`
+carries only the summed scalars and profile; per-species diagnostics (`torque_profile`, `records`,
+`panel_psis`, `resonance_psis`) are not aggregated and are left at their defaults.
+"""
+function combine_species_states(states::AbstractVector{KineticForcesState})
+    combined = KineticForcesState()
+    isempty(states) && return combined
+    method_names = sort(unique(Iterators.flatten(keys(s.method_results) for s in states)))
+    for mname in method_names
+        results = [s.method_results[mname] for s in states if haskey(s.method_results, mname)]
+        isempty(results) && continue
+        grid = sort(unique(reduce(vcat, (r.psi_grid for r in results); init=Float64[])))
+        # Sum each species' dT/dψ onto the union grid by linear interpolation of its own
+        # (psi_grid, dtdpsi) arrays (torque_profile interpolants are not populated here); zero
+        # outside a species' ψ range.
+        dtdpsi = zeros(ComplexF64, length(grid))
+        for r in results
+            length(r.psi_grid) >= 2 || continue
+            o = sortperm(r.psi_grid)
+            xs = r.psi_grid[o]
+            ys = r.dtdpsi[o]
+            for (k, ψ) in enumerate(grid)
+                (ψ < xs[1] || ψ > xs[end]) && continue
+                j = searchsortedlast(xs, ψ)
+                if j == length(xs) || xs[j+1] == xs[j]   # endpoint or duplicate node: no interpolation
+                    dtdpsi[k] += ys[j]
+                else
+                    t = (ψ - xs[j]) / (xs[j+1] - xs[j])
+                    dtdpsi[k] += (1 - t) * ys[j] + t * ys[j+1]
+                end
+            end
+        end
+        tcum = zeros(ComplexF64, length(grid))
+        for j in 2:length(grid)
+            tcum[j] = tcum[j-1] + 0.5 * (dtdpsi[j] + dtdpsi[j-1]) * (grid[j] - grid[j-1])
+        end
+        combined.method_results[mname] = MethodResult(;
+            method=mname, nn=first(results).nn,
+            total_torque=sum(r.total_torque for r in results),
+            total_energy=sum(r.total_energy for r in results),
+            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
+
 """
     compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal,
                                 ctrl::KineticForcesControl, equil, kinetic_profiles)
@@ -224,15 +283,16 @@ For multi-n calculations, loops over toroidal mode numbers and assembles
 block-diagonal kinetic matrices.
 
 # Arguments
-- `state::KineticForcesState`: Accumulates results for all methods
-- `intr::KineticForcesInternal`: Internal state with equilibrium data
-- `ctrl::KineticForcesControl`: Control parameters specifying which methods to run
-- `equil`: PlasmaEquilibrium with 2D interpolants
-- `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines
+
+  - `state::KineticForcesState`: Accumulates results for all methods
+  - `intr::KineticForcesInternal`: Internal state with equilibrium data
+  - `ctrl::KineticForcesControl`: Control parameters specifying which methods to run
+  - `equil`: PlasmaEquilibrium with 2D interpolants
+  - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines
 """
 function compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal,
-                                     ctrl::KineticForcesControl, equil,
-                                     kinetic_profiles::Equilibrium.KineticProfileSplines)
+    ctrl::KineticForcesControl, equil,
+    kinetic_profiles::Equilibrium.KineticProfileSplines)
 
     for entry in METHOD_REGISTRY
         getfield(ctrl, entry.flag) || continue
@@ -311,7 +371,7 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor
             t_cumulative=t_cum_out,
             psi_nsteps=psi_nsteps_total,
             panel_psis=panel_psis_out,
-            resonance_psis=resonance_psis_out,
+            resonance_psis=resonance_psis_out
         )
         state.method_results[method] = result_entry
 
diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl
index 7e533d8e9..5cca62539 100644
--- a/src/KineticForces/KineticForcesStructs.jl
+++ b/src/KineticForces/KineticForcesStructs.jl
@@ -3,15 +3,16 @@
 
 Single source of truth for the NTV calculation methods. Each entry is a NamedTuple
 `(name, flag, kind, doc)`:
-- `name`  — short method identifier used as the HDF5 group key and in `intr.method`
-- `flag`  — the `KineticForcesControl` field symbol that enables the method
-- `kind`  — dispatch routing tag consumed by `method_kind` / `Torque.jl`
-            (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the
-            three special-cased methods)
-- `doc`   — one-line description printed in verbose output
+
+  - `name`  — short method identifier used as the HDF5 group key and in `intr.method`
+  - `flag`  — the `KineticForcesControl` field symbol that enables the method
+  - `kind`  — dispatch routing tag consumed by `method_kind` / `Torque.jl`
+    (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the
+    three special-cased methods)
+  - `doc`   — one-line description printed in verbose output
 
 The method names/docs and the `Compute.jl` enable list are all derived from this
-tuple, and `Torque.jl` routes on `kind`, so the methods are enumerated in one place. 
+tuple, and `Torque.jl` routes on `kind`, so the methods are enumerated in one place.
 To add a method: append an entry here and add the matching `*_flag` field
 to `KineticForcesControl`.
 """
@@ -49,13 +50,53 @@ function method_kind(name::AbstractString)
 end
 
 
+"""
+    IonSpecies(; z, m, fraction=NaN, density="")
+
+One main-ion species in a multi-ion NTV run. `z`/`m` are the charge (e) and mass (proton
+masses). The density is given by **exactly one** of:
+
+  - `fraction`: this species' share of the total main-ion density (the kinetic file's `n_i`
+    column), e.g. 0.5 for each of 50/50 D-T. Convenient when species share a profile shape.
+  - `density`: the name of an explicit per-species density profile in the (HDF5) kinetic file,
+    for measured profiles with distinct shapes.
+
+Both resolve to a per-species density profile `n_s(ψ)` at load time (one downstream path).
+Constructed from a TOML `[[KineticForces.ion_species]]` table via `convert`.
+"""
+struct IonSpecies
+    z::Int
+    m::Int
+    fraction::Float64
+    density::String
+end
+IonSpecies(; z::Integer, m::Integer, fraction=NaN, density="") =
+    IonSpecies(Int(z), Int(m), Float64(fraction), String(density))
+IonSpecies(d::AbstractDict) = IonSpecies(; (Symbol(k) => v for (k, v) in d)...)
+# TOML.jl parses [[KineticForces.ion_species]] as a Vector{Any} of Dicts, so accept any
+# AbstractVector and build each element (pass IonSpecies through unchanged).
+Base.convert(::Type{Vector{IonSpecies}}, v::AbstractVector) =
+    IonSpecies[x isa IonSpecies ? x : IonSpecies(x) for x in v]
+
 """
     KineticForcesControl
 
 User-facing control parameters from the TOML `[KineticForces]` section.
 Configures which NTV methods to run, species parameters, tolerances, and output options.
 
+`ion_species` selects the plasma model: empty (default) ⇒ a single main ion from `zi`/`mi`
+(unchanged behaviour); non-empty ⇒ a multi-main-ion run where the NTV is computed per
+`IonSpecies` under one shared full-composition Zeff and summed. For a multi-ion run the kinetic
+file's `n_i` column is the **total** main-ion density, split across species by their `fraction`
+or explicit `density` profile.
+
+`electron = true` always adds the electron NTV **in addition to** the ion species, on both the
+single-ion and multi-ion paths (a single-ion + electron run resolves as the species pair
+{main ion, electron} plus any quasineutrality impurity, summed). The legacy Fortran PENTRC
+electron-*instead-of*-ions behaviour was an IO limitation and is not reproduced.
+
 Constructed via keyword arguments or from a TOML dict:
+
 ```julia
 ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...)
 ```
@@ -89,7 +130,8 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce
     mi::Int = 2                     # Ion mass (proton masses)
     zimp::Int = 6                   # Impurity charge
     mimp::Int = 12                  # Impurity mass
-    electron::Bool = false          # Include electron contribution
+    electron::Bool = false          # Add electron NTV in addition to the ion species (see docstring)
+    ion_species::Vector{IonSpecies} = IonSpecies[]   # multi-main-ion set (see docstring); empty ⇒ single ion
 
     # Mode numbers
     nn::Int = 1                     # Toroidal mode number
@@ -176,12 +218,13 @@ Internal working state for KineticForces calculations.
 Holds equilibrium-derived quantities, profile interpolants, and integration results.
 
 Fields replacing former module-level globals:
-- `ro`, `bo`, `chi1`: Equilibrium geometry parameters
-- `mthsurf`, `mfac`: Poloidal grid info
-- `dbob_m`, `divx_m`: Perturbation mode interpolants
-- `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as
-  panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on
-  Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection
+
+  - `ro`, `bo`, `chi1`: Equilibrium geometry parameters
+  - `mthsurf`, `mfac`: Poloidal grid info
+  - `dbob_m`, `divx_m`: Perturbation mode interpolants
+  - `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as
+    panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on
+    Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection
 
 Equilibrium and kinetic profile data are read directly from the
 `PlasmaEquilibrium` (`equil.profiles`, `equil.geometry`) and the
@@ -271,17 +314,17 @@ function KineticForcesInternal(equil; verbose::Bool=false)
     # Axis toroidal field F(0)/ro that normalizes λ = μ·bo/E; F_spline stores 2πF.
     bo_axis = abs(equil.profiles.F_spline(0.0)) / (2π * equil.ro)
     KineticForcesInternal(;
-        ro      = equil.ro,
-        bo      = bo_axis,
-        chi1    = 2π * equil.psio,
+        ro=equil.ro,
+        bo=bo_axis,
+        chi1=2π * equil.psio,
         mthsurf,
-        tpsi_xs       = collect(range(0.0, 1.0, length=nth)),
-        tpsi_B        = Vector{Float64}(undef, nth),
-        tpsi_dBdpsi   = Vector{Float64}(undef, nth),
-        tpsi_dBdtheta = Vector{Float64}(undef, nth),
-        tpsi_jac      = Vector{Float64}(undef, nth),
-        tpsi_djdpsi   = Vector{Float64}(undef, nth),
-        verbose,
+        tpsi_xs=collect(range(0.0, 1.0; length=nth)),
+        tpsi_B=Vector{Float64}(undef, nth),
+        tpsi_dBdpsi=Vector{Float64}(undef, nth),
+        tpsi_dBdtheta=Vector{Float64}(undef, nth),
+        tpsi_jac=Vector{Float64}(undef, nth),
+        tpsi_djdpsi=Vector{Float64}(undef, nth),
+        verbose
     )
 end
 
@@ -291,19 +334,21 @@ end
 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
+
+ 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
 
 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
+
+ 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
 """
 function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_intr,
-                                equil::Equilibrium.PlasmaEquilibrium,
-                                metric::ForceFreeStates.MetricData)
+    equil::Equilibrium.PlasmaEquilibrium,
+    metric::ForceFreeStates.MetricData)
     # Copy mode numbers from FFS
     kf_intr.mlow = ffs_intr.mlow
     kf_intr.mhigh = ffs_intr.mhigh
@@ -372,9 +417,9 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in
         psi = psi_grid[ipsi]
 
         # Get Clebsch displacement vectors at this ψ
-        xsp  = view(xi_modes.clebsch_psi,  ipsi, :)       # ξ^ψ [mpert]
+        xsp = view(xi_modes.clebsch_psi, ipsi, :)       # ξ^ψ [mpert]
         xmp1 = view(xi_modes.clebsch_psi1, ipsi, :)       # ∂ξ^ψ/∂ψ [mpert]
-        xms  = view(clebsch_alpha_mat, 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)
@@ -394,8 +439,8 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in
         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
+        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
@@ -422,9 +467,9 @@ 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})
+    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)
 
@@ -494,8 +539,8 @@ Accumulated results from all KineticForces computations.
 Written to gpec.h5 under the "kinetic_forces" group.
 """
 @kwdef mutable struct KineticForcesState
-    method_results::Dict{String, MethodResult} = Dict{String, MethodResult}()
+    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}}()
+    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 d1fef4e09..6fadf73a4 100644
--- a/src/KineticForces/Output.jl
+++ b/src/KineticForces/Output.jl
@@ -15,8 +15,8 @@ Write all KineticForces results to the "kinetic_forces" group in gpec.h5.
 - `h5file::HDF5.File`: Open HDF5 file handle
 - `state::KineticForcesState`: Accumulated computation results
 """
-function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState)
-    g = create_group(h5file, "kinetic_forces")
+function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; group_name::AbstractString="kinetic_forces")
+    g = create_group(h5file, group_name)
 
     for (method_name, result) in state.method_results
         mg = create_group(g, method_name)
diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl
index f9af86485..1cd40285e 100644
--- a/src/KineticForces/Torque.jl
+++ b/src/KineticForces/Torque.jl
@@ -153,6 +153,8 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int,
         dT_s_dpsi = kinetic_profiles.Te_deriv(psi)
         nu_s      = kinetic_profiles.nue_spline(psi)
     else
+        # ni_spline / nui_spline carry THIS species' resonant density and its
+        # full-composition collisionality (per-species views built by the loader).
         n_s       = kinetic_profiles.ni_spline(psi)
         T_s       = kinetic_profiles.Ti_spline(psi)
         dn_s_dpsi = kinetic_profiles.ni_deriv(psi)
diff --git a/test/runtests.jl b/test/runtests.jl
index 2036e610c..790c5a0d9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -45,6 +45,7 @@ else
     include("./runtests_dispersion_polish.jl")
     include("./runtests_slayer_runner.jl")
     include("./runtests_kinetic.jl")
+    include("./runtests_multiion.jl")
     include("./runtests_fullruns.jl")
     include("./runtests_coils.jl")
     include("./runtests_imas.jl")
diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl
new file mode 100644
index 000000000..a378ae02f
--- /dev/null
+++ b/test/runtests_multiion.jl
@@ -0,0 +1,201 @@
+using HDF5
+
+@testset "Multi-ion NTV" begin
+    EQ = GeneralizedPerturbedEquilibrium.Equilibrium
+    KF = GeneralizedPerturbedEquilibrium.KineticForces
+
+    @testset "multi_ion_composition" begin
+        ne = 1.24e20
+        ni = 1.13e20
+        zimp = 6
+        mimp = 12
+        # Reduces exactly to the single-ion formula Zeff = zimp - (ni/ne)*zi*(zimp-zi).
+        for zi in (1, 2)
+            zeff_old = zimp - (ni / ne) * zi * (zimp - zi)
+            zeff, zpitch, nmain, nimp = EQ.multi_ion_composition([zi], [ni], ne; zimp=zimp, mimp=mimp)
+            @test zeff ≈ zeff_old rtol = 1e-12
+            @test nmain ≈ ni
+        end
+        # 50/50 D-T with total ni == single ion at total ni; quasineutrality closes the impurity.
+        zdt, _, nmdt, nimpdt = EQ.multi_ion_composition([1, 1], [ni / 2, ni / 2], ne; zimp=zimp, mimp=mimp)
+        zsingle, _, _, _ = EQ.multi_ion_composition([1], [ni], ne; zimp=zimp, mimp=mimp)
+        @test zdt ≈ zsingle rtol = 1e-12
+        @test nmdt ≈ ni
+        @test ne ≈ ni + zimp * nimpdt rtol = 1e-12          # ne = Σ z_s n_s + z_imp n_imp
+    end
+
+    @testset "IonSpecies + TOML conversion" begin
+        # TOML.jl parses [[KineticForces.ion_species]] as a Vector{Any} of Dicts.
+        v = Any[Dict("z" => 1, "m" => 2, "fraction" => 0.5), Dict("z" => 1, "m" => 3, "fraction" => 0.5)]
+        c = KF.KineticForcesControl(; ion_species=v)
+        @test length(c.ion_species) == 2
+        @test (c.ion_species[1].z, c.ion_species[1].m, c.ion_species[1].fraction) == (1, 2, 0.5)
+        @test isempty(KF.KineticForcesControl().ion_species)     # default single-ion
+    end
+
+    @testset "resolve_ntv_species" begin
+        # Synthetic ASCII kinetic profile: psi, n_i(total), n_e, T_i, T_e, omega_E.
+        psi = collect(0.0:0.05:1.0)
+        Ni = @. 1.0e20 * (1 - 0.5psi)
+        ne = @. 1.1e20 * (1 - 0.5psi)
+        Ti = @. 2000.0 * (1 - 0.5psi)
+        Te = @. 2500.0 * (1 - 0.5psi)
+        wE = @. 1.0e4 * (1 - psi)
+        f = tempname() * ".txt"
+        open(f, "w") do io
+            println(io, "# psi ni ne Ti Te wexb")
+            for i in eachindex(psi)
+                println(io, join((psi[i], Ni[i], ne[i], Ti[i], Te[i], wE[i]), "  "))
+            end
+        end
+
+        # Single species [z=1,m=2,fraction=1] reproduces load_kinetic_profiles.
+        single = EQ.load_kinetic_profiles(f; zi=1, zimp=6, mi=2, mimp=12)
+        sp1 = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=1.0)]; electron=false, zimp=6, mimp=12)
+        main = sp1[1]
+        for ψ in (0.3, 0.6)
+            @test main.profiles.ni_spline(ψ) ≈ single.ni_spline(ψ) rtol = 1e-8
+            @test main.profiles.nui_spline(ψ) ≈ single.nui_spline(ψ) rtol = 1e-8
+            @test main.profiles.zeff_spline(ψ) ≈ single.zeff_spline(ψ) rtol = 1e-8
+        end
+
+        # D-T-e set = {D, T, impurity, electron}; ν scales as 1/√m at fixed z, T.
+        sp = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=0.5), KF.IonSpecies(; z=1, m=3, fraction=0.5)];
+            electron=true, zimp=6, mimp=12)
+        @test length(sp) == 4
+        @test count(s -> !s.electron && s.z == 1, sp) == 2   # D, T
+        @test any(s -> s.z == 6 && !s.electron, sp)          # impurity
+        @test any(s -> s.electron, sp)                       # electron
+        iD = findfirst(s -> s.m == 2, sp)
+        iT = findfirst(s -> s.m == 3, sp)
+        @test sp[iD].profiles.nui_spline(0.5) / sp[iT].profiles.nui_spline(0.5) ≈ sqrt(3 / 2) rtol = 1e-6
+
+        # Exactly one of {fraction, density} required.
+        @test_throws ErrorException EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2)]; zimp=6, mimp=12)
+
+        # Unified electron semantics: a single-ion + electron run resolves as
+        # {main ion, impurity, electron} — electrons in ADDITION to the ions.
+        spe = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=1.0)]; electron=true, zimp=6, mimp=12)
+        @test length(spe) == 3
+        @test !spe[1].electron && spe[1].z == 1                  # main ion unchanged
+        @test spe[1].profiles.ni_spline(0.5) ≈ main.profiles.ni_spline(0.5)
+        @test spe[end].electron                                  # electron appended
+        rm(f)
+    end
+
+    @testset "explicit per-species profiles (HDF5)" begin
+        psi = collect(0.0:0.05:1.0)
+        Ni = @. 1.0e20 * (1 - 0.5psi)
+        ne = @. 1.1e20 * (1 - 0.5psi)
+        Ti = @. 2000.0 * (1 - 0.5psi)
+        Te = @. 2500.0 * (1 - 0.5psi)
+        wE = @. 1.0e4 * (1 - psi)
+        h5 = tempname() * ".h5"
+        HDF5.h5open(h5, "w") do file
+            file["psi"] = psi
+            file["n_e"] = ne
+            file["T_i"] = Ti
+            file["T_e"] = Te
+            file["omega_E"] = wE
+            file["n_i"] = Ni
+            file["n_D"] = 0.6 .* Ni
+            file["n_T"] = 0.4 .* Ni     # unequal explicit shapes
+        end
+        sp = EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_D"), KF.IonSpecies(; z=1, m=3, density="n_T")];
+            electron=false, zimp=6, mimp=12)
+        iD = findfirst(s -> s.m == 2, sp)
+        iT = findfirst(s -> s.m == 3, sp)
+        @test sp[iD].profiles.ni_spline(0.5) ≈ 0.6 * 1.0e20 * (1 - 0.25) rtol = 1e-6
+        @test sp[iT].profiles.ni_spline(0.5) ≈ 0.4 * 1.0e20 * (1 - 0.25) rtol = 1e-6
+        # A missing named profile errors clearly.
+        @test_throws ErrorException EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_He")]; zimp=6, mimp=12)
+        rm(h5)
+    end
+
+    @testset "write_kinetic_h5 species-density roundtrip" begin
+        psi = collect(0.0:0.1:1.0)
+        mkprof(a) = a .* (1 .- 0.5 .* psi)
+        data = EQ.KineticProfileData(; psi=psi, n_i=mkprof(1.0e20), n_e=mkprof(1.1e20),
+            T_i=mkprof(2000.0), T_e=mkprof(2500.0), omega_E=mkprof(1.0e4),
+            species_densities=Dict("n_D" => mkprof(0.6e20), "n_T" => mkprof(0.4e20)),
+            provenance="roundtrip test")
+        h5 = tempname() * ".h5"
+        EQ.write_kinetic_h5(h5, data)
+        back = EQ.read_kinetic_file(h5)
+        @test back.species_densities !== nothing
+        @test sort(collect(keys(back.species_densities))) == ["n_D", "n_T"]
+        @test back.species_densities["n_D"] ≈ data.species_densities["n_D"]
+        @test back.species_densities["n_T"] ≈ data.species_densities["n_T"]
+        # The written file feeds the explicit-density multi-ion input directly.
+        sp = EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_D"), KF.IonSpecies(; z=1, m=3, density="n_T")];
+            zimp=6, mimp=12)
+        @test count(s -> !s.electron && s.z == 1, sp) == 2
+        # Non-`n_*` species names are rejected at write time (they would not round-trip).
+        bad = EQ.KineticProfileData(; psi=psi, n_e=mkprof(1.1e20),
+            species_densities=Dict("deuterium" => mkprof(0.6e20)))
+        @test_throws ErrorException EQ.write_kinetic_h5(tempname() * ".h5", bad)
+        rm(h5)
+    end
+
+    @testset "input-validation guards" begin
+        psi = collect(0.0:0.05:1.0)
+        Ni = @. 1.0e20 * (1 - 0.5psi)
+        ne = @. 1.1e20 * (1 - 0.5psi)
+        Ti = @. 2000.0 * (1 - 0.5psi)
+        Te = @. 2500.0 * (1 - 0.5psi)
+        wE = @. 1.0e4 * (1 - psi)
+        h5 = tempname() * ".h5"
+        HDF5.h5open(h5, "w") do file
+            file["psi"] = psi
+            file["n_i"] = Ni
+            file["n_e"] = ne
+            file["T_i"] = Ti
+            file["T_e"] = Te
+            file["omega_E"] = wE
+            file["n_T"] = 0.4 .* Ni
+            file["n_dilute"] = 0.05 .* Ni
+        end
+        DT(f1, f2) = [KF.IonSpecies(; z=1, m=2, fraction=f1), KF.IonSpecies(; z=1, m=3, fraction=f2)]
+        # An all-fraction list must account for the full n_i (the impurity is set by n_i/n_e, not a shortfall).
+        @test_throws ErrorException EQ.resolve_ntv_species(h5, DT(0.5, 0.4); zimp=6, mimp=12)
+        # Fractions are shares of n_i, so summing above 1 is fatal in any list, mixed included.
+        @test_throws ErrorException EQ.resolve_ntv_species(h5, DT(0.7, 0.5); zimp=6, mimp=12)
+        @test_throws ErrorException EQ.resolve_ntv_species(h5,
+            [KF.IonSpecies(; z=1, m=2, fraction=1.2), KF.IonSpecies(; z=1, m=3, density="n_T")]; zimp=6, mimp=12)
+        # A partial fraction alongside an explicit density is legal — the n_i/n_e deficit is the impurity.
+        sp = EQ.resolve_ntv_species(h5,
+            [KF.IonSpecies(; z=1, m=2, fraction=0.5), KF.IonSpecies(; z=1, m=3, density="n_T")]; zimp=6, mimp=12)
+        @test count(s -> !s.electron && s.z == 1, sp) == 2
+        # Main-ion z must stay below the impurity charge (zpitch closure pole at Zeff = zimp).
+        @test_throws ErrorException EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=6, m=12, fraction=1.0)]; zimp=6, mimp=12)
+        # A dilute main-ion mix drives Zeff toward zimp and warns about the zpitch pole.
+        @test_logs (:warn, r"zpitch closure near its pole") match_mode = :any EQ.resolve_ntv_species(
+            h5, [KF.IonSpecies(; z=1, m=2, density="n_dilute")]; zimp=6, mimp=12)
+        # Vacuum point: neutral composition, no impurity, finite zpitch.
+        @test EQ.multi_ion_composition([1], [1.0e19], 0.0; zimp=6, mimp=12) == (1.0, 1.0, 1.0e19, 0.0)
+        rm(h5)
+    end
+
+    @testset "combine_species_states" begin
+        # Combined profile = Σ species dT/dψ interpolated onto the union grid; total = Σ total.
+        mk(psi, dt) = (
+            st = KF.KineticForcesState();
+            st.method_results["fgar"] = KF.MethodResult(; method="fgar", nn=1,
+                total_torque=ComplexF64(sum(dt)), psi_grid=collect(Float64, psi),
+                dtdpsi=ComplexF64.(dt), t_cumulative=zeros(ComplexF64, length(dt)));
+            st
+        )
+        s1 = mk(0.0:0.1:1.0, fill(1.0, 11))
+        s2 = mk(0.05:0.1:0.95, fill(2.0, 10))     # different (offset) grid
+        c = KF.combine_species_states([s1, s2])
+        r = c.method_results["fgar"]
+        @test real(r.total_torque) ≈ 31.0         # 11·1 + 10·2
+        @test maximum(abs.(r.dtdpsi)) > 0          # not the all-zeros regression
+        overlap = [real(d) for (p, d) in zip(r.psi_grid, r.dtdpsi) if 0.1 <= p <= 0.9]
+        @test all(x -> isapprox(x, 3.0; atol=1e-9), overlap)   # 1 + 2 on the overlap
+        # A duplicate ψ node in one species' grid must not produce NaN (zero-width interpolation guard).
+        s3 = mk([0.0, 0.5, 0.5, 1.0], fill(1.0, 4))
+        c2 = KF.combine_species_states([s1, s3])
+        @test all(isfinite, real.(c2.method_results["fgar"].dtdpsi))
+    end
+end