diff --git a/adept/vfp1d/base.py b/adept/vfp1d/base.py
index 05270232..7f54f490 100644
--- a/adept/vfp1d/base.py
+++ b/adept/vfp1d/base.py
@@ -6,7 +6,7 @@
 from adept.normalization import UREG, laser_normalization
 from adept.utils import filter_scalars
 from adept.vfp1d.grid import Grid
-from adept.vfp1d.helpers import _initialize_total_distribution_, calc_logLambda
+from adept.vfp1d.helpers import _initialize_total_distribution_, calc_logLambda, load_profile_on_grid
 from adept.vfp1d.storage import get_save_quantities, post_process
 from adept.vfp1d.vector_field import OSHUN1D
 
@@ -169,8 +169,54 @@ def init_state_and_args(self) -> dict:
         for field in ["e", "b"]:
             state[field] = jnp.zeros(grid.nx + 1)
 
-        state["Z"] = jnp.ones(grid.nx)
-        state["ni"] = ne_prof / self.cfg["units"]["Z"]
+        ion_cfg = self.cfg["density"].get("ion", None)
+        if ion_cfg is not None:
+            # Spatially varying ion density (in units of n0) from a mass-density
+            # profile. Requires the ion mass number A (or a mixture that provides it).
+            # state["Z"] is relative to units.Z -- FLMCollisions' nuei_coeff already
+            # carries units.Z**2, and the runtime rate goes as (Z**2 * ni) * units.Z**2
+            rho_cfg = ion_cfg["mass_density"]
+            if rho_cfg.get("basis", "file") != "file":
+                raise NotImplementedError("density.ion.mass_density only supports basis: file")
+            rho_q = load_profile_on_grid(rho_cfg, grid.x, self.plasma_norm)
+
+            mixture = ion_cfg.get("mixture", None)
+            if mixture is not None:
+                # Multi-species atomic mixture represented as a single effective ion
+                # that preserves both quasineutrality (Z*ni = ne) and the e-i
+                # collision rate (Z^2*ni = sum_i n_i Z_i^2 = Zeff*ne), with
+                # Zeff = <Z^2>/<Z>. units.Z should be set to Zeff.
+                fracs = np.array([sp["fraction"] for sp in mixture])
+                if not np.isclose(fracs.sum(), 1.0):
+                    raise ValueError(f"Ion mixture fractions must sum to 1, got {fracs.sum()}")
+                Zs = np.array([sp["Z"] for sp in mixture])
+                zbar = float(np.sum(fracs * Zs))
+                zeff = float(np.sum(fracs * Zs**2)) / zbar
+                abar = float(ion_cfg.get("A", np.sum(fracs * np.array([sp["A"] for sp in mixture]))))
+
+                # cross-check: fully ionized quasineutrality between the mass-density
+                # and electron-density profiles
+                n_atom = (rho_q / (abar * UREG.Quantity(1.0, "amu")) / self.plasma_norm.n0).to("").magnitude
+                ne_from_rho = zbar * np.asarray(n_atom)
+                mismatch = float(np.max(np.abs(ne_from_rho - np.asarray(ne_prof))) / np.max(np.asarray(ne_prof)))
+                if mismatch > 0.05:
+                    print(
+                        f"WARNING: ne from the mass density (Zbar * rho / (A*amu)) deviates from the "
+                        f"electron-density profile by up to {mismatch:.1%} of its peak. "
+                        f"Check A, the mixture fractions, and the profile units."
+                    )
+
+                state["Z"] = zeff / self.cfg["units"]["Z"] * jnp.ones(grid.nx)
+                state["ni"] = jnp.asarray(ne_prof) / zeff
+            else:
+                # single species: ni from the mass density and pointwise Z = ne/ni
+                ni_phys = rho_q / (ion_cfg["A"] * UREG.Quantity(1.0, "amu"))
+                ni_prof = jnp.asarray((ni_phys / self.plasma_norm.n0).to("").magnitude)
+                state["Z"] = jnp.asarray(ne_prof) / ni_prof / self.cfg["units"]["Z"]
+                state["ni"] = ni_prof
+        else:
+            state["Z"] = jnp.ones(grid.nx)
+            state["ni"] = ne_prof / self.cfg["units"]["Z"]
 
         self.state = state
         self.args = {"drivers": self.cfg["drivers"]}
diff --git a/adept/vfp1d/grid.py b/adept/vfp1d/grid.py
index 22415485..92baba41 100644
--- a/adept/vfp1d/grid.py
+++ b/adept/vfp1d/grid.py
@@ -49,6 +49,7 @@ class Grid(eqx.Module):
     # Physics / mode parameters stored on the grid for convenience
     nl: int  # number of Legendre harmonics
     boundary: str  # "periodic" or "reflective"
