diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
new file mode 100644
index 00000000000..084dd40d3b3
--- /dev/null
+++ b/.pre-commit-config.yaml
@@ -0,0 +1,7 @@
+repos:
+  - repo: https://github.com/pre-commit/mirrors-clang-format
+    rev: v12.0.1  # Use the version of clang-format that your project uses
+    hooks:
+      - id: clang-format
+        files: \\.cpp$
+        args: [--style=file]
\ No newline at end of file
diff --git a/openmc/data/mf33_njoy.py b/openmc/data/mf33_njoy.py
new file mode 100644
index 00000000000..64a703742e6
--- /dev/null
+++ b/openmc/data/mf33_njoy.py
@@ -0,0 +1,362 @@
+#!/usr/bin/env python3
+"""
+Minimal NJOY/ERRORR driver to produce MF=33 (multigroup covariance) as text.
+Run a short NJOY chain and return dict with keys: mat, ek, tape33.
+
+Assumptions
+-----------
+- `energy_grid_ev` is a strictly increasing sequence of group boundaries in **eV** (length G+1).
+- Output is relative covariances (IRELCO=1) on an explicit group structure (IGN=1).
+"""
+
+from __future__ import annotations
+
+import os
+import shutil
+import subprocess as sp
+from pathlib import Path
+from tempfile import TemporaryDirectory
+from typing import Any, Dict, List, Optional, Sequence, Union
+
+
+NJOY_tolerance = 0.01
+NJOY_sig0 = [1e10]
+NJOY_thermr_emax = 10.0
+
+# ------------------------- small helpers -------------------------
+
+def _read_mat_from_endf(endf_path: Path) -> int:
+    """Infer MAT from an ENDF-6 evaluation (text)."""
+    first_nonzero = None
+    with endf_path.open("r", errors="ignore") as f:
+        for line in f:
+            if len(line) < 75:
+                continue
+            try:
+                mat = int(line[66:70])
+                mf = int(line[70:72])
+                mt = int(line[72:75])
+            except ValueError:
+                continue
+            if mat > 0 and first_nonzero is None:
+                first_nonzero = mat
+            if mat > 0 and mf == 1 and mt == 451:
+                return mat
+    if first_nonzero is None:
+        raise ValueError(f"Could not infer MAT from ENDF file: {endf_path}")
+    return first_nonzero
+
+
+def _validate_energy_grid_ev(ek: Sequence[float]) -> List[float]:
+    ek = [float(x) for x in ek]
+    if len(ek) < 2:
+        raise ValueError("Energy grid must have at least 2 boundaries (G+1).")
+
+    for i in range(1, len(ek)):
+        if not (ek[i] > ek[i-1]):
+            raise ValueError("Energy grid boundaries must be strictly increasing (in eV).")
+
+    return ek
+
+
+# ------------------------- NJOY deck builder -------------------------
+
+def _moder_input(nin: int, nout: int) -> str:
+    return f"moder\n{nin:d} {nout:d} /\n"
+
+def _reconr_input(endfin, pendfout, mat, err=NJOY_tolerance):
+    return (
+        "reconr\n"
+        f"{endfin:d} {pendfout:d} /\n"
+        f"'mf33'/\n"
+        f"{mat:d} 0 0 /\n"
+        f"{err:g} 0. /\n"
+        "0/\n"
+    )
+
+def _thermr_input(endfin, pendfin, pendfout, mat, temperature,
+                   err=NJOY_tolerance, emax=NJOY_thermr_emax):
+    return (
+        "thermr\n"
+        f"{endfin:d} {pendfin:d} {pendfout:d} /\n"
+        f"0 {mat:d} 20 1 1 0 0 1 221 0 /\n"
+        f"{temperature:.1f} /\n"
+        f"{err:g} {emax:g} /\n"
+    )
+
+def _broadr_input(endfin, pendfin, pendfout, mat, temperature, err=NJOY_tolerance):
+    return (
+        "broadr\n"
+        f"{endfin:d} {pendfin:d} {pendfout:d} /\n"
+        f"{mat:d} 1 0 0 0. /\n"
+        f"{err:g} /\n"
+        f"{temperature:.1f} /\n"
+        "0 /\n"
+    )
+
+def _unresr_input(endfin, pendfin, pendfout, mat, temperature,
+                   sig0=NJOY_sig0):
+    nsig0 = len(sig0)
+    return (
+        "unresr\n"
+        f"{endfin:d} {pendfin:d} {pendfout:d} /\n"
+        f"{mat:d} 1 {nsig0:d} 0 /\n"
+        f"{temperature:.1f} /\n"
+        + " ".join(f"{s:.2E}" for s in sig0) + " /\n"
+        "0 /\n"
+    )
+
+def _purr_input(endfin, pendfin, pendfout, mat, temperature,
+                 sig0=NJOY_sig0, bins=20, ladders=32):
+    nsig0 = len(sig0)
+    return (
+        "purr\n"
+        f"{endfin:d} {pendfin:d} {pendfout:d} /\n"
+        f"{mat:d} 1 {nsig0:d} {bins:d} {ladders:d} 0 /\n"
+        f"{temperature:.1f} /\n"
+        + " ".join(f"{s:.2E}" for s in sig0) + " /\n"
+        "0 /\n"
+    )
+
+
+def _errorr_mf33_input(
+    *,
+    endfin: int,
+    pendfin: int,
+    errorrout: int,
+    mat: int,
+    temperature: float,
+    ek: Sequence[float],
+    iwt: int = 2,
+    spectrum: Optional[Sequence[float]] = None,
+    relative: bool = True,
+    irespr: int = 1,
+    mt: Optional[Union[int, Sequence[int]]] = None,
+    iprint: bool = False,
+) -> str:
+    irelco = 1 if relative else 0
+    iread  = 1 if mt is not None else 0
+    iwt_   = 1 if spectrum is not None else iwt
+    pf     = int(iprint)
+    nk     = len(ek) - 1
+
+    lines = [
+        "errorr",
+        f"{endfin:d} {pendfin:d} 0 {errorrout:d} 0 /",
+        f"{mat:d} 1 {iwt_:d} {pf:d} {irelco} /",   # ign=1 (explicit grid)
+        f"{pf:d} {temperature:.1f} /",
+        f"{iread:d} 33 {irespr:d} /",
+    ]
+
+    # MT list
+    if iread == 1:
+        mtlist = [mt] if isinstance(mt, int) else list(mt)
+        lines.append(f"{len(mtlist):d} 0 /")
+        lines.append(" ".join(str(m) for m in mtlist) + " /")
+
+    # Explicit energy grid
+    lines.append(f"{nk} /")
+    lines.append(" ".join(f"{x:.5e}" for x in ek) + " /")
+
+    # User weight spectrum (TAB1, iwt=1)
+    if iwt_ == 1 and spectrum is not None:
+        n_pairs = len(spectrum) // 2
+        lines.append(f" 0.0 0.0 0 0 1 {n_pairs:d} /")
+        lines.append(f"{n_pairs:d} 1 /")
+        for k in range(n_pairs):
+            lines.append(f"{spectrum[2*k]:.6e} {spectrum[2*k+1]:.6e} /")
+        lines.append("/")
+
+    return "\n".join(lines) + "\n"
+
+
+# ---- deck assembly ---------------------------------------------------------
+
+def _build_deck(
+    mat: int,
+    ek: Sequence[float],
+    temperature: float,
+    *,
+    err: float = NJOY_tolerance,
+    thermr: bool = False,
+    unresr: bool = False,
+    purr: bool = False,
+    iwt: int = 2,
+    spectrum: Optional[Sequence[float]] = None,
+    relative: bool = True,
+    irespr: int = 1,
+    mt: Optional[Union[int, Sequence[int]]] = None,
+    iprint: bool = False,
+) -> str:
+    e = 21      # internal formatted ENDF
+    p = e + 1   # running PENDF tape number
+    deck = _moder_input(20, -e)
+
+    # RECONR
+    deck += _reconr_input(-e, -p, mat=mat, err=err)
+
+    # BROADR
+    if temperature > 0:
+        o = p + 1
+        deck += _broadr_input(-e, -p, -o, mat=mat,
+                               temperature=temperature, err=err)
+        p = o
+
+    # THERMR (free-gas)
+    if thermr:
+        o = p + 1
+        deck += _thermr_input(0, -p, -o, mat=mat,
+                               temperature=temperature, err=err)
+        p = o
+
+    # UNRESR
+    if unresr:
+        o = p + 1
+        deck += _unresr_input(-e, -p, -o, mat=mat,
+                               temperature=temperature)
+        p = o
+
+    # PURR
+    if purr:
+        o = p + 1
+        deck += _purr_input(-e, -p, -o, mat=mat,
+                             temperature=temperature)
+        p = o
+
+    # ERRORR (MF=33)
+    deck += _errorr_mf33_input(
+        endfin=-e, pendfin=-p, errorrout=33,
+        mat=mat, temperature=temperature, ek=ek,
+        iwt=iwt, spectrum=spectrum, relative=relative,
+        irespr=irespr, mt=mt, iprint=iprint,
+    )
+    deck += "stop\n"
+    return deck
+
+
+# ---- NJOY runner -----------------------------------------------------------
+
+def _run_njoy(deck: str, endf: Path, exe: Optional[str]) -> Dict[int, str]:
+    if exe is None:
+        exe = os.environ.get("NJOY")
+        if not exe:
+            raise ValueError(
+                "NJOY executable not provided and $NJOY env var is unset."
+            )
+
+    with TemporaryDirectory() as td:
+        tmpdir = Path(td)
+        shutil.copy(endf, tmpdir / "tape20")
+        (tmpdir / "input").write_text(deck)
+
+        proc = sp.Popen(
+            exe, shell=True, cwd=str(tmpdir),
+            stdin=sp.PIPE, stdout=sp.PIPE, stderr=sp.PIPE,
+        )
+        stdout, stderr = proc.communicate(input=deck.encode())
+
+        if proc.returncode != 0:
+            # Save work dir for debugging
+            dest = Path("njoy_failed_outputs")
+            if dest.exists():
+                shutil.rmtree(dest)
+            shutil.copytree(tmpdir, dest)
+            raise RuntimeError(
+                f"NJOY failed (rc={proc.returncode}). "
+                f"Work dir saved to '{dest}'.\n"
+                f"stderr: {stderr.decode(errors='replace')[:2000]}"
+            )
+
+        tp33 = tmpdir / "tape33"
+        if not tp33.exists():
+            raise RuntimeError("NJOY ran but did not produce tape33.")
+
+        return {33: tp33.read_text()}
+
+
+# ---- public API ------------------------------------------------------------
+
+def generate_errorr_mf33(
+    endf_path: str | Path,
+    energy_grid_ev: Sequence[float],
+    *,
+    njoy_exec: Optional[str] = None,
+    mat: Optional[int] = None,
+    temperature: float = 293.6,
+    # processing chain
+    err: float = NJOY_tolerance,
+    minimal_processing: bool = True,
+    thermr: bool = False,
+    unresr: bool = False,
+    purr: bool = False,
+    # ERRORR options
+    iwt: int = 2,
+    spectrum: Optional[Sequence[float]] = None,
+    relative: bool = True,
+    irespr: int = 1,
+    mt: Optional[Union[int, Sequence[int]]] = None,
+    iprint: bool = False,
+) -> Dict[str, Any]:
+    """Run NJOY/ERRORR and return MF=33 tape33 as text.
+
+    Parameters
+    ----------
+    endf_path
+        Path to ENDF-6 evaluation (text).
+    energy_grid_ev
+        Explicit group boundaries (G+1 values) in eV.
+    njoy_exec
+        NJOY executable/command.  Falls back to $NJOY.
+    mat
+        MAT number.  Inferred from the file when omitted.
+    temperature
+        Processing temperature in K.
+    err
+        Reconstruction tolerance for RECONR / BROADR.
+    minimal_processing
+        If True (default), force thermr/unresr/purr off.
+    thermr, unresr, purr
+        Individual module toggles.  All ignored when
+        ``minimal_processing=True``.
+    iwt
+        Weight function option (default 2 = constant).
+        Ignored when *spectrum* is given.
+    spectrum
+        User weight function as flat [E0, W0, E1, W1, …].
+        Forces iwt=1.
+    relative
+        Produce relative covariances (default True).
+    irespr
+        Resonance-parameter covariance processing
+        (0 = area sensitivity, 1 = 1%% sensitivity).
+    mt
+        Restrict ERRORR to specific reaction numbers.
+    iprint
+        NJOY print flag.
+
+    Returns
+    -------
+    dict
+        Keys: "mat", "ek" (validated eV grid), "tape33" (text).
+    """
+    endf_path = Path(endf_path).resolve()
+    if not endf_path.exists():
+        raise FileNotFoundError(endf_path)
+
+    ek = _validate_energy_grid_ev(energy_grid_ev)
+    if mat is None:
+        mat = _read_mat_from_endf(endf_path)
+
+    if minimal_processing:
+        thermr = unresr = purr = False
+
+    deck = _build_deck(
+        mat=mat, ek=ek, temperature=float(temperature),
+        err=err, thermr=thermr, unresr=unresr, purr=purr,
+        iwt=iwt, spectrum=spectrum, relative=relative,
+        irespr=irespr, mt=mt, iprint=iprint,
+    )
+    print("Running NJOY to produce MF33 covariance data")
+    tapes = _run_njoy(deck, endf_path, exe=njoy_exec)
+
+    return {"mat": int(mat), "ek": ek, "tape33": tapes[33]}
\ No newline at end of file
diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py
index 492cdd7f3b5..f0f4100e5cc 100644
--- a/openmc/data/neutron.py
+++ b/openmc/data/neutron.py
@@ -108,6 +108,7 @@ def __init__(self, name, atomic_number, mass_number, metastable,
         self.reactions = {}
         self._urr = {}
         self._resonances = None
+        self._mg_covariance = None
 
     def __contains__(self, mt):
         return mt in self.reactions
@@ -228,6 +229,30 @@ def urr(self, urr):
             cv.check_type('probability tables', value, ProbabilityTables)
         self._urr = urr
 
+    @property
+    def mg_covariance(self):
+        return self._mg_covariance
+
+    @mg_covariance.setter
+    def mg_covariance(self, cov):
+        if cov is None:
+            self._mg_covariance = None
+            return
+        if not hasattr(cov, 'write_mf33_group'):
+            raise TypeError(
+                "mg_covariance must provide write_mf33_group "
+                "(e.g. NeutronXSCovariances)")
+        self._mg_covariance = cov
+
+    # Optional alias if you already use `data.covariance = ...` somewhere
+    @property
+    def covariance(self):
+        return self._mg_covariance
+    
+    @covariance.setter
+    def covariance(self, cov):
+        self.mg_covariance = cov
+
     @property
     def temperatures(self):
         return [f"{int(round(kT / K_BOLTZMANN))}K" for kT in self.kTs]
@@ -428,7 +453,16 @@ def export_to_hdf5(self, path, mode='a', libver='earliest'):
             if self.fission_energy is not None:
                 fer_group = g.create_group('fission_energy_release')
                 self.fission_energy.to_hdf5(fer_group)
-
+            
+            # Write covariance data (per-reaction MF=33 schema)
+            if self._mg_covariance is not None:
+                cov = self._mg_covariance
+                cov_root = g.require_group("covariance")
+
+                # Prefer the shared writer if available (NeutronXSCovariances)
+                if hasattr(cov, "write_mf33_group"):
+                    cov.write_mf33_group(cov_root)
+        
     @classmethod
     def from_hdf5(cls, group_or_filename):
         """Generate continuous-energy neutron interaction data from HDF5 group
@@ -503,7 +537,13 @@ def from_hdf5(cls, group_or_filename):
         if 'fission_energy_release' in group:
             fer_group = group['fission_energy_release']
             data.fission_energy = FissionEnergyRelease.from_hdf5(fer_group)
-
+        
+        # Read covariance data (per-reaction MF=33 schema)
+        if 'covariance' in group and 'mf33' in group['covariance']:
+            from .xs_covariance_njoy import NeutronXSCovariances
+            data._mg_covariance = NeutronXSCovariances._read_mf33_group(
+                group['covariance']['mf33'], name=name,)
+            
         return data
 
     @classmethod
diff --git a/openmc/data/xs_covariance_njoy.py b/openmc/data/xs_covariance_njoy.py
new file mode 100644
index 00000000000..79ba67d176e
--- /dev/null
+++ b/openmc/data/xs_covariance_njoy.py
@@ -0,0 +1,460 @@
+"""openmc.data.xs_covariance
+
+Tiny model + HDF5 I/O for multigroup **cross section** covariances (ERRORR MF=33).
+
+This is intentionally narrow:
+- neutron MF=33 only
+- explicit energy grid (group edges in eV)
+- **relative** covariances (what ERRORR typically produces with IRELCO=1)
+
+Public API
+----------
+NeutronXSCovariances.from_endf(...)
+NeutronXSCovariances.to_hdf5(...)
+NeutronXSCovariances.from_hdf5(...)
+
+Everything else is internal parsing.
+"""
+
+from __future__ import annotations
+
+import logging
+from dataclasses import dataclass, field
+from pathlib import Path
+from typing import Any, Dict, List, Optional, Sequence, Tuple
+
+import numpy as np
+import scipy.linalg as la
+from .mf33_njoy import generate_errorr_mf33
+
+log = logging.getLogger(__name__)
+
+# -----------------------------------------------------------------------------
+# ERRORR tape33 (ENDF-6 text) parser for MF=33
+# -----------------------------------------------------------------------------
+
+from collections import namedtuple
+
+_CONT = namedtuple("CONT", "C1 C2 L1 L2 N1 N2")
+_LIST = namedtuple("LIST", "C1 C2 L1 L2 NPL N2 B")
+
+
+def _endf_int(field: str) -> int:
+    field = field.strip()
+    if not field:
+        return 0
+    try:
+        return int(field)
+    except ValueError:
+        return int(float(field))
+
+
+def _endf_float(field: str) -> float:
+    s = field.strip()
+    if not s:
+        return 0.0
+    if "e" in s.lower():
+        return float(s)
+    # ENDF-style: mantissa + exponent sign+digits, no 'e'
+    for i in range(len(s) - 1, 0, -1):
+        if s[i] in "+-" and s[i - 1].isdigit():
+            mant = s[:i]
+            exp = s[i:]
+            return float(f"{mant}e{exp}")
+    return float(s)
+
+
+def _parse_endf_ids(line: str) -> Tuple[int, int, int]:
+    if len(line) < 75:
+        return (0, 0, 0)
+    return (int(line[66:70]), int(line[70:72]), int(line[72:75]))
+
+
+def _read_cont_line(line: str) -> _CONT:
+    return _CONT(
+        _endf_float(line[0:11]),
+        _endf_float(line[11:22]),
+        _endf_int(line[22:33]),
+        _endf_int(line[33:44]),
+        _endf_int(line[44:55]),
+        _endf_int(line[55:66]),
+    )
+
+
+def _read_list_record(lines: List[str], ipos: int) -> Tuple[_LIST, int]:
+    C = _read_cont_line(lines[ipos])
+    npl = int(C.N1)
+    ipos += 1
+    vals: List[float] = []
+    nlines = int(np.ceil(npl / 6.0))
+    for _ in range(nlines):
+        ln = lines[ipos]
+        for j in range(6):
+            a = 11 * j
+            b = a + 11
+            if a >= 66:
+                break
+            vals.append(_endf_float(ln[a:b]))
+        ipos += 1
+    vals = vals[:npl]
+    return _LIST(C.C1, C.C2, C.L1, C.L2, npl, C.N2, vals), ipos
+
+
+def _parse_mf33_section_lines(section_lines: List[str], mat: int, mt: int) -> Dict[str, Any]:
+    i = 0
+    C0 = _read_cont_line(section_lines[i])
+    i += 1
+    out: Dict[str, Any] = {"MAT": mat, "MF": 33, "MT": mt, "ZA": C0.C1, "AWR": C0.C2}
+
+    reaction_pairs: Dict[int, np.ndarray] = {}
+    for _ in range(int(C0.N2)):  # number of reaction pairs
+        C = _read_cont_line(section_lines[i])
+        i += 1
+        mt1 = int(C.L2)
+        ng = int(C.N2)
+        M = np.zeros((ng, ng))
+        while True:
+            L, i = _read_list_record(section_lines, i)
+            ngcol = int(L.L1)
+            grow = int(L.N2)
+            gcol = int(L.L2)
+            # LIST rows are 1-based
+            M[grow - 1, gcol - 1 : gcol - 1 + ngcol] = np.asarray(L.B, dtype=float)
+            if grow >= ng:
+                break
+        reaction_pairs[mt1] = M
+
+    out["COVS"] = reaction_pairs
+    return out
+
+
+def parse_errorr_mf33_text(tape33_text: str, mat: Optional[int] = None) -> Dict[int, Dict[str, Any]]:
+    lines = tape33_text.splitlines(True)
+    reactions: Dict[int, Dict[str, Any]] = {}
+
+    current: List[str] = []
+    current_key: Optional[Tuple[int, int, int]] = None
+
+    def flush():
+        nonlocal current, current_key
+        if current_key is None or not current:
+            current = []
+            current_key = None
+            return
+        m, mf, mt = current_key
+        if mf == 33 and mt > 0 and (mat is None or m == mat):
+            reactions[mt] = _parse_mf33_section_lines(current, m, mt)
+        current = []
+        current_key = None
+
+    for ln in lines:
+        m, mf, mt = _parse_endf_ids(ln)
+
+        if m == 0 and mf == 0 and mt == 0:
+            flush()
+            continue
+        if mt == 0:
+            flush()
+            continue
+
+        key = (m, mf, mt)
+        if current_key is None:
+            current_key = key
+        elif key != current_key:
+            flush()
+            current_key = key
+
+        current.append(ln)
+
+    flush()
+    return reactions
+
+
+
+# -----------------------------------------------------------------------------
+# Covariance factor computation (eigendecomposition + QR)
+# -----------------------------------------------------------------------------
+
+@dataclass
+class CovFactorResult:
+    """Container for a covariance factorization and its diagnostics.
+
+    Attributes
+    ----------
+    L : np.ndarray
+        Lower-triangular factor, shape (G, r) where r = effective_rank.
+        Satisfies L @ L.T ≈ original covariance (after any regularization).
+    effective_rank : int
+        Number of positive eigenvalues retained.
+    full_size : int
+        Original matrix dimension G.
+    method : str
+        Which path was taken: "cholesky", "eigen_qr", or "zero_matrix".
+    """
+    L: np.ndarray
+    effective_rank: int
+    full_size: int
+    method: str = "eigen_qr"
+  
+
+def _enforce_symmetry(A: np.ndarray) -> np.ndarray:
+    """Force exact symmetry: A_sym = (A + A^T) / 2."""
+    return 0.5 * (A + A.T)
+
+
+def _strip_zero_variance(A: np.ndarray):
+    """
+    Remove rows/columns with zero diagonal (zero variance).
+    """
+    diag = np.diag(A)
+    nz = np.flatnonzero(diag > 0.0)
+    if len(nz) == 0:
+        return np.empty((0, 0), dtype=A.dtype), nz
+    return A[np.ix_(nz, nz)], nz
+
+def compute_covariance_factor(
+    cov_matrix: np.ndarray,
+    tol: float = 1e-10,
+
+) -> CovFactorResult:
+    """
+    Compute a lower-triangular factor L such that L @ L.T ≈ Sigma.
+    """
+    G = cov_matrix.shape[0]
+
+    # --- 0. Enforce exact symmetry ---
+    A = _enforce_symmetry(np.asarray(cov_matrix, dtype=np.float64))
+
+    # --- 1. Strip zero-variance rows/columns ---
+    Ar, nz = _strip_zero_variance(A)
+
+    if Ar.size == 0:
+        return CovFactorResult(
+            L=np.zeros((G, 0), dtype=np.float64),
+            effective_rank=0,
+            full_size=G,
+            method="zero_matrix",
+            condition_number=np.inf,
+            nonzero_variance_indices=nz,
+        )
+
+    Gr = Ar.shape[0]
+
+    # --- 2a. Fast path: standard Cholesky ---
+    try:
+        Lr = la.cholesky(Ar, lower=True)
+        eig_pos = np.diag(Lr) ** 2
+        cond = float(eig_pos.max() / eig_pos.min()) if eig_pos.min() > 0 else np.inf
+
+        L_full = np.zeros((G, Gr), dtype=np.float64)
+        L_full[nz, :] = Lr
+
+        return CovFactorResult(
+            L=L_full,
+            effective_rank=Gr,
+            full_size=G,
+            method="cholesky"
+        )
+    except la.LinAlgError:
+        pass
+
+    # --- 2b. Fallback: eigendecomposition + QR ---
+    eigenvalues, V = la.eigh(Ar)
+
+    max_eig = eigenvalues.max() if eigenvalues.max() > 0 else 1.0
+    eigenvalues[eigenvalues / max_eig < tol] = 0.0
+    eigenvalues[eigenvalues < 0.0] = 0.0
+
+    pos_mask = eigenvalues > 0.0
+    r = int(pos_mask.sum())
+
+    if r == 0:
+        return CovFactorResult(
+            L=np.zeros((G, 0), dtype=np.float64),
+            effective_rank=0,
+            full_size=G,
+            method="eigen_qr",
+        )
+
+    V_pos = V[:, pos_mask]
+    sqrt_lam = np.sqrt(eigenvalues[pos_mask])
+    Athin = V_pos * sqrt_lam[np.newaxis, :]
+
+    _, R = la.qr(Athin.T, mode="economic")
+    Lr = R.T
+
+    L_full = np.zeros((G, r), dtype=np.float64)
+    L_full[nz, :] = Lr
+
+    return CovFactorResult(
+        L=L_full,
+        effective_rank=r,
+        full_size=G,
+        method="eigen_qr",
+    )
+
+# -----------------------------------------------------------------------------
+# Data model + HDF5 I/O
+# -----------------------------------------------------------------------------
+
+@dataclass
+class NeutronXSCovariances:
+    """Neutron multigroup **relative** covariance data for cross sections (MF=33)."""
+
+    name: str
+    energy_grid_ev: np.ndarray              # shape (G+1,)
+    reactions: Dict[int, Dict[str, Any]]    # MT -> parsed dict from ERRORR
+    mat: Optional[int] = None
+    temperature_k: Optional[float] = None   # what we processed at (single T)
+    factor_results: Optional[Dict[int, Dict[int, CovFactorResult]]] = None 
+
+    @classmethod
+    def from_endf(
+        cls,
+        endf_path: str | Path,
+        energy_grid_ev: Sequence[float],
+        *,
+        njoy_exec: Optional[str] = None,
+        mat: Optional[int] = None,
+        temperature: float = 293.6,
+        name: Optional[str] = None,
+        compute_factors: bool = True,
+        eig_tol: float = 1e-10,
+    ) -> "NeutronXSCovariances":
+        res = generate_errorr_mf33(
+            endf_path,
+            energy_grid_ev,
+            njoy_exec=njoy_exec,
+            mat=mat,
+            temperature=temperature,
+        )
+        mat_used = int(res["mat"])
+        ek = np.asarray(res["ek"], dtype=float)
+
+        reactions = parse_errorr_mf33_text(res["tape33"], mat=mat_used)
+
+        # Pre-compute covariance factors for all sub-blocks
+        fcache: Optional[Dict[int, Dict[int, CovFactorResult]]] = None
+        if compute_factors:
+            fcache = {}
+            for mt_key, sec in reactions.items():
+                fcache[mt_key] = {}
+                for mt1, M in sec.get("COVS", {}).items():
+                    result = compute_covariance_factor(
+                        np.asarray(M, dtype=np.float64),
+                        tol=eig_tol,
+                    )
+                    fcache[mt_key][mt1] = result
+                    log.info(
+                        "  MT %s -> MT1 %s: method=%-9s  rank=%d/%d  "
+                        "cond=%.2e  neg_mass=%.2e%%",
+                        mt_key, mt1,
+                        result.method,
+                        result.effective_rank,
+                        result.full_size,
+                        result.condition_number,
+                        result.negative_eig_mass_pct,
+                    )
+
+        if name is None:
+            name = Path(endf_path).stem
+
+        return cls(
+            name=str(name),
+            energy_grid_ev=ek,
+            reactions=reactions,
+            mat=mat_used,
+            temperature_k=float(temperature),
+            factor_results=fcache,
+        )
+    
+    # -----------------------------------------------------------------
+    # HDF5 writing
+    # -----------------------------------------------------------------
+
+    def to_hdf5(self, filename: str | Path, store_raw_covariance: bool = True) -> None:
+        """
+        Write covariances to a standalone HDF5 file.
+        """
+        import h5py
+
+        filename = Path(filename)
+        filename.parent.mkdir(parents=True, exist_ok=True)
+
+        with h5py.File(filename, "w") as f:
+            self.write_mf33_group(f, store_raw_covariance=store_raw_covariance)
+
+    def write_mf33_group(self, h5_group, store_raw_covariance: bool = True,
+                         eig_tol: float = 1e-10) -> None:
+        """
+        Write the "mf33" sub-tree into an already-open HDF5 group.
+        """
+        if "mf33" in h5_group:
+            del h5_group["mf33"]
+        mf33 = h5_group.create_group("mf33")
+
+        mf33.attrs["format"] = np.bytes_("openmc.mf33.v1")
+        mf33.attrs["source"] = np.bytes_("njoy errorr")
+        mf33.attrs["relative"] = 1  # int flag – portable across languages
+        mf33.attrs["store_raw_covariance"] = int(store_raw_covariance)
+        #mf33.attrs["factorization"] = np.bytes_("eigen_qr_thin")
+        if self.mat is not None:
+            mf33.attrs["mat"] = int(self.mat)
+        if self.temperature_k is not None:
+            mf33.attrs["temperature_k"] = float(self.temperature_k)
+
+        mf33.create_dataset(
+            "energy_grid_ev",
+            data=np.asarray(self.energy_grid_ev, dtype=np.float64),
+        )
+        mf33.create_dataset(
+            "mts",
+            data=np.asarray(sorted(self.reactions.keys()), dtype=np.int32),
+        )
+
+        if store_raw_covariance:
+            greact = mf33.create_group("reactions")
+        gfact = mf33.create_group("factors")
+
+        cached = self.factor_results or {}
+
+        for mt, sec in self.reactions.items():
+            mt_str = str(int(mt))
+
+            if store_raw_covariance:
+                gmt = greact.create_group(mt_str)
+                for attr_name in ("ZA", "AWR"):
+                    if attr_name in sec:
+                        gmt.attrs[attr_name] = float(sec[attr_name])
+
+            gmt_fact = gfact.create_group(mt_str)
+            for attr_name in ("ZA", "AWR"):
+                if attr_name in sec:
+                    gmt_fact.attrs[attr_name] = float(sec[attr_name])
+
+            covs: Dict[int, np.ndarray] = sec.get("COVS", {})
+            for mt1, M in covs.items():
+                M_arr = np.asarray(M, dtype=np.float64)
+                ds_name = str(int(mt1))
+
+                # ---- raw covariance (optional) ----
+                if store_raw_covariance:
+                    gmt.create_dataset(
+                        ds_name,
+                        data=M_arr,
+                        compression="gzip",
+                        shuffle=True,
+                    )
+
+                # ---- Triangular factor (always stored) ----
+                result = cached.get(mt, {}).get(mt1, None)
+                if result is None:
+                    result = compute_covariance_factor(
+                        M_arr, tol=eig_tol,)
+
+                gmt_fact.create_dataset(
+                    ds_name,
+                    data=result.L,
+                    compression="gzip",
+                    shuffle=True,
+                )
\ No newline at end of file
diff --git a/src/simulation.cpp b/src/simulation.cpp
index 4fad196a604..1023b3d321e 100644
--- a/src/simulation.cpp
+++ b/src/simulation.cpp
@@ -160,8 +160,8 @@ int openmc_simulation_init()
   // load weight windows from file
   if (!settings::weight_windows_file.empty()) {
     openmc_weight_windows_import(settings::weight_windows_file.c_str());
-  }
-
+  } 
+  
   // Set flag indicating initialization is done
   simulation::initialized = true;
   return 0;