feat(analysis): bind CIF source identities to H5MD atom order

docs/source_order_binding.md

@@ -0,0 +1,94 @@
+# CIF / H5MD source-order binding, version 1
+
+The bundle saver binds source atom identities in its **final serialization
+order** when `source_atom_ids` is supplied. XPONGE resolves IDs by Atom identity
+after its preparation/reordering; XpongeCPP uses the native serialized
+`molecule.atoms` order. Exports without source IDs retain the previous hashes.
+
+## Hash contract
+
+`digest(value)` is `"sha256:" + SHA256(canonical_json(value)).hexdigest()`.
+Canonical JSON uses Python `json.dumps(sort_keys=True, separators=(",", ":"),
+ensure_ascii=True, allow_nan=False)`, encoded as UTF-8. IDs are unique nonempty
+strings. No Unicode normalization or whitespace trimming is performed.
+
+```text
+source_atom_order_hash = digest({
+  "schema": "sponge-source-atom-order-v1",
+  "source_atom_ids": [IDs in final simulation order]
+})
+
+atom_order_hash = digest({
+  "schema": "sponge-bound-atom-order-v1",
+  "base_atom_order_hash": original native atom_order_hash,
+  "source_atom_order_hash": source_atom_order_hash
+})
+```
+
+The original native order digest is retained as `base_atom_order_hash`. Thus
+permuting distinct IDs changes the bound order digest even if the swapped atoms
+have identical mass, charge and residue properties. The existing topology hash
+continues to identify the native physical topology payload; the added binding
+and source IDs are provenance metadata.
+
+Before staged bundle publication, the saver writes:
+
+- `/topology/source_atom_ids`: UTF-8 array in final simulation order.
+- `/topology/source_order_binding`: UTF-8 JSON object defined below.
+- `/topology/atom_order_hash`: new bound order digest.
+- Restart `/run/atom_order_hash`: the same digest.
+
+SPONGE already propagates this opaque topology order hash to
+`/parameters/sponge/topology_compatibility/atom_order_hash` in output H5MD,
+and restores global atom order before coordinate output. No per-frame IDs or
+SPONGE format change is required.
+
+## Mapping sidecar
+
+Mokda checks that the serialized source IDs exactly equal the final mapping's
+`external_id` sequence and that the binding agrees with current topology
+metadata. It then adds an optional `topology_binding` object to the existing
+`sponge-atom-order-mapping` version 1 JSON and save manifest:
+
+```json
+{
+  "schema": "sponge-source-order-binding",
+  "schema_version": 1,
+  "hash_algorithm": "sha256",
+  "atom_count": 2,
+  "base_atom_order_hash": "native exporter digest",
+  "source_atom_order_hash": "sha256:...",
+  "atom_order_hash": "sha256:...",
+  "topology_hash": "native topology digest"
+}
+```
+
+This is bound during native topology export, **never by copying values from a
+chosen trajectory**. Chemcore's ordered CIF and existing `mapping_hash` remain
+unchanged: they prove the mapping's complete identity rows correspond to the
+CIF atom order. The source-ID digest connects those rows to the native topology.
+
+## Analysis validation
+
+`load_cif_h5md_universe` checks:
+
+1. CIF mapping order, IDs and digest, then equality with the sidecar.
+2. Binding schema and recomputed source-ID / bound-order digests.
+3. CIF/H5MD atom counts.
+4. Both available H5MD `atom_order_hash` and `topology_hash` against the binding.
+
+Conflicting available hashes or malformed bindings are rejected even with
+`strict=False`. Missing sidecars, old sidecars without a binding, and older
+H5MD files missing hashes remain readable with a warning and unverified status.
+A custom particle stream cannot be proven by SPONGE's global hashes and is
+reported as unverified. Raw exports currently have no source-order binding.
+
+Results expose `universe.cif_metadata["trajectory_order_validation"]` with
+`verified`, `method` (when verified), or `reason` (when unverified). Mokda copies
+this into result `meta.trajectoryOrderValidation` and
+`analysis_inputs.json.trajectory_order_validation`. This remains available when
+the original topology H5 has been removed; the CIF and sidecar suffice.
+
+This verifies export provenance and order consistency under the writer's
+coordinate-order contract. It is not a signature, nor does it detect coordinates
+edited after simulation while their provenance metadata is deliberately retained.

