diff --git a/README.md b/README.md
index 82a6af4..7ba5edf 100644
--- a/README.md
+++ b/README.md
@@ -24,6 +24,7 @@ python -m sphinx -b html docs docs/_build/html   # open docs/_build/html/index.h
 - Thermochemistry calculations for molecular systems
 - DFTB+ geometry optimization, Hessian, and normal-mode integration
 - Stepwise and overall reference-calibrated reduction potentials
+- Structure-aware redox datasets, grouped delta-model validation, and Pareto selection
 - Readers for DFTB+ `.gen`, XYZ, and vibrational frequency files
 - Runtime type checking for public API calls
 - Test coverage for parsing, thermochemistry, and optional DFTB+ execution paths
diff --git a/ThermoScreening/thermo/__init__.py b/ThermoScreening/thermo/__init__.py
index 1130f92..94552ed 100644
--- a/ThermoScreening/thermo/__init__.py
+++ b/ThermoScreening/thermo/__init__.py
@@ -16,3 +16,11 @@
 from .ensemble import boltzmann_weights, ensemble_free_energy, lowest_gibbs, EnsembleThermo
 from .kinetics import eyring_rate_constant, wigner_tunneling_correction
 from .pka import pKa, calibrate_proton_reference, PROTON_AQUEOUS_FREE_ENERGY_KCAL