+    geometry: str  # "cartesian" or "spherical"
 
     def __init__(
         self,
@@ -63,12 +64,27 @@ def __init__(
         vmax: float,
         nl: int,
         boundary: str = "periodic",
+        geometry: str = "cartesian",
     ):
+        if geometry not in ("cartesian", "spherical"):
+            raise ValueError(f"Unknown geometry '{geometry}', expected 'cartesian' or 'spherical'")
+        if geometry == "spherical":
+            if xmin < 0.0:
+                raise ValueError(
+                    "Spherical geometry requires xmin >= 0 (the coordinate is the radius r; "
+                    "xmin > 0 gives an annular domain with a reflecting inner wall)"
+                )
+            if boundary != "reflective":
+                raise ValueError(
+                    "Spherical geometry requires boundary = 'reflective' "
+                    "(f10 and E vanish at the inner boundary -- by symmetry when xmin = 0 -- "
+                    "and at the outer wall)"
+                )
         # -- Spatial ----------------------------------------------------------
         self.xmin = xmin
         self.xmax = xmax
         self.nx = nx
-        self.dx = xmax / nx
+        self.dx = (xmax - xmin) / nx
 
         self.x = jnp.linspace(xmin + self.dx / 2, xmax - self.dx / 2, nx)
         self.x_edge = jnp.linspace(xmin, xmax, nx + 1)
@@ -99,6 +115,7 @@ def __init__(
         # -- Physics ----------------------------------------------------------
         self.nl = nl
         self.boundary = boundary
+        self.geometry = geometry
 
     @staticmethod
     def from_config(cfg_grid: dict, norm: PlasmaNormalization) -> Grid:
@@ -112,6 +129,7 @@ def from_config(cfg_grid: dict, norm: PlasmaNormalization) -> Grid:
 
         xmax = normalize(cfg_grid["xmax"], norm, dim="x")
         xmin = normalize(cfg_grid["xmin"], norm, dim="x")
+        geometry = cfg_grid.get("geometry", "cartesian")
         tmax = normalize(cfg_grid["tmax"], norm, dim="t")
         dt = normalize(cfg_grid["dt"], norm, dim="t")
 
@@ -127,7 +145,8 @@ def from_config(cfg_grid: dict, norm: PlasmaNormalization) -> Grid:
             nv=cfg_grid["nv"],
             vmax=vmax,
             nl=cfg_grid["nl"],
-            boundary=cfg_grid.get("boundary", "periodic"),
+            boundary=cfg_grid.get("boundary", "reflective" if geometry == "spherical" else "periodic"),
+            geometry=geometry,
         )
 
     def as_dict(self) -> dict:
diff --git a/adept/vfp1d/heat_flux.py b/adept/vfp1d/heat_flux.py
new file mode 100644
index 00000000..c1e8a548
--- /dev/null
+++ b/adept/vfp1d/heat_flux.py
@@ -0,0 +1,255 @@
+"""Spitzer-Härm and SNB heat-flux diagnostics for VFP-1D.
+
+These are pure-numpy post-processing diagnostics that compute reduced-model heat
+fluxes from the fluid moments of a VFP-1D run so they can be compared against the
+kinetic heat flux carried by f10. Everything is evaluated in code-normalized units
+(v in c, t in 1/w0, n in n0) so the comparison with the kinetic ``q`` is direct.
+
+References:
+
+1. Schurtz, G. P., Nicolaï, Ph. D. & Busquet, M. A nonlocal electron conduction
+   model for multidimensional radiation hydrodynamics codes. Phys. Plasmas 7,
+   4238 (2000).
+2. Brodrick, J. P. et al. Testing nonlocal models of electron thermal conduction
+   for magnetic and inertial confinement fusion applications. Phys. Plasmas 24,
+   092309 (2017).  [The "separated" SNB variant with r = 2 is implemented here.]
+3. Epperlein, E. M. & Haines, M. G. Plasma transport coefficients in a magnetic
+   field by direct numerical solution of the Fokker-Planck equation.
+   Physics of Fluids 29, 1029 (1986).
+"""
+
+import os
+
+import numpy as np
+import xarray as xr
+from matplotlib import pyplot as plt
+from scipy.linalg import solve_banded
+
+from adept.vfp1d.storage import calc_EH
+
+DEFAULT_SNB_CFG = {"ngroups": 32, "beta_max": 30.0, "r": 2.0}
+
+
+def _geometry_factors(x: np.ndarray, dx: float, geometry: str) -> tuple[np.ndarray, np.ndarray]:
+    """Face areas (nx+1,) and cell volumes (nx,) for the 1D divergence operator."""
+    x_edge = np.concatenate([[x[0] - dx / 2], 0.5 * (x[1:] + x[:-1]), [x[-1] + dx / 2]])
+    if geometry == "spherical":
+        area = x_edge**2
+        vol = np.diff(x_edge**3) / 3.0
+    else:
+        area = np.ones_like(x_edge)
+        vol = np.full_like(x, dx)
+    return area, vol
+
+
+def _interp_c2e(f: np.ndarray) -> np.ndarray:
+    """Cell centers (nx,) to edges (nx+1,) with clamped ends."""
+    return np.concatenate([f[:1], 0.5 * (f[1:] + f[:-1]), f[-1:]])
+
+
+def _ddx_c2e(f: np.ndarray, dx: float) -> np.ndarray:
+    """Gradient at interior edges, zero at the boundary edges (reflective/symmetry)."""
+    return np.concatenate([[0.0], np.diff(f) / dx, [0.0]])
+
+
+def local_transport_coefficients(n: np.ndarray, T: np.ndarray, Z: np.ndarray, cfg: dict) -> dict:
+    """Local (x-dependent) transport quantities in code-normalized units.
+
+    The local e-i momentum-relaxation rate is scaled from the reference
+    ``nuei_epphaines_norm`` (evaluated at the reference n, T, Z) using
+    nu_ei ~ Z * n * T^(-3/2). Spatial variation of log-Lambda is neglected.
+
+    Args:
+        n: electron density in units of n0 (nx,)
+        T: electron temperature in code units, T = P/n with v in c (nx,)
+        Z: physical ionization state (nx,)
+        cfg: full config with cfg["units"]["derived"] populated
+    """
+    derived = cfg["units"]["derived"]
+    n_ref = float((derived["ne"] / derived["n0"]).to("").magnitude)
+    T_ref = float(derived["vth_norm"]) ** 2 / 2.0
+    Z_ref = cfg["units"]["Z"]
+    nuei_ref = float(derived["nuei_epphaines_norm"].magnitude)
+
+    nu_ei = nuei_ref * (Z / Z_ref) * (n / n_ref) * (T / T_ref) ** -1.5
+    return {"nu_ei": nu_ei, "kappa_eh": calc_EH(Z, 0.0)}
+
+
+def spitzer_harm_heat_flux(
+    x: np.ndarray, dx: float, n: np.ndarray, T: np.ndarray, Z: np.ndarray, cfg: dict
+) -> np.ndarray:
+    """Local Spitzer-Härm(-Epperlein-Haines) heat flux at cell edges (nx+1,).
+
+    q_SH = -kappa_EH(Z) * n * T / nu_ei(n, T, Z) * dT/dx, which by construction is
+    consistent with the ``kappa`` metric logged by the VFP-1D post-processing
+    (kappa = -q / (n T dT/dx) * nuei_epphaines_norm -> kappa_EH in the local limit).
+    """
+    coeffs = local_transport_coefficients(n, T, Z, cfg)
+    chi_c = coeffs["kappa_eh"] * n * T / coeffs["nu_ei"]
+    return -_interp_c2e(chi_c) * _ddx_c2e(T, dx)
+
+
+def snb_heat_flux(
+    x: np.ndarray,
+    dx: float,
+    n: np.ndarray,
+    T: np.ndarray,
+    Z: np.ndarray,
+    ni: np.ndarray,
+    cfg: dict,
+    geometry: str = "cartesian",
+    snb_cfg: dict | None = None,
+) -> dict:
+    """Multigroup SNB heat flux at cell edges (nx+1,).
+
+    Implements the "separated" SNB variant recommended by Brodrick et al. (2017):
+    for each energy group g (beta = v^2 / 2T, group weights from the Spitzer-Härm
+    spectrum W(beta) = beta^4 exp(-beta) / 24),
+
+        H_g / (lambda_g^ee / r) - div( (xi * lambda_g^ei / 3) grad H_g ) = -div(W_g q_SH)
+
+    with r = 2, xi(Z) = (Z + 0.24)/(Z + 4.2), and
+
+        q_SNB = q_SH - sum_g (xi * lambda_g^ei / 3) grad H_g.
+
+    The mean free paths use the same normalized collision rates as the kinetic
+    solver: nu_ei(v) = nuee_coeff * logLam_ratio * Z^2 * ni / v^3 (the l=1 rate in
+    FLMCollisions) and nu_ee(v) = nuee_coeff * n / v^3, so lambda = v^4 / (...).
+
+    Boundary conditions are zero flux at both ends, matching the reflective /
+    r=0 symmetry boundaries of the spherical solver.
+
+    Returns:
+        dict with "q_snb" (nx+1,), "q_sh" (nx+1,), and "H" (ngroups, nx)
+    """
+    snb = {**DEFAULT_SNB_CFG, **(snb_cfg or {})}
+    ngroups, beta_max, r_snb = int(snb["ngroups"]), float(snb["beta_max"]), float(snb["r"])
+
+    derived = cfg["units"]["derived"]
+    nuee_coeff = float(derived["nuee_coeff"])
+    logLam_ratio = derived["logLam_ratio"]
+    logLam_ratio = float(logLam_ratio.magnitude if hasattr(logLam_ratio, "magnitude") else logLam_ratio)
+
+    area, vol = _geometry_factors(x, dx, geometry)
+    q_sh = spitzer_harm_heat_flux(x, dx, n, T, Z, cfg)
+
+    # group edges in beta = v^2/(2T); logarithmic spacing resolves the low-energy
+    # groups while still covering the suprathermal tail that carries q_SH
+    beta_edges = np.geomspace(beta_max / 400.0, beta_max, ngroups + 1)
+    beta_edges[0] = 0.0
+
+    def _cdf(b):
+        # int_0^b beta^4 exp(-beta)/24 dbeta = 1 - exp(-b)(1 + b + b^2/2 + b^3/6 + b^4/24)
+        return 1.0 - np.exp(-b) * (1.0 + b + b**2 / 2.0 + b**3 / 6.0 + b**4 / 24.0)
+
+    weights = np.diff(_cdf(beta_edges))
+    beta_c = np.sqrt(beta_edges[:-1] * beta_edges[1:])
+    beta_c[0] = 0.5 * beta_edges[1]
+
+    xi = (Z + 0.24) / (Z + 4.2)
+    T_edge = _interp_c2e(T)
+
+    q_nonlocal = np.zeros_like(q_sh)
+    H_all = np.zeros((ngroups, x.size))
+
+    for g in range(ngroups):
+        v2_c = 2.0 * beta_c[g] * T  # v_g^2 at centers
+        v2_e = 2.0 * beta_c[g] * T_edge
+        lam_ee_c = v2_c**2 / (nuee_coeff * n)
+        lam_ei_e = v2_e**2 / (nuee_coeff * logLam_ratio * _interp_c2e(Z**2 * ni))
+        D_e = _interp_c2e(xi) * lam_ei_e / 3.0
+
+        # source: -div(W_g q_SH), finite volume with zero flux through boundary edges
+        q_g = weights[g] * q_sh
+        src = -np.diff(area * q_g) / vol
+
+        # tridiagonal finite-volume diffusion operator
+        w = area[1:-1] * D_e[1:-1] / dx  # interior face conductances (nx-1,)
+        lower = np.concatenate([-w / vol[1:], [0.0]])
+        upper = np.concatenate([[0.0], -w / vol[:-1]])
+        diag = r_snb / lam_ee_c
+        diag[:-1] += w / vol[:-1]
+        diag[1:] += w / vol[1:]
+
+        ab = np.zeros((3, x.size))
+        ab[0, :] = upper
+        ab[1, :] = diag
+        ab[2, :] = lower
+        H_g = solve_banded((1, 1), ab, src)
+        H_all[g] = H_g
+
+        dH_e = _ddx_c2e(H_g, dx)
+        q_nonlocal += D_e * dH_e
+
+    return {"q_snb": q_sh - q_nonlocal, "q_sh": q_sh, "H": H_all}
+
+
+def compare_heat_flux(fields_xr: xr.Dataset, cfg: dict, td: str) -> dict:
+    """Computes Spitzer-Härm and SNB heat fluxes from the saved fluid moments and
+    compares them against the kinetic heat flux for every saved timestep.
+
+    Writes a netCDF with the three fluxes and a comparison plot at the final saved
+    time. Returns a dict with the xr.Dataset and summary metrics.
+    """
+    x = np.asarray(cfg["grid"]["x"])
+    dx = float(cfg["grid"]["dx"])
+    geometry = cfg["grid"].get("geometry", "cartesian")
+    Z_ref = cfg["units"]["Z"]
+    snb_cfg = cfg.get("diagnostics", {}).get("snb", None)
+
+    n_t = np.asarray(fields_xr["fields-n n_c"].data)
+    T_t = np.asarray(fields_xr["fields-T a.u."].data)
+    q_kin_t = np.asarray(fields_xr["fields-q a.u."].data)
+    ni_t = np.asarray(fields_xr["fields-ni a.u."].data)
+    if "fields-Z a.u." in fields_xr:
+        Z_t = np.asarray(fields_xr["fields-Z a.u."].data) * Z_ref
+    else:
+        Z_t = Z_ref * np.ones_like(n_t)
+
+    nt = n_t.shape[0]
+    q_sh_t = np.zeros_like(q_kin_t)
+    q_snb_t = np.zeros_like(q_kin_t)
+    for it in range(nt):
+        res = snb_heat_flux(x, dx, n_t[it], T_t[it], Z_t[it], ni_t[it], cfg, geometry=geometry, snb_cfg=snb_cfg)
+        # fluxes live at edges; average to centers to match the saved kinetic q
+        q_sh_t[it] = 0.5 * (res["q_sh"][1:] + res["q_sh"][:-1])
+        q_snb_t[it] = 0.5 * (res["q_snb"][1:] + res["q_snb"][:-1])
+
+    tax = fields_xr.coords["t (ps)"].data
+    xax = fields_xr.coords["x (um)"].data
+    coords = (("t (ps)", tax), ("x (um)", xax))
+    ds = xr.Dataset(
+        {
+            "q_kinetic": xr.DataArray(q_kin_t, coords=coords),
+            "q_spitzer_harm": xr.DataArray(q_sh_t, coords=coords),
+            "q_snb": xr.DataArray(q_snb_t, coords=coords),
+        }
+    )
+    ds.to_netcdf(os.path.join(td, "binary", "heat-flux-comparison.nc"))
+
+    fig, ax = plt.subplots(1, 1, figsize=(8, 5), tight_layout=True)
+    ax.plot(xax, q_kin_t[-1], label="kinetic (VFP)")
+    ax.plot(xax, q_sh_t[-1], "--", label="Spitzer-Härm")
+    ax.plot(xax, q_snb_t[-1], "-.", label="SNB")
+    xlabel = "r (um)" if geometry == "spherical" else "x (um)"
+    ax.set_xlabel(xlabel)
+    ax.set_ylabel("q (code units)")
+    ax.set_title(f"Heat flux comparison at t = {tax[-1]:.3f} ps")
+    ax.legend()
+    ax.grid()
+    fig.savefig(os.path.join(td, "plots", "fields", "heat-flux-comparison.png"), bbox_inches="tight")
+    plt.close(fig)
+
+    peak = np.argmax(np.abs(q_kin_t[-1]))
+    q_sh_peak = q_sh_t[-1][peak]
+
+    def _safe_ratio(num):
+        # avoid inf metrics when there is no temperature gradient (q_SH = 0)
+        return float(num / q_sh_peak) if abs(q_sh_peak) > 1e-30 else 0.0
+
+    metrics = {
+        "q_ratio_kinetic_over_sh": _safe_ratio(q_kin_t[-1][peak]),
+        "q_ratio_snb_over_sh": _safe_ratio(q_snb_t[-1][peak]),
+    }
+
+    return {"heat_flux": ds, "metrics": metrics}
diff --git a/adept/vfp1d/helpers.py b/adept/vfp1d/helpers.py
index a8f128ab..5dbcc115 100644
--- a/adept/vfp1d/helpers.py
+++ b/adept/vfp1d/helpers.py
@@ -8,7 +8,7 @@
 from scipy.special import gamma
 
 from adept._base_ import get_envelope
-from adept.normalization import PlasmaNormalization, normalize
+from adept.normalization import UREG, PlasmaNormalization, normalize
 
 # ideally this should be passed as as an argument and not re-initialised
 from adept.vfp1d.vector_field import OSHUN1D
@@ -56,7 +56,35 @@ def calc_logLambda(
 
     """
     if isinstance(cfg["units"]["logLambda"], str):
-        if cfg["units"]["logLambda"].casefold() == "nrl":
+        if cfg["units"]["logLambda"].casefold() in ("lee-more", "lee_more", "leemore"):
+            # Lee & More 1984 (Phys. Fluids 27, 1273):
+            #   lnLambda = max(2, 0.5 * ln(1 + (bmax/bmin)^2))
+            # bmax: Debye-Hückel length with electron + ion screening, floored at
+            #       the ion-sphere radius R0
+            # bmin: max(classical closest approach, electron thermal de Broglie)
+            # With the Zeff convention (Z*ni = ne, Z^2*ni = sum_i n_i Z_i^2) the ion
+            # screening term Z^2*ni/Ti is exact for mixtures.
+            ne_cc = ne.to("1/cc").magnitude
+            Te_eV = Te.to("eV").magnitude
+            Ti_eV = UREG.Quantity(cfg["units"]["reference ion temperature"]).to("eV").magnitude
+            ni_cc = ne_cc / Z
+
+            e2 = 1.44e-7  # e^2 in eV cm (Gaussian)
+            hbar_c = 1.9732697e-5  # eV cm
+            me_c2 = 510998.95  # eV
+
+            inv_lD2 = 4 * np.pi * e2 * (ne_cc / Te_eV + Z**2 * ni_cc / Ti_eV)
+            R0 = (3.0 / (4 * np.pi * ni_cc)) ** (1.0 / 3.0)
+            b_max = max(inv_lD2**-0.5, R0)
+
+            b_classical = Z * e2 / (3.0 * Te_eV)
+            b_quantum = hbar_c / (2.0 * np.sqrt(3.0 * Te_eV * me_c2))
+            b_min = max(b_classical, b_quantum)
+
+            logLambda_ei = max(2.0, 0.5 * np.log(1.0 + (b_max / b_min) ** 2))
+            logLambda_ee = logLambda_ei
+
+        elif cfg["units"]["logLambda"].casefold() == "nrl":
             log_ne = np.log(ne.to("1/cc").magnitude)
             log_Te = np.log(Te.to("eV").magnitude)
             log_Z = np.log(Z)
@@ -82,6 +110,67 @@ def calc_logLambda(
     return logLambda_ei, logLambda_ee
 
 
+def _read_two_column_file(path: str) -> tuple[np.ndarray, np.ndarray]:
+    """
+    Reads a two-column (coordinate, value) text file.
+
+    Tries comma-delimited first, then whitespace-delimited. Header rows and any
+    rows that fail to parse as numbers are dropped. If the file has more than two
+    columns, the last two are used (a leading row-index column is thereby ignored).
+    Rows are sorted by coordinate.
+
+    :param path: path to the file
+    :return: (coordinate, value) arrays
+    """
+    data = None
+    for delimiter in (",", None):
+        try:
+            candidate = np.genfromtxt(path, delimiter=delimiter, comments="#")
+        except ValueError:
+            continue
+        if candidate.ndim == 2 and candidate.shape[1] >= 2:
+            data = candidate
+            break
+
+    if data is None:
+        raise ValueError(f"Could not parse two numeric columns from {path}")
+
+    data = data[:, -2:]
+    data = data[~np.isnan(data).any(axis=1)]
+    if data.shape[0] < 2:
+        raise ValueError(f"Fewer than 2 valid (coordinate, value) rows in {path}")
+
+    order = np.argsort(data[:, 0])
+    return data[order, 0], data[order, 1]
+
+
+def load_profile_on_grid(prof_cfg: dict, x_grid: Array, norm: PlasmaNormalization):
+    """
+    Loads a (coordinate, value) profile from a file and interpolates it onto the
+    simulation grid. Values outside the file's coordinate range are clamped to
+    the endpoint values.
+
+    Expected config::
+
+        basis: file
+        path: /path/to/profile.csv
+        units: 1/cm^3          # units of the value column
+        coordinate_units: um   # units of the coordinate column (default um)
+
+    :param prof_cfg: the profile config dict (must contain "path" and "units")
+    :param x_grid: normalized spatial grid to interpolate onto (nx,)
+    :param norm: plasma normalization (provides L0 for coordinate conversion)
+    :return: pint Quantity of interpolated values on the grid (nx,)
+    """
+    x_file, val_file = _read_two_column_file(prof_cfg["path"])
+
+    coord_scale = float((UREG.Quantity(1.0, prof_cfg.get("coordinate_units", "um")) / norm.L0).to("").magnitude)
+    x_file_norm = x_file * coord_scale
+
+    on_grid = np.interp(np.asarray(x_grid), x_file_norm, val_file)
+    return UREG.Quantity(on_grid, prof_cfg["units"])
+
+
 def _initialize_distribution_(
     nv: int,
     m: float = 2.0,
@@ -175,6 +264,16 @@ def _initialize_total_distribution_(
 
                     profs[k] = baseline * (1.0 + amp * jnp.cos(2 * jnp.pi / ll * grid.x))
 
+                elif species_params[k]["basis"] == "file":
+                    prof_q = load_profile_on_grid(species_params[k], grid.x, norm)
+                    if k == "n":
+                        # relative to the reference electron density -- the physical
+                        # scaling to n/n0 happens in init_state_and_args
+                        ref = UREG.Quantity(cfg["units"]["reference electron density"])
+                        profs[k] = np.asarray((prof_q / ref).to("").magnitude)
+                    else:  # k == "T", relative to the reference electron temperature
+                        profs[k] = np.asarray((prof_q / norm.T0).to("").magnitude)
+
                 else:
                     raise NotImplementedError
 
@@ -236,4 +335,11 @@ def _initialize_total_distribution_(
     if not species_found:
         raise ValueError("No species found! Check the config")
 
+    # P1 positivity bound: f = f0 + mu*f1 >= 0 requires |f10| <= f0. The local
+    # (big-dt) initialization overshoots into free-streaming in strongly nonlocal
+    # regions (mean free path >> gradient length), which would immediately drive
+    # f0 negative; clamp to the flux-limited bound and let the kinetics relax.
+    f0_at_edges = np.concatenate([f0[:1], 0.5 * (f0[1:] + f0[:-1]), f0[-1:]], axis=0)
+    f10 = np.clip(f10, -f0_at_edges, f0_at_edges)
+
     return f0, f10, prof_total["n"]
diff --git a/adept/vfp1d/storage.py b/adept/vfp1d/storage.py
index 00491839..be22081d 100644
--- a/adept/vfp1d/storage.py
+++ b/adept/vfp1d/storage.py
@@ -288,7 +288,14 @@ def post_process(soln: Solution, cfg: dict, td: str, args: dict | None = None) -
         "kappa_eh": float(round(calc_EH(cfg["units"]["Z"], 0.0), 4)),
     }
 
-    return {"fields": fields_xr, "dists": f_xr, "metrics": metrics}
+    # Spitzer-Härm / SNB / kinetic heat-flux comparison (lazy import: heat_flux
+    # imports calc_EH from this module)
+    from adept.vfp1d.heat_flux import compare_heat_flux
+
+    hf = compare_heat_flux(fields_xr, cfg, td)
+    metrics.update(hf["metrics"])
+
+    return {"fields": fields_xr, "dists": f_xr, "metrics": metrics, "heat_flux": hf["heat_flux"]}
 
 
 def get_field_save_func(cfg: dict, k: str) -> Callable:
@@ -325,6 +332,7 @@ def fields_save_func(t, y, args):
             temp["e"] = _interp_e2c_(y["e"])
             temp["b"] = _interp_e2c_(y["b"])
             temp["ni"] = y["ni"]
+            temp["Z"] = y["Z"]
 
             return temp
 
diff --git a/adept/vfp1d/vector_field.py b/adept/vfp1d/vector_field.py
index bd9477f9..613759cd 100644
--- a/adept/vfp1d/vector_field.py
+++ b/adept/vfp1d/vector_field.py
@@ -104,6 +104,21 @@ def ddx_e2c(self, f: Array) -> Array:
         """d/dx of edge field, evaluated at centers. (nx+1, ...) -> (nx, ...)"""
         return jnp.diff(f, axis=0) / self.grid.dx
 
+    def div_e2c(self, f: Array) -> Array:
+        """Divergence of an edge flux, evaluated at centers. (nx+1, ...) -> (nx, ...)
+
+        Cartesian: df/dx. Spherical: (1/r^2) d(r^2 f)/dr in conservative
+        finite-volume form, using the edge radii so the flux through r=0
+        vanishes identically.
+        """
+        if self.grid.geometry == "spherical":
+            r_e = self.grid.x_edge.reshape((-1,) + (1,) * (f.ndim - 1))
+            # exact finite-volume cell volume (r_+^3 - r_-^3)/3 so that the
+            # volume-weighted sum of the divergence telescopes (conservation)
+            dvol = jnp.diff(r_e**3, axis=0) / 3.0
+            return jnp.diff(r_e**2 * f, axis=0) / dvol
+        return self.ddx_e2c(f)
+
     def interp_e2c(self, f: Array) -> Array:
         """Average edge values to centers. (nx+1, ...) -> (nx, ...)"""
         return 0.5 * (f[:-1] + f[1:])
@@ -317,7 +332,10 @@ def _vdfdx_(self, t: float, y: dict, args: dict) -> dict:
         f0 = y["f0"]
         f10 = y["f10"]
 
-        df0dt_sa = -self.grid.v[None, :] / 3.0 * self.ddx_e2c(f10)  # (nx+1,nv) -> (nx,nv)
+        # In spherical geometry the f0 equation picks up the geometric 2*f10/r term
+        # via the radial divergence; the f10 equation keeps the plain gradient of f0
+        # (the l=1 angular coupling only involves f2, which is dropped for nl=1).
+        df0dt_sa = -self.grid.v[None, :] / 3.0 * self.div_e2c(f10)  # (nx+1,nv) -> (nx,nv)
         df10dt_sa = -self.grid.v[None, :] * self.ddx_c2e(f0)  # (nx,nv) -> (nx+1,nv)
 
         if self.grid.boundary == "reflective":
diff --git a/configs/vfp-1d/spherical-heatflow.yaml b/configs/vfp-1d/spherical-heatflow.yaml
new file mode 100644
index 00000000..1ff61b0d
--- /dev/null
+++ b/configs/vfp-1d/spherical-heatflow.yaml
@@ -0,0 +1,103 @@
+# Spherical (radial) VFP-1D heat-flow calculation from tabulated profiles.
+#
+# Loads ne(r), Te(r), and mass density rho(r) from two-column CSV files
+# (coordinate, value) and runs the kinetic VFP solver in spherical geometry.
+# Post-processing computes the Spitzer-Härm and SNB heat fluxes from the same
+# profiles and compares them with the kinetic heat flux
+# (plots/fields/heat-flux-comparison.png and binary/heat-flux-comparison.nc).
+#
+# Update the three `path` entries, the reference values in `units`, and the ion
+# mass number `A` for your material before running.
+
+# Material: fully ionized CD (40% C / 60% D atomic mixture) in the conduction
+# zone. Zbar = <Z> = 3 (quasineutrality), Zeff = <Z^2>/<Z> = 5 (e-i collisionality
+# and transport), A = 6.304. units.Z must be Zeff; Zbar enters via the ion mixture.
+
+units:
+  laser_wavelength: 351nm
+  reference electron temperature: 1000eV # choose a typical Te of the profile
+  reference ion temperature: 1000eV
+  reference electron density: 1.0e21/cm^3 # choose a typical ne of the profile
+  Z: 5 # Zeff = <Z^2>/<Z> for the CD mixture
+  Ion: CD
+  logLambda: lee-more
+
+density:
+  quasineutrality: true
+  species-background:
+    noise_seed: 420
+    noise_type: gaussian
+    noise_val: 0.0
+    v0: 0.0
+    T0: 1.0
+    m: 2.0
+    n:
+      basis: file
+      path: /path/to/electrondensity.csv
+      units: 1/cm^3
+      coordinate_units: um
+    T:
+      basis: file
+      path: /path/to/electronTemperature.csv
+      units: eV
+      coordinate_units: um
+  ion:
+    A: 6.304 # mean atomic weight of the CD mixture (amu)
+    mixture: # fully ionized CD: sets Zbar = 3, Zeff = 5
+      - { Z: 6, A: 12.011, fraction: 0.4 } # C
+      - { Z: 1, A: 2.0141, fraction: 0.6 } # D
+    mass_density: # used to cross-check quasineutrality against ne(r)
+      basis: file
+      path: /path/to/massdensity.csv
+      units: g/cm^3
+      coordinate_units: um
+
+grid:
+  dt: 1fs
+  nv: 512
+  nx: 256
+  tmin: 0.
+  tmax: 5ps
+  vmax: 8.0
+  xmin: 0.0um
+  xmax: 500um # set to the outer radius of the profiles
+  nl: 1
+  boundary: reflective
+  geometry: spherical
+
+save:
+  fields:
+    t:
+      tmin: 0.0ps
+      tmax: 5ps
+      nt: 11
+  electron:
+    t:
+      tmin: 0.0ps
+      tmax: 5ps
+      nt: 6
+
+solver: vfp-1d
+
+mlflow:
+  experiment: vfp1d-spherical
+  run: heatflow-profiles
+
+drivers:
+  ex: {}
+  ey: {}
+
+terms:
+  fokker_planck:
+    flm:
+      ee: true
+    f00:
+      model: CoulombianKernel
+      scheme: chang_cooper
+  e_solver: oshun
+
+diagnostics:
+  snb:
+    ngroups: 32
+    beta_max: 30.0
+    r: 2.0
diff --git a/docs/source/solvers/vfp1d/config.md b/docs/source/solvers/vfp1d/config.md
index 335e682a..a2651d4c 100644
--- a/docs/source/solvers/vfp1d/config.md
+++ b/docs/source/solvers/vfp1d/config.md
@@ -41,7 +41,7 @@ Physical unit normalizations for the simulation.
 | `reference electron density` | string | Reference electron density with unit, e.g., `"2.275e21/cm^3"` |
 | `Z` | int | Ionization state |
 | `Ion` | string | Ion species label, e.g., `"Au+"` |
-| `logLambda` | string or float | Coulomb logarithm: `"nrl"` for NRL formula, or a numeric value |
+| `logLambda` | string or float | Coulomb logarithm: `"nrl"` (NRL formulary), `"lee-more"` (Lee & More 1984: $\ln\Lambda = \max[2, \tfrac{1}{2}\ln(1 + b_\max^2/b_\min^2)]$ with Debye-Hückel screening floored at the ion-sphere radius and quantum/classical $b_\min$), or a numeric value |
 
 Example:
 ```yaml
@@ -124,6 +124,48 @@ T:
   width: 20um
 ```
 
+**`file`**: Tabulated profile from a two-column (coordinate, value) CSV or whitespace-delimited text file. Header rows are skipped automatically; the profile is interpolated onto the grid and clamped to the endpoint values outside the tabulated range. Density values are normalized by `units.reference electron density`, temperature values by `units.reference electron temperature`.
+```yaml
+n:
+  basis: file
+  path: /path/to/electrondensity.csv
+  units: 1/cm^3          # units of the value column
+  coordinate_units: um   # units of the coordinate column (default um)
+```
+
+### Ion block (optional)
+
+An optional `density.ion` block sets a spatially varying ion density and ionization state from a tabulated mass-density profile:
+
+```yaml
+density:
+  ion:
+    A: 9.0122            # ion mass number (amu)
+    mass_density:
+      basis: file        # only "file" is supported
+      path: /path/to/massdensity.csv
+      units: g/cm^3
+      coordinate_units: um
+```
+
+This sets $n_i(x) = \rho(x) / (A \, m_u)$ and $Z(x) = n_e(x)/n_i(x)$ pointwise. When the block is omitted, $Z$ is uniform (`units.Z`) and $n_i = n_e / Z$.
+
+For atomic **mixtures** (e.g. fully ionized CD), add a `mixture` list. The plasma is then represented by a single effective ion species that preserves both quasineutrality ($Z n_i = n_e$) and the electron-ion collision rate ($Z^2 n_i = \sum_i n_i Z_i^2 = Z_\text{eff}\, n_e$), with $Z_\text{eff} = \langle Z^2\rangle / \langle Z\rangle$. **Set `units.Z` to $Z_\text{eff}$** — it is the collisional/transport Z that enters the Epperlein-Haines and SNB coefficients. The mass-density profile is used to cross-check full-ionization quasineutrality ($\bar{Z}\rho/(A m_u)$ vs $n_e$); a warning is printed if they disagree by more than 5%.
+
+```yaml
+density:
+  ion:
+    A: 6.304                                # mean atomic weight (optional; default sum of fractions * A)
+    mixture:
+      - { Z: 6, A: 12.011, fraction: 0.4 }  # C
+      - { Z: 1, A: 2.0141, fraction: 0.6 }  # D
+    mass_density:
+      basis: file
+      path: /path/to/massdensity.csv
+      units: g/cm^3
+      coordinate_units: um
+```
+
 ## grid
 
 Simulation grid parameters. All dimensioned quantities use astropy-style unit strings.
@@ -140,6 +182,7 @@ Simulation grid parameters. All dimensioned quantities use astropy-style unit st
 | `xmax` | string | Domain maximum x with unit |
 | `nl` | int | Maximum spherical harmonic order (typically `1`) |
 | `boundary` | string | Spatial boundary condition: `"periodic"` (default) or `"reflective"` |
+| `geometry` | string | Spatial geometry: `"cartesian"` (default) or `"spherical"` |
 
 Example:
 ```yaml
@@ -161,6 +204,15 @@ grid:
 - **`periodic`**: Standard periodic wrapping. Suitable for full-wavelength perturbations.
 - **`reflective`**: Zero-flux boundaries. Vector quantities ($f_{10}$, $E$) are zero at domain edges. Scalar quantities ($f_0$) use mirror ghost cells. Suitable for half-wavelength (cosine) perturbations where symmetry allows simulating half the domain.
 
+### Geometry
+
+- **`cartesian`**: Planar 1D geometry (default). The spatial coordinate is $x$.
+- **`spherical`**: Radial 1D geometry with spherical symmetry. The spatial coordinate is the radius $r$. Requires `xmin >= 0` and `boundary: reflective` (when `xmin: 0` the flux through $r=0$ vanishes by symmetry; `xmin > 0` gives an annular domain with a reflecting inner wall; the outer edge is a reflecting wall). Velocity space is unchanged. For $n_l = 1$ the only modification to the equations is the $f_0$ advection term becoming a radial divergence,
+
+$$\partial_t f_0 + \frac{v}{3}\frac{1}{r^2}\partial_r \left(r^2 f_{10}\right) = C[f_0], \qquad \partial_t f_{10} + v\, \partial_r f_0 - E\, \partial_v f_0 = C[f_{10}]$$
+
+  (the angular coupling term $\propto (1-\mu^2)/r$ only feeds $f_2$ into the $f_{10}$ equation, which is dropped at $n_l = 1$). The divergence is discretized in conservative finite-volume form with exact cell volumes $(r_+^3 - r_-^3)/3$.
+
 ## save
 
 Configures what data to save and at what times.
@@ -312,6 +364,24 @@ terms:
       ee: true
 ```
 
+## diagnostics
+
+Optional post-processing diagnostics configuration. The Spitzer-Härm / SNB / kinetic heat-flux comparison always runs during post-processing (it writes `plots/fields/heat-flux-comparison.png` and `binary/heat-flux-comparison.nc`, and logs the `q_ratio_kinetic_over_sh` and `q_ratio_snb_over_sh` metrics). The `diagnostics.snb` block tunes the SNB model (Schurtz et al. 2000; the "separated" variant recommended by Brodrick et al. 2017):
+
+| Field | Type | Default | Description |
+|-------|------|---------|-------------|
+| `ngroups` | int | `32` | Number of energy groups in $\beta = v^2/(2T)$ |
+| `beta_max` | float | `30.0` | Upper edge of the group grid |
+| `r` | float | `2.0` | e-e absorption scaling of the group mfp, $\lambda_g^{ee}/r$ |
+
+```yaml
+diagnostics:
+  snb:
+    ngroups: 32
+    beta_max: 30.0
+    r: 2.0
+```
+
 ## Complete Examples
 
 ### Epperlein-Short Heat Transport
@@ -321,3 +391,7 @@ See `configs/vfp-1d/epp-short.yaml` - Epperlein-Short heat transport benchmark w
 ### Hotspot Relaxation
 
 See `configs/vfp-1d/hotspot.yaml` - Temperature hotspot relaxation with `tanh` temperature profile.
+
+### Spherical Heat Flow from Tabulated Profiles
+
+See `configs/vfp-1d/spherical-heatflow.yaml` - Radial (spherical) geometry with $n_e(r)$, $T_e(r)$, and mass-density profiles loaded from CSV files, set up for a kinetic vs SNB vs Spitzer-Härm heat-flow comparison.
diff --git a/tests/test_vfp1d/spherical-hotspot.yaml b/tests/test_vfp1d/spherical-hotspot.yaml
new file mode 100644
index 00000000..309fd26b
--- /dev/null
+++ b/tests/test_vfp1d/spherical-hotspot.yaml
@@ -0,0 +1,79 @@
+units:
+  laser_wavelength: 351nm
+  reference electron temperature: 300eV
+  reference ion temperature: 300eV
+  reference electron density: 1.5e21/cm^3
+  Z: 6
+  Ion: Au+
+  logLambda: nrl
+
+
+density:
+  quasineutrality: true
+  species-background:
+    noise_seed: 420
+    noise_type: gaussian
+    noise_val: 0.0
+    v0: 0.0
+    T0: 1.0
+    m: 2.0
+    n:
+      basis: uniform
+      baseline: 1.0
+      bump_or_trough: bump
+      center: 0m
+      rise: 1m
+      bump_height: 0.0
+      width: 1m
+    T:
+      basis: tanh
+      baseline: 1.0
+      bump_or_trough: bump
+      center: 0um
+      rise: 10um
+      bump_height: 0.1
+      width: 80um
+
+grid:
+  dt: 1fs
+  nv: 256
+  nx: 32
+  tmin: 0.ps
+  tmax: 0.3ps
+  vmax: 8.0
+  xmax: 125um
+  xmin: 0.0um
+  nl: 1
+  boundary: reflective
+  geometry: spherical
+
+save:
+  fields:
+    t:
+      tmin: 0.0ps
+      tmax: 0.3ps
+      nt: 6
+  electron:
+    t:
+      tmin: 0.0ps
+      tmax: 0.3ps
+      nt: 3
+
+solver: vfp-1d
+
+mlflow:
+  experiment: vfp1d-tests
+  run: spherical-hotspot
+
+drivers:
+  ex: {}
+  ey: {}
+
+terms:
+  fokker_planck:
+    flm:
+      ee: false
+    f00:
+      model: CoulombianKernel
+      scheme: chang_cooper
+  e_solver: oshun
diff --git a/tests/test_vfp1d/spherical-uniform.yaml b/tests/test_vfp1d/spherical-uniform.yaml
new file mode 100644
index 00000000..a42742e9
--- /dev/null
+++ b/tests/test_vfp1d/spherical-uniform.yaml
@@ -0,0 +1,79 @@
+units:
+  laser_wavelength: 351nm
+  reference electron temperature: 300eV
+  reference ion temperature: 300eV
+  reference electron density: 1.5e21/cm^3
+  Z: 6
+  Ion: Au+
+  logLambda: nrl
+
+
+density:
+  quasineutrality: true
+  species-background:
+    noise_seed: 420
+    noise_type: gaussian
+    noise_val: 0.0
+    v0: 0.0
+    T0: 1.0
+    m: 2.0
+    n:
+      basis: uniform
+      baseline: 1.0
+      bump_or_trough: bump
+      center: 0m
+      rise: 1m
+      bump_height: 0.0
+      width: 1m
+    T:
+      basis: uniform
+      baseline: 1.0
+      bump_or_trough: bump
+      center: 0m
+      rise: 1m
+      bump_height: 0.0
+      width: 1m
+
+grid:
+  dt: 1fs
+  nv: 256
+  nx: 16
+  tmin: 0.ps
+  tmax: 0.1ps
+  vmax: 8.0
+  xmax: 125um
+  xmin: 0.0um
+  nl: 1
+  boundary: reflective
+  geometry: spherical
+
+save:
+  fields:
+    t:
+      tmin: 0.0ps
+      tmax: 0.1ps
+      nt: 6
+  electron:
+    t:
+      tmin: 0.0ps
+      tmax: 0.1ps
+      nt: 3
+
+solver: vfp-1d
+
+mlflow:
+  experiment: vfp1d-tests
+  run: spherical-uniform
+
+drivers:
+  ex: {}
+  ey: {}
+
+terms:
+  fokker_planck:
+    flm:
+      ee: false
+    f00:
+      model: CoulombianKernel
+      scheme: chang_cooper
+  e_solver: oshun
diff --git a/tests/test_vfp1d/test_spherical.py b/tests/test_vfp1d/test_spherical.py
new file mode 100644
index 00000000..cc01ef92
--- /dev/null
+++ b/tests/test_vfp1d/test_spherical.py
@@ -0,0 +1,281 @@
+"""Tests for the spherical (radial) geometry of the VFP-1D solver, the file-based
+profile initialization, and the Spitzer-Härm / SNB heat-flux diagnostics."""
+
+import os
+
+import numpy as np
+import pytest
+import yaml
+
+from adept import ergoExo
+from adept.normalization import UREG, laser_normalization
+
+TEST_DIR = os.path.join(os.getcwd(), "tests", "test_vfp1d")
+
+
+def _load_cfg(config_name: str) -> dict:
+    with open(os.path.join(TEST_DIR, f"{config_name}.yaml")) as fi:
+        return yaml.safe_load(fi)
+
+
+def test_spherical_uniform_preservation():
+    """A uniform plasma in spherical geometry must stay uniform: the geometric
+    2*f10/r term in the f0 equation only acts on a nonzero radial flux."""
+    cfg = _load_cfg("spherical-uniform")
+
+    exo = ergoExo()
+    exo.setup(cfg)
+    _, datasets, _ = exo(None)
+
+    dataT = datasets["fields"]["fields-T keV"].data
+    datan = datasets["fields"]["fields-n n_c"].data
+
+    # spatially uniform at the end and unchanged from the start
+    np.testing.assert_allclose(dataT[-1], np.mean(dataT[-1]), rtol=1e-6)
+    np.testing.assert_allclose(np.mean(dataT[-1]), np.mean(dataT[0]), rtol=1e-6)
+    np.testing.assert_allclose(datan[-1], np.mean(datan[-1]), rtol=1e-6)
+
+
+def test_spherical_hotspot():
+    """A central hotspot in spherical geometry: heat must flow outward (q >= 0),
+    volume-weighted energy must be conserved, and the heat-flux comparison
+    diagnostic must produce all three fluxes."""
+    cfg = _load_cfg("spherical-hotspot")
+
+    exo = ergoExo()
+    exo.setup(cfg)
+    _, datasets, _ = exo(None)
+
+    fields = datasets["fields"]
+    r = fields.coords["x (um)"].data
+    U = fields["fields-U a.u."].data
+
+    # volume-weighted energy conservation (reflective walls, no drivers); use the
+    # exact finite-volume cell volumes so the divergence telescopes
+    dr = r[1] - r[0]
+    r_edge = np.concatenate([[0.0], 0.5 * (r[1:] + r[:-1]), [r[-1] + dr / 2]])
+    vols = np.diff(r_edge**3) / 3.0
+    energy = np.sum(U * vols[None, :], axis=1)
+    np.testing.assert_allclose(energy[-1], energy[0], rtol=1e-3)
+
+    # outward heat flow in the gradient region
+    q = fields["fields-q a.u."].data[-1]
+    assert q[np.argmax(np.abs(q))] > 0.0
+
+    hf = datasets["heat_flux"]
+    for k in ["q_kinetic", "q_spitzer_harm", "q_snb"]:
+        assert np.all(np.isfinite(hf[k].data))
+
+    # both model fluxes should peak outward as well and be within an order of
+    # magnitude of the kinetic flux at the kinetic peak
+    ipk = np.argmax(np.abs(hf["q_kinetic"].data[-1]))
+    q_kin = hf["q_kinetic"].data[-1][ipk]
+    for k in ["q_spitzer_harm", "q_snb"]:
+        ratio = hf[k].data[-1][ipk] / q_kin
+        assert 0.1 < ratio < 10.0, f"{k}/q_kinetic = {ratio}"
+
+
+def test_snb_local_limit():
+    """In the highly collisional (local) limit the SNB flux must reduce to the
+    Spitzer-Härm flux; in the long-mean-free-path limit it must be inhibited."""
+    from adept.vfp1d.heat_flux import snb_heat_flux
+
+    nx = 200
+    dx = 1.0
+    x = np.linspace(dx / 2, nx * dx - dx / 2, nx)
+
+    vth_norm = 0.05
+    T0 = vth_norm**2 / 2.0
+    n = np.ones(nx)
+    Z = np.ones(nx)
+    ni = np.ones(nx)
+    T = T0 * (1.0 + 0.01 * np.cos(2 * np.pi * x / (nx * dx)))
+
+    cfg = {
+        "units": {
+            "Z": 1,
+            "derived": {
+                "ne": UREG.Quantity(1.5e21, "1/cc"),
+                "n0": UREG.Quantity(1.5e21, "1/cc"),
+                "vth_norm": vth_norm,
+                "nuei_epphaines_norm": UREG.Quantity(1e-3, ""),
+                "logLam_ratio": 1.0,
+                "nuee_coeff": None,  # set per case below
+            },
+        }
+    }
+
+    # collisional: group mfps << both the cell size and the gradient length
+    cfg["units"]["derived"]["nuee_coeff"] = 1e3
+    res = snb_heat_flux(x, dx, n, T, Z, ni, cfg)
+    np.testing.assert_allclose(res["q_snb"], res["q_sh"], rtol=0, atol=2e-2 * np.max(np.abs(res["q_sh"])))
+
+    # nonlocal: mfps comparable to the gradient length -> flux inhibition
+    cfg["units"]["derived"]["nuee_coeff"] = 1e-4
+    res_nl = snb_heat_flux(x, dx, n, T, Z, ni, cfg)
+    ipk = np.argmax(np.abs(res_nl["q_sh"]))
+    assert abs(res_nl["q_snb"][ipk]) < abs(res_nl["q_sh"][ipk])
+
+
+def test_file_basis_profiles(tmp_path):
+    """Profiles loaded from CSV files must end up in the initial state with the
+    right normalizations, including the ion mass-density -> ni, Z conversion."""
+    from adept.vfp1d.base import BaseVFP1D
+
+    r_um = np.linspace(0.0, 125.0, 100)
+    ne_cc = 1.5e21 * (1.0 + 0.5 * np.exp(-((r_um / 50.0) ** 2)))
+    Te_eV = 300.0 * (1.0 - 0.3 * r_um / 125.0)
+    A_ion = 9.0122  # beryllium
+    Z_ion = 4.0
+    m_p_g = float(UREG.Quantity(1.0, "amu").to("g").magnitude)
+    rho_gcc = ne_cc / Z_ion * A_ion * m_p_g
+
+    files = {}
+    for name, vals in [("ne", ne_cc), ("Te", Te_eV), ("rho", rho_gcc)]:
+        path = tmp_path / f"{name}.csv"
+        np.savetxt(path, np.stack([r_um, vals], axis=1), delimiter=",", header="r_um,value")
+        files[name] = str(path)
+
+    cfg = _load_cfg("spherical-uniform")
+    cfg["units"]["Z"] = Z_ion
+    cfg["density"]["species-background"]["n"] = {
+        "basis": "file",
+        "path": files["ne"],
+        "units": "1/cm^3",
+        "coordinate_units": "um",
+    }
+    cfg["density"]["species-background"]["T"] = {
+        "basis": "file",
+        "path": files["Te"],
+        "units": "eV",
+        "coordinate_units": "um",
+    }
+    cfg["density"]["ion"] = {
+        "A": A_ion,
+        "mass_density": {
+            "basis": "file",
+            "path": files["rho"],
+            "units": "g/cm^3",
+            "coordinate_units": "um",
+        },
+    }
+
+    module = BaseVFP1D(cfg)
+    module.write_units()
+    module.init_state_and_args()
+
+    grid = module.grid
+    state = module.state
+
+    # density moment of f0 should follow ne(r) (in units of n0)
+    n0_cc = float(module.plasma_norm.n0.to("1/cc").magnitude)
+    n_from_f0 = 4 * np.pi * np.sum(np.asarray(state["f0"]) * np.asarray(grid.v) ** 2, axis=1) * grid.dv
+    r_grid_um = np.asarray(grid.x) * float(module.plasma_norm.L0.to("um").magnitude)
+    ne_expected = np.interp(r_grid_um, r_um, ne_cc) / n0_cc
+    np.testing.assert_allclose(n_from_f0, ne_expected, rtol=2e-2)
+
+    # rho was constructed for uniform Z_ion, so the relative Z must be ~1
+    np.testing.assert_allclose(np.asarray(state["Z"]), 1.0, rtol=1e-6)
+    np.testing.assert_allclose(np.asarray(state["ni"]), ne_expected / Z_ion, rtol=2e-2)
+
+
+def test_ion_mixture(tmp_path):
+    """A fully ionized CD mixture must be represented by an effective ion that
+    preserves quasineutrality (Z*ni = ne) and e-i collisionality (Z^2*ni = Zeff*ne)."""
+    from adept.vfp1d.base import BaseVFP1D
+
+    r_um = np.linspace(0.0, 125.0, 100)
+    ne_cc = 1.0e21 * np.ones_like(r_um)
+    zbar, zeff, abar = 3.0, 5.0, 6.304
+    m_u_g = float(UREG.Quantity(1.0, "amu").to("g").magnitude)
+    rho_gcc = ne_cc / zbar * abar * m_u_g
+
+    files = {}
+    for name, vals in [("ne", ne_cc), ("rho", rho_gcc)]:
+        path = tmp_path / f"{name}.csv"
+        np.savetxt(path, np.stack([r_um, vals], axis=1), delimiter=",", header="r_um,value")
+        files[name] = str(path)
+
+    cfg = _load_cfg("spherical-uniform")
+    cfg["units"]["Z"] = zeff
+    cfg["units"]["reference electron density"] = "1.0e21/cm^3"
+    cfg["units"]["logLambda"] = "lee-more"
+    cfg["density"]["species-background"]["n"] = {
+        "basis": "file",
+        "path": files["ne"],
+        "units": "1/cm^3",
+        "coordinate_units": "um",
+    }
+    cfg["density"]["ion"] = {
+        "A": abar,
+        "mixture": [
+            {"Z": 6, "A": 12.011, "fraction": 0.4},
+            {"Z": 1, "A": 2.0141, "fraction": 0.6},
+        ],
+        "mass_density": {
+            "basis": "file",
+            "path": files["rho"],
+            "units": "g/cm^3",
+            "coordinate_units": "um",
+        },
+    }
+
+    module = BaseVFP1D(cfg)
+    module.write_units()
+    module.init_state_and_args()
+    state = module.state
+
+    Z_phys = np.asarray(state["Z"]) * cfg["units"]["Z"]
+    np.testing.assert_allclose(Z_phys, zeff, rtol=1e-12)
+
+    # quasineutrality: Z * ni = ne
+    n0_cc = float(module.plasma_norm.n0.to("1/cc").magnitude)
+    ne_expected = ne_cc[0] / n0_cc
+    np.testing.assert_allclose(Z_phys * np.asarray(state["ni"]), ne_expected, rtol=1e-6)
+
+    # e-i collisionality: Z^2 * ni = Zeff * ne = sum_i n_i Z_i^2
+    np.testing.assert_allclose(Z_phys**2 * np.asarray(state["ni"]), zeff * ne_expected, rtol=1e-6)
+
+
+def test_lee_more_loglambda():
+    """Lee-More Coulomb logarithm: sensible values at conduction-zone conditions,
+    close to NRL there, and floored at 2 for cold dense plasma."""
+    from adept.vfp1d.helpers import calc_logLambda
+
+    def _cfg(loglam, Ti="1000eV"):
+        return {"units": {"logLambda": loglam, "reference ion temperature": Ti}}
+
+    ne = UREG.Quantity(1.0e21, "1/cc")
+    Te = UREG.Quantity(1000.0, "eV")
+    lm_ei, lm_ee = calc_logLambda(_cfg("lee-more"), ne, Te, Z=5, ion_species="CD")
+    assert lm_ee == lm_ei
+    assert 5.0 < lm_ei < 9.0
+
+    nrl_ei, _ = calc_logLambda(_cfg("nrl"), ne, Te, Z=5, ion_species="CD")
+    assert abs(lm_ei - nrl_ei) / nrl_ei < 0.25
+
+    # cold, dense -> hits the floor of 2
+    lm_cold, _ = calc_logLambda(
+        _cfg("lee-more", Ti="1eV"), UREG.Quantity(1.0e24, "1/cc"), UREG.Quantity(1.0, "eV"), Z=5, ion_species="CD"
+    )
+    assert lm_cold == 2.0
+
+
+def test_spherical_config_validation():
+    """Spherical geometry requires reflective boundaries; xmin > 0 is a valid
+    annular domain with dx = (xmax - xmin)/nx."""
+    from adept.vfp1d.grid import Grid
+
+    cfg = _load_cfg("spherical-uniform")
+    norm = laser_normalization(cfg["units"]["laser_wavelength"], cfg["units"]["reference electron temperature"])
+
+    cfg["grid"]["boundary"] = "periodic"
+    with pytest.raises(ValueError):
+        Grid.from_config(cfg["grid"], norm)
+
+    # annular domain
+    cfg["grid"]["boundary"] = "reflective"
+    cfg["grid"]["xmin"] = "25.0um"
+    grid = Grid.from_config(cfg["grid"], norm)
+    np.testing.assert_allclose(grid.dx, (grid.xmax - grid.xmin) / grid.nx)
+    np.testing.assert_allclose(np.asarray(grid.x_edge)[0], grid.xmin)