src/Xponge/io_bundle/source_order.py

@@ -0,0 +1,2 @@
+"""Compatibility import for source identity order bindings."""
+from XpongeCPP.io_bundle.source_order import *  # noqa: F401,F403

src/XpongeCPP/analysis/cif_mdanalysis.py

@@ -15,6 +15,7 @@
 from MDAnalysis.lib.util import openany
 from MDAnalysis.topology.base import TopologyReaderBase
 
+from ..io_bundle.source_order import validate_source_order_binding, validate_trajectory_source_order
 from ..io_bundle.errors import BundleTopologyError, BundleTrajectoryError
 
 
@@ -823,8 +824,17 @@ def _verify_mokda_mapping_file(mapping_file, cif_mapping):
         raise BundleTopologyError(
             "Mokda mapping file atoms do not match the CIF trajectory mapping"
         )
+    binding = None
+    if "topology_binding" in document:
+        try:
+            binding = validate_source_order_binding(
+                document["topology_binding"], [row["external_id"] for row in sidecar_mapping]
+            )
+        except ValueError as exc:
+            raise BundleTopologyError(f"invalid mapping topology_binding: {exc}") from exc
     return {
         "path": os.fspath(mapping_file),
+        "topology_binding": binding,
         "verified": True,
         "mapping_hash": cif_mapping_hash,
         "atom_count": len(cif_identity_mapping),
@@ -887,6 +897,15 @@ def load_cif_h5md_universe(
             f"{trajectory!r} has {trajectory_atom_count} atoms, CIF topology has "
             f"{cif_topology.n_atoms}"
         )
+    try:
+        order_validation = validate_trajectory_source_order(
+            trajectory,
+            cif_topology.cif_metadata["mapping_file_validation"].get("topology_binding"),
+            particle_stream=particle_stream,
+        )
+    except ValueError as exc:
+        raise BundleTrajectoryError(str(exc)) from exc
+    cif_topology.cif_metadata["trajectory_order_validation"] = order_validation
     if cif_topology.cif_metadata["mokda_trajectory_mapping"]:
         mapping_file_verified = cif_topology.cif_metadata["mapping_file_validation"][
             "verified"
@@ -906,7 +925,8 @@ def load_cif_h5md_universe(
             "CIF/H5MD compatibility is checked by atom count only; ensure the CIF "
             "atom order matches the trajectory atom order."
         )
-    warnings.warn(warning, RuntimeWarning, stacklevel=2)
+    if not order_validation["verified"]:
+        warnings.warn(warning, RuntimeWarning, stacklevel=2)
 
     universe = mda.Universe(
         cif_topology,

src/XpongeCPP/io_bundle/saver.py

@@ -16,6 +16,7 @@
 )
 from .errors import BundlePathError
 from .protocol import SpongeProtocol, add_protocol_to_bundle
+from .source_order import bind_bundle_source_order
 
 
 def save_sponge_input_bundle(
@@ -65,6 +66,29 @@ def save_sponge_input_bundle(
                 "ResidueType, or registered template-like object"
             )
 
+    values = None
+    if return_mapping and source_atom_ids is None:
+        raise ValueError("return_mapping=True requires source_atom_ids")
+    if source_atom_ids is not None:
+        atoms = list(target.atoms)
+        if isinstance(source_atom_ids, dict):
+            by_index = {
+                int(atom.index): str(value) for atom, value in source_atom_ids.items()
+            }
+            if set(by_index) != {int(atom.index) for atom in atoms}:
+                raise ValueError(
+                    "source_atom_ids mapping must cover every input Atom exactly once"
+                )
+            values = tuple(by_index[int(atom.index)] for atom in atoms)
+        else:
+            values = tuple(str(value) for value in source_atom_ids)
+            if len(values) != len(atoms):
+                raise ValueError(
+                    "source_atom_ids must contain one ID for every input Atom"
+                )
+        if len(set(values)) != len(values):
+            raise ValueError("source_atom_ids must be unique")
+
     output_root = Path(dirname).resolve()
     output_root.mkdir(parents=True, exist_ok=True)
     normalized_prefix = str(prefix or getattr(molecule, "name", "system"))
@@ -79,6 +103,8 @@ def save_sponge_input_bundle(
             target, staged_prefix, str(staging_root), protocol=None
         )
         staged_paths = _prefixed_bundle_paths(staging_root, staged_prefix)
+        if values is not None:
+            bind_bundle_source_order(staged_paths, values)
         _apply_protocol(staged_paths, protocol)
         for source, destination in (
             (staged_paths.topology, final_paths.topology),
@@ -90,26 +116,6 @@ def save_sponge_input_bundle(
     del prepared
     if not return_mapping:
         return target
-    if source_atom_ids is None:
-        raise ValueError("return_mapping=True requires source_atom_ids")
-    atoms = list(target.atoms)
-    if isinstance(source_atom_ids, dict):
-        by_index = {
-            int(atom.index): str(value) for atom, value in source_atom_ids.items()
-        }
-        if set(by_index) != {int(atom.index) for atom in atoms}:
-            raise ValueError(
-                "source_atom_ids mapping must cover every input Atom exactly once"
-            )
-        values = tuple(by_index[int(atom.index)] for atom in atoms)
-    else:
-        values = tuple(str(value) for value in source_atom_ids)
-        if len(values) != len(atoms):
-            raise ValueError(
-                "source_atom_ids must contain one ID for every input Atom"
-            )
-    if len(set(values)) != len(values):
-        raise ValueError("source_atom_ids must be unique")
     return target, tuple(
         {"simulation_index": index, "source_atom_id": source_id}
         for index, source_id in enumerate(values)

src/XpongeCPP/io_bundle/source_order.py

@@ -0,0 +1,132 @@
+"""Versioned binding of serialized source identities to native atom order."""
+
+from __future__ import annotations
+
+import hashlib
+import json
+
+import h5py
+import numpy as np
+
+
+def _digest(value):
+    payload = json.dumps(value, sort_keys=True, separators=(",", ":"), allow_nan=False)
+    return "sha256:" + hashlib.sha256(payload.encode("utf-8")).hexdigest()
+
+
+def make_source_order_binding(source_atom_ids, base_atom_order_hash, topology_hash):
+    """Bind the actual serialized identity sequence to its native order hash."""
+    identities = list(source_atom_ids)
+    if not identities or any(not isinstance(value, str) or not value for value in identities):
+        raise ValueError("source atom IDs must be nonempty strings")
+    if len(set(identities)) != len(identities):
+        raise ValueError("source atom IDs must be unique")
+    if not all(isinstance(value, str) and value for value in (base_atom_order_hash, topology_hash)):
+        raise ValueError("source order binding requires native topology and atom-order hashes")
+    source_hash = _digest({"schema": "sponge-source-atom-order-v1", "source_atom_ids": identities})
+    order_hash = _digest({
+        "schema": "sponge-bound-atom-order-v1",
+        "base_atom_order_hash": base_atom_order_hash,
+        "source_atom_order_hash": source_hash,
+    })
+    return {
+        "schema": "sponge-source-order-binding",
+        "schema_version": 1,
+        "hash_algorithm": "sha256",
+        "atom_count": len(identities),
+        "base_atom_order_hash": base_atom_order_hash,
+        "source_atom_order_hash": source_hash,
+        "atom_order_hash": order_hash,
+        "topology_hash": topology_hash,
+    }
+
+
+def validate_source_order_binding(binding, source_atom_ids):
+    """Reject malformed bindings and identity sequences different from export."""
+    if not isinstance(binding, dict):
+        raise ValueError("source order binding must be an object")
+    expected = make_source_order_binding(
+        source_atom_ids, binding.get("base_atom_order_hash"), binding.get("topology_hash")
+    )
+    if binding != expected:
+        raise ValueError("source order binding does not match the atom identity sequence or schema")
+    return expected
+
+
+def _text(value):
+    return value.decode("utf-8") if isinstance(value, bytes) else str(value)
+
+
+def _write_text(handle, path, value):
+    if path in handle:
+        del handle[path]
+    handle.create_dataset(path, data=value, dtype=h5py.string_dtype("utf-8"))
+
+
+def bind_bundle_source_order(paths, source_atom_ids):
+    """Bind a staged bundle before publication; keep restart lineage consistent.
+
+    Only the exporter may call this with IDs in its actual serialization order.
+    No coordinates or force-field datasets are changed.
+    """
+    identities = list(source_atom_ids)
+    with h5py.File(paths.topology, "r+") as top, h5py.File(paths.restart, "r+") as restart:
+        if "/topology/source_order_binding" in top:
+            raise ValueError("bundle source atom order is already bound")
+        count = int(np.asarray(top["/topology/atom_count"][()]).reshape(-1)[0])
+        if len(identities) != count:
+            raise ValueError("source atom IDs must cover every serialized atom")
+        binding = make_source_order_binding(
+            identities,
+            _text(top["/topology/atom_order_hash"][()]),
+            _text(top["/topology/topology_hash"][()]),
+        )
+        top.create_dataset("/topology/source_atom_ids", data=identities, dtype=h5py.string_dtype("utf-8"))
+        _write_text(top, "/topology/source_order_binding", json.dumps(binding, sort_keys=True))
+        _write_text(top, "/topology/atom_order_hash", binding["atom_order_hash"])
+        _write_text(restart, "/run/atom_order_hash", binding["atom_order_hash"])
+    return binding
+
+
+def read_source_order_binding(topology, source_atom_ids):
+    """Read an export binding, checking its stored IDs and current H5 metadata."""
+    with h5py.File(topology, "r") as top:
+        if "/topology/source_order_binding" not in top:
+            return None
+        binding = validate_source_order_binding(
+            json.loads(_text(top["/topology/source_order_binding"][()])), source_atom_ids
+        )
+        if top["/topology/source_atom_ids"].asstr()[...].tolist() != list(source_atom_ids):
+            raise ValueError("topology source atom IDs do not match the final mapping")
+        for name in ("atom_order_hash", "topology_hash"):
+            if _text(top[f"/topology/{name}"][()]) != binding[name]:
+                raise ValueError(f"topology {name} does not match its source order binding")
+        count = int(np.asarray(top["/topology/atom_count"][()]).reshape(-1)[0])
+        if count != binding["atom_count"]:
+            raise ValueError("topology atom count does not match its source order binding")
+    return binding
+
+
+def validate_trajectory_source_order(trajectory, binding, *, particle_stream=None):
+    """Check available trajectory hashes; legacy missing hashes stay unverified.
+
+    SPONGE compatibility hashes describe the global ``all`` particle stream.
+    They cannot establish the identity order of an arbitrary custom stream.
+    """
+    if binding is None:
+        return {"verified": False, "reason": "missing_source_order_binding"}
+    missing = []
+    with h5py.File(trajectory, "r") as handle:
+        streams = list(handle.get("particles", {}))
+        selected = particle_stream or (streams[0] if len(streams) == 1 else None)
+        if selected != "all":
+            return {"verified": False, "reason": "unbound_particle_stream"}
+        for name in ("atom_order_hash", "topology_hash"):
+            path = f"/parameters/sponge/topology_compatibility/{name}"
+            if path not in handle or not _text(handle[path][()]):
+                missing.append(name)
+            elif _text(handle[path][()]) != binding[name]:
+                raise ValueError(f"H5MD {name} does not match the CIF mapping source order binding")
+    if missing:
+        return {"verified": False, "reason": "missing_trajectory_hash", "missing": missing}
+    return {"verified": True, "method": "source_order_binding_v1", **binding}