+from .redox_screening import (
+    DeltaRedoxModel,
+    audit_redox_dataset,
+    balanced_group_folds,
+    canonical_structure_identity,
+    grouped_delta_validation,
+    pareto_front,
+)
diff --git a/ThermoScreening/thermo/redox_screening.py b/ThermoScreening/thermo/redox_screening.py
new file mode 100644
index 0000000..c39d114
--- /dev/null
+++ b/ThermoScreening/thermo/redox_screening.py
@@ -0,0 +1,557 @@
+"""Dataset-level tools for calibrated redox screening."""
+
+from collections import Counter
+from collections.abc import Iterable, Mapping
+from dataclasses import dataclass, replace
+import math
+
+import numpy as np
+from rdkit import Chem, rdBase
+from rdkit.Chem import inchi
+
+from ThermoScreening.exceptions import TSValueError
+
+
+@dataclass(frozen=True)
+class MoleculeIdentity:
+    """Canonical structure identity derived from a SMILES string."""
+
+    canonical_smiles: str
+    structure_id: str
+
+
+@dataclass(frozen=True)
+class DuplicateConflict:
+    """Duplicate structure whose reported values disagree beyond tolerance."""
+
+    structure_id: str
+    indices: tuple[int, ...]
+    spreads: tuple[tuple[str, float], ...]
+
+
+@dataclass(frozen=True)
+class DatasetAudit:
+    """Canonical identities and duplicate findings for a redox dataset."""
+
+    identities: tuple[MoleculeIdentity, ...]
+    duplicate_groups: tuple[tuple[int, ...], ...]
+    conflicts: tuple[DuplicateConflict, ...]
+
+    @property
+    def unique_structures(self):
+        """Number of unique molecular structures in the dataset."""
+
+        return len({identity.structure_id for identity in self.identities})
+
+
+@dataclass(frozen=True)
+class DeltaPrediction:
+    """Corrected potentials and their model-domain diagnostics."""
+
+    corrected_potential: np.ndarray
+    correction: np.ndarray
+    uncertainty: np.ndarray
+    extrapolated: np.ndarray
+    unknown_features: tuple[tuple[str, ...], ...]
+
+
+@dataclass(frozen=True)
+class CrossValidationResult:
+    """Out-of-fold diagnostics from grouped delta-model validation."""
+
+    predictions: np.ndarray
+    errors: np.ndarray
+    folds: tuple[object, ...]
+    mae: float
+    rmse: float
+    bias: float
+
+
+def canonical_structure_identity(smiles):
+    """Return canonical isomeric SMILES and an InChIKey for ``smiles``."""
+
+    if not isinstance(smiles, str) or not smiles.strip():
+        raise TSValueError("SMILES must be a non-empty string.")
+    molecule = Chem.MolFromSmiles(smiles)  # pylint: disable=no-member
+    if molecule is None:
+        raise TSValueError(f"Invalid SMILES: {smiles!r}.")
+
+    canonical = Chem.MolToSmiles(  # pylint: disable=no-member
+        molecule, canonical=True, isomericSmiles=True
+    )
+    with rdBase.BlockLogs():  # pylint: disable=c-extension-no-member
+        fixed_h_inchi = inchi.MolToInchi(molecule, options="/FixedH")
+    structure_id = inchi.InchiToInchiKey(fixed_h_inchi)
+    if not structure_id:
+        raise TSValueError(f"Could not derive an InChIKey from SMILES: {smiles!r}.")
+    return MoleculeIdentity(canonical, structure_id)
+
+
+def audit_redox_dataset(smiles, potentials=None, tolerance=0.02):
+    """
+    Canonicalize structures and identify duplicate measurement conflicts.
+
+    ``potentials`` maps labels such as ``"E1"`` to sequences in volts. Missing
+    values may be ``None`` or ``NaN``. Duplicate structures are always reported;
+    a conflict is reported when any potential spans more than ``tolerance``.
+    """
+
+    smiles = tuple(smiles)
+    if not math.isfinite(tolerance) or tolerance < 0.0:
+        raise TSValueError("tolerance must be a finite, non-negative value.")
+
+    values = {}
+    for name, column in (potentials or {}).items():
+        column = tuple(column)
+        if len(column) != len(smiles):
+            raise TSValueError(
+                f"Potential column {name!r} has {len(column)} rows; expected {len(smiles)}."
+            )
+        parsed = []
+        for index, value in enumerate(column):
+            if value is None or value == "":
+                parsed.append(np.nan)
+                continue
+            try:
+                number = float(value)
+            except (TypeError, ValueError) as exc:
+                raise TSValueError(
+                    f"Potential {name!r} at row {index} is not numeric: {value!r}."
+                ) from exc
+            if math.isinf(number):
+                raise TSValueError(
+                    f"Potential {name!r} at row {index} must not be infinite."
+                )
+            parsed.append(number)
+        values[str(name)] = tuple(parsed)
+
+    identities = tuple(canonical_structure_identity(value) for value in smiles)
+    grouped = {}
+    for index, identity in enumerate(identities):
+        grouped.setdefault(identity.structure_id, []).append(index)
+
+    duplicate_groups = tuple(
+        tuple(indices) for indices in grouped.values() if len(indices) > 1
+    )
+    conflicts = []
+    for indices in duplicate_groups:
+        spreads = []
+        for name, column in values.items():
+            finite = [column[index] for index in indices if math.isfinite(column[index])]
+            if len(finite) > 1:
+                spread = max(finite) - min(finite)
+                if spread > tolerance:
+                    spreads.append((name, spread))
+        if spreads:
+            conflicts.append(
+                DuplicateConflict(
+                    identities[indices[0]].structure_id,
+                    indices,
+                    tuple(spreads),
+                )
+            )
+
+    return DatasetAudit(identities, duplicate_groups, tuple(conflicts))
+
+
+def _as_finite_vector(values, name):
+    vector = np.asarray(values, dtype=float)
+    if vector.ndim != 1:
+        raise TSValueError(f"{name} must be one-dimensional.")
+    if not np.all(np.isfinite(vector)):
+        raise TSValueError(f"{name} must contain only finite values.")
+    return vector
+
+
+def _group_counts(groups):
+    if isinstance(groups, Mapping):
+        items = groups.items()
+    else:
+        if isinstance(groups, (str, bytes)) or not isinstance(groups, Iterable):
+            raise TSValueError(
+                "Each functional-group descriptor must be a mapping or an iterable of names."
+            )
+        items = Counter(groups).items()
+
+    counts = {}
+    for name, value in items:
+        if not isinstance(name, str) or not name:
+            raise TSValueError("Functional-group names must be non-empty strings.")
+        try:
+            count = float(value)
+        except (TypeError, ValueError) as exc:
+            raise TSValueError(f"Count for functional group {name!r} is not numeric.") from exc
+        if not math.isfinite(count) or count < 0.0 or not count.is_integer():
+            raise TSValueError(
+                f"Count for functional group {name!r} must be a non-negative integer."
+            )
+        if count:
+            counts[name] = int(count)
+    return counts
+
+
+def _normalise_descriptors(descriptors, expected_length):
+    descriptors = tuple(_group_counts(groups) for groups in descriptors)
+    if len(descriptors) != expected_length:
+        raise TSValueError(
+            f"functional_groups has {len(descriptors)} rows; expected {expected_length}."
+        )
+    return descriptors
+
+
+def _pair_value(counts, first, second):
+    if first == second:
+        count = counts.get(first, 0.0)
+        return count * (count - 1.0) / 2.0
+    return counts.get(first, 0.0) * counts.get(second, 0.0)
+
+
+def _available_interactions(descriptors, minimum_count):
+    group_names = sorted({name for row in descriptors for name in row})
+    interactions = []
+    for first_index, first in enumerate(group_names):
+        for second in group_names[first_index:]:
+            support = sum(_pair_value(row, first, second) for row in descriptors)
+            if support >= minimum_count:
+                interactions.append((first, second))
+    return tuple(interactions)
+
+
+def _design_matrix(
+    approximate,
+    descriptors,
+    approximate_center,
+    group_names,
+    interactions,
+):
+    matrix = np.ones(
+        (len(approximate), 2 + len(group_names) + len(interactions)), dtype=float
+    )
+    matrix[:, 1] = approximate - approximate_center
+    for column, name in enumerate(group_names, start=2):
+        matrix[:, column] = [row.get(name, 0.0) for row in descriptors]
+    offset = 2 + len(group_names)
+    for column, (first, second) in enumerate(interactions, start=offset):
+        matrix[:, column] = [
+            _pair_value(row, first, second) for row in descriptors
+        ]
+    return matrix
+
+
+@dataclass(frozen=True)
+class DeltaRedoxModel:
+    """Transparent linear correction from approximate to reference potentials."""
+
+    approximate_center: float
+    group_names: tuple[str, ...]
+    interactions: tuple[tuple[str, str], ...]
+    coefficients: np.ndarray
+    covariance_factor: np.ndarray
+    residual_std: float
+    approximate_range: tuple[float, float]
+    ridge: float
+    include_interactions: bool
+    validation: CrossValidationResult | None = None
+
+    @classmethod
+    def fit(
+        cls,
+        approximate,
+        reference,
+        functional_groups,
+        *,
+        include_interactions=False,
+        min_interaction_count=5,
+        ridge=1.0e-8,
+        validation_groups=None,
+    ):
+        """
+        Fit ``reference - approximate`` using group counts and optional pairs.
+
+        ``validation_groups`` assigns each row to a held-out fold. Rows with the
+        same label are never split between training and validation.
+        """
+
+        approximate = _as_finite_vector(approximate, "approximate")
+        reference = _as_finite_vector(reference, "reference")
+        if len(reference) != len(approximate):
+            raise TSValueError("approximate and reference must have the same length.")
+        descriptors = _normalise_descriptors(functional_groups, len(approximate))
+        model = cls._fit(
+            approximate,
+            reference,
+            descriptors,
+            include_interactions=include_interactions,
+            min_interaction_count=min_interaction_count,
+            ridge=ridge,
+        )
+        if validation_groups is None:
+            return model
+        validation = grouped_delta_validation(
+            approximate,
+            reference,
+            descriptors,
+            validation_groups,
+            include_interactions=include_interactions,
+            min_interaction_count=min_interaction_count,
+            ridge=ridge,
+        )
+        return replace(model, validation=validation)
+
+    @classmethod
+    def _fit(
+        cls,
+        approximate,
+        reference,
+        descriptors,
+        *,
+        include_interactions,
+        min_interaction_count,
+        ridge,
+    ):
+        if len(approximate) < 2:
+            raise TSValueError("At least two matched rows are required to fit a delta model.")
+        if not math.isfinite(ridge) or ridge < 0.0:
+            raise TSValueError("ridge must be a finite, non-negative value.")
+        if (
+            isinstance(min_interaction_count, bool)
+            or not isinstance(min_interaction_count, int)
+            or min_interaction_count < 1
+        ):
+            raise TSValueError("min_interaction_count must be at least 1.")
+
+        group_names = tuple(sorted({name for row in descriptors for name in row}))
+        interactions = (
+            _available_interactions(descriptors, min_interaction_count)
+            if include_interactions
+            else ()
+        )
+        center = float(np.mean(approximate))
+        matrix = _design_matrix(
+            approximate, descriptors, center, group_names, interactions
+        )
+        target = reference - approximate
+        penalty = np.eye(matrix.shape[1]) * ridge
+        penalty[0, 0] = 0.0
+        cross_product = matrix.T @ matrix
+        normal = cross_product + penalty
+        inverse = np.linalg.pinv(normal)
+        coefficients = inverse @ matrix.T @ target
+        residual = target - matrix @ coefficients
+        degrees_of_freedom = max(len(target) - np.linalg.matrix_rank(matrix), 1)
+        residual_std = float(np.sqrt(np.sum(residual**2) / degrees_of_freedom))
+        return cls(
+            center,
+            group_names,
+            interactions,
+            coefficients,
+            inverse @ cross_product @ inverse,
+            residual_std,
+            (float(np.min(approximate)), float(np.max(approximate))),
+            float(ridge),
+            bool(include_interactions),
+        )
+
+    @property
+    def correction_coefficients(self):
+        """Model coefficients with an uncentered approximate-potential term."""
+
+        result = {
+            "intercept": float(
+                self.coefficients[0]
+                - self.coefficients[1] * self.approximate_center
+            ),
+            "approximate_potential": float(self.coefficients[1]),
+        }
+        offset = 2
+        for name, value in zip(self.group_names, self.coefficients[offset:]):
+            result[f"group:{name}"] = float(value)
+        offset += len(self.group_names)
+        for pair, value in zip(self.interactions, self.coefficients[offset:]):
+            result[f"pair:{pair[0]}|{pair[1]}"] = float(value)
+        return result
+
+    def predict(self, approximate, functional_groups):
+        """Correct approximate potentials and flag model-domain extrapolation."""
+
+        approximate = _as_finite_vector(approximate, "approximate")
+        descriptors = _normalise_descriptors(functional_groups, len(approximate))
+        matrix = _design_matrix(
+            approximate,
+            descriptors,
+            self.approximate_center,
+            self.group_names,
+            self.interactions,
+        )
+        correction = matrix @ self.coefficients
+        leverage = np.einsum(
+            "ij,jk,ik->i", matrix, self.covariance_factor, matrix
+        )
+        uncertainty_floor = self.residual_std
+        if self.validation is not None:
+            uncertainty_floor = max(uncertainty_floor, self.validation.rmse)
+        uncertainty = uncertainty_floor * np.sqrt(np.maximum(1.0 + leverage, 1.0))
+
+        known_groups = set(self.group_names)
+        known_pairs = set(self.interactions)
+        unknown = []
+        for potential, row in zip(approximate, descriptors):
+            missing = {f"group:{name}" for name in row if name not in known_groups}
+            if self.include_interactions:
+                row_pairs = _available_interactions((row,), 1)
+                missing.update(
+                    f"pair:{first}|{second}"
+                    for first, second in row_pairs
+                    if (first, second) not in known_pairs
+                )
+            if potential < self.approximate_range[0] or potential > self.approximate_range[1]:
+                missing.add("approximate_potential_range")
+            unknown.append(tuple(sorted(missing)))
+
+        extrapolated = np.asarray([bool(items) for items in unknown], dtype=bool)
+        uncertainty[extrapolated] = np.nan
+        return DeltaPrediction(
+            approximate + correction,
+            correction,
+            uncertainty,
+            extrapolated,
+            tuple(unknown),
+        )
+
+
+def grouped_delta_validation(
+    approximate,
+    reference,
+    functional_groups,
+    validation_groups,
+    *,
+    include_interactions=False,
+    min_interaction_count=5,
+    ridge=1.0e-8,
+):
+    """Validate a delta model while keeping each labelled group in one fold."""
+
+    approximate = _as_finite_vector(approximate, "approximate")
+    reference = _as_finite_vector(reference, "reference")
+    if len(reference) != len(approximate):
+        raise TSValueError("approximate and reference must have the same length.")
+    descriptors = _normalise_descriptors(functional_groups, len(approximate))
+    labels = tuple(validation_groups)
+    if len(labels) != len(approximate):
+        raise TSValueError(
+            f"validation_groups has {len(labels)} rows; expected {len(approximate)}."
+        )
+    try:
+        unique_labels = tuple(dict.fromkeys(labels))
+    except TypeError as exc:
+        raise TSValueError("validation_groups must contain hashable labels.") from exc
+    if len(unique_labels) < 2:
+        raise TSValueError("validation_groups must contain at least two distinct labels.")
+
+    predictions = np.empty(len(approximate), dtype=float)
+    for label in unique_labels:
+        test = np.fromiter((candidate == label for candidate in labels), dtype=bool)
+        train = ~test
+        if np.count_nonzero(train) < 2:
+            raise TSValueError(
+                f"Validation fold {label!r} leaves fewer than two training rows."
+            )
+        model = DeltaRedoxModel._fit(
+            approximate[train],
+            reference[train],
+            tuple(row for row, keep in zip(descriptors, train) if keep),
+            include_interactions=include_interactions,
+            min_interaction_count=min_interaction_count,
+            ridge=ridge,
+        )
+        prediction = model.predict(
+            approximate[test],
+            tuple(row for row, keep in zip(descriptors, test) if keep),
+        )
+        predictions[test] = prediction.corrected_potential
+
+    errors = predictions - reference
+    return CrossValidationResult(
+        predictions,
+        errors,
+        labels,
+        float(np.mean(np.abs(errors))),
+        float(np.sqrt(np.mean(errors**2))),
+        float(np.mean(errors)),
+    )
+
+
+def balanced_group_folds(group_ids, n_folds=5):
+    """Assign complete groups to deterministic folds with similar row counts."""
+
+    group_ids = tuple(group_ids)
+    if isinstance(n_folds, bool) or not isinstance(n_folds, int) or n_folds < 2:
+        raise TSValueError("n_folds must be an integer of at least 2.")
+    try:
+        counts = Counter(group_ids)
+    except TypeError as exc:
+        raise TSValueError("group_ids must contain hashable values.") from exc
+    if len(counts) < n_folds:
+        raise TSValueError(
+            f"Cannot build {n_folds} folds from only {len(counts)} distinct groups."
+        )
+
+    loads = [0] * n_folds
+    assignment = {}
+    ordered = sorted(counts.items(), key=lambda item: (-item[1], repr(item[0])))
+    for group_id, count in ordered:
+        fold = min(range(n_folds), key=lambda index: (loads[index], index))
+        assignment[group_id] = fold
+        loads[fold] += count
+    return tuple(assignment[group_id] for group_id in group_ids)
+
+
+def pareto_front(records, objectives):
+    """
+    Return the non-dominated records for named minimization/maximization goals.
+
+    ``objectives`` maps record keys to ``"min"`` or ``"max"``. Equal objective
+    vectors are retained, and the input order is preserved.
+    """
+
+    records = list(records)
+    if not objectives:
+        raise TSValueError("At least one Pareto objective is required.")
+    if not records:
+        return []
+
+    columns = []
+    for name, direction in objectives.items():
+        if direction not in {"min", "max"}:
+            raise TSValueError(
+                f"Objective {name!r} direction must be 'min' or 'max', got {direction!r}."
+            )
+        try:
+            column = np.asarray([record[name] for record in records], dtype=float)
+        except (KeyError, TypeError, ValueError) as exc:
+            raise TSValueError(f"Objective {name!r} is missing or non-numeric.") from exc
+        if not np.all(np.isfinite(column)):
+            raise TSValueError(f"Objective {name!r} must contain only finite values.")
+        columns.append(column if direction == "min" else -column)
+
+    points = np.column_stack(columns)
+    unique_points, inverse = np.unique(points, axis=0, return_inverse=True)
+    frontier = []
+    for index, point in enumerate(unique_points):
+        if not frontier:
+            frontier.append(index)
+            continue
+        current = unique_points[frontier]
+        dominated = np.any(
+            np.all(current <= point, axis=1) & np.any(current < point, axis=1)
+        )
+        if dominated:
+            continue
+        removes = np.all(point <= current, axis=1) & np.any(point < current, axis=1)
+        frontier = [
+            kept for kept, remove in zip(frontier, removes) if not remove
+        ]
+        frontier.append(index)
+
+    selected = set(frontier)
+    return [record for record, point_index in zip(records, inverse) if point_index in selected]
diff --git a/docs/api.rst b/docs/api.rst
index 905226f..c73c98c 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -36,6 +36,17 @@ Reactions and redox
 .. autofunction:: ThermoScreening.thermo.reactions.reduction_potential
 .. autofunction:: ThermoScreening.thermo.reactions.calibrate_reduction_reference
 
