diff --git a/packages/peptide-property-ng/.gitignore b/packages/peptide-property-ng/.gitignore
new file mode 100644
index 000000000..015345982
--- /dev/null
+++ b/packages/peptide-property-ng/.gitignore
@@ -0,0 +1,11 @@
+# training / evaluation outputs
+runs/
+
+# python build / cache
+__pycache__/
+*.py[cod]
+*.egg-info/
+build/
+dist/
+.pytest_cache/
+.venv/
diff --git a/packages/peptide-property-ng/README.md b/packages/peptide-property-ng/README.md
new file mode 100644
index 000000000..162a6b118
--- /dev/null
+++ b/packages/peptide-property-ng/README.md
@@ -0,0 +1,40 @@
+# peptide-property-ng
+
+A next-generation, clean-slate neural network for **peptide property prediction**
+on timsTOF proteomics data — fragment intensity, ion mobility / CCS, retention
+time and precursor charge from **one shared encoder**.
+
+It is a research-track sibling of `imspy-predictors`, not a replacement for the
+production model. See `docs/nn-architecture-exploration/` in the
+`claudius-proteomics` project for the design rationale.
+
+## What is different from the production `UnifiedPeptideModel`
+
+| | production `imspy_predictors` | this model |
+|---|---|---|
+| Modifications | opaque `[UNIMOD:N]` tokens (no chemistry) | **hybrid**: per-mod token **+** atomic-composition feature |
+| Unseen mods | embedding stays ~random | generalize via chemistry; open-vocab |
+| Conditioning | dormant instrument embedding, never used | instrument / acq-mode in the encoder; charge / m-z / CE at the heads (keeps the charge head leak-free) |
+| Intensity output | fixed 174-vector → 30-residue cap | **fragment-indexed** → no length cap |
+| Backbone | bespoke transformer | built on **Depthcharge** primitives |
+
+## Layout
+
+```
+src/peptide_property_ng/
+  modifications/   UNIMOD atomic-composition table (from sagepy)
+  model/           hybrid embedding, Depthcharge-based encoder, task heads
+  data/            Sage-parquet dataset, fragment targets, splits, collate
+  train/           multi-task training
+  eval/            per-task + OOD + unseen-modification evaluation
+```
+
+## Quick start
+
+```bash
+# 1. build the modification-composition table (one-off)
+python -m peptide_property_ng.modifications.build_table
+
+# 2. train the first-prototype model on processed Sage results
+python -m peptide_property_ng.train.train --datasets-glob '/scratch/claudius-proteomics/*'
+```
diff --git a/packages/peptide-property-ng/pyproject.toml b/packages/peptide-property-ng/pyproject.toml
new file mode 100644
index 000000000..2c77fac91
--- /dev/null
+++ b/packages/peptide-property-ng/pyproject.toml
@@ -0,0 +1,38 @@
+[build-system]
+requires = ["hatchling"]
+build-backend = "hatchling.build"
+
+[project]
+name = "peptide-property-ng"
+version = "0.0.1"
+description = "Next-generation unified neural predictor for peptide properties (fragment intensity, ion mobility/CCS, retention time, charge) on timsTOF proteomics data."
+readme = "README.md"
+requires-python = ">=3.11"
+license = "MIT"
+
+# Core, pip-installable dependencies.
+dependencies = [
+    "torch>=2.0",
+    "depthcharge-ms>=0.4.9",
+    "numpy>=1.24",
+    "pyarrow>=15",
+    "scipy>=1.10",
+]
+
+# NOTE: two further dependencies are Rust-built monorepo artifacts, not on PyPI,
+# and must be provided by the environment:
+#   - imspy_connector  -> the Rust-backed `ProformaTokenizer`
+#   - sagepy / sagepy_connector -> the UNIMOD atomic-composition tables
+# During development we run against the claudius-proteomics `.venv`, which has both.
+
+[project.optional-dependencies]
+dev = ["pytest>=7"]
+
+[tool.hatch.build.targets.wheel]
+packages = ["src/peptide_property_ng"]
+
+[tool.hatch.build.targets.wheel.force-include]
+"src/peptide_property_ng/modifications/data" = "peptide_property_ng/modifications/data"
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/__init__.py b/packages/peptide-property-ng/src/peptide_property_ng/__init__.py
new file mode 100644
index 000000000..c23002225
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/__init__.py
@@ -0,0 +1,8 @@
+"""peptide-property-ng — unified neural peptide property predictor.
+
+A clean-slate, Depthcharge-based model that predicts fragment intensity, ion
+mobility / CCS, retention time and precursor charge from one shared encoder,
+with a hybrid (token + atomic-composition) modification representation.
+"""
+
+__version__ = "0.0.1"
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/__init__.py b/packages/peptide-property-ng/src/peptide_property_ng/data/__init__.py
new file mode 100644
index 000000000..6e5ae62bb
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/__init__.py
@@ -0,0 +1 @@
+"""Data — Sage-parquet dataset, fragment-intensity targets, splits, collation."""
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/ccs_pretrain.py b/packages/peptide-property-ng/src/peptide_property_ng/data/ccs_pretrain.py
new file mode 100644
index 000000000..329914934
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/ccs_pretrain.py
@@ -0,0 +1,137 @@
+"""ionmob CCS datasets -> ion-mobility-pretraining examples.
+
+Uses the deduplicated, UNIMOD-annotated ionmob CCS parquets
+(`*_unique_unimod*.parquet`, Meier / Tenzer / Chang / ...): schema
+``mz, charge, sequence-tokenized, ccs, name``.
+
+Target-unit note: these provide **CCS** (Angstrom^2); the Sage fine-tuning data
+provides **inverse ion mobility** (1/K0). They are deliberately *not* converted
+into each other — that physics conversion is error-prone, and it is unnecessary:
+pretraining is single-task per stage, so the CCS head simply warms on CCS here
+and the unit/scale shift to 1/K0 is absorbed during the campaign fine-tune —
+exactly as the RT head pretrains on Chronologer ``HI`` and fine-tunes on Sage
+``aligned_rt``.
+"""
+from __future__ import annotations
+
+import glob
+import re
+
+import numpy as np
+import pyarrow.parquet as pq
+
+from peptide_property_ng.data.sage_dataset import SagePropertyDataset
+
+DEFAULT_CCS_GLOB = (
+    "/home/administrator/Documents/promotion/rust/rustims-submission/"
+    "re-scoring/models/unimod/**/*_unique_unimod*.parquet"
+)
+_UNIMOD_RE = re.compile(r"\[UNIMOD:\d+\]")
+
+# CCS (Angstrom^2) is min-max normalised to ~[0,1] with fixed bounds so the CCS
+# head sees the *same target scale* in pretraining (CCS) and fine-tuning (Sage
+# 1/K0, normalised the same way in sage_dataset.py) — this avoids the negative
+# transfer of feeding a ~400-scale target into a [0,1]-scale head.
+CCS_MIN, CCS_MAX = 200.0, 1200.0
+
+
+def normalize_ccs(ccs: float) -> float:
+    """Map a CCS value (Angstrom^2) to [0,1], clipped at the fixed bounds."""
+    return float(np.clip((float(ccs) - CCS_MIN) / (CCS_MAX - CCS_MIN), 0.0, 1.0))
+
+
+def _tokens_to_proforma(tokens: list[str]) -> str | None:
+    """ionmob ``sequence-tokenized`` list -> a ProForma string, or ``None`` to skip.
+
+    A ``<START>`` / ``<END>`` marker carrying a terminal modification (e.g.
+    ``<START>[UNIMOD:1]``) flags a terminal mod -> skip, consistent with the
+    other loaders.
+    """
+    if len(tokens) < 5 or tokens[0] != "<START>" or tokens[-1] != "<END>":
+        return None
+    return "".join(tokens[1:-1])
+
+
+def prepare_ccs_examples(
+    data_glob: str = DEFAULT_CCS_GLOB,
+    *,
+    cap: int | None = None,
+    instrument: int = 0,
+    acq_mode: int = 0,
+    max_len: int = 64,
+    seed: int = 0,
+) -> list[dict]:
+    """Load the ionmob CCS parquets into ion-mobility-pretraining example dicts."""
+    from imspy_predictors.utilities.tokenizers import ProformaTokenizer
+
+    files = sorted(set(glob.glob(data_glob, recursive=True)))
+    if not files:
+        return []
+    seqs: list[list[str]] = []
+    charge: list[int] = []
+    mz: list[float] = []
+    ccs: list[float] = []
+    for path in files:
+        d = pq.read_table(path, columns=["mz", "charge", "sequence-tokenized", "ccs"]).to_pydict()
+        seqs += d["sequence-tokenized"]
+        charge += d["charge"]
+        mz += d["mz"]
+        ccs += d["ccs"]
+
+    idx = np.arange(len(seqs))
+    if cap is not None and len(idx) > cap:
+        idx = np.sort(np.random.RandomState(seed).choice(len(idx), cap, replace=False))
+
+    tok = ProformaTokenizer.with_defaults()
+    proformas: list[str] = []
+    src: list[int] = []
+    for i in idx:
+        pf = _tokens_to_proforma(seqs[i])
+        if pf is not None:
+            proformas.append(pf)
+            src.append(int(i))
+
+    examples: list[dict] = []
+    chunk = 20000  # tokenise in chunks (batch tokenisation right-pads to the chunk max)
+    for start in range(0, len(proformas), chunk):
+        block = proformas[start : start + chunk]
+        block_src = src[start : start + chunk]
+        enc = tok(block)
+        ids, attn = enc["input_ids"], enc["attention_mask"]
+        for j, pf in enumerate(block):
+            i = block_src[j]
+            n_real = int(sum(attn[j]))
+            residue_ids = ids[j][1 : n_real - 1]
+            length = len(residue_ids)
+            stripped = _UNIMOD_RE.sub("", pf)
+            if length != len(stripped) or not 3 <= length <= max_len:
+                continue
+            ccs_val, mz_val, z = float(ccs[i]), float(mz[i]), int(charge[i])
+            if not (np.isfinite(ccs_val) and ccs_val > 0
+                    and np.isfinite(mz_val) and mz_val > 0 and z >= 1):
+                continue  # skip malformed rows rather than emit a NaN target
+            examples.append(
+                {
+                    "accession": "ionmob",
+                    "stripped": stripped,
+                    "modseq": pf,
+                    "tokens": np.asarray(residue_ids, dtype=np.int64),
+                    "charge": z,
+                    "precursor_mz": mz_val,
+                    "collision_energy": 0.0,
+                    "instrument": int(instrument),
+                    "acq_mode": int(acq_mode),
+                    # masked placeholder — CCS pretraining runs only the CCS task
+                    "intensity_target": np.full((length - 1, 6), -1.0, dtype=np.float32),
+                    "ccs_target": normalize_ccs(ccs_val),
+                    "ccs_valid": True,
+                    "rt_target": float("nan"),
+                    "rt_valid": False,
+                }
+            )
+    return examples
+
+
+def build_ccs_dataset(data_glob: str = DEFAULT_CCS_GLOB, *, cap: int | None = None, **kwargs):
+    """Prepared ionmob CCS examples wrapped as a Dataset."""
+    return SagePropertyDataset(prepare_ccs_examples(data_glob, cap=cap, **kwargs))
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/chronologer_rt.py b/packages/peptide-property-ng/src/peptide_property_ng/data/chronologer_rt.py
new file mode 100644
index 000000000..038ed8fe2
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/chronologer_rt.py
@@ -0,0 +1,106 @@
+"""Chronologer retention-time database -> RT-pretraining examples.
+
+The Chronologer DB (``Chronologer_DB_220308.gz``, ~2.64 M peptide-RT records
+harmonised across 11 community datasets) is the RT pretraining corpus. Its
+``HI`` (hydrophobicity index) column is the cross-dataset-harmonised RT target
+-- so it sidesteps the per-run RT-scale problem -- and is min-max normalised to
+``[0, 1]`` to match the model's RT-head output range and the Sage ``aligned_rt``
+fine-tuning target.
+
+Peptides use Sage-style ``[+mass]`` delta notation; they are converted to UNIMOD
+form with the same ``sagepy_rescore`` converter the Sage loader uses.
+"""
+from __future__ import annotations
+
+import re
+from pathlib import Path
+
+import numpy as np
+
+from peptide_property_ng.data.sage_dataset import SagePropertyDataset
+
+DEFAULT_CHRONOLOGER_DB = Path(
+    "/home/administrator/Documents/promotion/chronologer/data/Chronologer_DB_220308.gz"
+)
+# HI range across the full DB — fixed so normalisation is reproducible for any cap.
+HI_MIN, HI_MAX = -1.01, 31.07
+_UNIMOD_RE = re.compile(r"\[UNIMOD:\d+\]")
+
+
+def _normalise_hi(hi: float) -> float:
+    """HI -> ~[0,1], clipped (a few records sit just outside the fixed bounds).
+
+    A non-finite HI propagates as NaN through ``np.clip`` and is filtered by the
+    caller.
+    """
+    return float(np.clip((float(hi) - HI_MIN) / (HI_MAX - HI_MIN), 0.0, 1.0))
+
+
+def prepare_chronologer_examples(
+    db_path: str | Path = DEFAULT_CHRONOLOGER_DB,
+    *,
+    cap: int | None = None,
+    instrument: int = 0,
+    acq_mode: int = 0,
+    max_len: int = 64,
+    seed: int = 0,
+) -> list[dict]:
+    """Load the Chronologer DB into RT-pretraining example dicts."""
+    import pandas as pd
+    from imspy_predictors.utilities.tokenizers import ProformaTokenizer
+
+    from peptide_property_ng.data.mass_to_unimod import parse_delta_mass_peptide
+
+    df = pd.read_csv(db_path, sep="\t", compression="gzip", usecols=["PeptideModSeq", "HI"])
+    if cap is not None and len(df) > cap:
+        df = df.sample(cap, random_state=seed)
+    peptides = df["PeptideModSeq"].tolist()
+    hi = df["HI"].to_numpy(dtype=np.float32)
+
+    tok = ProformaTokenizer.with_defaults()
+    parsed = [parse_delta_mass_peptide(p) for p in peptides]  # [+mass] -> UNIMOD
+
+    examples: list[dict] = []
+    chunk = 20000  # tokenise in chunks (batch tokenisation right-pads to the chunk max)
+    for start in range(0, len(parsed), chunk):
+        block = parsed[start : start + chunk]
+        enc = tok([m for _, m in block])
+        ids, attn = enc["input_ids"], enc["attention_mask"]
+        for j, (stripped, modseq) in enumerate(block):
+            if "[+" in modseq or "[-" in modseq:
+                continue  # an unconverted delta mass — skip
+            n_real = int(sum(attn[j]))
+            residue_ids = ids[j][1 : n_real - 1]  # strip [CLS]/[SEP] + right-padding
+            length = len(residue_ids)
+            if length != len(stripped) or not 3 <= length <= max_len:
+                continue
+            rt = _normalise_hi(hi[start + j])
+            if not np.isfinite(rt):
+                continue  # non-finite HI -> useless as an RT example
+            examples.append(
+                {
+                    "accession": "Chronologer",
+                    "stripped": stripped,
+                    "modseq": modseq,
+                    "tokens": np.asarray(residue_ids, dtype=np.int64),
+                    "charge": 2,  # dummy — RT is charge-independent; the RT head ignores charge
+                    "precursor_mz": 0.0,
+                    "collision_energy": 0.0,
+                    "instrument": int(instrument),
+                    "acq_mode": int(acq_mode),
+                    # masked placeholder — RT pretraining runs only the RT task
+                    "intensity_target": np.full((length - 1, 6), -1.0, dtype=np.float32),
+                    "ccs_target": float("nan"),
+                    "ccs_valid": False,
+                    "rt_target": rt,
+                    "rt_valid": True,
+                }
+            )
+    return examples
+
+
+def build_chronologer_dataset(
+    db_path: str | Path = DEFAULT_CHRONOLOGER_DB, *, cap: int | None = None, **kwargs
+) -> SagePropertyDataset:
+    """Prepared Chronologer RT examples wrapped as a Dataset."""
+    return SagePropertyDataset(prepare_chronologer_examples(db_path, cap=cap, **kwargs))
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/collate.py b/packages/peptide-property-ng/src/peptide_property_ng/data/collate.py
new file mode 100644
index 000000000..2dc3fd511
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/collate.py
@@ -0,0 +1,61 @@
+"""Variable-length collation of prepared examples into model batches."""
+from __future__ import annotations
+
+import numpy as np
+import torch
+
+from peptide_property_ng.data.fragment_targets import N_ION_CHANNELS
+
+
+def collate_examples(
+    examples: list[dict],
+    *,
+    pad_token_id: int = 61,
+    max_charge: int = 8,
+) -> dict[str, torch.Tensor]:
+    """Collate a list of example dicts into a padded batch.
+
+    Token sequences pad to the batch-max length with ``pad_token_id``; intensity
+    targets pad to ``(batch_max_len - 1, n_ion_channels)`` with ``-1`` (the
+    masked-out marker the spectral-angle loss ignores).
+    """
+    batch_size = len(examples)
+    max_len = max(len(e["tokens"]) for e in examples)
+
+    tokens = np.full((batch_size, max_len), pad_token_id, dtype=np.int64)
+    intensity = np.full((batch_size, max_len - 1, N_ION_CHANNELS), -1.0, dtype=np.float32)
+    for i, e in enumerate(examples):
+        length = len(e["tokens"])
+        tokens[i, :length] = e["tokens"]
+        tgt = e["intensity_target"]
+        intensity[i, : tgt.shape[0], :] = tgt
+
+    def _col(key: str, dtype) -> torch.Tensor:
+        return torch.tensor([e[key] for e in examples], dtype=dtype)
+
+    charge = _col("charge", torch.long)
+    return {
+        "tokens": torch.from_numpy(tokens),
+        "charge": charge,
+        "precursor_mz": _col("precursor_mz", torch.float32),
+        "collision_energy": _col("collision_energy", torch.float32),
+        "instrument": _col("instrument", torch.long),
+        "acq_mode": _col("acq_mode", torch.long),
+        # targets
+        "intensity_target": torch.from_numpy(intensity),
+        "ccs_target": _col("ccs_target", torch.float32),
+        "ccs_valid": _col("ccs_valid", torch.bool),
+        "rt_target": _col("rt_target", torch.float32),
+        "rt_valid": _col("rt_valid", torch.bool),
+        # charge head target: class j == charge j+1
+        "charge_target": charge.clamp(1, max_charge) - 1,
+    }
+
+
+def make_collate_fn(pad_token_id: int = 61, max_charge: int = 8):
+    """Return a ``collate_fn`` bound to the given padding / charge settings."""
+
+    def _fn(examples: list[dict]) -> dict[str, torch.Tensor]:
+        return collate_examples(examples, pad_token_id=pad_token_id, max_charge=max_charge)
+
+    return _fn
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/fragment_targets.py b/packages/peptide-property-ng/src/peptide_property_ng/data/fragment_targets.py
new file mode 100644
index 000000000..dc4f58177
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/fragment_targets.py
@@ -0,0 +1,67 @@
+"""Fragment-intensity targets.
+
+The MS2 intensity target the model trains against is **cleavage-site-indexed**:
+``(L-1, 6)`` — one row per cleavage site, channels ``[y+1, y+2, y+3, b+1, b+2, b+3]``.
+
+Both data sources reach it through the *same single conversion*,
+``prosit174_to_sites``:
+  - Public Wilhelmlab HF datasets store intensities natively as the Prosit
+    174-vector.
+  - Sage search results are turned into that 174-vector by the proven
+    ``imspy_predictors`` encoder ``observed_fragments_to_intensity_target``
+    (see ``data/sage_dataset.py``) — the fragment→vector encoding is *not*
+    re-implemented here.
+
+So there is exactly one place where the ordinal→site remap happens, and it is
+this file.
+
+Convention (Prosit / imspy compatible): ``-1`` = impossible/masked channel,
+``0`` = possible but unobserved, ``0..1`` = observed (base-peak-normalised).
+"""
+from __future__ import annotations
+
+import numpy as np
+
+# Channel order — verified against imspy_predictors/intensity/utility.py:211
+# ("Ion types in order: y+1, y+2, y+3, b+1, b+2, b+3").
+ION_CHANNELS: tuple[str, ...] = ("y+1", "y+2", "y+3", "b+1", "b+2", "b+3")
+N_ION_CHANNELS = len(ION_CHANNELS)
+
+# The Prosit 174-vector covers fragment ordinals 1..29 (peptides up to 30 aa).
+PROSIT_MAX_ORDINAL = 29
+PROSIT_VECTOR_LEN = PROSIT_MAX_ORDINAL * N_ION_CHANNELS  # 174
+
+
+def prosit174_to_sites(intensities_174, peptide_len: int) -> np.ndarray:
+    """Convert a Prosit 174-element intensity vector to our site-indexed target.
+
+    The Prosit vector is **ordinal-indexed**: reshaped to (29, 6), row ``r`` holds
+    the ions of fragment ordinal ``r+1`` — ``y_{r+1}`` in columns 0:3, ``b_{r+1}``
+    in columns 3:6.
+
+    Our ``(L-1, 6)`` target is **cleavage-site-indexed**: site ``i`` holds the two
+    ions produced by cleaving bond ``i`` — the b ion ``b_{i+1}`` and its
+    *complementary* y ion ``y_{L-1-i}``. The b and y of one cleavage have
+    different ordinals (``b_k`` pairs with ``y_{L-k}``), so this is a genuine
+    remap, not a reshape.
+
+    Mapping (``L`` = peptide_len, ``n_sites`` = L-1):
+        site i, b channels (3:6)  <-  ordinal i+1    ->  Prosit row i
+        site i, y channels (0:3)  <-  ordinal L-1-i  ->  Prosit row (L-1-i)-1 = L-2-i
+
+    The ``-1`` masked-channel markers in the Prosit vector carry through.
+    """
+    arr = np.asarray(intensities_174, dtype=np.float32).reshape(-1)
+    if arr.size != PROSIT_VECTOR_LEN:
+        raise ValueError(f"expected {PROSIT_VECTOR_LEN} intensity values, got {arr.size}")
+    prosit = arr.reshape(PROSIT_MAX_ORDINAL, N_ION_CHANNELS)  # row r -> fragment ordinal r+1
+
+    n_sites = peptide_len - 1
+    if not 1 <= n_sites <= PROSIT_MAX_ORDINAL:
+        raise ValueError(f"peptide_len {peptide_len} outside the 174-vector range [2, 30]")
+
+    target = np.full((n_sites, N_ION_CHANNELS), -1.0, dtype=np.float32)
+    for i in range(n_sites):
+        target[i, 3:6] = prosit[i, 3:6]                  # b_{i+1}
+        target[i, 0:3] = prosit[n_sites - 1 - i, 0:3]    # y_{L-1-i}  (Prosit row L-2-i)
+    return target
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/hf_corpus_dataset.py b/packages/peptide-property-ng/src/peptide_property_ng/data/hf_corpus_dataset.py
new file mode 100644
index 000000000..1cc18b852
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/hf_corpus_dataset.py
@@ -0,0 +1,188 @@
+"""HF-corpus -> multi-task training examples (validation of the published corpus).
+
+Mirrors ``sage_dataset.prepare_examples`` but sources from the aggregated HF
+corpus parquets (TO_HF_CORPUS.md): Tier-1 (per-precursor RT/IM/charge/CE) +
+Tier-3 (per-precursor matched b/y fragments). Produces the SAME example dicts
+``SagePropertyDataset`` wraps, so it drops into the existing collate + train loop.
+
+Differences vs the Sage-parquet path, by design:
+- RT: Tier-1 ships ``rt_seconds`` (raw apex), not Sage ``aligned_rt`` -> we
+  **per-accession min-max normalise** RT so the head sees a comparable scale.
+- CE: Tier-1 ships real per-precursor ``collision_energy_mean_v`` (volts) ->
+  ``CE = clip(volts/100, 0, 1)`` (the production convention), not a flat 0.26.
+- The split is already baked into the parquet (``split`` column).
+"""
+from __future__ import annotations
+
+from collections import defaultdict
+from pathlib import Path
+from types import SimpleNamespace
+
+import numpy as np
+import pyarrow.dataset as pds
+
+from peptide_property_ng.data.fragment_targets import prosit174_to_sites
+from peptide_property_ng.data.sage_dataset import (SagePropertyDataset,
+                                                   normalize_inverse_mobility)
+
+_PROTON = 1.007276
+
+_T1_COLS = ["accession", "raw_file", "precursor_id", "sage_psm_id",
+            "modified_sequence", "sequence", "charge", "mz", "rt_seconds",
+            "rt_aligned", "mobility", "collision_energy_mean_v"]
+
+
+def load_rt_aligned_lookup(path):
+    """{(accession, int sage_psm_id): aligned_rt} from a lookup parquet."""
+    import pyarrow.parquet as pq
+    t = pq.read_table(path, columns=["accession", "psm_id", "aligned_rt"])
+    accs = t["accession"].to_pylist()
+    pids = t["psm_id"].to_pylist()
+    rts = t["aligned_rt"].to_pylist()
+    return {(a, int(p)): r for a, p, r in zip(accs, pids, rts) if p is not None}
+_T3_COLS = ["accession", "raw_file", "precursor_id", "fragment_type",
+            "fragment_ordinal", "fragment_charge", "fragment_intensity"]
+
+
+def _key(acc, raw, pid):
+    return (acc, raw, int(pid))
+
+
+def prepare_examples_hf(corpus_dir, split, *, cap=4000, max_len=30, seed=0,
+                        accessions=None, rt_lookup=None):
+    """Build example dicts from ``corpus_dir/{tier1,tier3}/{split}.parquet``.
+
+    rt_lookup: optional {(accession, int sage_psm_id): aligned_rt} (or a path)
+    — when a precursor's aligned RT is found it is used (cross-run consistent,
+    ~[0,1]); otherwise we fall back to per-accession-normalised raw rt_seconds.
+    """
+    from imspy_predictors.intensity.predictors import observed_fragments_to_intensity_target
+    from imspy_predictors.utilities.tokenizers import ProformaTokenizer
+
+    from peptide_property_ng.data.mass_to_unimod import parse_delta_mass_peptide
+
+    if isinstance(rt_lookup, (str, Path)):
+        rt_lookup = load_rt_aligned_lookup(rt_lookup)
+
+    corpus_dir = Path(corpus_dir)
+    afilter = (pds.field("accession").isin(list(accessions))
+               if accessions else None)
+
+    # Tier-3 is ~391M rows -> NEVER materialise it. Pass 1: stream Tier-1
+    # (pushdown + batches), pick up to `cap` precursors/accession + props +
+    # per-accession RT range. Pass 2: stream Tier-3, keep only the selected
+    # precursors' fragments. Memory stays bounded by the selection (~50k prec).
+    t1set = pds.dataset(corpus_dir / "tier1" / f"{split}.parquet")
+    cols = [c for c in _T1_COLS if c in t1set.schema.names]   # rt_aligned only in v0.2+
+    sel: dict = {}                       # key -> props
+    per_acc_count: dict = defaultdict(int)
+    rt_lo, rt_hi = {}, {}
+    for rb in t1set.scanner(columns=cols, filter=afilter,
+                            batch_size=131072).to_batches():
+        d = rb.to_pydict()
+        for i in range(len(d["precursor_id"])):
+            acc = d["accession"][i]
+            rt = d["rt_seconds"][i]
+            if rt is not None and np.isfinite(rt):
+                rt_lo[acc] = rt if acc not in rt_lo else min(rt_lo[acc], rt)
+                rt_hi[acc] = rt if acc not in rt_hi else max(rt_hi[acc], rt)
+            if per_acc_count[acc] >= cap or not d["modified_sequence"][i]:
+                continue
+            key = _key(acc, d["raw_file"][i], d["precursor_id"][i])
+            sel[key] = {"acc": acc, "modseq_raw": d["modified_sequence"][i],
+                        "charge": d["charge"][i], "mz": d["mz"][i],
+                        "rt": rt, "mobility": d["mobility"][i],
+                        "ce_v": d["collision_energy_mean_v"][i],
+                        "sage_psm_id": d["sage_psm_id"][i],
+                        "rt_aligned": d["rt_aligned"][i] if "rt_aligned" in d else None}
+            per_acc_count[acc] += 1
+
+    by_prec: dict = defaultdict(lambda: ([], [], [], []))
+    t3set = pds.dataset(corpus_dir / "tier3" / f"{split}.parquet")
+    for rb in t3set.scanner(columns=_T3_COLS, filter=afilter,
+                            batch_size=262144).to_batches():
+        fd = rb.to_pydict()
+        for acc, raw, pid, ft, fo, fc, fi in zip(
+                fd["accession"], fd["raw_file"], fd["precursor_id"],
+                fd["fragment_type"], fd["fragment_ordinal"],
+                fd["fragment_charge"], fd["fragment_intensity"]):
+            key = _key(acc, raw, pid)
+            if key not in sel:
+                continue
+            cols = by_prec[key]
+            cols[0].append(ft); cols[1].append(int(fo))
+            cols[2].append(int(fc)); cols[3].append(float(fi))
+
+    tok = ProformaTokenizer.with_defaults()
+    examples: list[dict] = []
+    for key, p in sel.items():
+        acc = p["acc"]
+        cols = by_prec.get(key)
+        if not cols or not cols[0]:
+            continue  # no matched fragments -> no intensity target
+        modseq_raw = p["modseq_raw"]
+        stripped, modseq = parse_delta_mass_peptide(modseq_raw)
+        if "[+" in modseq or "[-" in modseq:
+            continue  # unconverted delta mass
+        enc = tok([modseq])
+        ids, mask = enc["input_ids"][0], enc["attention_mask"][0]
+        n_real = int(sum(mask))
+        residue_ids = ids[1:n_real - 1]
+        length = len(residue_ids)
+        if length < 3 or length > max_len:
+            continue
+
+        charge = int(p["charge"])
+        mz = float(p["mz"]) if p["mz"] is not None else 0.0
+        frag = SimpleNamespace(ion_types=cols[0], fragment_ordinals=cols[1],
+                               charges=cols[2], intensities=cols[3])
+        prosit174 = observed_fragments_to_intensity_target(stripped, charge, frag)
+        target = prosit174_to_sites(prosit174, length)
+
+        im = p["mobility"]
+        im = float(im) if im is not None else float("nan")
+        # RT target: prefer the native rt_aligned column (v0.2+), else an
+        # external aligned_rt lookup, else per-accession-normalised rt_seconds.
+        al = None
+        ra = p.get("rt_aligned")
+        if ra is not None and np.isfinite(ra):
+            al = ra
+        elif rt_lookup is not None and p["sage_psm_id"] is not None:
+            al = rt_lookup.get((acc, int(p["sage_psm_id"])))
+        if al is not None and np.isfinite(al):
+            rt_norm = float(min(1.0, max(0.0, al)))
+            rt_valid = True
+        else:
+            rt = float(p["rt"]) if p["rt"] is not None else float("nan")
+            lo, hi = rt_lo.get(acc), rt_hi.get(acc)
+            if np.isfinite(rt) and lo is not None and hi is not None and hi > lo:
+                rt_norm = (rt - lo) / (hi - lo)
+                rt_valid = True
+            else:
+                rt_norm, rt_valid = float("nan"), False
+
+        ce = float(np.clip((p["ce_v"] or 26.0) / 100.0, 0.0, 1.0))
+
+        examples.append({
+            "accession": acc, "psm_id": int(key[2]),
+            "stripped": stripped, "modseq": modseq,
+            "tokens": np.asarray(residue_ids, dtype=np.int64),
+            "charge": charge, "precursor_mz": mz, "collision_energy": ce,
+            "instrument": 0, "acq_mode": 0,
+            "intensity_target": target,
+            "ccs_target": normalize_inverse_mobility(im),
+            "ccs_valid": bool(np.isfinite(im) and im > 0.0),
+            "rt_target": rt_norm, "rt_valid": rt_valid,
+        })
+    return examples
+
+
+def build_split_datasets_hf(corpus_dir, *, cap=4000, seed=0, accessions=None,
+                            rt_lookup=None):
+    """{split: SagePropertyDataset} from the HF corpus parquets."""
+    if isinstance(rt_lookup, (str, Path)):
+        rt_lookup = load_rt_aligned_lookup(rt_lookup)   # load once, reuse per split
+    return {sp: SagePropertyDataset(
+                prepare_examples_hf(corpus_dir, sp, cap=cap, seed=seed,
+                                    accessions=accessions, rt_lookup=rt_lookup))
+            for sp in ("train", "val", "test")}
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/hf_intensity.py b/packages/peptide-property-ng/src/peptide_property_ng/data/hf_intensity.py
new file mode 100644
index 000000000..8f380a969
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/hf_intensity.py
@@ -0,0 +1,117 @@
+"""Wilhelmlab HuggingFace MS2 datasets -> intensity-pretraining examples.
+
+The three Wilhelmlab MS2 datasets — `timsTOF-ms2`, `prospect-ptms-ms2`,
+`Prosit-2025-lac-ms2` — share one schema:
+  - `intensities_raw`         : the Prosit 174-vector (`-1` = masked)
+  - `modified_sequence`       : `[]-...-[]` terminal placeholders + inline `[UNIMOD:N]`
+  - `precursor_charge_onehot` : 6-way one-hot charge
+  - `collision_energy_aligned_normed` : normalised collision energy
+
+Each row becomes an example dict identical in shape to the Sage loader's, so the
+same collate function and model consume it. CCS / RT are absent here (masked).
+The intensity target is produced by the *same* `prosit174_to_sites` conversion
+the Sage path uses.
+"""
+from __future__ import annotations
+
+import itertools
+import re
+
+import numpy as np
+
+from peptide_property_ng.data.fragment_targets import prosit174_to_sites
+from peptide_property_ng.data.sage_dataset import SagePropertyDataset
+
+_UNIMOD_RE = re.compile(r"\[UNIMOD:\d+\]")
+
+
+def _strip_terminal_placeholders(modseq: str) -> str:
+    """Drop empty `[]-` / `-[]` terminal-modification placeholders.
+
+    A non-empty terminal bracket (a real terminal mod) is left in place; such a
+    peptide then fails the token-count check and is skipped — same policy as the
+    Sage loader.
+    """
+    modseq = re.sub(r"^\[\]-", "", modseq)
+    modseq = re.sub(r"-\[\]$", "", modseq)
+    return modseq
+
+
+def prepare_hf_intensity_examples(
+    dataset_name: str,
+    *,
+    split: str = "train",
+    cap: int | None = None,
+    instrument: int = 0,
+    acq_mode: int = 0,
+    max_len: int = 30,
+    batch: int = 2048,
+) -> list[dict]:
+    """Stream a Wilhelmlab HF MS2 dataset into intensity-pretraining examples."""
+    from datasets import load_dataset
+    from imspy_predictors.utilities.tokenizers import ProformaTokenizer
+
+    tok = ProformaTokenizer.with_defaults()
+    rows = load_dataset(dataset_name, split=split, streaming=True)
+    if cap is not None:
+        rows = itertools.islice(rows, cap)
+
+    examples: list[dict] = []
+    buf: list[dict] = []
+
+    def _flush() -> None:
+        if not buf:
+            return
+        modseqs = [_strip_terminal_placeholders(r["modified_sequence"]) for r in buf]
+        enc = tok(modseqs)
+        ids, attn = enc["input_ids"], enc["attention_mask"]
+        for row, modseq, tok_ids, am in zip(buf, modseqs, ids, attn):
+            if "-" in modseq or "[+" in modseq:
+                continue  # an unresolved terminal mod / delta mass — skip
+            n_real = int(sum(am))
+            residue_ids = tok_ids[1 : n_real - 1]  # strip [CLS]/[SEP] + right-padding
+            length = len(residue_ids)
+            stripped = _UNIMOD_RE.sub("", modseq)
+            if length != len(stripped) or not 3 <= length <= max_len:
+                continue
+            onehot = row["precursor_charge_onehot"]
+            if sum(onehot) != 1:
+                continue  # malformed / all-zero one-hot -> skip rather than default to charge 1
+            charge = int(np.argmax(onehot)) + 1
+            try:
+                target = prosit174_to_sites(row["intensities_raw"], length)
+            except ValueError:
+                continue
+            examples.append(
+                {
+                    "accession": dataset_name,
+                    "stripped": stripped,
+                    "modseq": modseq,
+                    "tokens": np.asarray(residue_ids, dtype=np.int64),
+                    "charge": charge,
+                    "precursor_mz": 0.0,  # absent here; unused by the intensity task
+                    "collision_energy": float(row["collision_energy_aligned_normed"]),
+                    "instrument": int(instrument),
+                    "acq_mode": int(acq_mode),
+                    "intensity_target": target,
+                    "ccs_target": float("nan"), "ccs_valid": False,
+                    "rt_target": float("nan"), "rt_valid": False,
+                }
+            )
+        buf.clear()
+
+    for row in rows:
+        buf.append(row)
+        if len(buf) >= batch:
+            _flush()
+    _flush()
+    return examples
+
+
+def build_hf_intensity_dataset(
+    dataset_name: str, *, split: str = "train", cap: int | None = None, **kwargs
+) -> SagePropertyDataset:
+    """Prepared HF intensity examples wrapped as a Dataset."""
+    return SagePropertyDataset(
+        prepare_hf_intensity_examples(dataset_name, split=split, cap=cap, **kwargs)
+    )
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/mass_to_unimod.py b/packages/peptide-property-ng/src/peptide_property_ng/data/mass_to_unimod.py
new file mode 100644
index 000000000..c06167d26
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/mass_to_unimod.py
@@ -0,0 +1,94 @@
+"""Delta-mass peptide -> UNIMOD conversion.
+
+Sage and the Chronologer DB write peptides with delta-mass brackets
+(``MC[+57.0216]PEPTIDE``); the tokenizer needs UNIMOD form
+(``MC[UNIMOD:4]PEPTIDE``). The mass->UNIMOD step is residue-specificity aware
+(the residue immediately before the bracket), with a small hardcoded table for
+the common search modifications and the quantized sagepy table as the general
+case.
+
+This logic is **vendored verbatim** from ``sagepy_rescore`` (``sage_loader.
+_parse_sage_peptide`` + ``pipeline._mass_to_unimod`` and its tables) so this
+package does not import ``sagepy_rescore.pipeline`` — a heavy module whose
+unrelated ``sagepy.core.fdr`` imports are version-fragile across sagepy
+releases. The conversion itself is unchanged; it needs only ``sagepy.core.unimod``
+and ``sagepy_connector.py_unimod``.
+"""
+from __future__ import annotations
+
+import re
+
+from sagepy.core.unimod import unimod_to_mass
+from sagepy_connector import py_unimod
+
+# Sage / Chronologer write only positive delta masses, e.g. "[+57.0216]".
+_BRACKET_MOD = re.compile(r"\[\+([0-9.]+)\]")
+
+
+def _build_mass_to_unimod() -> dict[float, str]:
+    out: dict[float, str] = {}
+    for unimod_id, mass in unimod_to_mass().items():
+        out[round(float(mass), 4)] = unimod_id
+    return out
+
+
+_MASS_TO_UNIMOD = _build_mass_to_unimod()
+
+# Residue-specificity-disambiguated entries for the common search modifications.
+_COMMON_SAGE_MASS_TO_UNIMOD = {
+    ("C", py_unimod.quanzied_mass(57.0216)): "[UNIMOD:4]",
+    ("M", py_unimod.quanzied_mass(15.9949)): "[UNIMOD:35]",
+    ("[", py_unimod.quanzied_mass(42.0)): "[UNIMOD:1]",
+}
+
+
+def _mass_to_unimod(value, specificity: str | None = None):
+    """Resolve a delta mass (+ optional residue specificity) to a UNIMOD token."""
+    if isinstance(value, list):
+        return [_mass_to_unimod(v, specificity=specificity) for v in value]
+    if isinstance(value, str):
+        return value
+    if isinstance(value, int):
+        return f"[UNIMOD:{value}]"
+
+    quantized = py_unimod.quanzied_mass(float(value))
+    common = _COMMON_SAGE_MASS_TO_UNIMOD.get((specificity, quantized))
+    if common is not None:
+        return common
+
+    candidates = py_unimod.quantized_mass_to_unimod_candidates().get(quantized, [])
+    if candidates:
+        return candidates[0]
+
+    key = round(float(value), 4)
+    if key in _MASS_TO_UNIMOD:
+        return _MASS_TO_UNIMOD[key]
+    raise ValueError(
+        f"Mass {value!r} for specificity {specificity!r} has no UNIMOD entry "
+        "in sagepy's quantized table."
+    )
+
+
+def parse_delta_mass_peptide(peptide: str) -> tuple[str, str]:
+    """Convert a delta-mass peptide string to ``(unmodified, UNIMOD-modified)``.
+
+    Example: ``IIPGFMC[+57.0216]QGGDFTR`` ->
+    ``("IIPGFMCQGGDFTR", "IIPGFMC[UNIMOD:4]QGGDFTR")``. A delta mass with no
+    UNIMOD match is left as the raw ``[+mass]`` bracket (the caller skips it).
+    """
+    out_parts: list[str] = []
+    last_end = 0
+    for m in _BRACKET_MOD.finditer(peptide):
+        start = m.start()
+        out_parts.append(peptide[last_end:start])
+        specificity = "[" if start == 0 else peptide[start - 1]
+        try:
+            unimod = _mass_to_unimod(float(m.group(1)), specificity=specificity)
+        except ValueError:
+            unimod = m.group(0)  # unknown mass -> leave the raw bracket
+        out_parts.append(unimod)
+        last_end = m.end()
+    out_parts.append(peptide[last_end:])
+    modified = "".join(out_parts).replace("-", "")
+    unmodified = _BRACKET_MOD.sub("", peptide).replace("-", "")
+    return unmodified, modified
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/sage_dataset.py b/packages/peptide-property-ng/src/peptide_property_ng/data/sage_dataset.py
new file mode 100644
index 000000000..c1540a0db
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/sage_dataset.py
@@ -0,0 +1,292 @@
+"""Sage-parquet -> multi-task training examples.
+
+Reads a dataset's ``results.sage.parquet`` + ``matched_fragments.sage.parquet``,
+filters to confident target PSMs, converts the Sage delta-mass peptide strings
+to UNIMOD form (via ``sagepy_rescore``'s residue-specificity-aware converter),
+tokenises, and builds per-task targets.
+"""
+from __future__ import annotations
+
+import csv
+import glob
+from collections import defaultdict
+from pathlib import Path
+from types import SimpleNamespace
+
+import numpy as np
+import pyarrow as pa
+import pyarrow.compute as pc
+import pyarrow.parquet as pq
+from torch.utils.data import Dataset
+
+from peptide_property_ng.data.fragment_targets import prosit174_to_sites
+from peptide_property_ng.data.splits import peptide_split
+
+_PROTON = 1.007276
+
+# Inverse ion mobility (1/K0) is min-max normalised to ~[0,1] with fixed bounds,
+# matching the CCS normalisation in ccs_pretrain.py so the CCS/IM head sees one
+# consistent target scale across pretraining (CCS) and fine-tuning (1/K0).
+IM_MIN, IM_MAX = 0.5, 1.9
+
+
+def normalize_inverse_mobility(one_over_k0: float) -> float:
+    """Map an inverse ion mobility (1/K0) value to [0,1], clipped at the fixed bounds."""
+    return float(np.clip((float(one_over_k0) - IM_MIN) / (IM_MAX - IM_MIN), 0.0, 1.0))
+
+_RES_COLS = [
+    "psm_id", "peptide", "stripped_peptide", "charge", "calcmass",
+    "aligned_rt", "ion_mobility", "spectrum_q", "is_decoy", "rank",
+    "matched_peaks", "peptide_len",
+]
+_FRAG_COLS = [
+    "psm_id", "fragment_type", "fragment_ordinals", "fragment_charge", "fragment_intensity",
+]
+
+
+def discover_sage_dirs(glob_pattern: str) -> list[Path]:
+    """Find non-empty ``.../processed/sage`` directories under a dataset-root glob.
+
+    Datasets whose ``results.sage.parquet`` has zero rows (e.g. the legacy
+    ``TimsCompressionType=1`` deposits Sage cannot read) are skipped.
+    """
+    out: list[Path] = []
+    for root in sorted(glob.glob(glob_pattern)):
+        res = Path(root) / "processed" / "sage" / "results.sage.parquet"
+        if not res.exists():
+            continue
+        try:
+            if pq.read_metadata(res).num_rows == 0:
+                continue
+        except Exception:
+            continue
+        out.append(res.parent)
+    return out
+
+
+def _accession(sage_dir: Path) -> str:
+    # .../<ACCESSION>/processed/sage  ->  <ACCESSION>
+    return sage_dir.parts[-3] if len(sage_dir.parts) >= 3 else sage_dir.name
+
+
+def prepare_examples(
+    sage_dir: str | Path,
+    *,
+    cap: int = 8000,
+    q_max: float = 0.01,
+    min_peaks: int = 6,
+    max_len: int = 30,  # the Prosit 174-vector intensity encoding caps peptides at 30 aa
+    seed: int = 0,
+    instrument: int = 0,
+    acq_mode: int = 0,
+    # Sage doesn't surface per-PSM CE for timsTOF diaPASEF; the intensity head's CE
+    # FiLM was tuned by pretraining around the timsTOF-ms2 median (~0.26 normalised).
+    # Default to that here so the head is fed an in-distribution CE rather than 0.
+    default_ce: float = 0.26,
+    ce_calibration: dict[tuple[str, int], float] | None = None,
+) -> list[dict]:
+    """Load and prepare one dataset's PSMs into a list of example dicts."""
+    from imspy_predictors.intensity.predictors import observed_fragments_to_intensity_target
+    from imspy_predictors.utilities.tokenizers import ProformaTokenizer
+
+    from peptide_property_ng.data.mass_to_unimod import parse_delta_mass_peptide
+
+    sage_dir = Path(sage_dir)
+    res_path = sage_dir / "results.sage.parquet"
+    frag_path = sage_dir / "matched_fragments.sage.parquet"
+    if not res_path.exists():
+        return []
+    res = pq.read_table(res_path, columns=_RES_COLS)
+    if res.num_rows == 0:
+        return []
+
+    keep = pc.and_(
+        pc.and_(pc.equal(res["is_decoy"], False), pc.equal(res["rank"], 1)),
+        pc.and_(
+            pc.less_equal(res["spectrum_q"], q_max),
+            pc.and_(
+                pc.greater_equal(res["matched_peaks"], min_peaks),
+                pc.less_equal(res["peptide_len"], max_len),
+            ),
+        ),
+    )
+    res = res.filter(keep)
+    if res.num_rows == 0:
+        return []
+    if res.num_rows > cap:
+        idx = np.sort(np.random.RandomState(seed).choice(res.num_rows, cap, replace=False))
+        res = res.take(pa.array(idx))
+    df = res.to_pydict()
+    n = len(df["psm_id"])
+
+    # Fragments for the kept PSMs, grouped by psm_id.
+    frag = pq.read_table(frag_path, columns=_FRAG_COLS)
+    frag = frag.filter(pc.is_in(frag["psm_id"], value_set=pa.array(list(set(df["psm_id"])))))
+    fd = frag.to_pydict()
+    by_psm: dict[int, tuple[list, list, list, list]] = defaultdict(lambda: ([], [], [], []))
+    for pid, ft, fo, fc, fi in zip(
+        fd["psm_id"], fd["fragment_type"], fd["fragment_ordinals"],
+        fd["fragment_charge"], fd["fragment_intensity"],
+    ):
+        cols = by_psm[pid]
+        cols[0].append(ft); cols[1].append(fo); cols[2].append(fc); cols[3].append(fi)
+
+    # Sage emits delta-mass peptides; convert to UNIMOD, then tokenise in one batch.
+    # Batch tokenisation right-pads to the batch-max length, so the attention
+    # mask is needed to recover each peptide's true token count.
+    parsed = [parse_delta_mass_peptide(p) for p in df["peptide"]]
+    tok = ProformaTokenizer.with_defaults()
+    _encoded = tok([m for _, m in parsed])
+    token_ids, attn_mask = _encoded["input_ids"], _encoded["attention_mask"]
+
+    accession = _accession(sage_dir)
+    examples: list[dict] = []
+    for k in range(n):
+        stripped, modseq = parsed[k]
+        if "[+" in modseq or "[-" in modseq:
+            continue  # an unconverted delta mass — skip rather than train on a dropped mod
+        n_real = int(sum(attn_mask[k]))  # real tokens incl. [CLS]/[SEP], excl. padding
+        residue_ids = token_ids[k][1 : n_real - 1]  # strip [CLS]/[SEP] and right-padding
+        length = len(residue_ids)
+        # Terminal mods become standalone tokens -> token count != residue count;
+        # skip those rather than mis-align the fragment indexing.
+        if length != int(df["peptide_len"][k]) or length < 3 or length > max_len:
+            continue
+
+        charge = int(df["charge"][k])
+        mz = (float(df["calcmass"][k]) + charge * _PROTON) / charge
+        psm_id = int(df["psm_id"][k])
+        cols = by_psm.get(psm_id)
+        if not cols or not cols[0]:
+            continue  # no matched fragments -> no intensity target
+        # Encode fragments with the proven imspy encoder (Sage fragments ->
+        # Prosit 174-vector), then the single shared 174 -> site conversion.
+        frag = SimpleNamespace(
+            ion_types=cols[0], fragment_ordinals=cols[1],
+            charges=cols[2], intensities=cols[3],
+        )
+        prosit174 = observed_fragments_to_intensity_target(stripped, charge, frag)
+        target = prosit174_to_sites(prosit174, length)
+
+        im = df["ion_mobility"][k]
+        rt = df["aligned_rt"][k]
+        im = float(im) if im is not None else float("nan")
+        rt = float(rt) if rt is not None else float("nan")
+        ce_val = float(default_ce)
+        if ce_calibration is not None:
+            cal = ce_calibration.get((accession, psm_id))
+            if cal is not None:
+                ce_val = float(cal)
+        examples.append(
+            {
+                "accession": accession,
+                "psm_id": psm_id,
+                "stripped": stripped,
+                "modseq": modseq,
+                "tokens": np.asarray(residue_ids, dtype=np.int64),
+                "charge": charge,
+                "precursor_mz": mz,
+                "collision_energy": ce_val,
+                "instrument": int(instrument),
+                "acq_mode": int(acq_mode),
+                "intensity_target": target,
+                "ccs_target": normalize_inverse_mobility(im),
+                "ccs_valid": bool(np.isfinite(im) and im > 0.0),
+                "rt_target": rt,
+                "rt_valid": bool(np.isfinite(rt)),
+            }
+        )
+    return examples
+
+
+class SagePropertyDataset(Dataset):
+    """A thin wrapper over a list of prepared example dicts."""
+
+    def __init__(self, examples: list[dict]):
+        self.examples = examples
+
+    def __len__(self) -> int:
+        return len(self.examples)
+
+    def __getitem__(self, i: int) -> dict:
+        return self.examples[i]
+
+
+def load_catalog_metadata(catalog_path: str | Path) -> dict[str, tuple[int, int]]:
+    """Read ``timstof_catalog.tsv`` -> ``{accession: (instrument_id, acq_id)}``.
+
+    Unknown / missing values map to id 0 via ``instrument_id`` / ``acquisition_id``.
+    """
+    from peptide_property_ng.model.config import acquisition_id, instrument_id
+
+    out: dict[str, tuple[int, int]] = {}
+    with open(catalog_path) as f:
+        reader = csv.DictReader(f, delimiter="\t")
+        for row in reader:
+            acc = (row.get("accession") or "").strip()
+            if not acc:
+                continue
+            out[acc] = (
+                instrument_id(row.get("instrument_model")),
+                acquisition_id(row.get("acquisition_mode")),
+            )
+    return out
+
+
+def load_ce_calibration(parquet_path: str | Path) -> dict[tuple[str, int], float]:
+    """Read a CE calibration parquet -> ``{(accession, psm_id): calibrated_ce}``.
+
+    The parquet must have columns ``accession`` (string), ``psm_id`` (int64),
+    ``calibrated_ce`` (float). Missing keys at lookup time fall back to the
+    dataset's ``default_ce`` — so partial calibrations are safe.
+    """
+    table = pq.read_table(parquet_path)
+    accs = table["accession"].to_pylist()
+    pids = table["psm_id"].to_pylist()
+    ces = table["calibrated_ce"].to_pylist()
+    return {(a, int(p)): float(c) for a, p, c in zip(accs, pids, ces)}
+
+
+def build_split_datasets(
+    sage_dirs: list[str | Path],
+    *,
+    cap: int = 8000,
+    seed: int = 0,
+    val_frac: float = 0.1,
+    test_frac: float = 0.1,
+    catalog_path: str | Path | None = None,
+    ce_calibration: dict[tuple[str, int], float] | None = None,
+    **prepare_kwargs,
+) -> dict[str, SagePropertyDataset]:
+    """Prepare every dataset and split peptide-level into train / val / test.
+
+    If ``catalog_path`` is given (a ``timstof_catalog.tsv``), per-dataset
+    ``instrument`` and ``acq_mode`` ids are looked up by accession and passed to
+    ``prepare_examples`` — so the encoder's metadata conditioning actually sees
+    the right instrument instead of falling back to ``unknown``. Otherwise the
+    loader's defaults (0/unknown) are used.
+
+    If ``ce_calibration`` is given, each PSM's collision energy is overridden by
+    the calibrated value for that ``(accession, psm_id)``; PSMs not in the map
+    fall back to ``default_ce``.
+    """
+    metadata: dict[str, tuple[int, int]] = (
+        load_catalog_metadata(catalog_path) if catalog_path else {}
+    )
+    examples: list[dict] = []
+    for sage_dir in sage_dirs:
+        sage_dir_p = Path(sage_dir)
+        accession = _accession(sage_dir_p)
+        kw = dict(prepare_kwargs)
+        if accession in metadata:
+            instr, acq = metadata[accession]
+            kw.setdefault("instrument", instr)
+            kw.setdefault("acq_mode", acq)
+        examples.extend(prepare_examples(sage_dir_p, cap=cap, seed=seed,
+                                         ce_calibration=ce_calibration, **kw))
+
+    buckets: dict[str, list[dict]] = {"train": [], "val": [], "test": []}
+    for ex in examples:
+        split = peptide_split(ex["stripped"], val_frac=val_frac, test_frac=test_frac, seed=seed)
+        buckets[split].append(ex)
+    return {name: SagePropertyDataset(rows) for name, rows in buckets.items()}
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/data/splits.py b/packages/peptide-property-ng/src/peptide_property_ng/data/splits.py
new file mode 100644
index 000000000..c85b05b9f
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/data/splits.py
@@ -0,0 +1,48 @@
+"""Leakage-free dataset splits.
+
+Peptide-level split — a peptide (its *stripped*, modification-free sequence) is
+hashed to a fixed bucket, so the same peptide never appears in two splits
+regardless of which dataset, charge or run it came from. This is essential:
+peptides recur heavily across the deposits.
+
+Also provides the unseen-modification split used by the key generalisation
+probe — hold every PSM bearing a chosen modification out of train/val.
+"""
+from __future__ import annotations
+
+import hashlib
+import re
+
+_UNIMOD_RE = re.compile(r"\[UNIMOD:(\d+)\]")
+
+
+def _bucket(key: str, seed: int) -> float:
+    """Stable hash of ``key`` into ``[0, 1)``."""
+    digest = hashlib.sha1(f"{seed}:{key}".encode()).hexdigest()
+    return (int(digest[:8], 16) % 1_000_000) / 1_000_000.0
+
+
+def peptide_split(
+    stripped_peptide: str,
+    *,
+    val_frac: float = 0.1,
+    test_frac: float = 0.1,
+    seed: int = 0,
+) -> str:
+    """Assign a peptide to ``"train"`` / ``"val"`` / ``"test"`` by stable hash."""
+    f = _bucket(stripped_peptide, seed)
+    if f < test_frac:
+        return "test"
+    if f < test_frac + val_frac:
+        return "val"
+    return "train"
+
+
+def modifications_in(modified_peptide: str) -> set[int]:
+    """UNIMOD ids present in a UNIMOD-format peptide string."""
+    return {int(x) for x in _UNIMOD_RE.findall(modified_peptide)}
+
+
+def has_modification(modified_peptide: str, unimod_id: int) -> bool:
+    """True if the peptide carries the given UNIMOD modification."""
+    return unimod_id in modifications_in(modified_peptide)
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/eval/__init__.py b/packages/peptide-property-ng/src/peptide_property_ng/eval/__init__.py
new file mode 100644
index 000000000..a1e61af4e
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/eval/__init__.py
@@ -0,0 +1 @@
+"""Evaluation — per-task metrics, OOD probes, the unseen-modification split."""
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/eval/evaluate.py b/packages/peptide-property-ng/src/peptide_property_ng/eval/evaluate.py
new file mode 100644
index 000000000..0b32fa81f
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/eval/evaluate.py
@@ -0,0 +1,52 @@
+"""Evaluate a trained checkpoint on the held-out test split.
+
+    python -m peptide_property_ng.eval.evaluate --checkpoint runs/prototype/best.pt
+"""
+from __future__ import annotations
+
+import argparse
+import json
+
+import torch
+from torch.utils.data import DataLoader
+
+from peptide_property_ng.data.collate import make_collate_fn
+from peptide_property_ng.data.sage_dataset import build_split_datasets, discover_sage_dirs
+from peptide_property_ng.eval.metrics import evaluate_split
+from peptide_property_ng.model.config import PRESETS
+from peptide_property_ng.model.multitask import UnifiedPeptidePropertyModel
+from peptide_property_ng.modifications.composition import CompositionTable
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description="Evaluate a trained checkpoint.")
+    ap.add_argument("--checkpoint", required=True)
+    ap.add_argument("--datasets-glob", default="/scratch/claudius-proteomics/*")
+    ap.add_argument("--cap", type=int, default=4000)
+    ap.add_argument("--max-datasets", type=int, default=0)
+    ap.add_argument("--batch-size", type=int, default=64)
+    ap.add_argument("--seed", type=int, default=0)
+    ap.add_argument("--split", default="test", choices=["train", "val", "test"])
+    ap.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu")
+    args = ap.parse_args()
+
+    ckpt = torch.load(args.checkpoint, map_location=args.device)
+    cfg = PRESETS[ckpt.get("preset", "small")]
+    model = UnifiedPeptidePropertyModel(cfg, CompositionTable.load()).to(args.device)
+    model.load_state_dict(ckpt["model_state_dict"])
+
+    sage_dirs = discover_sage_dirs(args.datasets_glob)
+    if args.max_datasets:
+        sage_dirs = sage_dirs[: args.max_datasets]
+    splits = build_split_datasets(sage_dirs, cap=args.cap, seed=args.seed)
+    loader = DataLoader(
+        splits[args.split],
+        batch_size=args.batch_size,
+        collate_fn=make_collate_fn(cfg.pad_token_id, cfg.max_charge),
+    )
+    metrics = evaluate_split(model, loader, args.device)
+    print(json.dumps({"split": args.split, "checkpoint": args.checkpoint, "metrics": metrics}, indent=2))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/eval/evaluate_optimal_nce.py b/packages/peptide-property-ng/src/peptide_property_ng/eval/evaluate_optimal_nce.py
new file mode 100644
index 000000000..a38b72daa
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/eval/evaluate_optimal_nce.py
@@ -0,0 +1,164 @@
+"""Per-PSM optimal NCE calibration — v4's inference recipe applied to ppng.
+
+For each PSM, sweep the model over a grid of collision energies and pick the CE
+that maximises spectral-angle similarity to the observed spectrum. Reports the
+mean optimal-CE SA, the mean fixed-CE SA at the original (per-PSM dataset/catalog)
+value, and the optimal-CE distribution. No retraining — this isolates "how good
+is the model when the right CE is used per spectrum" from "how good is the
+flat CE we used at fine-tune time."
+
+Example:
+    python -m peptide_property_ng.eval.evaluate_optimal_nce \\
+        --checkpoint runs/finetune-canonicalSA/best.pt \\
+        --datasets-glob '/scratch/claudius-proteomics/*' \\
+        --catalog .../timstof_catalog.tsv
+"""
+from __future__ import annotations
+
+import argparse
+from pathlib import Path
+
+import torch
+from torch.utils.data import DataLoader
+
+from peptide_property_ng.data.collate import make_collate_fn
+from peptide_property_ng.data.sage_dataset import (
+    build_split_datasets, discover_sage_dirs, load_ce_calibration,
+)
+from peptide_property_ng.losses import intensity_signal_mask, masked_spectral_angle
+from peptide_property_ng.model.config import PRESETS
+from peptide_property_ng.model.multitask import UnifiedPeptidePropertyModel
+from peptide_property_ng.modifications.composition import CompositionTable
+
+
+@torch.no_grad()
+def evaluate_optimal_nce(
+    model,
+    loader,
+    device: str,
+    ce_grid: torch.Tensor,
+):
+    """Return per-PSM (opt_sa, opt_ce, fixed_sa) tensors over the signal subset.
+
+    ``fixed_sa`` uses the loader's original (per-PSM dataset/catalog) CE — i.e.
+    the same CE the model saw during fine-tune; ``opt_sa`` is the per-PSM max
+    over the supplied grid.
+    """
+    model.eval()
+    ce_grid = ce_grid.to(device)
+    n_ce = ce_grid.numel()
+
+    opt_sa_list: list[torch.Tensor] = []
+    opt_ce_list: list[torch.Tensor] = []
+    fixed_sa_list: list[torch.Tensor] = []
+
+    for batch in loader:
+        batch = {k: v.to(device) for k, v in batch.items()}
+        target = batch["intensity_target"]
+        signal = intensity_signal_mask(target)
+        if not signal.any():
+            continue
+
+        B = target.shape[0]
+        ce_orig = batch["collision_energy"].clone()
+
+        # Per-CE SA similarity (B, n_ce).
+        sa_matrix = torch.zeros(B, n_ce, device=device)
+        for j, ce in enumerate(ce_grid):
+            batch["collision_energy"] = ce.expand(B).float()
+            pred = model(batch, tasks=["intensity"])["intensity"]
+            sa_matrix[:, j] = 1.0 - masked_spectral_angle(pred, target)
+
+        opt_sa, opt_idx = sa_matrix.max(dim=1)
+        opt_ce = ce_grid[opt_idx]
+
+        # Fixed reference: model's original CE (reproduces the training test SA).
+        batch["collision_energy"] = ce_orig
+        pred = model(batch, tasks=["intensity"])["intensity"]
+        fixed_sa = 1.0 - masked_spectral_angle(pred, target)
+
+        opt_sa_list.append(opt_sa[signal].cpu())
+        opt_ce_list.append(opt_ce[signal].cpu())
+        fixed_sa_list.append(fixed_sa[signal].cpu())
+
+    return (
+        torch.cat(opt_sa_list),
+        torch.cat(opt_ce_list),
+        torch.cat(fixed_sa_list),
+    )
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--checkpoint", required=True, help="path to a best.pt")
+    ap.add_argument("--datasets-glob", default="/scratch/claudius-proteomics/*")
+    ap.add_argument("--catalog", default=None)
+    ap.add_argument("--default-ce", type=float, default=0.26)
+    ap.add_argument("--ce-calibration", default=None,
+                    help="parquet of per-PSM calibrated CEs; if supplied, the 'fixed' baseline "
+                         "uses these (matches a fine-tune that was trained with the same calibration)")
+    ap.add_argument("--cap", type=int, default=4000)
+    ap.add_argument("--max-datasets", type=int, default=0)
+    ap.add_argument("--batch-size", type=int, default=64)
+    ap.add_argument("--seed", type=int, default=0)
+    ap.add_argument("--ce-min", type=float, default=0.10)
+    ap.add_argument("--ce-max", type=float, default=0.40)
+    ap.add_argument("--ce-step", type=float, default=0.02)
+    ap.add_argument("--split", default="test", choices=["train", "val", "test"])
+    ap.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu")
+    args = ap.parse_args()
+
+    print(f"checkpoint: {args.checkpoint}")
+    sage_dirs = discover_sage_dirs(args.datasets_glob)
+    if args.max_datasets:
+        sage_dirs = sage_dirs[: args.max_datasets]
+    print(f"  {len(sage_dirs)} datasets, cap {args.cap}/dataset")
+    ce_cal = load_ce_calibration(args.ce_calibration) if args.ce_calibration else None
+    splits = build_split_datasets(sage_dirs, cap=args.cap, seed=args.seed,
+                                  catalog_path=args.catalog, default_ce=args.default_ce,
+                                  ce_calibration=ce_cal)
+    ds = splits[args.split]
+
+    import dataclasses
+    ckpt = torch.load(args.checkpoint, map_location=args.device, weights_only=False)
+    preset = ckpt.get("preset", "small")
+    cfg = PRESETS[preset]
+    intensity_head = ckpt.get("intensity_head")
+    if intensity_head is None:
+        keys = ckpt["model_state_dict"].keys()
+        intensity_head = "pooled" if any("heads.intensity.attention." in k for k in keys) else "site"
+    cfg = dataclasses.replace(cfg, intensity_head=intensity_head)
+    collate = make_collate_fn(cfg.pad_token_id, cfg.max_charge)
+    loader = DataLoader(ds, batch_size=args.batch_size, shuffle=False, collate_fn=collate)
+    model = UnifiedPeptidePropertyModel(cfg, CompositionTable.load()).to(args.device)
+    model.load_state_dict(ckpt["model_state_dict"])
+
+    n_pts = int(round((args.ce_max - args.ce_min) / args.ce_step)) + 1
+    ce_grid = torch.linspace(args.ce_min, args.ce_max, n_pts)
+    print(f"split={args.split}, |dataset|={len(ds):,}, preset={preset}")
+    print(f"CE grid ({n_pts} pts): {args.ce_min:.3f} .. {args.ce_max:.3f} step {args.ce_step:.3f}")
+    print()
+
+    opt_sa, opt_ce, fixed_sa = evaluate_optimal_nce(model, loader, args.device, ce_grid)
+    n = opt_sa.numel()
+    lift = float(opt_sa.mean() - fixed_sa.mean())
+    print(f"PSMs scored:                       {n:,}")
+    print(f"fixed (original per-PSM CE)  SA:   {float(fixed_sa.mean()):.4f}")
+    print(f"per-PSM optimal CE           SA:   {float(opt_sa.mean()):.4f}   (lift +{lift:.4f})")
+    print()
+    print("optimal CE distribution (per-PSM):")
+    for q in [0.05, 0.25, 0.5, 0.75, 0.95]:
+        print(f"  q{int(q*100):02d}   {float(opt_ce.quantile(q)):.3f}")
+    print(f"  mean   {float(opt_ce.mean()):.3f}")
+    # Mode by histogram (CE grid is discrete, so this is exact).
+    vals, counts = torch.unique(opt_ce, return_counts=True)
+    print(f"  mode   {float(vals[counts.argmax()]):.3f}  ({int(counts.max())} / {n} PSMs)")
+    print()
+    print("share of PSMs at each CE grid point:")
+    for v, c in zip(vals.tolist(), counts.tolist()):
+        bar = "#" * int(40 * c / n)
+        print(f"  {v:.3f}  {c:>5d}  {c/n:5.1%}  {bar}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/eval/metrics.py b/packages/peptide-property-ng/src/peptide_property_ng/eval/metrics.py
new file mode 100644
index 000000000..fc2f0b96c
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/eval/metrics.py
@@ -0,0 +1,73 @@
+"""Per-task evaluation metrics."""
+from __future__ import annotations
+
+import torch
+
+from peptide_property_ng.losses import intensity_signal_mask, masked_spectral_angle
+
+
+def _pearson(x: torch.Tensor, y: torch.Tensor) -> float:
+    if x.numel() < 2:
+        return float("nan")
+    x = x - x.mean()
+    y = y - y.mean()
+    denom = x.norm() * y.norm()
+    return float((x * y).sum() / denom) if denom > 0 else float("nan")
+
+
+@torch.no_grad()
+def evaluate_split(model, loader, device: str = "cpu",
+                   tasks: list[str] | None = None) -> dict[str, float]:
+    """Run the model over a DataLoader and return per-task metrics.
+
+    ``tasks`` restricts the forward to the named tasks (e.g. ``["intensity"]``);
+    metrics for skipped tasks come back as ``nan``.
+
+    Returns: intensity spectral angle (similarity, higher better), CCS / RT
+    median absolute error, RT Pearson r, charge accuracy.
+    """
+    model.eval()
+    sa_sum = sa_n = 0.0
+    ccs_abs: list[torch.Tensor] = []
+    rt_pred: list[torch.Tensor] = []
+    rt_true: list[torch.Tensor] = []
+    charge_correct = charge_n = 0.0
+
+    for batch in loader:
+        batch = {k: v.to(device) for k, v in batch.items()}
+        out = model(batch, tasks=tasks)
+
+        if "intensity" in out:
+            sa = 1.0 - masked_spectral_angle(out["intensity"], batch["intensity_target"])
+            signal = intensity_signal_mask(batch["intensity_target"])
+            sa_sum += float(sa[signal].sum())
+            sa_n += int(signal.sum())
+
+        if "ccs" in out:
+            mean, _ = out["ccs"]
+            v = batch["ccs_valid"]
+            if v.any():
+                ccs_abs.append((mean[v] - batch["ccs_target"][v]).abs().cpu())
+
+        if "rt" in out:
+            v = batch["rt_valid"]
+            if v.any():
+                rt_pred.append(out["rt"][v].cpu())
+                rt_true.append(batch["rt_target"][v].cpu())
+
+        if "charge" in out:
+            pred = out["charge"].argmax(dim=-1)
+            charge_correct += float((pred == batch["charge_target"]).sum())
+            charge_n += pred.numel()
+
+    ccs = torch.cat(ccs_abs) if ccs_abs else torch.tensor([])
+    rtp = torch.cat(rt_pred) if rt_pred else torch.tensor([])
+    rtt = torch.cat(rt_true) if rt_true else torch.tensor([])
+    return {
+        "intensity_sa": sa_sum / sa_n if sa_n else float("nan"),
+        "ccs_mae": float(ccs.mean()) if ccs.numel() else float("nan"),
+        "ccs_median_ae": float(ccs.median()) if ccs.numel() else float("nan"),
+        "rt_mae": float((rtp - rtt).abs().mean()) if rtp.numel() else float("nan"),
+        "rt_pearson": _pearson(rtp, rtt),
+        "charge_acc": charge_correct / charge_n if charge_n else float("nan"),
+    }
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/losses.py b/packages/peptide-property-ng/src/peptide_property_ng/losses.py
new file mode 100644
index 000000000..b2389b04e
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/losses.py
@@ -0,0 +1,131 @@
+"""Multi-task losses for the unified peptide property model."""
+from __future__ import annotations
+
+import math
+
+import torch
+from torch.nn import functional as F
+
+
+def masked_spectral_angle(
+    pred: torch.Tensor,
+    target: torch.Tensor,
+    norm_eps: float = 1e-8,
+    clamp_eps: float = 1e-4,
+) -> torch.Tensor:
+    """Per-sample spectral-angle distance in ``[0, 1]`` (0 = identical spectra).
+
+    Canonical Prosit / Sage-fine-tune loss: mask ``target > 0`` (only observed
+    peaks contribute). Zeros are treated as "unmatched / unknown" rather than
+    real "no peak" labels — which is the correct interpretation for Sage
+    ``matched_fragments`` (Sage only reports peaks it matched) and matches the
+    production v4 / imspy_predictors `masked_spectral_distance` convention.
+
+    The two eps are deliberately decoupled:
+      - ``norm_eps`` (1e-8) keeps ``F.normalize`` accurate for non-degenerate
+        magnitudes;
+      - ``clamp_eps`` (1e-4) bounds ``arccos``'s near-boundary derivative
+        (≈ -1/sqrt(2·clamp_eps), so ~70 instead of ~7000 at 1e-8). With 1e-8
+        the SA-loss backward NaN-s on sparse one-positive-fragment targets
+        whose normalize-to-one-hot forces cos exactly to ±1; with 1e-4 those
+        same samples are absorbed by the clamp (grad 0 above 1-clamp_eps).
+
+    Pretraining on PROSPECT-style densely-annotated targets is unaffected in
+    practice — observable b/y positions almost all carry positive intensity,
+    so the > 0 vs >= 0 masks coincide and cos rarely hits the boundary.
+    """
+    # Pooled head outputs a fixed (B, 29, n_ion) grid; site head outputs
+    # (B, L-1, n_ion). Slice pred to the batch's target sites so both shapes
+    # work transparently here.
+    if pred.shape[1] > target.shape[1]:
+        pred = pred[:, : target.shape[1], :]
+    mask = (target > 0).float()
+    p = (pred * mask).reshape(pred.shape[0], -1)
+    t = (target * mask).reshape(target.shape[0], -1)
+    p = F.normalize(p, dim=1, eps=norm_eps)
+    t = F.normalize(t, dim=1, eps=norm_eps)
+    cos = (p * t).sum(dim=1).clamp(-1.0 + clamp_eps, 1.0 - clamp_eps)
+    return 2.0 * torch.arccos(cos) / math.pi
+
+
+def intensity_signal_mask(target: torch.Tensor) -> torch.Tensor:
+    """``(B,)`` bool — True where a target has at least one observed (positive) peak.
+
+    A target with no positive b/y charge-1-3 peak is degenerate: the spectral
+    angle against it is a constant with zero gradient, so such samples are
+    excluded from the intensity loss and metric rather than diluting them.
+    """
+    return (target > 0).flatten(1).any(dim=1)
+
+
+class MultiTaskLoss:
+    """Weighted sum of per-task losses; tasks with no valid targets are skipped.
+
+    - intensity : masked spectral-angle distance
+    - ccs       : L1 on the mean 1/K0
+    - im_sigma  : L1 on the ccs head's std vs the measured IM Gaussian width
+                  (timsim shape target; masked to peptides with a valid label)
+    - rt        : L1 on normalised retention time
+    - rt_sigma  : L1 on the gradient-normalised RT EMG peak-width (timsim shape)
+    - charge    : cross-entropy
+
+    The two shape terms (``im_sigma``, ``rt_sigma``) are *additive and optional*:
+    they contribute only when the batch carries the corresponding target+mask,
+    so existing (shape-free) training is unchanged.
+    """
+
+    DEFAULT_WEIGHTS = {"intensity": 1.0, "ccs": 1.0, "im_sigma": 0.5,
+                       "rt": 0.5, "rt_sigma": 0.5, "charge": 0.3}
+
+    def __init__(self, weights: dict[str, float] | None = None):
+        self.weights = dict(self.DEFAULT_WEIGHTS if weights is None else weights)
+
+    def __call__(
+        self,
+        outputs: dict[str, object],
+        batch: dict[str, torch.Tensor],
+    ) -> tuple[torch.Tensor, dict[str, float]]:
+        parts: dict[str, torch.Tensor] = {}
+
+        if "intensity" in outputs:
+            sa = masked_spectral_angle(outputs["intensity"], batch["intensity_target"])
+            signal = intensity_signal_mask(batch["intensity_target"])
+            if signal.any():
+                parts["intensity"] = sa[signal].mean()
+
+        if "ccs" in outputs:
+            mean, std = outputs["ccs"]
+            valid = batch["ccs_valid"]
+            if valid.any():
+                # L1 on the mean — stable for the prototype. Gaussian-NLL on
+                # (mean, std) is the eventual target but is unbounded-below and
+                # destabilises the fixed-weight multi-task sum; revisit together
+                # with uncertainty-based task weighting.
+                parts["ccs"] = F.l1_loss(mean[valid], batch["ccs_target"][valid])
+            # IM peak-width (timsim Gaussian sigma): supervise the head's std on
+            # the subset of peptides that carry a measured IM-sigma label. L1
+            # (not Gaussian-NLL) per the design — the std is the timsim shape param.
+            if "im_sigma_target" in batch and "im_sigma_valid" in batch:
+                sv = batch["im_sigma_valid"]
+                if sv.any():
+                    parts["im_sigma"] = F.l1_loss(std[sv], batch["im_sigma_target"][sv])
+
+        if "rt" in outputs:
+            valid = batch["rt_valid"]
+            if valid.any():
+                parts["rt"] = F.l1_loss(outputs["rt"][valid], batch["rt_target"][valid])
+
+        if "rt_sigma" in outputs and "rt_sigma_target" in batch:
+            sv = batch["rt_sigma_valid"]
+            if sv.any():
+                parts["rt_sigma"] = F.l1_loss(
+                    outputs["rt_sigma"][sv], batch["rt_sigma_target"][sv]
+                )
+
+        if "charge" in outputs:
+            parts["charge"] = F.cross_entropy(outputs["charge"], batch["charge_target"])
+
+        if not parts:
+            raise ValueError("no task produced a loss for this batch")
+        total = sum(self.weights.get(name, 1.0) * value for name, value in parts.items())
+        return total, {name: float(value.detach()) for name, value in parts.items()}
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/__init__.py b/packages/peptide-property-ng/src/peptide_property_ng/model/__init__.py
new file mode 100644
index 000000000..d3956ded4
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/__init__.py
@@ -0,0 +1 @@
+"""Model — hybrid residue embedding, Depthcharge-based encoder, multi-task heads."""
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/config.py b/packages/peptide-property-ng/src/peptide_property_ng/model/config.py
new file mode 100644
index 000000000..1582cde9c
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/config.py
@@ -0,0 +1,127 @@
+"""Model configuration — presets, instrument / acquisition vocabularies.
+
+Conditioning policy (deliberate, see the charge-leakage note):
+  - Instrument model and acquisition mode are conditioned *inside* the encoder
+    via the prepended global token — they cannot leak any task label.
+  - Charge, precursor m/z and collision energy are conditioned at the *heads*
+    that legitimately use them (intensity, CCS). They are deliberately kept out
+    of the encoder so the charge head — which reads the shared encoder output —
+    never sees its own label. One encoder pass, no leakage.
+"""
+from __future__ import annotations
+
+from dataclasses import dataclass, field
+
+# Instrument vocabulary — timsTOF-focused (the target domain), plus `orbitrap`
+# so Orbitrap-HCD pretraining corpora (PROSPECT) carry a distinct id and the
+# instrument embedding can learn the HCD<->PASEF domain shift.
+INSTRUMENTS: tuple[str, ...] = (
+    "unknown",
+    "timsTOF",
+    "timsTOF Pro",
+    "timsTOF Pro 2",
+    "timsTOF HT",
+    "timsTOF SCP",
+    "timsTOF Ultra",
+    "timsTOF fleX",
+    "orbitrap",
+)
+ACQUISITION_MODES: tuple[str, ...] = (
+    "unknown",
+    "DDA",
+    "DDA-PASEF",
+    "DIA",
+    "diaPASEF",
+    "PRM",
+    "MALDI",
+    "PASEF",        # catalog's dominant value (716 deposits)
+    "single-cell",  # catalog: 70 deposits
+)
+INSTRUMENT_TO_ID: dict[str, int] = {n: i for i, n in enumerate(INSTRUMENTS)}
+ACQ_TO_ID: dict[str, int] = {n: i for i, n in enumerate(ACQUISITION_MODES)}
+
+
+def instrument_id(name: str | None) -> int:
+    """Map an instrument-model string to an id, tolerating unknowns/aliases."""
+    if not name:
+        return 0
+    n = " ".join(str(name).split())  # normalise whitespace
+    if n in INSTRUMENT_TO_ID:
+        return INSTRUMENT_TO_ID[n]
+    low = n.lower()
+    for cand, idx in INSTRUMENT_TO_ID.items():
+        if cand.lower() == low:
+            return idx
+    return 0  # "unknown"
+
+
+def acquisition_id(name: str | None) -> int:
+    """Map an acquisition-mode string to an id, tolerating unknowns."""
+    if not name:
+        return 0
+    n = str(name).strip()
+    for cand, idx in ACQ_TO_ID.items():
+        if cand.lower() == n.lower():
+            return idx
+    return 0
+
+
+@dataclass
+class ModelConfig:
+    """Architecture + vocabulary configuration for the unified predictor."""
+
+    # encoder
+    d_model: int = 256
+    n_layers: int = 6
+    n_heads: int = 8
+    dim_feedforward: int = 1024
+    dropout: float = 0.1
+    max_seq_len: int = 64  # peptide-length cap — far above the old 174-vector's 30
+
+    # tokenizer / modification vocabulary
+    vocab_size: int = 2175
+    pad_token_id: int = 61
+    n_elements: int = 31  # atomic-composition channels (see modifications.composition)
+
+    # hybrid embedding fusion: "add" | "gate" | "token_only" | "composition_only"
+    # ("token_only" / "composition_only" are ablation modes)
+    comp_fusion: str = "add"
+
+    # conditioning
+    max_charge: int = 8
+    # Embedding *capacity* for instrument / acquisition-mode ids, deliberately
+    # larger than the current vocabularies. A new instrument or mode appended to
+    # INSTRUMENTS / ACQUISITION_MODES later gets the next free id without
+    # resizing the embedding -- so existing checkpoints stay loadable and the
+    # new id simply starts from an untrained row, ready to fine-tune.
+    n_instruments: int = 64  # 9 ids in use (see INSTRUMENTS); the rest reserved
+    n_acq_modes: int = 32    # 7 ids in use (see ACQUISITION_MODES); the rest reserved
+
+    # heads
+    n_ion_channels: int = 6  # b/y ions x fragment charges 1-3 (Prosit-compatible)
+    tasks: tuple[str, ...] = ("intensity", "ccs", "rt", "charge")
+    # intensity head topology: "site" (local FiLM, per-cleavage-site) or
+    # "pooled" (v4-style attention-pooled Prosit-174); see heads/intensity_pooled.py
+    intensity_head: str = "site"
+
+    def __post_init__(self) -> None:
+        if self.d_model % self.n_heads != 0:
+            raise ValueError(f"d_model {self.d_model} not divisible by n_heads {self.n_heads}")
+        if len(INSTRUMENTS) > self.n_instruments:
+            raise ValueError(
+                f"{len(INSTRUMENTS)} instruments exceed embedding capacity {self.n_instruments}"
+            )
+        if len(ACQUISITION_MODES) > self.n_acq_modes:
+            raise ValueError(
+                f"{len(ACQUISITION_MODES)} acquisition modes exceed capacity {self.n_acq_modes}"
+            )
+
+
+# First-prototype preset — roughly on par with the production BASE encoder so
+# iteration is fast on the ~19 datasets available now.
+SMALL = ModelConfig(d_model=256, n_layers=6, n_heads=8, dim_feedforward=1024)
+
+# Research preset — larger; for the full ~138-dataset corpus.
+RESEARCH = ModelConfig(d_model=384, n_layers=8, n_heads=8, dim_feedforward=1536)
+
+PRESETS: dict[str, ModelConfig] = {"small": SMALL, "research": RESEARCH}
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/embedding.py b/packages/peptide-property-ng/src/peptide_property_ng/model/embedding.py
new file mode 100644
index 000000000..38e28b6a5
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/embedding.py
@@ -0,0 +1,102 @@
+"""Hybrid residue embedding — learnable token embedding + atomic-composition.
+
+The central design decision: every residue is represented by BOTH
+  (a) a learnable per-token embedding (great for common modifications, which
+      have plenty of training data), and
+  (b) a chemistry feature derived from the modification's atomic composition
+      (lets rare / unseen modifications generalise — they inherit meaning from
+      chemically similar mods, and the scheme is open-vocab).
+
+A bare residue or special token has an all-zero composition delta; with the
+bias-free CompositionEncoder this maps to an exact zero contribution, so
+unmodified residues fall back to the pure token embedding.
+"""
+from __future__ import annotations
+
+import torch
+from torch import nn
+
+from peptide_property_ng.modifications.composition import CompositionTable
+
+
+class CompositionEncoder(nn.Module):
+    """Maps a signed atomic-composition vector to a d_model chemistry embedding.
+
+    Bias-free so that an all-zero composition (a bare residue) maps to an exact
+    zero vector — the additive-fusion semantics then mean "unmodified residue =
+    pure token embedding".
+    """
+
+    def __init__(self, n_elements: int, d_model: int, hidden: int | None = None):
+        super().__init__()
+        hidden = hidden or d_model // 2
+        self.net = nn.Sequential(
+            nn.Linear(n_elements, hidden, bias=False),
+            nn.GELU(),
+            nn.Linear(hidden, d_model, bias=False),
+        )
+
+    def forward(self, comp: torch.Tensor) -> torch.Tensor:  # (..., n_elements) -> (..., d_model)
+        return self.net(comp)
+
+
+class HybridResidueEmbedding(nn.Module):
+    """Token embedding + composition embedding, fused per residue.
+
+    Parameters
+    ----------
+    vocab_size, d_model, pad_token_id
+        Tokenizer / model dimensions.
+    composition_table
+        Per-token atomic-composition deltas (see modifications.composition).
+    fusion
+        ``"add"`` — ``token + composition`` (default).
+        ``"gate"`` — ``token + sigmoid(g) * composition`` with a learnable
+        per-channel gate ``g`` (init 0 → gate 0.5); the recommended upgrade if
+        the unseen-modification probe shows the composition branch is ignored.
+        ``"token_only"`` / ``"composition_only"`` — ablation modes that drop one
+        branch, for measuring each branch's contribution (e.g. the
+        unseen-modification probe's control is ``token_only``).
+    """
+
+    _FUSIONS = ("add", "gate", "token_only", "composition_only")
+
+    def __init__(
+        self,
+        vocab_size: int,
+        d_model: int,
+        composition_table: CompositionTable,
+        pad_token_id: int,
+        fusion: str = "add",
+    ):
+        super().__init__()
+        if fusion not in self._FUSIONS:
+            raise ValueError(f"unknown fusion '{fusion}'; expected one of {self._FUSIONS}")
+        self.d_model = d_model
+        self.pad_token_id = pad_token_id
+        self.fusion = fusion
+
+        # +1 row to match depthcharge's nn.Embedding(n_tokens + 1, ...) convention.
+        self.token_emb = nn.Embedding(vocab_size + 1, d_model, padding_idx=pad_token_id)
+        self.comp_encoder = CompositionEncoder(composition_table.n_elements, d_model)
+
+        # Composition lookup as a frozen buffer, indexed by token id; signed-log1p
+        # normalised. One extra zero row to align with token_emb's vocab_size + 1.
+        comp = composition_table.as_tensor(normalize=True)  # (vocab_size, n_elements)
+        comp = torch.cat([comp, torch.zeros(1, comp.shape[1])], dim=0)
+        self.register_buffer("composition", comp, persistent=False)
+
+        if fusion == "gate":
+            self.gate = nn.Parameter(torch.zeros(d_model))
+
+    def forward(self, tokens: torch.Tensor) -> torch.Tensor:
+        """tokens: (batch, seq_len) long  ->  (batch, seq_len, d_model)."""
+        if self.fusion == "token_only":
+            return self.token_emb(tokens)
+        comp_vec = self.comp_encoder(self.composition[tokens])
+        if self.fusion == "composition_only":
+            return comp_vec
+        token_vec = self.token_emb(tokens)
+        if self.fusion == "gate":
+            comp_vec = torch.sigmoid(self.gate) * comp_vec
+        return token_vec + comp_vec
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/encoder.py b/packages/peptide-property-ng/src/peptide_property_ng/model/encoder.py
new file mode 100644
index 000000000..227870bc9
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/encoder.py
@@ -0,0 +1,128 @@
+"""Depthcharge-based shared encoder for the unified peptide property predictor."""
+from __future__ import annotations
+
+import torch
+from depthcharge.transformers import AnalyteTransformerEncoder
+from torch import nn
+
+from peptide_property_ng.model.config import ModelConfig
+from peptide_property_ng.model.embedding import HybridResidueEmbedding
+from peptide_property_ng.modifications.composition import CompositionTable
+
+
+class PeptidePropertyEncoder(AnalyteTransformerEncoder):
+    """Shared transformer encoder.
+
+    Built on Depthcharge's ``AnalyteTransformerEncoder``, but:
+      - the plain token embedding is replaced with the hybrid (token +
+        atomic-composition) residue embedding, and
+      - conditioning on instrument model + acquisition mode happens via the
+        prepended global token. Charge / precursor m-z / collision energy are
+        deliberately *not* fed here — they are conditioned at the heads that use
+        them, so the charge head (which reads this shared output) never sees its
+        own label. One encoder pass, no label leakage.
+
+    ``forward(tokens, instrument, acq_mode) -> (latent, padding_mask)``
+      - ``latent``       ``(batch, 1 + seq_len, d_model)`` — index 0 is the global token
+      - ``padding_mask`` ``(batch, 1 + seq_len)`` bool — ``True`` where padded
+    """
+
+    def __init__(self, cfg: ModelConfig, composition_table: CompositionTable):
+        super().__init__(
+            n_tokens=cfg.vocab_size,
+            d_model=cfg.d_model,
+            nhead=cfg.n_heads,
+            dim_feedforward=cfg.dim_feedforward,
+            n_layers=cfg.n_layers,
+            dropout=cfg.dropout,
+            positional_encoder=True,
+            padding_int=cfg.pad_token_id,
+        )
+        self.cfg = cfg
+
+        # Fail loudly on a composition-table / config mismatch rather than
+        # silently indexing the wrong chemistry vectors.
+        if composition_table.vocab_size != cfg.vocab_size:
+            raise ValueError(
+                f"composition table vocab {composition_table.vocab_size} != "
+                f"cfg.vocab_size {cfg.vocab_size} — rebuild the table or fix the config"
+            )
+        if composition_table.n_elements != cfg.n_elements:
+            raise ValueError(
+                f"composition table has {composition_table.n_elements} elements != "
+                f"cfg.n_elements {cfg.n_elements}"
+            )
+
+        # The base class built a plain nn.Embedding token encoder; drop it and
+        # use the hybrid residue embedding instead.
+        del self.token_encoder
+        self.hybrid_embedding = HybridResidueEmbedding(
+            vocab_size=cfg.vocab_size,
+            d_model=cfg.d_model,
+            composition_table=composition_table,
+            pad_token_id=cfg.pad_token_id,
+            fusion=cfg.comp_fusion,
+        )
+
+        # Leak-free conditioning (instrument + acquisition mode only).
+        self.instrument_emb = nn.Embedding(cfg.n_instruments, cfg.d_model)
+        self.acq_emb = nn.Embedding(cfg.n_acq_modes, cfg.d_model)
+
+        # PyTorch's TransformerEncoder converts padded batches to a NestedTensor
+        # for the SDPA fast path; the backward pass on that path returns NaN
+        # gradients on certain shape combinations (silently observed as the
+        # `_nested_tensor_from_mask` warning we kept seeing — actually the smoke).
+        # Disable the fast path: slightly slower forward, finite gradients.
+        if hasattr(self.transformer_encoder, "enable_nested_tensor"):
+            self.transformer_encoder.enable_nested_tensor = False
+
+    def global_token_hook(  # noqa: D102 (documented on the class)
+        self,
+        tokens: torch.Tensor,
+        *,
+        instrument: torch.Tensor | None = None,
+        acq_mode: torch.Tensor | None = None,
+        **_: object,
+    ) -> torch.Tensor:
+        g = torch.zeros(
+            tokens.shape[0],
+            self.d_model,
+            device=tokens.device,
+            dtype=self.instrument_emb.weight.dtype,
+        )
+        if instrument is not None:
+            g = g + self.instrument_emb(instrument)
+        if acq_mode is not None:
+            g = g + self.acq_emb(acq_mode)
+        return g
+
+    def forward(  # type: ignore[override]
+        self,
+        tokens: torch.Tensor,
+        *,
+        instrument: torch.Tensor | None = None,
+        acq_mode: torch.Tensor | None = None,
+        attn_mask: torch.Tensor | None = None,
+    ) -> tuple[torch.Tensor, torch.Tensor]:
+        encoded = self.hybrid_embedding(tokens)  # (B, L, d)
+        global_token = self.global_token_hook(
+            tokens, instrument=instrument, acq_mode=acq_mode
+        )  # (B, d)
+        encoded = torch.cat([global_token[:, None, :], encoded], dim=1)  # (B, 1+L, d)
+
+        # Explicit padding mask from the token ids — robust, and the global
+        # token at position 0 is never padding.
+        pad = tokens == self.cfg.pad_token_id  # (B, L)
+        padding_mask = torch.cat(
+            [
+                torch.zeros(tokens.shape[0], 1, dtype=torch.bool, device=tokens.device),
+                pad,
+            ],
+            dim=1,
+        )  # (B, 1+L)
+
+        encoded = self.positional_encoder(encoded)
+        latent = self.transformer_encoder(
+            encoded, mask=attn_mask, src_key_padding_mask=padding_mask
+        )
+        return latent, padding_mask
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/heads/__init__.py b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/__init__.py
new file mode 100644
index 000000000..627540428
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/__init__.py
@@ -0,0 +1 @@
+"""Task heads — fragment-indexed intensity, CCS/IM, retention time, charge."""
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/heads/ccs.py b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/ccs.py
new file mode 100644
index 000000000..97c6002fa
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/ccs.py
@@ -0,0 +1,67 @@
+"""Ion-mobility head — a per-charge physics prior plus a learned correction.
+
+Predicts a min-max-normalised ion-mobility target in ~[0,1]. Both data sources
+are normalised to that same scale with fixed bounds (campaign data = Sage
+``ion_mobility`` 1/K0; pretraining data = ionmob CCS — see the data loaders),
+so the head sees one consistent target scale and no Mason-Schamp CCS<->1/K0
+conversion is needed.
+
+The physics prior ``1/K0 ~ slope[z]*sqrt(m/z) + intercept[z]`` is ported from
+the production ``SquareRootProjectionLayer`` (imspy_predictors); slopes and
+intercepts are learnable. A Gaussian-NLL ``(mean, std)`` output keeps a
+calibrated uncertainty, useful downstream for the pep-centric SNR work.
+"""
+from __future__ import annotations
+
+import torch
+from torch import nn
+from torch.nn import functional as F
+
+from peptide_property_ng.model.config import ModelConfig
+
+
+class SqrtMzProjection(nn.Module):
+    """Per-charge physics prior: ``property ~ slope[z]*sqrt(m/z) + intercept[z]``."""
+
+    def __init__(self, max_charge: int, slope_init: float = 0.0, intercept_init: float = 0.5):
+        super().__init__()
+        # Index by charge state directly (row 0 unused). Init at the centre of the
+        # normalised [0,1] target range; the per-charge slopes are learnable.
+        self.slopes = nn.Parameter(torch.full((max_charge + 1,), float(slope_init)))
+        self.intercepts = nn.Parameter(torch.full((max_charge + 1,), float(intercept_init)))
+
+    def forward(self, mz: torch.Tensor, charge: torch.Tensor) -> torch.Tensor:
+        z = charge.clamp(0, self.slopes.shape[0] - 1)
+        return self.slopes[z] * torch.sqrt(mz.clamp(min=0.0)) + self.intercepts[z]
+
+
+class IonMobilityHead(nn.Module):
+    """Predict inverse ion mobility (1/K0) as ``(mean, std)``."""
+
+    def __init__(self, cfg: ModelConfig):
+        super().__init__()
+        d = cfg.d_model
+        self.max_charge = cfg.max_charge
+        self.physics = SqrtMzProjection(cfg.max_charge)
+        self.charge_emb = nn.Embedding(cfg.max_charge + 1, d)
+        self.net = nn.Sequential(
+            nn.Linear(2 * d, d),
+            nn.LayerNorm(d),
+            nn.GELU(),
+            nn.Linear(d, d // 2),
+            nn.GELU(),
+            nn.Linear(d // 2, 2),  # correction to the mean, log-std
+        )
+
+    def forward(
+        self,
+        latent: torch.Tensor,        # (B, 1+L, d)
+        charge: torch.Tensor,        # (B,) long
+        precursor_mz: torch.Tensor,  # (B,) float
+    ) -> tuple[torch.Tensor, torch.Tensor]:
+        glob = latent[:, 0]                                       # (B, d)
+        h = torch.cat([glob, self.charge_emb(charge.clamp(0, self.max_charge))], dim=-1)
+        correction, log_std = self.net(h).unbind(dim=-1)          # (B,), (B,)
+        mean = self.physics(precursor_mz, charge) + correction
+        std = F.softplus(log_std) + 1e-4
+        return mean, std
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/heads/intensity.py b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/intensity.py
new file mode 100644
index 000000000..484329a0a
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/intensity.py
@@ -0,0 +1,59 @@
+"""Fragment-indexed intensity head.
+
+Unlike Prosit's fixed 174-value vector (29 positions x 6 channels -> a hard
+30-residue cap), this head predicts one intensity per *cleavage site*. A peptide
+of length L has L-1 sites; the site representation is built from the two
+flanking residue outputs, so the output is ``(batch, L-1, n_ion_channels)`` and
+the peptide length is bounded only by the encoder's ``max_seq_len``.
+
+Charge and collision energy are conditioned here (head-level) via FiLM — they
+are deliberately kept out of the shared encoder (see config.py).
+"""
+from __future__ import annotations
+
+import torch
+from torch import nn
+
+from peptide_property_ng.model.config import ModelConfig
+
+
+class IntensityHead(nn.Module):
+    """Predict per-cleavage-site fragment intensities in ``[0, 1]``."""
+
+    def __init__(self, cfg: ModelConfig):
+        super().__init__()
+        d = cfg.d_model
+        self.charge_emb = nn.Embedding(cfg.max_charge + 1, d)
+        self.ce_proj = nn.Linear(1, d)
+
+        self.site_in = nn.Linear(2 * d, d)
+        # FiLM modulation from (charge, CE); zero-init -> starts as identity.
+        self.film = nn.Linear(d, 2 * d)
+        nn.init.zeros_(self.film.weight)
+        nn.init.zeros_(self.film.bias)
+
+        self.mlp = nn.Sequential(
+            nn.LayerNorm(d),
+            nn.GELU(),
+            nn.Linear(d, d),
+            nn.GELU(),
+            nn.Linear(d, cfg.n_ion_channels),
+        )
+        self.max_charge = cfg.max_charge
+
+    def forward(
+        self,
+        latent: torch.Tensor,        # (B, 1+L, d) — index 0 is the global token
+        charge: torch.Tensor,        # (B,) long
+        collision_energy: torch.Tensor,  # (B,) float (0.0 when unknown)
+    ) -> torch.Tensor:
+        residues = latent[:, 1:]                       # (B, L, d)
+        sites = torch.cat([residues[:, :-1], residues[:, 1:]], dim=-1)  # (B, L-1, 2d)
+        h = self.site_in(sites)                        # (B, L-1, d)
+
+        cond = self.charge_emb(charge.clamp(0, self.max_charge))
+        cond = cond + self.ce_proj(collision_energy.unsqueeze(-1).float())  # (B, d)
+        scale, shift = self.film(cond).chunk(2, dim=-1)                      # (B, d) each
+        h = h * (1.0 + scale[:, None, :]) + shift[:, None, :]
+
+        return torch.sigmoid(self.mlp(h))              # (B, L-1, n_ion_channels)
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/heads/intensity_pooled.py b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/intensity_pooled.py
new file mode 100644
index 000000000..71d376733
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/intensity_pooled.py
@@ -0,0 +1,101 @@
+"""V4-compatible pooled intensity head — Prosit-174 ablation.
+
+Mirrors ``imspy_predictors.models.heads.IntensityHead``: attention-pool the
+residue latents, concat with charge / CE / instrument side-embeddings, MLP →
+sigmoid → reshape to ``(B, max_sites, n_ion_channels)``.
+
+The hypothesis (codex review): a global head that predicts the *whole* Prosit-174
+vector can model cross-fragment competition ("this peptide's base peak should
+dominate, suppress the rest"), which the per-cleavage-site FiLM head in
+``intensity.py`` cannot express except through encoder context. If swapping
+heads on the same pretrained encoder closes the remaining SA gap to v4 (~0.83),
+the gap is head topology, not capacity.
+"""
+from __future__ import annotations
+
+import torch
+from torch import nn
+from torch.nn import functional as F
+
+from peptide_property_ng.model.config import ModelConfig
+
+
+# Prosit-174 ceiling — fragments per peptide cap (29 = 30 aa - 1 cleavage site).
+# This is the public-corpus ground truth shape; even though ``cfg.max_seq_len``
+# allows longer peptides in the encoder, intensity supervision is bounded here.
+_MAX_SITES = 29
+
+
+class PooledIntensityHead(nn.Module):
+    """Predict ``(B, 29, n_ion_channels)`` via attention-pool + side-conditioned MLP."""
+
+    def __init__(self, cfg: ModelConfig):
+        super().__init__()
+        d = cfg.d_model
+        self.max_sites = _MAX_SITES
+        self.n_ion = cfg.n_ion_channels
+        self.num_outputs = self.max_sites * self.n_ion
+        self.max_charge = cfg.max_charge
+
+        # Single-head self-attention over the residue latents.
+        self.attention = nn.MultiheadAttention(
+            embed_dim=d, num_heads=1, batch_first=True, dropout=cfg.dropout,
+        )
+
+        # Side embeddings — match v4's d_model // 4 sizing.
+        side = max(d // 4, 16)
+        self.charge_embed = nn.Linear(cfg.max_charge + 1, side)
+        self.ce_embed = nn.Linear(1, side)
+        # Instrument conditioning lives in the head (v4 parity) on top of the
+        # encoder's global-token conditioning, not as a replacement.
+        self.instrument_embed = nn.Embedding(cfg.n_instruments, side)
+
+        in_dim = d + 3 * side
+        hidden = max(d * 2, 512)
+        self.fc = nn.Sequential(
+            nn.Linear(in_dim, hidden),
+            nn.LayerNorm(hidden),
+            nn.GELU(),
+            nn.Dropout(cfg.dropout),
+            nn.Linear(hidden, hidden),
+            nn.LayerNorm(hidden),
+            nn.GELU(),
+            nn.Dropout(cfg.dropout),
+            nn.Linear(hidden, self.num_outputs),
+        )
+
+    def forward(
+        self,
+        latent: torch.Tensor,            # (B, 1+L, d) — index 0 is the global token
+        charge: torch.Tensor,            # (B,) long
+        collision_energy: torch.Tensor,  # (B,) float
+        instrument: torch.Tensor | None = None,    # (B,) long
+        padding_mask: torch.Tensor | None = None,  # (B, L) bool, True at pad
+    ) -> torch.Tensor:
+        residues = latent[:, 1:]                       # (B, L, d)
+        B = residues.size(0)
+        device = residues.device
+
+        # Self-attention pool, then mean over valid residues.
+        pooled, _ = self.attention(
+            residues, residues, residues, key_padding_mask=padding_mask,
+        )
+        if padding_mask is not None:
+            valid = (~padding_mask).unsqueeze(-1).float()
+            denom = valid.sum(dim=1).clamp(min=1.0)
+            pooled = (pooled * valid).sum(dim=1) / denom    # (B, d)
+        else:
+            pooled = pooled.mean(dim=1)
+
+        # Side-info embeddings (v4 layout).
+        chg = charge.view(-1).long().clamp(0, self.max_charge)
+        chg_onehot = F.one_hot(chg, num_classes=self.max_charge + 1).float()
+        chg_emb = self.charge_embed(chg_onehot)
+        ce_emb = self.ce_embed(collision_energy.view(-1, 1).float())
+        if instrument is None:
+            instrument = torch.zeros(B, dtype=torch.long, device=device)
+        inst_emb = self.instrument_embed(instrument.view(-1).long())
+
+        features = torch.cat([pooled, chg_emb, ce_emb, inst_emb], dim=-1)
+        out = torch.sigmoid(self.fc(features))           # (B, max_sites * n_ion)
+        return out.view(B, self.max_sites, self.n_ion)
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/heads/scalar.py b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/scalar.py
new file mode 100644
index 000000000..3e735cc7f
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/heads/scalar.py
@@ -0,0 +1,67 @@
+"""Scalar heads — retention time and precursor charge.
+
+Both read only the encoder's global token (index 0). That token aggregates the
+peptide sequence plus instrument / acquisition-mode conditioning — but *not*
+charge, m/z or collision energy (those are kept out of the encoder). So the
+charge head never sees its own label: no leakage, one encoder pass.
+"""
+from __future__ import annotations
+
+import torch
+from torch import nn
+
+from peptide_property_ng.model.config import ModelConfig
+
+
+def _mlp(d_in: int, d_hidden: int, d_out: int) -> nn.Sequential:
+    return nn.Sequential(
+        nn.Linear(d_in, d_hidden),
+        nn.LayerNorm(d_hidden),
+        nn.GELU(),
+        nn.Linear(d_hidden, d_hidden // 2),
+        nn.GELU(),
+        nn.Linear(d_hidden // 2, d_out),
+    )
+
+
+class RetentionTimeHead(nn.Module):
+    """Predict (per-dataset min-max normalised) retention time."""
+
+    def __init__(self, cfg: ModelConfig):
+        super().__init__()
+        self.net = _mlp(cfg.d_model, cfg.d_model, 1)
+
+    def forward(self, latent: torch.Tensor) -> torch.Tensor:
+        return self.net(latent[:, 0]).squeeze(-1)  # (B,)
+
+
+class RTSigmaHead(nn.Module):
+    """Predict the RT elution-peak *shape* width (gradient-normalised EMG sigma).
+
+    A positive scalar per peptide, read from the encoder global token only (same
+    leakage-free footing as the RT-apex head). Output is the gradient-normalised
+    sigma; timsim rescales by the target gradient at simulation time (see the
+    SIM shape-predictor design). ``lambda`` (EMG tail) is deliberately *not*
+    predicted — it is unidentifiable at timsTOF MS1 frame sampling, so it stays
+    sampled by timsim. softplus keeps sigma strictly positive.
+    """
+
+    def __init__(self, cfg: ModelConfig):
+        super().__init__()
+        self.net = _mlp(cfg.d_model, cfg.d_model, 1)
+
+    def forward(self, latent: torch.Tensor) -> torch.Tensor:
+        raw = self.net(latent[:, 0]).squeeze(-1)  # (B,)
+        return nn.functional.softplus(raw) + 1e-4
+
+
+class ChargeHead(nn.Module):
+    """Predict precursor-charge logits — class ``j`` corresponds to charge ``j + 1``."""
+
+    def __init__(self, cfg: ModelConfig):
+        super().__init__()
+        self.max_charge = cfg.max_charge
+        self.net = _mlp(cfg.d_model, cfg.d_model, cfg.max_charge)
+
+    def forward(self, latent: torch.Tensor) -> torch.Tensor:
+        return self.net(latent[:, 0])  # (B, max_charge) logits
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/model/multitask.py b/packages/peptide-property-ng/src/peptide_property_ng/model/multitask.py
new file mode 100644
index 000000000..b70155d3e
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/model/multitask.py
@@ -0,0 +1,95 @@
+"""Unified multi-task peptide property model — one encoder, four heads."""
+from __future__ import annotations
+
+import torch
+from torch import nn
+
+from peptide_property_ng.model.config import SMALL, ModelConfig
+from peptide_property_ng.model.encoder import PeptidePropertyEncoder
+from peptide_property_ng.model.heads.ccs import IonMobilityHead
+from peptide_property_ng.model.heads.intensity import IntensityHead
+from peptide_property_ng.model.heads.intensity_pooled import PooledIntensityHead
+from peptide_property_ng.model.heads.scalar import ChargeHead, RetentionTimeHead, RTSigmaHead
+from peptide_property_ng.modifications.composition import CompositionTable
+
+
+class UnifiedPeptidePropertyModel(nn.Module):
+    """Shared Depthcharge-based encoder driving all property heads.
+
+    ``forward(batch)`` returns a dict over the active tasks:
+      - ``intensity`` : ``(B, L-1, n_ion_channels)`` in ``[0, 1]``
+      - ``ccs``       : ``(mean (B,), std (B,))`` — inverse ion mobility, 1/K0
+      - ``rt``        : ``(B,)`` — normalised retention time
+      - ``rt_sigma``  : ``(B,)`` — gradient-normalised RT EMG peak-width (timsim shape)
+      - ``charge``    : ``(B, max_charge)`` logits (class j -> charge j+1)
+
+    The input ``batch`` dict carries: ``tokens`` (B, L) long, ``charge`` (B,)
+    long, ``precursor_mz`` (B,) float, ``collision_energy`` (B,) float, and
+    optionally ``instrument`` / ``acq_mode`` (B,) long.
+    """
+
+    HEAD_TYPES: dict[str, type[nn.Module]] = {
+        "intensity": IntensityHead,
+        "ccs": IonMobilityHead,
+        "rt": RetentionTimeHead,
+        "rt_sigma": RTSigmaHead,
+        "charge": ChargeHead,
+    }
+
+    def __init__(
+        self,
+        cfg: ModelConfig = SMALL,
+        composition_table: CompositionTable | None = None,
+    ):
+        super().__init__()
+        self.cfg = cfg
+        if composition_table is None:
+            composition_table = CompositionTable.load()
+        self.encoder = PeptidePropertyEncoder(cfg, composition_table)
+        # Intensity head is dispatched on cfg.intensity_head ("site" | "pooled").
+        # Other heads always use the canonical class. Stored under one key
+        # ("intensity") regardless, so checkpoint code keys are stable.
+        head_types: dict[str, type[nn.Module]] = dict(self.HEAD_TYPES)
+        if cfg.intensity_head == "pooled":
+            head_types["intensity"] = PooledIntensityHead
+        elif cfg.intensity_head != "site":
+            raise ValueError(f"unknown intensity_head '{cfg.intensity_head}'")
+        self.heads = nn.ModuleDict({t: head_types[t](cfg) for t in cfg.tasks})
+
+    def forward(
+        self,
+        batch: dict[str, torch.Tensor],
+        tasks: list[str] | None = None,
+    ) -> dict[str, object]:
+        active = list(self.heads) if tasks is None else [t for t in tasks if t in self.heads]
+        latent, _ = self.encoder(
+            batch["tokens"],
+            instrument=batch.get("instrument"),
+            acq_mode=batch.get("acq_mode"),
+        )
+        out: dict[str, object] = {}
+        if "intensity" in active:
+            head = self.heads["intensity"]
+            if isinstance(head, PooledIntensityHead):
+                pad_mask = batch["tokens"] == self.cfg.pad_token_id
+                out["intensity"] = head(
+                    latent, batch["charge"], batch["collision_energy"],
+                    instrument=batch.get("instrument"),
+                    padding_mask=pad_mask,
+                )
+            else:
+                out["intensity"] = head(
+                    latent, batch["charge"], batch["collision_energy"]
+                )
+        if "ccs" in active:
+            out["ccs"] = self.heads["ccs"](latent, batch["charge"], batch["precursor_mz"])
+        if "rt" in active:
+            out["rt"] = self.heads["rt"](latent)
+        if "rt_sigma" in active:
+            out["rt_sigma"] = self.heads["rt_sigma"](latent)
+        if "charge" in active:
+            out["charge"] = self.heads["charge"](latent)
+        return out
+
+    def num_parameters(self) -> int:
+        return sum(p.numel() for p in self.parameters() if p.requires_grad)
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/modifications/__init__.py b/packages/peptide-property-ng/src/peptide_property_ng/modifications/__init__.py
new file mode 100644
index 000000000..8892dff88
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/modifications/__init__.py
@@ -0,0 +1 @@
+"""UNIMOD modification chemistry — atomic-composition features for the hybrid encoding."""
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/modifications/build_table.py b/packages/peptide-property-ng/src/peptide_property_ng/modifications/build_table.py
new file mode 100644
index 000000000..9509ab81e
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/modifications/build_table.py
@@ -0,0 +1,95 @@
+"""Build the UNIMOD atomic-composition lookup table for the hybrid encoding.
+
+Maps every ProformaTokenizer vocabulary token to a signed element-count vector
+and writes ``data/mod_composition_table.npz``:
+  - ``counts``   (vocab_size, n_elements) int16 — signed element deltas
+  - ``elements`` (n_elements,) str          — element order (== composition.ELEMENTS)
+  - ``tokens``   (vocab_size,) str          — id-indexed token strings (provenance)
+
+Run once; re-run when UNIMOD or the tokenizer vocabulary changes — no model
+vocab re-mint is needed, only a fresh table:
+
+    python -m peptide_property_ng.modifications.build_table
+"""
+from __future__ import annotations
+
+import re
+from pathlib import Path
+
+import numpy as np
+
+from peptide_property_ng.modifications.composition import DATA_DIR, DEFAULT_TABLE, ELEMENTS
+
+# Matches every UNIMOD reference inside a token, e.g. "C[UNIMOD:4]" -> "4".
+_UNIMOD_RE = re.compile(r"\[UNIMOD:(\d+)\]")
+
+
+def build(out_path: Path = DEFAULT_TABLE) -> dict:
+    """Build the composition table and write it to ``out_path``. Returns a summary."""
+    from imspy_predictors.utilities.tokenizers import ProformaTokenizer
+    from sagepy.core.unimod import modification_atomic_composition
+
+    tok = ProformaTokenizer.with_defaults()
+    vocab = list(tok.get_vocab())             # id-indexed list of token strings
+    comp = modification_atomic_composition()  # {'[UNIMOD:N]': {element: count}}
+
+    # Every element sagepy reports must be in our canonical alphabet.
+    seen = {e for c in comp.values() for e in c}
+    unknown = seen - set(ELEMENTS)
+    if unknown:
+        raise RuntimeError(
+            f"UNIMOD composition uses elements absent from ELEMENTS: {sorted(unknown)} "
+            "— add them to composition.ELEMENTS and rebuild."
+        )
+
+    elem_index = {e: i for i, e in enumerate(ELEMENTS)}
+    counts = np.zeros((len(vocab), len(ELEMENTS)), dtype=np.int16)
+
+    n_mod_tokens = n_resolved = 0
+    unresolved: set[str] = set()
+    for tid, token in enumerate(vocab):
+        ids = _UNIMOD_RE.findall(token)
+        if not ids:
+            continue  # bare residue or special token -> zero delta
+        n_mod_tokens += 1
+        resolved_any = False
+        for uid in ids:  # sum all mods on the token (terminal + residue, etc.)
+            key = f"[UNIMOD:{uid}]"
+            mod = comp.get(key)
+            if mod is None:
+                unresolved.add(key)
+                continue
+            resolved_any = True
+            for elem, c in mod.items():
+                counts[tid, elem_index[elem]] += c
+        n_resolved += int(resolved_any)
+
+    DATA_DIR.mkdir(parents=True, exist_ok=True)
+    np.savez_compressed(
+        out_path,
+        counts=counts,
+        elements=np.array(ELEMENTS),
+        tokens=np.array(vocab),
+    )
+
+    summary = {
+        "vocab_size": len(vocab),
+        "mod_tokens": n_mod_tokens,
+        "resolved": n_resolved,
+        "unresolved": len(unresolved),
+        "n_elements": len(ELEMENTS),
+        "nonzero_rows": int((counts != 0).any(axis=1).sum()),
+        "path": str(out_path),
+    }
+    return summary
+
+
+if __name__ == "__main__":
+    s = build()
+    print(f"vocab size          : {s['vocab_size']}")
+    print(f"modification tokens : {s['mod_tokens']}")
+    print(f"  with composition  : {s['resolved']}")
+    print(f"  unresolved UNIMOD : {s['unresolved']}  (id absent from sagepy table)")
+    print(f"elements            : {s['n_elements']}")
+    print(f"nonzero rows        : {s['nonzero_rows']}")
+    print(f"written             : {s['path']}")
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/modifications/composition.py b/packages/peptide-property-ng/src/peptide_property_ng/modifications/composition.py
new file mode 100644
index 000000000..4b39bad0f
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/modifications/composition.py
@@ -0,0 +1,86 @@
+"""UNIMOD atomic-composition features for the hybrid modification encoding.
+
+Each vocabulary token gets a fixed vector of element-count *deltas* relative to
+the bare residue: a modified-residue token (e.g. ``C[UNIMOD:4]``) carries the
+modification's atomic composition; a bare residue or a special token carries
+zeros. Deltas are *signed* — UNIMOD neutral-loss modifications have negative
+counts.
+
+Source: ``sagepy.core.unimod.modification_atomic_composition()`` — the same
+UNIMOD view the Sage search engine uses, so it is consistent with the training
+labels produced by Sage.
+"""
+from __future__ import annotations
+
+from pathlib import Path
+
+import numpy as np
+import torch
+
+# Canonical element order. Verified against sagepy's UNIMOD composition tables
+# (31 elements incl. stable isotopes and metals); ``build_table.py`` asserts the
+# live data is a subset of this set and fails loudly if UNIMOD ever adds one.
+ELEMENTS: tuple[str, ...] = (
+    "C", "H", "N", "O", "S", "P", "Se",             # common organic
+    "2H", "13C", "15N", "18O",                       # stable isotopes (labels)
+    "Ag", "Al", "As", "B", "Br", "Ca", "Cl", "Cu",   # metals / halogens / ...
+    "F", "Fe", "Hg", "I", "K", "Li", "Mg", "Mo",
+    "Na", "Ni", "Si", "Zn",
+)
+N_ELEMENTS = len(ELEMENTS)
+
+DATA_DIR = Path(__file__).parent / "data"
+DEFAULT_TABLE = DATA_DIR / "mod_composition_table.npz"
+
+
+def signed_log1p(x):
+    """Sign-preserving log compression: ``sign(x) * log1p(|x|)``.
+
+    Plain ``log1p`` is invalid for the ~430 UNIMOD modifications with negative
+    (neutral-loss) atom counts; this keeps the sign and compresses magnitude so
+    a 120-atom glycan does not dwarf a single-atom methylation.
+    """
+    if isinstance(x, torch.Tensor):
+        return torch.sign(x) * torch.log1p(torch.abs(x))
+    x = np.asarray(x, dtype=np.float32)
+    return np.sign(x) * np.log1p(np.abs(x))
+
+
+class CompositionTable:
+    """Per-token atomic-composition delta vectors, indexed by tokenizer id."""
+
+    def __init__(self, counts: np.ndarray, elements: tuple[str, ...]):
+        if counts.ndim != 2 or counts.shape[1] != len(elements):
+            raise ValueError(
+                f"counts shape {counts.shape} inconsistent with {len(elements)} elements"
+            )
+        self.counts = counts.astype(np.float32)  # (vocab_size, n_elements), signed
+        self.elements = tuple(elements)
+
+    @classmethod
+    def load(cls, path: str | Path = DEFAULT_TABLE) -> "CompositionTable":
+        path = Path(path)
+        if not path.exists():
+            raise FileNotFoundError(
+                f"composition table not found at {path}; run "
+                "`python -m peptide_property_ng.modifications.build_table` first."
+            )
+        with np.load(path, allow_pickle=False) as z:
+            return cls(z["counts"], tuple(str(e) for e in z["elements"]))
+
+    @property
+    def vocab_size(self) -> int:
+        return self.counts.shape[0]
+
+    @property
+    def n_elements(self) -> int:
+        return self.counts.shape[1]
+
+    def as_tensor(self, normalize: bool = True) -> torch.Tensor:
+        """Return the ``(vocab_size, n_elements)`` table as a float tensor.
+
+        With ``normalize`` (default) the signed-log1p compression is applied —
+        this is the form the model's ``CompositionEncoder`` consumes.
+        """
+        t = torch.from_numpy(self.counts).float()
+        return signed_log1p(t) if normalize else t
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/modifications/data/.gitkeep b/packages/peptide-property-ng/src/peptide_property_ng/modifications/data/.gitkeep
new file mode 100644
index 000000000..1758240c3
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/modifications/data/.gitkeep
@@ -0,0 +1 @@
+# Generated artifacts (mod_composition_table.npz) land here — see build_table.py
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/modifications/data/mod_composition_table.npz b/packages/peptide-property-ng/src/peptide_property_ng/modifications/data/mod_composition_table.npz
new file mode 100644
index 000000000..6f31aea44
Binary files /dev/null and b/packages/peptide-property-ng/src/peptide_property_ng/modifications/data/mod_composition_table.npz differ
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/train/__init__.py b/packages/peptide-property-ng/src/peptide_property_ng/train/__init__.py
new file mode 100644
index 000000000..99224b64a
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/train/__init__.py
@@ -0,0 +1 @@
+"""Training — multi-task optimisation."""
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/train/calibrate_ce.py b/packages/peptide-property-ng/src/peptide_property_ng/train/calibrate_ce.py
new file mode 100644
index 000000000..06ca203c5
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/train/calibrate_ce.py
@@ -0,0 +1,194 @@
+"""Compute per-PSM optimal CE with a frozen checkpoint and write it to a parquet.
+
+For each prepared PSM, sweep collision energy over a grid and pick the value
+that maximises spectral-angle similarity to the observed spectrum. The output
+is a parquet of ``(accession, psm_id, calibrated_ce)`` that ``train.train``
+can load via ``--ce-calibration`` to fine-tune with per-PSM CE instead of a
+flat default.
+
+Recipe: calibrate the *pretrained* (or other freeze-point) checkpoint once
+over the full prepared corpus; then fine-tune from that same checkpoint with
+the calibration applied. v4's effective recipe.
+
+Example:
+    python -m peptide_property_ng.train.calibrate_ce \\
+        --checkpoint /scratch/.../small-pretrained.pt \\
+        --datasets-glob '/scratch/claudius-proteomics/*' \\
+        --catalog .../timstof_catalog.tsv \\
+        --out /scratch/.../ce_calibration_pretrained_small.parquet
+"""
+from __future__ import annotations
+
+import argparse
+import time
+from pathlib import Path
+
+import pyarrow as pa
+import pyarrow.parquet as pq
+import torch
+from torch.utils.data import DataLoader
+
+from peptide_property_ng.data.collate import make_collate_fn
+from peptide_property_ng.data.sage_dataset import (
+    SagePropertyDataset,
+    build_split_datasets,
+    discover_sage_dirs,
+)
+from peptide_property_ng.losses import intensity_signal_mask, masked_spectral_angle
+from peptide_property_ng.model.config import PRESETS
+from peptide_property_ng.model.multitask import UnifiedPeptidePropertyModel
+from peptide_property_ng.modifications.composition import CompositionTable
+
+
+@torch.no_grad()
+def calibrate(model, loader, device: str, ce_grid: torch.Tensor) -> tuple[list, list, list, list]:
+    """Sweep the grid per batch, return per-example ``(acc, psm, ce, fixed_sa, opt_sa)``."""
+    model.eval()
+    ce_grid = ce_grid.to(device)
+    n_ce = ce_grid.numel()
+
+    accs: list[str] = []
+    pids: list[int] = []
+    opt_ces: list[float] = []
+    fixed_sas: list[float] = []
+    opt_sas: list[float] = []
+
+    for batch in loader:
+        meta_acc = batch.pop("_accession")
+        meta_pid = batch.pop("_psm_id")
+        batch = {k: v.to(device) for k, v in batch.items()}
+        target = batch["intensity_target"]
+        signal = intensity_signal_mask(target)
+        B = target.shape[0]
+
+        ce_orig = batch["collision_energy"].clone()
+        sa_matrix = torch.full((B, n_ce), float("nan"), device=device)
+        for j, ce in enumerate(ce_grid):
+            batch["collision_energy"] = ce.expand(B).float()
+            pred = model(batch, tasks=["intensity"])["intensity"]
+            sa_matrix[:, j] = 1.0 - masked_spectral_angle(pred, target)
+
+        # PSMs with no signal: keep their CE at the loader's default (no calibration);
+        # also mark fixed_sa / opt_sa as NaN so the parquet records "uncalibrated".
+        opt_sa, opt_idx = sa_matrix.max(dim=1)
+        opt_ce = ce_grid[opt_idx]
+        opt_ce = torch.where(signal, opt_ce, ce_orig)
+
+        batch["collision_energy"] = ce_orig
+        pred = model(batch, tasks=["intensity"])["intensity"]
+        fixed_sa = 1.0 - masked_spectral_angle(pred, target)
+
+        for i in range(B):
+            accs.append(meta_acc[i])
+            pids.append(int(meta_pid[i].item()))
+            opt_ces.append(float(opt_ce[i].item()))
+            if bool(signal[i]):
+                fixed_sas.append(float(fixed_sa[i].item()))
+                opt_sas.append(float(opt_sa[i].item()))
+            else:
+                fixed_sas.append(float("nan"))
+                opt_sas.append(float("nan"))
+    return accs, pids, opt_ces, fixed_sas, opt_sas
+
+
+def _collate_with_meta(pad_token_id: int, max_charge: int):
+    """Same as ``make_collate_fn`` but also returns per-example accession + psm_id."""
+    inner = make_collate_fn(pad_token_id, max_charge)
+
+    def _fn(examples: list[dict]) -> dict:
+        out = inner(examples)
+        out["_accession"] = [e["accession"] for e in examples]
+        out["_psm_id"] = torch.tensor([e["psm_id"] for e in examples], dtype=torch.int64)
+        return out
+
+    return _fn
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description=__doc__,
+                                 formatter_class=argparse.RawDescriptionHelpFormatter)
+    ap.add_argument("--checkpoint", required=True, help="checkpoint to calibrate against")
+    ap.add_argument("--datasets-glob", default="/scratch/claudius-proteomics/*")
+    ap.add_argument("--catalog", default=None)
+    ap.add_argument("--default-ce", type=float, default=0.26)
+    ap.add_argument("--cap", type=int, default=4000)
+    ap.add_argument("--max-datasets", type=int, default=0)
+    ap.add_argument("--batch-size", type=int, default=64)
+    ap.add_argument("--seed", type=int, default=0,
+                    help="MUST match what train.train uses (controls per-dataset PSM sampling)")
+    ap.add_argument("--ce-min", type=float, default=0.00)
+    ap.add_argument("--ce-max", type=float, default=0.60)
+    ap.add_argument("--ce-step", type=float, default=0.02)
+    ap.add_argument("--out", required=True, help="output parquet path")
+    ap.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu")
+    args = ap.parse_args()
+
+    out_path = Path(args.out)
+    out_path.parent.mkdir(parents=True, exist_ok=True)
+
+    sage_dirs = discover_sage_dirs(args.datasets_glob)
+    if args.max_datasets:
+        sage_dirs = sage_dirs[: args.max_datasets]
+    print(f"[{time.strftime('%H:%M:%S')}] {len(sage_dirs)} datasets (cap {args.cap}/dataset)")
+    splits = build_split_datasets(sage_dirs, cap=args.cap, seed=args.seed,
+                                  catalog_path=args.catalog, default_ce=args.default_ce)
+    examples = splits["train"].examples + splits["val"].examples + splits["test"].examples
+    print(f"  total {len(examples):,} prepared examples to calibrate "
+          f"({len(splits['train']):,} train + {len(splits['val']):,} val + {len(splits['test']):,} test)")
+
+    import dataclasses
+    ckpt = torch.load(args.checkpoint, map_location=args.device, weights_only=False)
+    preset = ckpt.get("preset", "small")
+    cfg = PRESETS[preset]
+    intensity_head = ckpt.get("intensity_head")
+    if intensity_head is None:
+        # Legacy checkpoint without metadata — infer from the state-dict keys.
+        keys = ckpt["model_state_dict"].keys()
+        if any("heads.intensity.attention." in k for k in keys):
+            intensity_head = "pooled"
+        else:
+            intensity_head = "site"
+    cfg = dataclasses.replace(cfg, intensity_head=intensity_head)
+    model = UnifiedPeptidePropertyModel(cfg, CompositionTable.load()).to(args.device)
+    model.load_state_dict(ckpt["model_state_dict"])
+    print(f"  checkpoint preset={preset}, intensity_head={intensity_head}, "
+          f"{model.num_parameters():,} parameters")
+
+    n_pts = int(round((args.ce_max - args.ce_min) / args.ce_step)) + 1
+    ce_grid = torch.linspace(args.ce_min, args.ce_max, n_pts)
+    print(f"  CE grid ({n_pts} pts): {args.ce_min:.3f} .. {args.ce_max:.3f} step {args.ce_step:.3f}")
+
+    loader = DataLoader(
+        SagePropertyDataset(examples), batch_size=args.batch_size, shuffle=False,
+        collate_fn=_collate_with_meta(cfg.pad_token_id, cfg.max_charge),
+    )
+    t0 = time.time()
+    accs, pids, opt_ces, fixed_sas, opt_sas = calibrate(model, loader, args.device, ce_grid)
+    print(f"[{time.strftime('%H:%M:%S')}] calibrated {len(accs):,} PSMs in {time.time() - t0:.0f}s")
+
+    table = pa.table({
+        "accession": pa.array(accs, type=pa.string()),
+        "psm_id": pa.array(pids, type=pa.int64()),
+        "calibrated_ce": pa.array(opt_ces, type=pa.float32()),
+        "fixed_sa": pa.array(fixed_sas, type=pa.float32()),
+        "opt_sa": pa.array(opt_sas, type=pa.float32()),
+    })
+    pq.write_table(table, out_path)
+    print(f"  wrote {out_path}  ({out_path.stat().st_size / 1e6:.2f} MB)")
+
+    # Quick summary of the calibration so the user can see the lift.
+    import numpy as np
+    fa = np.array(fixed_sas, dtype=np.float64)
+    oa = np.array(opt_sas, dtype=np.float64)
+    valid = ~np.isnan(fa) & ~np.isnan(oa)
+    print(f"\nover the {int(valid.sum()):,} signal PSMs:")
+    print(f"  fixed (per-PSM default/catalog CE)  mean SA  {fa[valid].mean():.4f}")
+    print(f"  per-PSM optimal CE                  mean SA  {oa[valid].mean():.4f}  "
+          f"(lift +{(oa[valid] - fa[valid]).mean():.4f})")
+    ce_arr = np.array(opt_ces, dtype=np.float64)[valid]
+    print(f"  CE: mean {ce_arr.mean():.3f}  median {np.median(ce_arr):.3f}  "
+          f"q05 {np.quantile(ce_arr, 0.05):.3f}  q95 {np.quantile(ce_arr, 0.95):.3f}")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/train/pretrain.py b/packages/peptide-property-ng/src/peptide_property_ng/train/pretrain.py
new file mode 100644
index 000000000..c9d106964
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/train/pretrain.py
@@ -0,0 +1,199 @@
+"""Staged pretraining — warm the shared encoder + heads on large public corpora,
+then hand off to ``train.py`` for the multi-task campaign fine-tune.
+
+Default curriculum (each stage is single-task; the shared encoder accumulates):
+  1. intensity  — Wilhelmlab/prospect-ptms-ms2  (Orbitrap-HCD, broad)
+  2. intensity  — Wilhelmlab/timsTOF-ms2         (timsTOF — instrument domain)
+  3. ccs        — ionmob CCS parquets            (timsTOF)
+  4. rt         — Chronologer DB                 (harmonised retention time)
+
+    python -m peptide_property_ng.train.pretrain --cap 50000 --epochs 3
+    python -m peptide_property_ng.train.train --init-from runs/pretrain/pretrained.pt \\
+        --datasets-glob '/scratch/claudius-proteomics/*'
+"""
+from __future__ import annotations
+
+import argparse
+import dataclasses
+import json
+import time
+from pathlib import Path
+
+import torch
+from torch.utils.data import DataLoader
+
+from peptide_property_ng.data.ccs_pretrain import prepare_ccs_examples
+from peptide_property_ng.data.chronologer_rt import prepare_chronologer_examples
+from peptide_property_ng.data.collate import make_collate_fn
+from peptide_property_ng.data.hf_intensity import prepare_hf_intensity_examples
+from peptide_property_ng.data.sage_dataset import SagePropertyDataset
+from peptide_property_ng.data.splits import peptide_split
+from peptide_property_ng.eval.metrics import evaluate_split
+from peptide_property_ng.losses import MultiTaskLoss
+from peptide_property_ng.model.config import PRESETS, instrument_id
+from peptide_property_ng.model.multitask import UnifiedPeptidePropertyModel
+from peptide_property_ng.modifications.composition import CompositionTable
+
+_STAGE_METRIC = {"intensity": "intensity_sa", "ccs": "ccs_mae", "rt": "rt_mae"}
+
+
+def _default_stages(
+    chronologer_db: str | None = None, ccs_glob: str | None = None
+) -> list[dict]:
+    """Curriculum: each stage = name, task, and a builder taking a per-stage cap.
+
+    The order is deliberate: the auxiliary tasks (RT, CCS) run first and
+    intensity runs last, so the shared encoder ends tuned for the priority task
+    rather than left drifted toward RT — a catastrophic-forgetting mitigation.
+    Within intensity, the broad Orbitrap corpora precede the timsTOF one.
+    """
+    rt_kw = {} if chronologer_db is None else {"db_path": chronologer_db}
+    ccs_kw = {} if ccs_glob is None else {"data_glob": ccs_glob}
+    return [
+        {
+            "name": "rt-chronologer", "task": "rt",
+            "build": lambda cap: prepare_chronologer_examples(cap=cap, **rt_kw),
+        },
+        {
+            "name": "ccs-ionmob", "task": "ccs",
+            "build": lambda cap: prepare_ccs_examples(
+                cap=cap, instrument=instrument_id("timsTOF Pro"), **ccs_kw),
+        },
+        {
+            "name": "intensity-prospect", "task": "intensity",
+            "build": lambda cap: prepare_hf_intensity_examples(
+                "Wilhelmlab/prospect-ptms-ms2", split="train", cap=cap,
+                instrument=instrument_id("orbitrap")),
+        },
+        {
+            "name": "intensity-prosit2025", "task": "intensity",
+            "build": lambda cap: prepare_hf_intensity_examples(
+                "Wilhelmlab/Prosit-2025-lac-ms2", split="train", cap=cap,
+                instrument=instrument_id("orbitrap")),
+        },
+        {
+            "name": "intensity-timstof", "task": "intensity",
+            "build": lambda cap: prepare_hf_intensity_examples(
+                "Wilhelmlab/timsTOF-ms2", split="train", cap=cap,
+                instrument=instrument_id("timsTOF Pro")),
+        },
+    ]
+
+
+def _train_task_epoch(model, loader, task, loss_fn, optimizer, device) -> float:
+    model.train()
+    total, n_batches = 0.0, 0
+    for batch in loader:
+        batch = {k: v.to(device) for k, v in batch.items()}
+        optimizer.zero_grad()
+        loss, _ = loss_fn(model(batch, tasks=[task]), batch)
+        loss.backward()
+        torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)
+        optimizer.step()
+        total += float(loss.detach())
+        n_batches += 1
+    return total / max(n_batches, 1)
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description="Staged pretraining for the unified model.")
+    ap.add_argument("--preset", default="small", choices=sorted(PRESETS))
+    ap.add_argument("--intensity-head", default=None,
+                    choices=["site", "pooled"],
+                    help="intensity head topology (default: preset value); "
+                         "use 'pooled' to pretrain the v4-style Prosit-174 head")
+    ap.add_argument("--cap", type=int, default=50000, help="max examples per stage (0 = all)")
+    ap.add_argument("--epochs", type=int, default=3, help="epochs per stage")
+    ap.add_argument("--batch-size", type=int, default=128)
+    ap.add_argument("--lr", type=float, default=3e-4)
+    ap.add_argument("--val-frac", type=float, default=0.03)
+    ap.add_argument("--seed", type=int, default=0)
+    ap.add_argument("--stages", default="all", help="comma-separated stage names, or 'all'")
+    ap.add_argument("--chronologer-db", default=None,
+                    help="path to Chronologer_DB_*.gz (default: the package's known local path)")
+    ap.add_argument("--ccs-glob", default=None,
+                    help="glob for ionmob *_unique_unimod*.parquet (default: known local path)")
+    ap.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu")
+    ap.add_argument("--out", default="runs/pretrain")
+    args = ap.parse_args()
+
+    torch.manual_seed(args.seed)
+    out_dir = Path(args.out)
+    out_dir.mkdir(parents=True, exist_ok=True)
+    cfg = PRESETS[args.preset]
+    if args.intensity_head:
+        cfg = dataclasses.replace(cfg, intensity_head=args.intensity_head)
+    device = args.device
+    cap = None if args.cap == 0 else args.cap
+
+    model = UnifiedPeptidePropertyModel(cfg, CompositionTable.load()).to(device)
+    loss_fn = MultiTaskLoss()
+    collate = make_collate_fn(cfg.pad_token_id, cfg.max_charge)
+    print(f"model: '{args.preset}' preset, {model.num_parameters():,} parameters, device={device}")
+
+    stages = _default_stages(args.chronologer_db, args.ccs_glob)
+    if args.stages != "all":
+        wanted = {s.strip() for s in args.stages.split(",")}
+        stages = [s for s in stages if s["name"] in wanted]
+
+    history: list[dict] = []
+    for stage in stages:
+        name, task = stage["name"], stage["task"]
+        t0 = time.time()
+        print(f"\n=== stage '{name}' (task={task}) ===", flush=True)
+        examples = stage["build"](cap)
+        if len(examples) < args.batch_size * 2:
+            print(f"  only {len(examples)} examples — skipping", flush=True)
+            continue
+
+        # peptide-level hold-out — a peptide never crosses train/val within a
+        # stage (stage metrics are pretraining diagnostics, not the final eval).
+        train_ex, val_ex = [], []
+        for ex in examples:
+            bucket = peptide_split(ex["stripped"], val_frac=args.val_frac,
+                                   test_frac=0.0, seed=args.seed)
+            (val_ex if bucket == "val" else train_ex).append(ex)
+        print(f"  {len(train_ex):,} train / {len(val_ex):,} val  [prepared in {time.time()-t0:.0f}s]",
+              flush=True)
+
+        train_loader = DataLoader(SagePropertyDataset(train_ex), batch_size=args.batch_size,
+                                  shuffle=True, collate_fn=collate, drop_last=True)
+        val_loader = DataLoader(SagePropertyDataset(val_ex), batch_size=args.batch_size,
+                                collate_fn=collate)
+        optimizer = torch.optim.AdamW(model.parameters(), lr=args.lr, weight_decay=0.01)
+        metric = _STAGE_METRIC[task]
+        for epoch in range(1, args.epochs + 1):
+            te = time.time()
+            train_loss = _train_task_epoch(model, train_loader, task, loss_fn, optimizer, device)
+            # Restrict eval to the current stage's task: other heads' outputs are
+            # at shapes the current corpus's targets/placeholders don't share
+            # (eg. pooled-head intensity (B,29,6) vs Chronologer's (B,L-1,6) for
+            # L up to max_seq_len), and per-stage metrics are diagnostic anyway.
+            val = evaluate_split(model, val_loader, device, tasks=[task])
+            print(f"  epoch {epoch}: train_loss={train_loss:.4f}  "
+                  f"val[{metric}]={val[metric]:.4f}  [{time.time()-te:.0f}s]", flush=True)
+            history.append({"stage": name, "epoch": epoch,
+                             "train_loss": train_loss, "val": val})
+
+        # Save a per-stage checkpoint (so the campaign fine-tune can pick the
+        # handoff point, not just the post-RT state) and the running pretrained.pt.
+        ckpt = {
+            "model_state_dict": model.state_dict(),
+            "preset": args.preset,
+            "intensity_head": cfg.intensity_head,
+            "stage": name,
+        }
+        torch.save(ckpt, out_dir / f"after-{name}.pt")
+        torch.save(ckpt, out_dir / "pretrained.pt")
+        print(f"  saved -> {out_dir}/after-{name}.pt  (+ pretrained.pt)", flush=True)
+
+    (out_dir / "pretrain_history.json").write_text(json.dumps(history, indent=2))
+    print(
+        f"\npretraining complete. fine-tune on the campaign data with:\n"
+        f"  python -m peptide_property_ng.train.train "
+        f"--init-from {out_dir}/pretrained.pt --datasets-glob '/scratch/claudius-proteomics/*'"
+    )
+
+
+if __name__ == "__main__":
+    main()
diff --git a/packages/peptide-property-ng/src/peptide_property_ng/train/train.py b/packages/peptide-property-ng/src/peptide_property_ng/train/train.py
new file mode 100644
index 000000000..d1ec08c6d
--- /dev/null
+++ b/packages/peptide-property-ng/src/peptide_property_ng/train/train.py
@@ -0,0 +1,234 @@
+"""Train the unified peptide property model — first prototype.
+
+Example:
+    python -m peptide_property_ng.train.train \\
+        --datasets-glob '/scratch/claudius-proteomics/*' --cap 4000 --epochs 12
+"""
+from __future__ import annotations
+
+import argparse
+import dataclasses
+import json
+import time
+from pathlib import Path
+
+import torch
+from torch.utils.data import DataLoader
+
+from peptide_property_ng.data.collate import make_collate_fn
+from peptide_property_ng.data.sage_dataset import (
+    build_split_datasets, discover_sage_dirs, load_ce_calibration,
+)
+from peptide_property_ng.eval.metrics import evaluate_split
+from peptide_property_ng.losses import MultiTaskLoss
+from peptide_property_ng.model.config import PRESETS
+from peptide_property_ng.model.multitask import UnifiedPeptidePropertyModel
+from peptide_property_ng.modifications.composition import CompositionTable
+
+
+def _to_device(batch: dict, device: str) -> dict:
+    return {k: v.to(device) for k, v in batch.items()}
+
+
+def train_epoch(model, loader, loss_fn, optimizer, device, tasks=None) -> dict[str, float]:
+    model.train()
+    totals: dict[str, float] = {}
+    n_batches = 0
+    for batch in loader:
+        batch = _to_device(batch, device)
+        optimizer.zero_grad()
+        loss, parts = loss_fn(model(batch, tasks=tasks), batch)
+        loss.backward()
+        torch.nn.utils.clip_grad_norm_(model.parameters(), 1.0)
+        optimizer.step()
+        for name, value in parts.items():
+            totals[name] = totals.get(name, 0.0) + value
+        totals["_total"] = totals.get("_total", 0.0) + float(loss.detach())
+        n_batches += 1
+    return {k: v / max(n_batches, 1) for k, v in totals.items()}
+
+
+def main() -> None:
+    ap = argparse.ArgumentParser(description="Train the unified peptide property model.")
+    ap.add_argument("--datasets-glob", default="/scratch/claudius-proteomics/*")
+    ap.add_argument("--hf-corpus", default=None,
+                    help="path to aggregated HF corpus dir (tier1/ + tier3/ with "
+                         "{train,val,test}.parquet); uses the HF loader instead of Sage")
+    ap.add_argument("--hf-max-datasets", type=int, default=0,
+                    help="limit the HF run to the first N accessions (0 = all)")
+    ap.add_argument("--hf-rt-lookup", default=None,
+                    help="aligned_rt lookup parquet (accession,psm_id,aligned_rt); "
+                         "uses Sage cross-run-aligned RT instead of raw rt_seconds")
+    ap.add_argument("--catalog", default=None,
+                    help="path to timstof_catalog.tsv for per-dataset instrument/acq conditioning")
+    ap.add_argument("--default-ce", type=float, default=0.26,
+                    help="default normalized collision energy for Sage PSMs that lack one "
+                         "(~0.26 = timsTOF-ms2 median, in-distribution for the pretrained CE FiLM)")
+    ap.add_argument("--ce-calibration", default=None,
+                    help="parquet of per-PSM calibrated CEs (from calibrate_ce.py); "
+                         "overrides default_ce per (accession, psm_id) where present")
+    ap.add_argument("--preset", default="small", choices=sorted(PRESETS))
+    ap.add_argument("--comp-fusion", default=None,
+                    choices=["add", "gate", "token_only", "composition_only"],
+                    help="modification-encoding fusion / ablation (default: preset value)")
+    ap.add_argument("--intensity-head", default=None,
+                    choices=["site", "pooled"],
+                    help="intensity head topology: 'site' (local FiLM, default) or "
+                         "'pooled' (v4-style Prosit-174 attention-pooled ablation)")
+    ap.add_argument("--cap", type=int, default=4000, help="max PSMs sampled per dataset")
+    ap.add_argument("--max-datasets", type=int, default=0, help="limit datasets (0 = all)")
+    ap.add_argument("--epochs", type=int, default=12)
+    ap.add_argument("--batch-size", type=int, default=64)
+    ap.add_argument("--lr", type=float, default=3e-4)
+    ap.add_argument("--weight-decay", type=float, default=0.01)
+    ap.add_argument("--seed", type=int, default=0)
+    ap.add_argument("--patience", type=int, default=4)
+    ap.add_argument("--init-from", default=None,
+                    help="checkpoint to initialise weights from (e.g. a pretrained.pt)")
+    ap.add_argument("--tasks", default="intensity,ccs,rt,charge",
+                    help="comma-separated tasks to train + evaluate; restricts forward and loss "
+                         "(e.g. 'intensity' for a specialized fine-tune)")
+    ap.add_argument("--freeze-encoder", action="store_true",
+                    help="freeze the shared encoder (only the task heads train) -- linear-probe / "
+                         "feature-extraction mode on the pretrained encoder")
+    ap.add_argument("--device", default="cuda" if torch.cuda.is_available() else "cpu")
+    ap.add_argument("--out", default="runs/prototype")
+    args = ap.parse_args()
+
+    torch.manual_seed(args.seed)
+    out_dir = Path(args.out)
+    out_dir.mkdir(parents=True, exist_ok=True)
+
+    print(f"[{time.strftime('%H:%M:%S')}] discovering datasets ...", flush=True)
+    t0 = time.time()
+    if args.hf_corpus:
+        import pyarrow.compute as _pc
+        import pyarrow.dataset as _pds
+        from peptide_property_ng.data.hf_corpus_dataset import build_split_datasets_hf
+        accs = None
+        if args.hf_max_datasets:
+            all_accs = sorted(_pc.unique(_pds.dataset(
+                f"{args.hf_corpus}/tier1/val.parquet").to_table(
+                columns=["accession"])["accession"]).to_pylist())
+            accs = all_accs[: args.hf_max_datasets]
+        print(f"[{time.strftime('%H:%M:%S')}] HF corpus {args.hf_corpus} "
+              f"(cap {args.cap}/dataset, {len(accs) if accs else 'all'} datasets) ...", flush=True)
+        splits = build_split_datasets_hf(args.hf_corpus, cap=args.cap,
+                                         seed=args.seed, accessions=accs,
+                                         rt_lookup=args.hf_rt_lookup)
+    else:
+        sage_dirs = discover_sage_dirs(args.datasets_glob)
+        if args.max_datasets:
+            sage_dirs = sage_dirs[: args.max_datasets]
+        print(f"  {len(sage_dirs)} datasets")
+        print(f"[{time.strftime('%H:%M:%S')}] preparing examples (cap {args.cap}/dataset) ...", flush=True)
+        ce_cal = load_ce_calibration(args.ce_calibration) if args.ce_calibration else None
+        if ce_cal is not None:
+            print(f"  CE calibration loaded: {len(ce_cal):,} keys from {args.ce_calibration}")
+        splits = build_split_datasets(sage_dirs, cap=args.cap, seed=args.seed,
+                                      catalog_path=args.catalog, default_ce=args.default_ce,
+                                      ce_calibration=ce_cal)
+    for name, ds in splits.items():
+        print(f"  {name}: {len(ds):,} examples")
+    print(f"  prepared in {time.time() - t0:.0f}s")
+    if len(splits["train"]) == 0:
+        raise SystemExit("no training examples — check --datasets-glob")
+
+    cfg = PRESETS[args.preset]
+    if args.comp_fusion:
+        cfg = dataclasses.replace(cfg, comp_fusion=args.comp_fusion)
+    if args.intensity_head:
+        cfg = dataclasses.replace(cfg, intensity_head=args.intensity_head)
+    collate = make_collate_fn(cfg.pad_token_id, cfg.max_charge)
+    loaders = {
+        name: DataLoader(
+            ds, batch_size=args.batch_size, shuffle=(name == "train"),
+            collate_fn=collate, drop_last=(name == "train"),
+        )
+        for name, ds in splits.items()
+    }
+
+    device = args.device
+    model = UnifiedPeptidePropertyModel(cfg, CompositionTable.load()).to(device)
+    if args.init_from:
+        ckpt = torch.load(args.init_from, map_location=device)
+        ckpt_preset = ckpt.get("preset")
+        if ckpt_preset is not None and ckpt_preset != args.preset:
+            raise SystemExit(
+                f"--init-from checkpoint is preset '{ckpt_preset}', but --preset is "
+                f"'{args.preset}'; they must match (the architectures differ)."
+            )
+        missing, unexpected = model.load_state_dict(ckpt["model_state_dict"], strict=False)
+        print(f"initialised from {args.init_from} "
+              f"(missing={len(missing)}, unexpected={len(unexpected)})")
+
+    task_list = [t.strip() for t in args.tasks.split(",") if t.strip()]
+    if args.freeze_encoder:
+        for p in model.encoder.parameters():
+            p.requires_grad = False
+        n_frozen = sum(p.numel() for p in model.encoder.parameters())
+        print(f"encoder frozen: {n_frozen:,} params no longer training")
+
+    n_train = sum(p.numel() for p in model.parameters() if p.requires_grad)
+    print(f"model: '{args.preset}' preset, {model.num_parameters():,} total parameters, "
+          f"{n_train:,} trainable, device={device}, tasks={task_list}")
+    optimizer = torch.optim.AdamW([p for p in model.parameters() if p.requires_grad],
+                                  lr=args.lr, weight_decay=args.weight_decay)
+    loss_fn = MultiTaskLoss()
+
+    best_sa, best_epoch, bad = -1.0, -1, 0
+    history: list[dict] = []
+    for epoch in range(1, args.epochs + 1):
+        t0 = time.time()
+        tr = train_epoch(model, loaders["train"], loss_fn, optimizer, device, tasks=task_list)
+        val = evaluate_split(model, loaders["val"], device, tasks=task_list)
+        history.append({"epoch": epoch, "train": tr, "val": val})
+        print(
+            f"[{time.strftime('%H:%M:%S')}] epoch {epoch:2d}  "
+            f"train_loss={tr.get('_total', 0):.4f} "
+            f"(int={tr.get('intensity', 0):.3f} ccs={tr.get('ccs', 0):.3f} "
+            f"rt={tr.get('rt', 0):.3f} chg={tr.get('charge', 0):.3f})  "
+            f"val: SA={val['intensity_sa']:.4f} ccs_mae={val['ccs_mae']:.4f} "
+            f"rt_mae={val['rt_mae']:.4f} chg_acc={val['charge_acc']:.3f}  "
+            f"[{time.time() - t0:.0f}s]",
+            flush=True,
+        )
+        if val["intensity_sa"] > best_sa:
+            best_sa, best_epoch, bad = val["intensity_sa"], epoch, 0
+            torch.save(
+                {"model_state_dict": model.state_dict(), "preset": args.preset,
+                 "intensity_head": cfg.intensity_head,
+                 "epoch": epoch, "val": val},
+                out_dir / "best.pt",
+            )
+        else:
+            bad += 1
+            if bad >= args.patience:
+                print(f"early stop (no val SA gain for {args.patience} epochs)")
+                break
+
+    if best_epoch == -1:
+        print("warning: no epoch improved val SA; saving the final model as best.pt")
+        torch.save(
+            {"model_state_dict": model.state_dict(), "preset": args.preset,
+             "epoch": epoch, "val": val},
+            out_dir / "best.pt",
+        )
+        best_epoch = epoch
+    ckpt = torch.load(out_dir / "best.pt", map_location=device)
+    model.load_state_dict(ckpt["model_state_dict"])
+    test = evaluate_split(model, loaders["test"], device, tasks=task_list)
+    print(f"\n=== test (best epoch {best_epoch}) ===")
+    for k, v in test.items():
+        print(f"  {k:16s} {v:.4f}")
+    (out_dir / "metrics.json").write_text(
+        json.dumps(
+            {"args": vars(args), "best_epoch": best_epoch, "history": history, "test": test},
+            indent=2,
+        )
+    )
+    print(f"\nsaved -> {out_dir}/best.pt , {out_dir}/metrics.json")
+
+
+if __name__ == "__main__":
+    main()
diff --git a/packages/peptide-property-ng/tests/__init__.py b/packages/peptide-property-ng/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/packages/peptide-property-ng/tests/test_composition.py b/packages/peptide-property-ng/tests/test_composition.py
new file mode 100644
index 000000000..f0e98a67c
--- /dev/null
+++ b/packages/peptide-property-ng/tests/test_composition.py
@@ -0,0 +1,59 @@
+"""Tests for the UNIMOD atomic-composition table."""
+import numpy as np
+import torch
+
+from peptide_property_ng.modifications.composition import (
+    DEFAULT_TABLE,
+    ELEMENTS,
+    CompositionTable,
+    signed_log1p,
+)
+
+
+def _tokens() -> list[str]:
+    """Id-indexed token strings, read straight from the built table."""
+    with np.load(DEFAULT_TABLE, allow_pickle=False) as z:
+        return [str(t) for t in z["tokens"]]
+
+
+def test_table_loads():
+    t = CompositionTable.load()
+    assert t.vocab_size == 2175
+    assert t.n_elements == len(ELEMENTS) == 31
+    assert t.elements == ELEMENTS
+
+
+def test_known_modifications():
+    """Carbamidomethyl (UNIMOD:4, +C2H3NO) and Oxidation (UNIMOD:35, +O)."""
+    t = CompositionTable.load()
+    vocab = _tokens()
+    ei = {e: i for i, e in enumerate(t.elements)}
+
+    cam = next(i for i, v in enumerate(vocab) if "[UNIMOD:4]" in v)
+    row = t.counts[cam]
+    assert (row[ei["C"]], row[ei["H"]], row[ei["N"]], row[ei["O"]]) == (2, 3, 1, 1)
+
+    ox = next(i for i, v in enumerate(vocab) if "[UNIMOD:35]" in v)
+    assert t.counts[ox][ei["O"]] == 1
+
+
+def test_bare_residue_is_zero():
+    """A bare amino-acid token carries an all-zero composition delta."""
+    t = CompositionTable.load()
+    ala = _tokens().index("A")
+    assert not t.counts[ala].any()
+
+
+def test_signed_log1p_handles_negative_losses():
+    """Neutral-loss modifications have negative counts; signed_log1p must stay finite."""
+    x = torch.tensor([-40.0, -1.0, 0.0, 3.0, 120.0])
+    out = signed_log1p(x)
+    assert torch.isfinite(out).all()
+    assert out[0] < 0 and out[2] == 0 and out[4] > 0
+    # numpy path agrees
+    assert np.allclose(signed_log1p(x.numpy()), out.numpy())
+
+
+def test_table_has_negative_entries():
+    """The table must preserve neutral-loss (negative) compositions."""
+    assert (CompositionTable.load().counts < 0).any()
diff --git a/packages/peptide-property-ng/tests/test_encoder.py b/packages/peptide-property-ng/tests/test_encoder.py
new file mode 100644
index 000000000..f8c181430
--- /dev/null
+++ b/packages/peptide-property-ng/tests/test_encoder.py
@@ -0,0 +1,69 @@
+"""Tests for the hybrid embedding and the Depthcharge-based encoder."""
+import torch
+
+from peptide_property_ng.model.config import SMALL
+from peptide_property_ng.model.encoder import PeptidePropertyEncoder
+from peptide_property_ng.modifications.composition import CompositionTable
+
+# Random token ids in 1..59 — safely valid and never the pad id (61).
+_LO, _HI = 1, 60
+
+
+def _encoder() -> PeptidePropertyEncoder:
+    return PeptidePropertyEncoder(SMALL, CompositionTable.load())
+
+
+def test_forward_shape():
+    enc = _encoder()
+    b, length = 4, 12
+    tokens = torch.randint(_LO, _HI, (b, length))
+    latent, mask = enc(tokens)
+    assert latent.shape == (b, 1 + length, SMALL.d_model)
+    assert mask.shape == (b, 1 + length)
+    assert not mask[:, 0].any()  # global token is never masked
+
+
+def test_padding_mask():
+    enc = _encoder()
+    tokens = torch.randint(_LO, _HI, (2, 8))
+    tokens[0, 5:] = SMALL.pad_token_id  # pad the last 3 residues of sample 0
+    _, mask = enc(tokens)
+    assert mask[0, 6:].all()       # +1 offset for the global token
+    assert not mask[0, :6].any()
+    assert not mask[1].any()
+
+
+def test_instrument_conditioning_changes_output():
+    enc = _encoder().eval()
+    tokens = torch.randint(_LO, _HI, (3, 10))
+    base, _ = enc(tokens)
+    cond, _ = enc(tokens, instrument=torch.tensor([2, 2, 2]))
+    assert not torch.allclose(base, cond)
+
+
+def test_long_peptide_runs():
+    """A 50-residue peptide must encode fine — no fixed-length cap."""
+    enc = _encoder()
+    tokens = torch.randint(_LO, _HI, (1, 50))
+    latent, _ = enc(tokens)
+    assert latent.shape == (1, 51, SMALL.d_model)
+
+
+def test_unmodified_residue_is_pure_token_embedding():
+    """Bias-free composition encoder => a zero-composition token gets no chemistry term."""
+    enc = _encoder()
+    he = enc.hybrid_embedding
+    zero_rows = (~he.composition.bool().any(dim=1)).nonzero().flatten()
+    tid = int(zero_rows[zero_rows > 0][0])
+    tokens = torch.tensor([[tid]])
+    assert torch.allclose(he(tokens), he.token_emb(tokens), atol=1e-6)
+
+
+def test_modified_residue_differs_from_token_only():
+    """A modified-residue token's hybrid embedding includes a chemistry term."""
+    enc = _encoder()
+    he = enc.hybrid_embedding
+    mod_rows = he.composition.bool().any(dim=1).nonzero().flatten()
+    tid = int(mod_rows[0])
+    tokens = torch.tensor([[tid]])
+    assert not torch.allclose(he(tokens), he.token_emb(tokens), atol=1e-6)
diff --git a/packages/peptide-property-ng/tests/test_fragment_targets.py b/packages/peptide-property-ng/tests/test_fragment_targets.py
new file mode 100644
index 000000000..2cb063e90
--- /dev/null
+++ b/packages/peptide-property-ng/tests/test_fragment_targets.py
@@ -0,0 +1,74 @@
+"""Tests for the Prosit-174 -> site-indexed intensity-target conversion.
+
+This conversion is the single place the ordinal->site remap happens, and it is
+easy to get wrong — hence the explicit b/y placement tests below.
+"""
+import numpy as np
+import pytest
+
+from peptide_property_ng.data.fragment_targets import (
+    N_ION_CHANNELS,
+    PROSIT_MAX_ORDINAL,
+    PROSIT_VECTOR_LEN,
+    prosit174_to_sites,
+)
+
+
+def _prosit_vector(entries: dict[tuple[int, int], float]) -> np.ndarray:
+    """Build a 174-vector; ``entries`` maps (fragment_ordinal, channel) -> value."""
+    v = np.zeros((PROSIT_MAX_ORDINAL, N_ION_CHANNELS), np.float32)
+    for (ordinal, channel), value in entries.items():
+        v[ordinal - 1, channel] = value
+    return v.reshape(-1)
+
+
+def test_b_ion_site_mapping():
+    """b_k (Prosit ordinal k) lands at site k-1, b channels (3:6)."""
+    length = 5  # 4 cleavage sites
+    t = prosit174_to_sites(_prosit_vector({(1, 3): 0.11, (4, 3): 0.44}), length)
+    assert t.shape == (length - 1, N_ION_CHANNELS)
+    assert t[0, 3] == 0.11   # b1 -> site 0
+    assert t[3, 3] == 0.44   # b4 -> site 3
+
+
+def test_y_ion_site_mapping():
+    """y_o (Prosit ordinal o) lands at site L-1-o, y channels (0:3)."""
+    length = 5
+    t = prosit174_to_sites(_prosit_vector({(1, 0): 0.81, (4, 0): 0.84}), length)
+    assert t[3, 0] == 0.81   # y1 -> site L-1-1 = 3
+    assert t[0, 0] == 0.84   # y4 -> site L-1-4 = 0
+
+
+def test_complementary_ions_share_a_site():
+    """Site i must hold b_{i+1} and y_{L-1-i} — the complementary pair of one cleavage."""
+    length = 6  # site 2 should carry b3 and y3
+    t = prosit174_to_sites(_prosit_vector({(3, 3): 0.3, (3, 0): 0.9}), length)
+    assert t[2, 3] == 0.3   # b3 at site 2
+    assert t[2, 0] == 0.9   # y3 at site 2
+
+
+def test_all_six_channels_pinned():
+    """Pin every channel — a charge permutation within a y/b triplet must fail."""
+    length = 5  # 4 sites; ordinal 2 -> b2 at site 1, y2 at site L-1-2 = 2
+    v = _prosit_vector({
+        (2, 0): 0.1, (2, 1): 0.2, (2, 2): 0.3,   # y2 charges 1,2,3
+        (2, 3): 0.4, (2, 4): 0.5, (2, 5): 0.6,   # b2 charges 1,2,3
+    })
+    t = prosit174_to_sites(v, length)
+    assert list(t[1, 3:6]) == [0.4, 0.5, 0.6]    # b2 -> site 1, b+1/b+2/b+3
+    assert list(t[2, 0:3]) == [0.1, 0.2, 0.3]    # y2 -> site 2, y+1/y+2/y+3
+
+
+def test_minus_one_mask_carries_through():
+    t = prosit174_to_sites(np.full(PROSIT_VECTOR_LEN, -1.0, np.float32), 10)
+    assert (t == -1.0).all()
+
+
+def test_shape_and_bounds():
+    v = np.zeros(PROSIT_VECTOR_LEN, np.float32)
+    assert prosit174_to_sites(v, 30).shape == (29, N_ION_CHANNELS)  # max peptide length
+    assert prosit174_to_sites(v, 3).shape == (2, N_ION_CHANNELS)
+    with pytest.raises(ValueError):
+        prosit174_to_sites(v, 31)               # peptide longer than the 174-vector
+    with pytest.raises(ValueError):
+        prosit174_to_sites(np.zeros(100), 10)   # wrong vector length
diff --git a/packages/peptide-property-ng/tests/test_losses.py b/packages/peptide-property-ng/tests/test_losses.py
new file mode 100644
index 000000000..3fa536551
--- /dev/null
+++ b/packages/peptide-property-ng/tests/test_losses.py
@@ -0,0 +1,33 @@
+"""Tests for the multi-task losses and degenerate-target handling."""
+import torch
+
+from peptide_property_ng.losses import MultiTaskLoss, intensity_signal_mask, masked_spectral_angle
+
+
+def test_spectral_angle_identical_is_zero():
+    x = torch.rand(4, 6, 6)
+    sa = masked_spectral_angle(x, x.clone())
+    assert torch.allclose(sa, torch.zeros(4), atol=1e-3)
+
+
+def test_spectral_angle_fully_masked_is_finite():
+    """A target with every entry masked (-1) must not produce NaN."""
+    sa = masked_spectral_angle(torch.rand(2, 5, 6), torch.full((2, 5, 6), -1.0))
+    assert torch.isfinite(sa).all()
+
+
+def test_signal_mask():
+    target = torch.zeros(3, 4, 6)
+    target[0, 0, 0] = 0.5   # sample 0 has an observed peak
+    target[2] = -1.0        # sample 2 fully masked -> degenerate
+    assert intensity_signal_mask(target).tolist() == [True, False, False]
+
+
+def test_loss_skips_degenerate_intensity_samples():
+    """A batch mixing real and degenerate intensity targets yields a finite loss."""
+    target = torch.zeros(4, 5, 6)
+    target[0, 1, 0] = 1.0  # only sample 0 carries signal
+    total, parts = MultiTaskLoss()({"intensity": torch.rand(4, 5, 6)},
+                                   {"intensity_target": target})
+    assert "intensity" in parts
+    assert torch.isfinite(torch.tensor(parts["intensity"]))
diff --git a/packages/peptide-property-ng/tests/test_multitask.py b/packages/peptide-property-ng/tests/test_multitask.py
new file mode 100644
index 000000000..0dd3b3f4f
--- /dev/null
+++ b/packages/peptide-property-ng/tests/test_multitask.py
@@ -0,0 +1,69 @@
+"""Tests for the unified multi-task model — output shapes and the leakage guard."""
+import torch
+
+from peptide_property_ng.model.config import SMALL
+from peptide_property_ng.model.multitask import UnifiedPeptidePropertyModel
+
+
+def _batch(b: int = 5, length: int = 14) -> dict[str, torch.Tensor]:
+    return {
+        "tokens": torch.randint(1, 60, (b, length)),
+        "charge": torch.randint(1, 5, (b,)),
+        "precursor_mz": torch.rand(b) * 800 + 300,
+        "collision_energy": torch.zeros(b),
+        "instrument": torch.randint(0, SMALL.n_instruments, (b,)),
+        "acq_mode": torch.randint(0, SMALL.n_acq_modes, (b,)),
+    }
+
+
+def test_forward_all_heads():
+    model = UnifiedPeptidePropertyModel(SMALL).eval()
+    b, length = 5, 14
+    out = model(_batch(b, length))
+
+    assert set(out) == {"intensity", "ccs", "rt", "charge"}
+    assert out["intensity"].shape == (b, length - 1, SMALL.n_ion_channels)
+    assert (out["intensity"] >= 0).all() and (out["intensity"] <= 1).all()
+
+    mean, std = out["ccs"]
+    assert mean.shape == (b,) and std.shape == (b,)
+    assert (std > 0).all()
+
+    assert out["rt"].shape == (b,)
+    assert out["charge"].shape == (b, SMALL.max_charge)
+
+
+def test_task_subset():
+    model = UnifiedPeptidePropertyModel(SMALL).eval()
+    out = model(_batch(), tasks=["intensity"])
+    assert set(out) == {"intensity"}
+
+
+def test_charge_head_does_not_see_charge():
+    """The charge head must be invariant to the precursor charge input
+    (charge is conditioned only at the intensity/CCS heads, never the encoder)."""
+    model = UnifiedPeptidePropertyModel(SMALL).eval()
+    batch = _batch()
+    with torch.no_grad():
+        logits_a = model(batch, tasks=["charge"])["charge"]
+        batch2 = dict(batch, charge=torch.full_like(batch["charge"], 5))
+        logits_b = model(batch2, tasks=["charge"])["charge"]
+    assert torch.allclose(logits_a, logits_b), "charge leaked into the charge head"
+
+
+def test_backward_runs():
+    model = UnifiedPeptidePropertyModel(SMALL)
+    out = model(_batch())
+    loss = (
+        out["intensity"].mean()
+        + out["ccs"][0].mean()
+        + out["rt"].mean()
+        + out["charge"].mean()
+    )
+    loss.backward()
+    assert any(p.grad is not None and p.grad.abs().sum() > 0 for p in model.parameters())
+
+
+def test_parameter_count_reported():
+    model = UnifiedPeptidePropertyModel(SMALL)
+    assert model.num_parameters() > 1_000_000
diff --git a/packages/peptide-property-ng/tests/test_shape_heads.py b/packages/peptide-property-ng/tests/test_shape_heads.py
new file mode 100644
index 000000000..9f3e41e4a
--- /dev/null
+++ b/packages/peptide-property-ng/tests/test_shape_heads.py
@@ -0,0 +1,77 @@
+"""Tests for the timsim peak-shape supervision: RT-sigma head + IM-sigma loss.
+
+Covers the additive wiring: a new ``rt_sigma`` head, and an ``im_sigma`` loss
+term riding the existing ``ccs`` head's ``std``. Both must be optional — absent
+targets ⇒ no term, existing training unchanged.
+"""
+import dataclasses
+
+import torch
+
+from peptide_property_ng.model.config import SMALL
+from peptide_property_ng.model.multitask import UnifiedPeptidePropertyModel
+from peptide_property_ng.losses import MultiTaskLoss
+
+
+def _cfg_with_shape():
+    return dataclasses.replace(SMALL, tasks=tuple(SMALL.tasks) + ("rt_sigma",))
+
+
+def _batch(b: int = 6, length: int = 14) -> dict[str, torch.Tensor]:
+    return {
+        "tokens": torch.randint(1, 60, (b, length)),
+        "charge": torch.randint(1, 5, (b,)),
+        "precursor_mz": torch.rand(b) * 800 + 300,
+        "collision_energy": torch.zeros(b),
+        "instrument": torch.randint(0, SMALL.n_instruments, (b,)),
+        "acq_mode": torch.randint(0, SMALL.n_acq_modes, (b,)),
+    }
+
+
+def test_rt_sigma_head_emits_positive_scalar():
+    model = UnifiedPeptidePropertyModel(_cfg_with_shape()).eval()
+    out = model(_batch(6))
+    assert "rt_sigma" in out
+    assert out["rt_sigma"].shape == (6,)
+    assert (out["rt_sigma"] > 0).all()  # softplus keeps sigma strictly positive
+
+
+def test_shape_loss_terms_present_and_masked():
+    b = 6
+    model = UnifiedPeptidePropertyModel(_cfg_with_shape())
+    out = model(_batch(b))
+    # only a subset of peptides carry a measured shape label
+    rt_valid = torch.tensor([True, True, False, True, False, False])
+    im_valid = torch.tensor([True, False, True, False, False, True])
+    batch = {
+        "rt_sigma_target": torch.rand(b),
+        "rt_sigma_valid": rt_valid,
+        "im_sigma_target": torch.rand(b) * 0.02,
+        "im_sigma_valid": im_valid,
+        # CCS-mean term inactive here (no ccs_valid set true)
+        "ccs_valid": torch.zeros(b, dtype=torch.bool),
+        "rt_valid": torch.zeros(b, dtype=torch.bool),
+        "charge_target": torch.randint(0, SMALL.max_charge, (b,)),
+        "intensity_target": torch.zeros(b, 13, SMALL.n_ion_channels),
+    }
+    total, logd = MultiTaskLoss()(out, batch)
+    assert "rt_sigma" in logd and "im_sigma" in logd
+    total.backward()
+    assert any(p.grad is not None and p.grad.abs().sum() > 0 for p in model.parameters())
+
+
+def test_shape_terms_absent_when_no_target():
+    """No shape targets in the batch ⇒ no shape loss terms (backward-compatible)."""
+    b = 5
+    model = UnifiedPeptidePropertyModel(SMALL)  # default tasks, no rt_sigma head
+    out = model(_batch(b))
+    batch = {
+        "ccs_valid": torch.ones(b, dtype=torch.bool),
+        "ccs_target": torch.rand(b),
+        "rt_valid": torch.ones(b, dtype=torch.bool),
+        "rt_target": torch.rand(b),
+        "charge_target": torch.randint(0, SMALL.max_charge, (b,)),
+        "intensity_target": torch.zeros(b, 13, SMALL.n_ion_channels),
+    }
+    _total, logd = MultiTaskLoss()(out, batch)
+    assert "rt_sigma" not in logd and "im_sigma" not in logd