feat(cv): support inline H5 RMSD references with end-to-end coverage

SPONGE/collective_variable/collective_variable.cpp

@@ -208,9 +208,11 @@ bool Load_H5_CV_Config(CONTROLLER* controller,
         constexpr const char* cv_root = "/cv/config";
         constexpr const char* restraint_root = "/restraint/config";
         constexpr const char* restraint_cv_root = "/restraint/cv/config";
+        constexpr const char* steer_root = "/steer/config";
         const bool has_cv = file->exist(cv_root);
         const bool has_restraint = file->exist(restraint_root);
         const bool has_restraint_cv = file->exist(restraint_cv_root);
+        const bool has_steer = file->exist(steer_root);
         SpongeH5MD::ProtocolCVH5Reader cv_reader;
         std::vector<SpongeH5MD::ProtocolCVDefinition> typed_cvs;
         std::vector<SpongeH5MD::ProtocolVirtualAtomDefinition>
@@ -261,7 +263,7 @@ bool Load_H5_CV_Config(CONTROLLER* controller,
         {
             throw std::runtime_error(steering_reader.Last_Error());
         }
-        if (!has_cv && !has_restraint && !has_restraint_cv &&
+        if (!has_cv && !has_restraint && !has_restraint_cv && !has_steer &&
             typed_cvs.empty() && typed_virtual_atoms.empty() &&
             typed_restraints.empty() && !has_typed_metadynamics &&
             !has_typed_steering)
@@ -271,14 +273,16 @@ bool Load_H5_CV_Config(CONTROLLER* controller,
         const bool has_legacy_cv = controller->Command_Exist("cv_in_file");
         const bool has_legacy_restraint =
             controller->Command_Exist("restrain_in_file") ||
-            controller->Command_Exist("restrain_cv_in_file");
+            controller->Command_Exist("restrain_cv_in_file") ||
+            controller->Command_Exist("steer_cv_in_file");
         if (has_legacy_cv || has_legacy_restraint)
         {
             return false;
         }
 
         std::vector<CVConfigSection> sections;
-        for (const auto& root : {cv_root, restraint_root, restraint_cv_root})
+        for (const auto& root :
+             {cv_root, restraint_root, restraint_cv_root, steer_root})
         {
             if (file->exist(root))
             {
@@ -416,7 +420,8 @@ void COLLECTIVE_VARIABLE_CONTROLLER::Initial(
     const bool has_h5_cv = Load_H5_CV_Config(controller, this);
     if (has_h5_cv || controller->Command_Exist("cv_in_file") ||
         controller->Command_Exist("restrain_in_file") ||
-        controller->Command_Exist("restrain_cv_in_file"))
+        controller->Command_Exist("restrain_cv_in_file") ||
+        controller->Command_Exist("steer_cv_in_file"))
     {
         int CV_numbers = 0;
         Commands_From_In_File(controller);
@@ -522,8 +527,8 @@ static void Set_CV_Config_Command(COLLECTIVE_VARIABLE_CONTROLLER* manager,
 void COLLECTIVE_VARIABLE_CONTROLLER::Commands_From_In_File(
     CONTROLLER* controller)
 {
-    for (const char* input_key :
-         {"cv_in_file", "restrain_in_file", "restrain_cv_in_file"})
+    for (const char* input_key : {"cv_in_file", "restrain_in_file",
+                                  "restrain_cv_in_file", "steer_cv_in_file"})
     {
         if (!controller->Command_Exist(input_key)) continue;
         const std::string cv_path = controller->Command(input_key);

SPONGE/utils/h5md/protocol_cv_h5.hpp

@@ -156,7 +156,7 @@ class ProtocolCVH5Reader
                 const std::size_t selected_atom_count =
                     atom_indices.empty() ? atom_refs.size()
                                          : atom_indices.size();
-                Read_Restart_Reference(root, selected_atom_count, &definition);
+                Read_Reference(root, selected_atom_count, &definition);
                 Validate_Current_Runtime_Shape(selected_atom_count, definition);
                 definitions->push_back(std::move(definition));
             }
@@ -429,44 +429,57 @@ class ProtocolCVH5Reader
         }
     }
 
-    void Read_Restart_Reference(const std::string& root,
-                                std::size_t selected_atom_count,
-                                ProtocolCVDefinition* definition)
+    void Read_Reference(const std::string& root,
+                        std::size_t selected_atom_count,
+                        ProtocolCVDefinition* definition)
     {
+        const std::string inline_path = root + "/coordinate";
+        if (protocol_->exist(inline_path))
+        {
+            Read_Reference_Dataset(*protocol_, inline_path, selected_atom_count,
+                                   definition);
+        }
         const std::string path = "/parameters/restart/references/cv/" +
                                  definition->name + "/coordinate";
-        if (restart_ == nullptr || !restart_->exist(path))
+        if (restart_ != nullptr && restart_->exist(path))
         {
-            if (definition->type == "rmsd" &&
-                !Has_Runtime_Parameter(*definition, "coordinate"))
-            {
-                throw std::runtime_error(path +
-                                         " is required for a native rmsd CV");
-            }
-            return;
+            Read_Reference_Dataset(*restart_, path, selected_atom_count,
+                                   definition);
+        }
+        if (definition->type == "rmsd" &&
+            !Has_Runtime_Parameter(*definition, "coordinate"))
+        {
+            throw std::runtime_error(inline_path + " or restart " + path +
+                                     " is required for a native rmsd CV");
         }
+    }
+
+    void Read_Reference_Dataset(HighFive::File& file, const std::string& path,
+                                std::size_t selected_atom_count,
+                                ProtocolCVDefinition* definition)
+    {
         if (definition->type != "rmsd")
         {
             throw std::runtime_error(path +
                                      " is only supported for rmsd CV objects");
         }
-        const auto dims = restart_->getDataSet(path).getSpace().getDimensions();
+        const auto dims = file.getDataSet(path).getSpace().getDimensions();
         if (dims != std::vector<std::size_t>{selected_atom_count, 3})
         {
-            throw std::runtime_error(path +
-                                     " must have shape [atom_indices,3]");
+            throw std::runtime_error(
+                path + " must have shape [selected_atom_count,3]");
         }
         std::vector<float> values(selected_atom_count * 3);
-        auto dataset = restart_->getDataSet(path);
+        auto dataset = file.getDataSet(path);
         if (H5Dread(dataset.getId(), H5T_NATIVE_FLOAT, H5S_ALL, H5S_ALL,
                     H5P_DEFAULT, values.data()) < 0)
         {
             throw std::runtime_error("failed to read " + path);
         }
         Validate_Finite(values, path);
-        definition->reference_coordinates = values;
         Add_Runtime_Parameter(definition, "coordinate", Join_Values(values),
                               path);
+        definition->reference_coordinates = values;
     }
 
     void Validate_Current_Runtime_Shape(std::size_t selected_atom_count,

docs/input-reference/collective-variables.md

@@ -119,6 +119,19 @@ rotate = true
 
 `rotate = true` enables optimal rotational alignment before RMSD evaluation.
 
+For native H5 input, set `/cv/<name>/type` to `rmsd` and store reference
+coordinates in the protocol dataset `/cv/<name>/coordinate`, with shape
+`[selected_atom_count, 3]`. Rows follow the order of `atom_indices` or
+`atom_refs`; coordinates must be finite and use the same units as system
+coordinates. Xponge and XpongeCPP write this dataset from
+`ProtocolCollectiveVariable.reference_coordinates`.
+
+The legacy restart dataset
+`/parameters/restart/references/cv/<name>/coordinate` remains supported.
+If both datasets are supplied, their coordinates must agree; conflicting
+references are rejected. A native RMSD CV requires one of these references.
+The `coordinate` dataset is only supported for RMSD CVs.
+
 Parameters:
 
 | Parameter | Type | Description |

tests/h5_bundle/README.md

@@ -22,6 +22,44 @@ pixi run -e dev-cpu ctest --test-dir build-h5-tests --output-on-failure
 
 ## Test targets
 
+### RMSD CV execution across Xponge and XpongeCPP
+
+`test_rmsd_cv_e2e.py` launches the actual SPONGE executable, using each
+producer's own Python environment to generate its input. It is a separate,
+opt-in pytest suite requiring `pytest`, `numpy`, and `h5py` in the test runner
+and the corresponding Xponge package in each producer environment.
+
+From the repository root, select a freshly built CPU or GPU SPONGE executable:
+
+```bash
+export SPONGE_EXECUTABLE="$PWD/build-dev-cpu/SPONGE"
+export XPONGE_PYTHON="/path/to/XPONGE/.venv/bin/python"
+export XPONGE_CPP_PYTHON="/path/to/XpongeCPP/.pixi/envs/default/bin/python"
+# Optional: test source checkouts instead of installed package versions.
+export XPONGE_SOURCE="/path/to/XPONGE"
+export XPONGE_CPP_SOURCE="/path/to/XpongeCPP/src"
+python -m pytest -c /dev/null --confcutdir="$PWD" \
+  tests/h5_bundle/test_rmsd_cv_e2e.py -q -p no:cacheprovider
+```
+
+All eight cases must pass without skips to cover both producers. Each covers
+`rotate=false` or `rotate=true` and either:
+
+- Native inline RMSD reference and converted legacy input evaluated by SPONGE,
+  with RMSD checked against an independent NumPy/Kabsch calculation and bias
+  forces checked against finite differences after subtracting a zero-bias run.
+- Four uninterrupted NVE steps versus two steps followed by two restarted
+  steps, checking RMSD, forces, coordinates, velocities, physical time, and the
+  original reference preserved in the generated restart.
+
+The fixture keeps the peptide inside the periodic box and uses an asymmetric
+reference with an unsorted atom selection. Scalar tolerances account for the
+runtime's printed precision; force trajectories retain float32 precision.
+Temporary inputs, H5 outputs, and subprocess logs remain in pytest's temporary
+directory. CPU success does not imply GPU execution coverage.
+
+### CTest targets
+
 | Target | Scope |
 |---|---|
 | `test_h5_output_plan` | Parser-visible H5 output keys, defaults, suffix helpers, helper null/empty-key behavior, empty H5 path handling, full legacy sidecar resolution matrix, explicit legacy sidecar provenance collection, VDS chunk size, repair policy validation. |

tests/h5_bundle/rmsd_cv_e2e_producer.py

@@ -0,0 +1,103 @@
+"""Create native and converted RMSD fixtures using a producer's own Python."""
+
+import importlib
+from pathlib import Path
+import sys
+
+import h5py
+import numpy as np
+
+
+def main():
+    package, output, rotate_text = sys.argv[1:]
+    xponge = importlib.import_module(package)
+    importlib.import_module(package + ".forcefield.amber.ff14sb")
+    converter = importlib.import_module(package + ".io_bundle")
+    root = Path(output)
+    molecule = xponge.get_peptide_from_sequence("AA")
+    if package == "XpongeCPP":
+        molecule.set_box_padding(15.0)
+    else:
+        molecule.box_length = [50.0, 50.0, 50.0]
+        # Keep the whole peptide inside the box so periodic wrapping is not
+        # conflated with reference serialization or restart behavior.
+        center = np.mean(
+            [(atom.x, atom.y, atom.z) for atom in molecule.atoms], axis=0
+        )
+        shift = 25.0 - center
+        for atom in molecule.atoms:
+            atom.x += shift[0]
+            atom.y += shift[1]
+            atom.z += shift[2]
+    xponge.save_sponge_input_bundle(molecule, "system", root / "seed")
+    with h5py.File(root / "seed/system_restart.spgr.h5") as handle:
+        positions = handle["/particles/all/position/value"][0]
+        box = handle["/particles/all/box/edges/value"][0]
+    # Deliberately noncontiguous and unsorted, with an asymmetric deformation.
+    selection = np.asarray([7, 1, 11, 4])
+    reference = positions[selection].copy()
+    reference += np.asarray(
+        [[0.2, 0.4, -0.2], [-0.3, 0.1, 0.5], [0.1, -0.6, -0.2], [0.3, 0.2, 0.1]]
+    )
+    reference += np.asarray([1.5, -2.0, 0.75])
+    protocol = xponge.SpongeProtocol(
+        collective_variables=(
+            xponge.ProtocolCollectiveVariable(
+                name="rmsd_cv",
+                type="rmsd",
+                atom_indices=tuple(map(int, selection)),
+                reference_coordinates=tuple(map(tuple, reference)),
+                rotate=rotate_text == "true",
+            ),
+        )
+    )
+    for name, weight in (("native", 2.0), ("baseline", 0.0)):
+        case = root / name
+        xponge.save_sponge_input_bundle(
+            molecule, "system", case, protocol=protocol
+        )
+        # Print and bias configuration goes through the existing /cv/config
+        # route, while the RMSD definition and reference remain fully native.
+        with h5py.File(case / "system_protocol.spgp.h5", "a") as handle:
+            config = handle.require_group("/cv/config")
+            text = h5py.string_dtype()
+            config.create_dataset(
+                "section/name", data=["print", "restrain"], dtype=text
+            )
+            config.create_dataset("section/key_offset", data=[0, 1, 4])
+            config.create_dataset("section/count", data=2)
+            config.create_dataset(
+                "key", data=["CV", "CV", "weight", "reference"], dtype=text
+            )
+            config.create_dataset(
+                "value",
+                data=["rmsd_cv", "rmsd_cv", str(weight), "0.2"],
+                dtype=text,
+            )
+            assert handle["/cv/rmsd_cv/coordinate"].shape == (4, 3)
+        with h5py.File(case / "system_restart.spgr.h5") as handle:
+            assert (
+                "/parameters/restart/references/cv/rmsd_cv/coordinate"
+                not in handle
+            )
+        (case / "mdin.bundled.spg.toml").write_text(
+            'mode = "minimization"\ncutoff = 8.0\n'
+            'input_h5_topology_path = "system_topology.spgt.h5"\n'
+            'input_h5_protocol_path = "system_protocol.spgp.h5"\n'
+            'input_h5_restart_path = "system_restart.spgr.h5"\n'
+            'input_h5_restart_load = "structural"\n'
+        )
+    converter.convert_bundle_to_legacy(
+        root / "native", root / "legacy", prefix="system"
+    )
+    np.savez(
+        root / "oracle.npz",
+        positions=positions,
+        box=box,
+        reference=reference,
+        selection=selection,
+    )
+
+
+if __name__ == "__main__":
+    main()