+Dataset redox screening
+-----------------------
+
+.. autoclass:: ThermoScreening.thermo.redox_screening.DeltaRedoxModel
+   :members: fit, predict, correction_coefficients
+.. autofunction:: ThermoScreening.thermo.redox_screening.canonical_structure_identity
+.. autofunction:: ThermoScreening.thermo.redox_screening.audit_redox_dataset
+.. autofunction:: ThermoScreening.thermo.redox_screening.balanced_group_folds
+.. autofunction:: ThermoScreening.thermo.redox_screening.grouped_delta_validation
+.. autofunction:: ThermoScreening.thermo.redox_screening.pareto_front
+
 Acid dissociation (pKa)
 ------------------------
 
diff --git a/docs/index.rst b/docs/index.rst
index 2b59452..d6ca86c 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -73,8 +73,8 @@ Features
       :link: usage
       :link-type: doc
 
-      Reaction free energies and one-electron reduction potentials
-      (vs. SHE or a calibrated reference) from any two ``Thermo`` objects.
+      Reaction free energies, reference-calibrated reduction potentials, and
+      validated DFTB-to-DFT correction models for screening series.
 
    .. grid-item-card:: :octicon:`flame` Kinetics
       :link: usage
@@ -105,6 +105,7 @@ Features
 
    installation
    usage
