diff --git a/SPONGE/collective_variable/collective_variable.cpp b/SPONGE/collective_variable/collective_variable.cpp index 867ced1a..43099ed5 100644 --- a/SPONGE/collective_variable/collective_variable.cpp +++ b/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 typed_cvs; std::vector @@ -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 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); diff --git a/SPONGE/utils/h5md/protocol_cv_h5.hpp b/SPONGE/utils/h5md/protocol_cv_h5.hpp index 90e2df6d..c571a1d3 100644 --- a/SPONGE/utils/h5md/protocol_cv_h5.hpp +++ b/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{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 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, diff --git a/docs/input-reference/collective-variables.md b/docs/input-reference/collective-variables.md index fd9af895..5b9f7d0a 100644 --- a/docs/input-reference/collective-variables.md +++ b/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//type` to `rmsd` and store reference +coordinates in the protocol dataset `/cv//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//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 | diff --git a/tests/h5_bundle/README.md b/tests/h5_bundle/README.md index f3812a7b..d6426e40 100644 --- a/tests/h5_bundle/README.md +++ b/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. | diff --git a/tests/h5_bundle/rmsd_cv_e2e_producer.py b/tests/h5_bundle/rmsd_cv_e2e_producer.py new file mode 100644 index 00000000..6bc35fce --- /dev/null +++ b/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() diff --git a/tests/h5_bundle/test_h5_input_validation.cpp b/tests/h5_bundle/test_h5_input_validation.cpp index 296021b0..3d7bb0ac 100644 --- a/tests/h5_bundle/test_h5_input_validation.cpp +++ b/tests/h5_bundle/test_h5_input_validation.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include #include @@ -228,6 +229,107 @@ static void Test_Protocol_Reader_Loads_Native_CV_Objects() std::filesystem::remove_all(dir); } +static void Test_Protocol_Reader_Inline_RMSD_Reference() +{ + const std::vector reference = {1, 2, 3, 4, 5, 6}; + struct Case + { + std::string type; + bool inline_reference; + bool restart_reference; + bool conflict; + int rows; + int columns; + bool nonfinite; + bool atom_refs; + std::string error; + }; + const std::vector cases = { + {"rmsd", true, false, false, 2, 3, false, false, ""}, + {"rmsd", true, false, false, 2, 3, false, true, ""}, + {"rmsd", true, true, false, 2, 3, false, false, ""}, + {"rmsd", false, true, false, 2, 3, false, false, ""}, + {"rmsd", true, true, true, 2, 3, false, false, "conflicts"}, + {"rmsd", false, false, false, 2, 3, false, false, "is required"}, + {"rmsd", true, false, false, 1, 3, false, false, "must have shape"}, + {"rmsd", true, false, false, 3, 2, false, false, "must have shape"}, + {"rmsd", true, false, false, 2, 3, true, false, "finite"}, + {"distance", true, false, false, 2, 3, false, false, + "only supported for rmsd"}, + }; + for (const auto& item : cases) + { + const auto dir = Unique_Temp_Path("protocol_inline_rmsd"); + std::filesystem::create_directories(dir); + const auto protocol = dir / "protocol.spgp.h5"; + const auto restart = dir / "restart.spgr.h5"; + { + HighFive::File file(protocol.string(), HighFive::File::Overwrite); + Write_Scalar(file, "/cv/backbone/type", item.type); + if (item.atom_refs) + Write_String_Vector(file, "/cv/backbone/atom_refs", {"0", "1"}); + else + Write_Int64_Vector(file, "/cv/backbone/atom_indices", {0, 1}); + if (item.inline_reference) + { + auto values = reference; + values.resize(item.rows * item.columns); + if (item.nonfinite) + values[0] = std::numeric_limits::quiet_NaN(); + Write_Matrix(file, "/cv/backbone/coordinate", values, item.rows, + item.columns); + } + } + if (item.restart_reference) + { + auto values = reference; + if (item.conflict) values[0] += 1; + HighFiveBackend backend; + RestartH5Writer writer(&backend); + SpongeH5OutputPlan::ResolvedOutputPlan plan; + plan.restart.enabled = true; + plan.restart.path = restart.string(); + REQUIRE_TRUE(writer.Open(plan, SpongeH5MD::kInputSchemaVersion)); + REQUIRE_TRUE(writer.Write_Lineage("top", "atoms", "protocol")); + REQUIRE_TRUE(writer.Define_Structural_State(2, false)); + const std::array box = {10, 0, 0, 0, 10, 0, 0, 0, 10}; + REQUIRE_TRUE(writer.Write_Structural_State(0, 0.0, reference.data(), + box.data(), nullptr)); + REQUIRE_TRUE( + writer.Write_CV_Reference("backbone", values.data(), 2)); + REQUIRE_TRUE(writer.Finalize()); + REQUIRE_TRUE(writer.Close()); + } + { + ProtocolCVH5Reader reader; + REQUIRE_TRUE(reader.Open_Protocol(protocol.string())); + if (item.restart_reference) + REQUIRE_TRUE(reader.Open_Restart(restart.string())); + std::vector definitions; + const bool ok = reader.Read_Definitions(2, &definitions); + if (item.error.empty()) + { + REQUIRE_TRUE(ok); + REQUIRE_EQ(definitions.size(), static_cast(1)); + REQUIRE_EQ(definitions[0].reference_coordinates, reference); + REQUIRE_TRUE( + std::find(definitions[0].runtime_parameters.begin(), + definitions[0].runtime_parameters.end(), + std::make_pair(std::string("coordinate"), + std::string("1 2 3 4 5 6"))) != + definitions[0].runtime_parameters.end()); + } + else + { + REQUIRE_TRUE(!ok); + REQUIRE_TRUE(reader.Last_Error().find(item.error) != + std::string::npos); + } + } + std::filesystem::remove_all(dir); + } +} + static void Test_Protocol_Reader_Loads_Typed_Virtual_Atoms() { const auto dir = Unique_Temp_Path("protocol_virtual_atom"); @@ -1494,6 +1596,7 @@ int main() { Test_Protocol_Reader_Loads_Typed_CV_Restraint(); Test_Protocol_Reader_Loads_Native_CV_Objects(); + Test_Protocol_Reader_Inline_RMSD_Reference(); Test_Protocol_Reader_Loads_Typed_Virtual_Atoms(); Test_Protocol_Reader_Loads_Native_Metadynamics_Object(); Test_Protocol_Reader_Loads_Native_Steering_Object(); diff --git a/tests/h5_bundle/test_rmsd_cv_e2e.py b/tests/h5_bundle/test_rmsd_cv_e2e.py new file mode 100644 index 00000000..ea644982 --- /dev/null +++ b/tests/h5_bundle/test_rmsd_cv_e2e.py @@ -0,0 +1,268 @@ +"""Real SPONGE execution gates for both RMSD H5 producers and conversion. + +Set SPONGE_EXECUTABLE, XPONGE_PYTHON, and XPONGE_CPP_PYTHON to run both +producers. Optional XPONGE_SOURCE / XPONGE_CPP_SOURCE select source roots +(the latter points to XpongeCPP/src). All subprocess output stays in tmp_path. +""" + +import json +import os +from pathlib import Path +import shutil +import subprocess +import tomllib + +import h5py +import numpy as np +import pytest + + +@pytest.fixture(params=[("Xponge", "XPONGE"), ("XpongeCPP", "XPONGE_CPP")]) +def producer(request): + package, prefix = request.param + python = os.environ.get(prefix + "_PYTHON") + executable = os.environ.get("SPONGE_EXECUTABLE") + if not python or not executable: + pytest.skip( + f"set {prefix}_PYTHON and SPONGE_EXECUTABLE for real RMSD execution" + ) + for path in (python, executable): + assert Path(path).is_file() and os.access(path, os.X_OK), path + env = dict( + os.environ, + OMP_NUM_THREADS="1", + MKL_NUM_THREADS="1", + PYTHONDONTWRITEBYTECODE="1", + ) + source = os.environ.get(prefix + "_SOURCE") + if source: + env["PYTHONPATH"] = source + return package, python, str(Path(executable).resolve()), env + + +def _checked(command, cwd, env): + result = subprocess.run( + command, + cwd=cwd, + env=env, + capture_output=True, + text=True, + timeout=120, + check=False, + ) + (cwd / "execution.log").write_text(result.stdout + "\n" + result.stderr) + assert result.returncode == 0, ( + f"{command}\n{result.stdout}\n{result.stderr}" + ) + + +def _generate(producer, root, rotate): + package, python, _, env = producer + script = Path(__file__).with_name("rmsd_cv_e2e_producer.py") + _checked( + [python, str(script), package, str(root), str(rotate).lower()], + root, + env, + ) + with np.load(root / "oracle.npz") as data: + return {key: data[key] for key in data.files} + + +def _run(producer, case, *, rerun=True, steps=1, restart=None): + _, _, executable, env = producer + mdin = ( + "mdin.legacy.spg.toml" + if case.name == "legacy" + else "mdin.bundled.spg.toml" + ) + command = [ + executable, + "-mdin", + mdin, + "-mode", + "rerun" if rerun else "nve", + "-step_limit", + str(steps + 1 if rerun else steps), + "-dt", + "0.000001", + "-print_zeroth_frame", + "1", + "-write_information_interval", + "1", + "-write_mdout_interval", + "1", + "-write_trajectory_interval", + "1", + "-write_restart_file_interval", + str(steps), + "-mdout", + "mdout.txt", + "-mdinfo", + "mdinfo.txt", + "-frc", + "forces.dat", + "-output_h5_trajectory_path", + "trajectory.h5", + "-output_h5_trajectory_vds", + "0", + ] + if rerun: + command += [ + "-crd", + "frames.dat", + "-box", + "frames.box", + "-rerun_frame_limit", + "2", + "-rerun_start", + "0", + "-rerun_strip", + "0", + ] + else: + command += [ + "-output_h5_restart_path", + "final.spgr.h5", + "-vel", + "velocity.dat", + ] + if restart is not None: + command += ["-input_h5_restart_path", str(restart)] + # SPONGE rejects duplicate mdin/CLI keys; write a single resolved launch + # file instead of relying on command-line override semantics. + settings = tomllib.loads((case / mdin).read_text()) + settings.update( + { + key.removeprefix("-"): value + for key, value in zip(command[3::2], command[4::2]) + } + ) + (case / "run.spg.toml").write_text( + "\n".join( + f"{key} = {json.dumps(value)}" for key, value in settings.items() + ) + + "\n" + ) + _checked([executable, "-mdin", "run.spg.toml"], case, env) + with h5py.File(case / "trajectory.h5") as handle: + return { + "rmsd": handle["/observables/all/rmsd_cv/value"][...], + "bias": handle["/observables/all/restrain_cv/value"][...], + "force": handle["/particles/all/force/value"][...], + "position": handle["/particles/all/position/value"][...], + } + + +def _rmsd(positions, reference, rotate): + points = np.asarray(positions, dtype=np.float64) + target = np.asarray(reference, dtype=np.float64) + points = points - points.mean(axis=0) + target = target - target.mean(axis=0) + if rotate: + left, _, right = np.linalg.svd(target.T @ points) + correction = np.diag([1.0, 1.0, np.linalg.det(left @ right)]) + target = target @ (left @ correction @ right) + return np.sqrt(np.sum((points - target) ** 2) / len(points)) + + +def _bias_force(positions, reference, rotate): + # Finite differences of an independent NumPy/Kabsch energy oracle. + points = np.asarray(positions, dtype=np.float64).copy() + force = np.zeros_like(points) + delta = 1e-5 + for index in np.ndindex(points.shape): + original = points[index] + points[index] = original + delta + plus = 2.0 * (_rmsd(points, reference, rotate) - 0.2) ** 2 + points[index] = original - delta + minus = 2.0 * (_rmsd(points, reference, rotate) - 0.2) ** 2 + points[index] = original + force[index] = -(plus - minus) / (2 * delta) + return force + + +@pytest.mark.parametrize("rotate", [False, True]) +def test_rmsd_native_and_converted_legacy_value_and_force( + producer, tmp_path, rotate +): + data = _generate(producer, tmp_path, rotate) + results = {} + for name in ("native", "baseline", "legacy"): + case = tmp_path / name + # SPONGE excludes the final rerun frame from trajectory output. Two + # identical frames give one force snapshot at the exact oracle input. + np.stack([data["positions"]] * 2).astype(np.float32).tofile( + case / "frames.dat" + ) + np.savetxt( + case / "frames.box", [[*np.diag(data["box"]), 90, 90, 90]] * 2 + ) + results[name] = _run(producer, case) + selected = data["positions"][data["selection"]] + expected = _rmsd(selected, data["reference"], rotate) + for name in ("native", "baseline", "legacy"): + np.testing.assert_allclose( + results[name]["rmsd"], [expected] * 2, atol=6e-5, rtol=0 + ) + assert results[name]["force"].shape == (1, len(data["positions"]), 3) + # Observables use the runtime's printed precision (4 decimals for CV, + # 2 for bias energy); force trajectories retain float32 precision. + np.testing.assert_allclose( + results["native"]["bias"], [2 * (expected - 0.2) ** 2] * 2, atol=0.006 + ) + expected_force = np.zeros_like(data["positions"]) + expected_force[data["selection"]] = _bias_force( + selected, data["reference"], rotate + ) + difference = results["native"]["force"][0] - results["baseline"]["force"][0] + np.testing.assert_allclose(difference, expected_force, atol=2e-3, rtol=2e-3) + np.testing.assert_allclose( + results["legacy"]["force"], + results["native"]["force"], + atol=2e-2, + rtol=2e-5, + ) + + +@pytest.mark.parametrize("rotate", [False, True]) +def test_rmsd_restart_matches_uninterrupted_nve(producer, tmp_path, rotate): + data = _generate(producer, tmp_path, rotate) + for name in ("continuous", "first", "resumed"): + shutil.copytree(tmp_path / "native", tmp_path / name) + continuous = _run(producer, tmp_path / "continuous", rerun=False, steps=4) + _run(producer, tmp_path / "first", rerun=False, steps=2) + restart = tmp_path / "first/final.spgr.h5" + with h5py.File(restart) as handle: + np.testing.assert_array_equal( + handle["/parameters/restart/references/cv/rmsd_cv/coordinate"][...], + data["reference"], + ) + assert "/particles/all/velocity/value" in handle + resumed = _run( + producer, tmp_path / "resumed", rerun=False, steps=2, restart=restart + ) + np.testing.assert_allclose( + resumed["rmsd"][-1], continuous["rmsd"][-1], atol=1e-4, rtol=0 + ) + for field in ("position", "force"): + assert len(resumed[field]) > 0 and len(continuous[field]) > 0 + np.testing.assert_allclose( + resumed[field][-1], continuous[field][-1], atol=2e-3, rtol=2e-5 + ) + with ( + h5py.File(tmp_path / "continuous/final.spgr.h5") as full, + h5py.File(tmp_path / "resumed/final.spgr.h5") as split, + ): + np.testing.assert_allclose( + split["/particles/all/position/time"][...], + full["/particles/all/position/time"][...], + atol=1e-12, + rtol=0, + ) + for field in ("position", "velocity"): + np.testing.assert_allclose( + split[f"/particles/all/{field}/value"][...], + full[f"/particles/all/{field}/value"][...], + atol=2e-5, + rtol=2e-5, + )