+   redox_screening
    configuration
    api
 
diff --git a/docs/redox_screening.rst b/docs/redox_screening.rst
new file mode 100644
index 0000000..a000fb2
--- /dev/null
+++ b/docs/redox_screening.rst
@@ -0,0 +1,205 @@
+Calibrated redox screening
+==========================
+
+A single reference compound removes a constant offset from computed reduction
+potentials. Across a chemically diverse series, the remaining error can still
+depend on the approximate potential, functional groups, substitution positions,
+and interactions between substituents. ``DeltaRedoxModel`` fits that residual:
+
+.. math::
+
+   \Delta E = E_{\mathrm{reference}} - E_{\mathrm{approximate}}
+
+The model is a regularized linear regression. Its coefficients remain available
+through ``correction_coefficients``. It can include integer functional-group
+counts and selected pair terms without requiring a machine-learning dependency.
+
+Data audit
+----------
+
+Canonicalize structures before splitting or fitting. The identity function uses
+canonical isomeric SMILES and a fixed-hydrogen InChIKey, so stereoisomers,
+protonation states, and fixed-hydrogen tautomers remain distinct.
+
+.. code-block:: python
+
+   from ThermoScreening.thermo import audit_redox_dataset
+
+   audit = audit_redox_dataset(
+       smiles,
+       {"E1_DFT": dft_e1, "E1_DFTB": dftb_e1},
+       tolerance=0.02,  # volts
+   )
+
+   print(audit.unique_structures)
+   print(audit.duplicate_groups)
+   print(audit.conflicts)
+
+Every duplicate is reported. ``conflicts`` identifies structures whose values
+differ by more than the tolerance. Inspect or exclude conflicting calculations;
+do not silently average them. Keep the original calculation provenance,
+convergence state, geometry, charge, spin, solvent, and method alongside each
+record.
+
+Grouped validation
+------------------
+
+Random row splits can place symmetry-equivalent structures or closely related
+positional isomers in both training and validation. Define a group identifier
+that represents the generalization claim, then assign complete groups to folds.
+For a substituted scaffold, a sorted substituent-composition tuple is a useful
+default because all positional isomers remain together.
+
+.. code-block:: python
+
+   from ThermoScreening.thermo import balanced_group_folds, DeltaRedoxModel
+
+   # One mapping per molecule. Include position labels when their effects matter.
+   descriptors = [
+       {"nitro": 1, "nitro@beta": 1},
+       {"hydroxy": 2, "hydroxy@alpha": 1, "hydroxy@beta": 1},
+       {"amino": 1, "amino@alpha": 1},
+   ]
+   composition_groups = [
+       ("nitro",),
+       ("hydroxy", "hydroxy"),
+       ("amino",),
+   ]
+   folds = balanced_group_folds(composition_groups, n_folds=3)
+
+   model = DeltaRedoxModel.fit(
+       dftb_e1,
+       dft_e1,
+       descriptors,
+       ridge=1.0e-4,
+       validation_groups=folds,
+   )
+
+   print(model.validation.mae)
+   print(model.validation.rmse)
+   print(model.correction_coefficients)
+
+Fit separate models for the first and second reductions. They describe different
+charge-state transitions and generally have different residual errors. Select
+the ridge strength and descriptor scheme without consulting the final
+experimental test set.
+
+Prediction and domain checks
+----------------------------
+
+.. code-block:: python
+
+   prediction = model.predict(candidate_dftb_e1, candidate_descriptors)
+
+   corrected = prediction.corrected_potential
+   uncertainty = prediction.uncertainty
+   extrapolated = prediction.extrapolated
+   reasons = prediction.unknown_features
+
+The uncertainty combines the regression residual with leverage and, when
+grouped validation was supplied, uses its RMSE as a lower bound. It is a model
+diagnostic rather than a calibrated probabilistic confidence interval. An
+unknown functional group, unsupported pair interaction, or potential outside
+the training range sets ``extrapolated`` and leaves the uncertainty as ``NaN``.
+Such candidates require higher-level calculations instead of automatic ranking.
+
+Pair interactions are optional:
+
+.. code-block:: python
+
+   interaction_model = DeltaRedoxModel.fit(
+       dftb_e1,
+       dft_e1,
+       descriptors,
+       include_interactions=True,
+       min_interaction_count=10,
+       ridge=1.0e-3,
+       validation_groups=folds,
+   )
+
+Use pair terms only when they improve grouped validation. A minimum support
+avoids assigning unconstrained coefficients to rare combinations.
+
+Candidate selection
+-------------------
+
+Do not replace electrochemical objectives with potential divided by molecular
+mass. That ratio is not specific energy. Keep competing properties separate and
+select the non-dominated candidates:
+
+.. code-block:: python
+
+   from ThermoScreening.thermo import pareto_front
+
+   candidates = [
+       {"name": "A", "E1": -0.65, "mass": 240.2, "uncertainty": 0.04},
+       {"name": "B", "E1": -0.58, "mass": 278.3, "uncertainty": 0.03},
+   ]
+   selected = pareto_front(
+       candidates,
+       {"E1": "max", "mass": "min", "uncertainty": "min"},
+   )
+
+Add solubility, stability, reversibility, or a physically defined cell voltage
+when those data are available. Pareto selection exposes the trade-offs; it does
+not decide which trade-off is appropriate for an application.
+
+Final validation
+----------------
+
+The corrected DFTB+ result remains a screening estimate. Recalculate Pareto-front
+candidates, extrapolations, and any case where ``E2 >= E1`` with the higher-level
+method. Check multiple conformers, confirm that optimized structures are minima,
+and inspect spin state and charge localization for the radical anion and dianion.
+Keep solvent, temperature, reference electrode, and standard-state conventions
+consistent across computed and experimental values.
+
+Reserve experimental measurements as a final test set. Do not use the same
+measurements to select descriptors, tune regularization, calibrate the reference,
+and report final accuracy.
+
+Anthraquinone benchmark
+-----------------------
+
+The workflow was checked against the matched DFTB+/DFT table from the
+`Anthraquinone-screening repository
+<https://github.com/jolint1/Anthraquinone-screening>`_. Every row belonging to
+one of the 22 duplicated canonical structures was excluded, leaving 5,035
+singleton structures before filtering missing potentials. Five balanced folds
+kept all positional isomers of a substituent composition together. The models
+used a ridge value of ``1e-4`` and no pair interactions.
+
+.. list-table:: Grouped five-fold errors against DFT
+   :header-rows: 1
+   :widths: 22 13 13 13 13
+
+   * - Model
+     - E1 MAE
+     - E1 RMSE
+     - E2 MAE
+     - E2 RMSE
+   * - Uncorrected DFTB+
+     - 82.9 mV
+     - 117.7 mV
+     - 92.6 mV
+     - 135.2 mV
+   * - Global intercept and slope
+     - 58.7 mV
+     - 78.2 mV
+     - 64.4 mV
+     - 92.6 mV
+   * - Substituent counts
+     - 31.7 mV
+     - 46.5 mV
+     - 44.7 mV
+     - 71.3 mV
+   * - Substituent and alpha/beta position
+     - 29.2 mV
+     - 42.7 mV
+     - 35.4 mV
+     - 59.6 mV
+
+The E1 and E2 comparisons contain 5,022 and 5,013 structures, respectively.
+These numbers measure agreement with the available DFT calculations, not with
+experiment. Experimental potentials measured under consistent conditions must
+remain an independent final validation set.
diff --git a/tests/thermo/test_redox_screening.py b/tests/thermo/test_redox_screening.py
new file mode 100644
index 0000000..2aa7e0b
--- /dev/null
+++ b/tests/thermo/test_redox_screening.py
@@ -0,0 +1,305 @@
+import numpy as np
+import pytest
+
+from ThermoScreening.exceptions import TSValueError
+from ThermoScreening.thermo.redox_screening import (
+    DeltaRedoxModel,
+    audit_redox_dataset,
+    balanced_group_folds,
+    canonical_structure_identity,
+    grouped_delta_validation,
+    pareto_front,
+)
+
+
+def _matched_series():
+    approximate = np.array([-1.15, -0.82, -1.31, -0.73, -1.02, -0.91, -1.24, -0.67])
+    groups = [
+        {"amino": 1},
+        {"nitro": 1},
+        {"amino": 1, "hydroxy": 1},
+        {"nitro": 1, "hydroxy": 1},
+        {"hydroxy": 1},
+        {"amino": 1, "nitro": 1},
+        {"amino": 2},
+        {"nitro": 2},
+    ]
+    correction = np.array(
+        [
+            0.04 + 0.12 * value + 0.03 * row.get("amino", 0) - 0.05 * row.get("nitro", 0)
+            + 0.02 * row.get("hydroxy", 0)
+            for value, row in zip(approximate, groups)
+        ]
+    )
+    return approximate, approximate + correction, groups
+
+
+def test_canonical_structure_identity_matches_equivalent_smiles():
+    first = canonical_structure_identity("CCO")
+    second = canonical_structure_identity("OCC")
+
+    assert first.canonical_smiles == "CCO"
+    assert first == second
+
+
+def test_canonical_structure_identity_preserves_stereochemistry():
+    clockwise = canonical_structure_identity("F[C@H](Cl)Br")
+    anticlockwise = canonical_structure_identity("F[C@@H](Cl)Br")
+
+    assert clockwise.structure_id != anticlockwise.structure_id
+
+
+def test_canonical_structure_identity_distinguishes_fixed_h_tautomers():
+    pyridone = canonical_structure_identity("O=c1cc[nH]cc1")
+    hydroxypyridine = canonical_structure_identity("Oc1ccncc1")
+
+    assert pyridone.structure_id != hydroxypyridine.structure_id
+
+
+@pytest.mark.parametrize("smiles", ["", None, "not a smiles"])
+def test_canonical_structure_identity_rejects_invalid_smiles(smiles):
+    with pytest.raises(TSValueError, match="SMILES|Invalid"):
+        canonical_structure_identity(smiles)
+
+
+def test_canonical_structure_identity_rejects_missing_inchi_key(monkeypatch):
+    monkeypatch.setattr(
+        "ThermoScreening.thermo.redox_screening.inchi.InchiToInchiKey",
+        lambda value: "",
+    )
+
+    with pytest.raises(TSValueError, match="derive an InChIKey"):
+        canonical_structure_identity("CCO")
+
+
+def test_audit_redox_dataset_reports_duplicates_and_conflicts():
+    audit = audit_redox_dataset(
+        ["CCO", "OCC", "CCN", "NCC"],
+        {"E1": [-0.5, -0.49, -0.8, -0.72], "E2": [None, np.nan, -1.1, -1.1]},
+        tolerance=0.02,
+    )
+
+    assert audit.unique_structures == 2
+    assert audit.duplicate_groups == ((0, 1), (2, 3))
+    assert len(audit.conflicts) == 1
+    assert audit.conflicts[0].indices == (2, 3)
+    assert audit.conflicts[0].spreads[0][0] == "E1"
+    assert audit.conflicts[0].spreads[0][1] == pytest.approx(0.08)
+
+
+def test_audit_redox_dataset_validates_columns_and_values():
+    with pytest.raises(TSValueError, match="expected 1"):
+        audit_redox_dataset(["CCO"], {"E1": [0.1, 0.2]})
+    with pytest.raises(TSValueError, match="not numeric"):
+        audit_redox_dataset(["CCO", "OCC"], {"E1": [0.1, "bad"]})
+    with pytest.raises(TSValueError, match="non-negative"):
+        audit_redox_dataset(["CCO"], tolerance=-1.0)
+    with pytest.raises(TSValueError, match="infinite"):
+        audit_redox_dataset(["CCO"], {"E1": [np.inf]})
+
+
+def test_delta_model_recovers_a_transparent_linear_correction():
+    approximate, reference, groups = _matched_series()
+
+    model = DeltaRedoxModel.fit(approximate, reference, groups, ridge=0.0)
+    prediction = model.predict(approximate, groups)
+
+    assert prediction.corrected_potential == pytest.approx(reference, abs=1.0e-10)
+    assert not prediction.extrapolated.any()
+    assert np.all(np.isfinite(prediction.uncertainty))
+    assert model.correction_coefficients["approximate_potential"] == pytest.approx(0.12)
+    assert model.correction_coefficients["group:amino"] == pytest.approx(0.03)
+    assert model.correction_coefficients["group:nitro"] == pytest.approx(-0.05)
+
+
+def test_delta_model_accepts_repeated_group_names_as_counts():
+    approximate, reference, groups = _matched_series()
+    names = [[name for name, count in row.items() for _ in range(int(count))] for row in groups]
+
+    model = DeltaRedoxModel.fit(approximate, reference, names, ridge=0.0)
+
+    assert model.predict(approximate, groups).corrected_potential == pytest.approx(reference)
+
+
+def test_delta_model_fits_supported_pair_interactions():
+    approximate = np.array([-1.2, -0.7, -1.0, -0.8, -1.3, -0.6, -1.1, -0.9])
+    groups = [
+        {},
+        {"a": 1},
+        {"b": 1},
+        {"a": 1, "b": 1},
+        {"a": 2},
+        {"b": 2},
+        {"a": 2, "b": 1},
+        {"a": 1, "b": 2},
+    ]
+    correction = []
+    for value, row in zip(approximate, groups):
+        a_count = row.get("a", 0)
+        b_count = row.get("b", 0)
+        correction.append(
+            0.02
+            + 0.05 * value
+            + 0.01 * a_count
+            - 0.03 * b_count
+            + 0.08 * a_count * b_count
+            + 0.04 * a_count * (a_count - 1) / 2
+            - 0.02 * b_count * (b_count - 1) / 2
+        )
+    reference = approximate + correction
+
+    model = DeltaRedoxModel.fit(
+        approximate,
+        reference,
+        groups,
+        include_interactions=True,
+        min_interaction_count=1,
+        ridge=0.0,
+    )
+
+    assert model.predict(approximate, groups).corrected_potential == pytest.approx(
+        reference, abs=1.0e-10
+    )
+    assert "pair:a|b" in model.correction_coefficients
+
+
+def test_delta_model_flags_unknown_groups_pairs_and_potential_range():
+    approximate = [-1.1, -0.9, -0.7, -0.8]
+    reference = [-1.0, -0.8, -0.6, -0.7]
+    groups = [{"a": 1}, {"b": 1}, {"c": 1}, {"a": 1, "b": 1}]
+    model = DeltaRedoxModel.fit(
+        approximate,
+        reference,
+        groups,
+        include_interactions=True,
+        min_interaction_count=1,
+    )
+
+    prediction = model.predict([-0.8, -1.3], [{"a": 1, "c": 1}, {"unknown": 1}])
+
+    assert prediction.extrapolated.tolist() == [True, True]
+    assert "pair:a|c" in prediction.unknown_features[0]
+    assert set(prediction.unknown_features[1]) == {
+        "approximate_potential_range",
+        "group:unknown",
+    }
+    assert np.isnan(prediction.uncertainty).all()
+
+
+def test_grouped_delta_validation_keeps_labels_in_single_folds():
+    approximate, reference, groups = _matched_series()
+    labels = ["one", "one", "two", "two", "three", "three", "four", "four"]
+
+    validation = grouped_delta_validation(
+        approximate, reference, groups, labels, ridge=1.0e-6
+    )
+    fitted = DeltaRedoxModel.fit(
+        approximate,
+        reference,
+        groups,
+        validation_groups=labels,
+        ridge=1.0e-6,
+    )
+
+    assert validation.predictions.shape == approximate.shape
+    assert validation.folds == tuple(labels)
+    assert validation.rmse >= 0.0
+    assert fitted.validation is not None
+    assert fitted.validation.rmse == pytest.approx(validation.rmse)
+    assert np.all(np.isfinite(fitted.predict(approximate, groups).uncertainty))
+
+
+@pytest.mark.parametrize(
+    ("approximate", "reference", "groups", "message"),
+    [
+        ([1.0], [1.0], [{}], "At least two"),
+        ([1.0, np.nan], [1.0, 2.0], [{}, {}], "finite"),
+        ([1.0, 2.0], [1.0], [{}, {}], "same length"),
+        ([1.0, 2.0], [1.0, 2.0], [{}], "expected 2"),
+        ([1.0, 2.0], [1.0, 2.0], [{"nitro": -1}, {}], "non-negative"),
+        ([1.0, 2.0], [1.0, 2.0], [{"nitro": 0.5}, {}], "integer"),
+        ([[1.0, 2.0]], [1.0], [{}], "one-dimensional"),
+        ([1.0, 2.0], [1.0, 2.0], ["nitro", []], "mapping or an iterable"),
+        ([1.0, 2.0], [1.0, 2.0], [{1: 1}, {}], "non-empty strings"),
+        ([1.0, 2.0], [1.0, 2.0], [{"nitro": "many"}, {}], "not numeric"),
+    ],
+)
+def test_delta_model_rejects_invalid_training_data(
+    approximate, reference, groups, message
+):
+    with pytest.raises(TSValueError, match=message):
+        DeltaRedoxModel.fit(approximate, reference, groups)
+
+
+def test_delta_model_rejects_invalid_regularization_options():
+    with pytest.raises(TSValueError, match="ridge"):
+        DeltaRedoxModel.fit([1, 2], [1, 2], [{}, {}], ridge=-1)
+    with pytest.raises(TSValueError, match="min_interaction_count"):
+        DeltaRedoxModel.fit(
+            [1, 2],
+            [1, 2],
+            [{}, {}],
+            include_interactions=True,
+            min_interaction_count=0,
+        )
+
+
+def test_grouped_delta_validation_rejects_invalid_folds():
+    with pytest.raises(TSValueError, match="at least two distinct"):
+        grouped_delta_validation([1, 2, 3], [1, 2, 3], [{}, {}, {}], ["same"] * 3)
+    with pytest.raises(TSValueError, match="fewer than two training"):
+        grouped_delta_validation([1, 2, 3], [1, 2, 3], [{}, {}, {}], ["a", "a", "b"])
+    with pytest.raises(TSValueError, match="same length"):
+        grouped_delta_validation([1, 2], [1], [{}, {}], ["a", "b"])
+    with pytest.raises(TSValueError, match="expected 2"):
+        grouped_delta_validation([1, 2], [1, 2], [{}, {}], ["a"])
+    with pytest.raises(TSValueError, match="hashable"):
+        grouped_delta_validation([1, 2], [1, 2], [{}, {}], [["a"], ["b"]])
+
+
+def test_balanced_group_folds_never_splits_a_group():
+    groups = ["large"] * 5 + ["medium"] * 3 + ["small-a"] * 2 + ["small-b"]
+
+    folds = balanced_group_folds(groups, n_folds=3)
+
+    assert len(folds) == len(groups)
+    assert all(len({fold for fold, value in zip(folds, groups) if value == group}) == 1 for group in set(groups))
+    loads = [folds.count(index) for index in range(3)]
+    assert max(loads) - min(loads) <= 2
+
+
+def test_balanced_group_folds_validates_inputs():
+    with pytest.raises(TSValueError, match="at least 2"):
+        balanced_group_folds(["a", "b"], n_folds=1)
+    with pytest.raises(TSValueError, match="only 2 distinct"):
+        balanced_group_folds(["a", "b"], n_folds=3)
+    with pytest.raises(TSValueError, match="hashable"):
+        balanced_group_folds([["a"], ["b"]], n_folds=2)
+
+
+def test_pareto_front_preserves_tradeoffs_and_equal_points():
+    records = [
+        {"name": "a", "potential": -0.8, "mass": 200, "uncertainty": 0.05},
+        {"name": "b", "potential": -0.7, "mass": 210, "uncertainty": 0.04},
+        {"name": "c", "potential": -0.9, "mass": 220, "uncertainty": 0.10},
+        {"name": "d", "potential": -0.8, "mass": 200, "uncertainty": 0.05},
+    ]
+
+    front = pareto_front(
+        records,
+        {"potential": "max", "mass": "min", "uncertainty": "min"},
+    )
+
+    assert [record["name"] for record in front] == ["a", "b", "d"]
+
+
+def test_pareto_front_handles_empty_input_and_validates_objectives():
+    assert pareto_front([], {"potential": "max"}) == []
+    with pytest.raises(TSValueError, match="At least one"):
+        pareto_front([{"potential": 1.0}], {})
+    with pytest.raises(TSValueError, match="'min' or 'max'"):
+        pareto_front([{"potential": 1.0}], {"potential": "up"})
+    with pytest.raises(TSValueError, match="finite"):
+        pareto_front([{"potential": np.nan}], {"potential": "max"})
+    with pytest.raises(TSValueError, match="missing or non-numeric"):
+        pareto_front([{"mass": 10.0}], {"potential": "max"})