Add support for dynamics crash reports

doc/source/changelog.rst

@@ -31,6 +31,8 @@ organisation on `GitHub <https://github.com/openbiosim/sire>`__.
 
 * Add support for ``openmm.MonteCarloMembraneBarostat`` in Sire-to-OpenMM conversion.
 
+* Added support for saving crash reports during Sire dynamics runs.
+
 `2025.2.0 <https://github.com/openbiosim/sire/compare/2025.1.0...2025.2.0>`__ - October 2025
 --------------------------------------------------------------------------------------------
 

src/sire/mol/_dynamics.py

@@ -972,6 +972,45 @@ def _rebuild_and_minimise(self):
             self._gcmc_sampler._set_water_state(self._omm_mols)
             self._gcmc_sampler.pop()
 
+        if self._save_crash_report:
+            import openmm
+            import numpy as np
+            from copy import deepcopy
+            from uuid import uuid4
+
+            # Create a unique identifier for this crash report.
+            crash_id = str(uuid4())[:8]
+
+            # Get the current context and system.
+            context = self._omm_mols
+            system = deepcopy(context.getSystem())
+
+            # Add each force to a unique group.
+            for i, f in enumerate(system.getForces()):
+                f.setForceGroup(i)
+
+            # Create a new context.
+            new_context = openmm.Context(system, deepcopy(context.getIntegrator()))
+            new_context.setPositions(context.getState(getPositions=True).getPositions())
+
+            # Write the  energies for each force group.
+            with open(f"crash_{crash_id}.log", "w") as f:
+                f.write(f"Current lambda: {str(self.get_lambda())}\n")
+                for i, force in enumerate(system.getForces()):
+                    state = new_context.getState(getEnergy=True, groups={i})
+                    f.write(f"{force.getName()}, {state.getPotentialEnergy()}\n")
+
+            # Save the serialised system.
+            with open(f"system_{crash_id}.xml", "w") as f:
+                f.write(openmm.XmlSerializer.serialize(system))
+
+            # Save the positions.
+            positions = (
+                new_context.getState(getPositions=True).getPositions(asNumpy=True)
+                / openmm.unit.nanometer
+            )
+            np.savetxt(f"positions_{crash_id}.txt", positions)
+
         self.run_minimisation()
 
     def run(
@@ -990,6 +1029,7 @@ def run(
         null_energy: str = None,
         excess_chemical_potential: float = None,
         num_waters: int = None,
+        save_crash_report: bool = False,
     ):
         if self.is_null():
             return
@@ -1009,6 +1049,7 @@ def run(
             "null_energy": null_energy,
             "excess_chemical_potential": excess_chemical_potential,
             "num_waters": num_waters,
+            "save_crash_report": save_crash_report,
         }
 
         from concurrent.futures import ThreadPoolExecutor
@@ -1048,6 +1089,13 @@ def run(
             except:
                 raise ValueError("'num_waters' must be an integer")
 
+        if save_crash_report is not None:
+            if not isinstance(save_crash_report, bool):
+                raise ValueError("'save_crash_report' must be True or False")
+            self._save_crash_report = save_crash_report
+        else:
+            self._save_crash_report = False
+
         try:
             steps_to_run = int(time.to(picosecond) / self.timestep().to(picosecond))
         except Exception:
@@ -1649,6 +1697,7 @@ def run(
         null_energy: str = None,
         excess_chemical_potential: str = None,
         num_waters: int = None,
+        save_crash_report: bool = False,
     ):
         """
         Perform dynamics on the molecules.
@@ -1722,6 +1771,13 @@ def run(
             Whether to save the energy on exit, regardless of whether
             the energy frequency has been reached.
 
+        save_crash_report: bool
+            Whether to save a crash report if the dynamics fails due to an
+            instability. This will save a named log file containing the energy
+            for each force, an XML file containing the OpenMM system at the
+            start of the dynamics block, and a NumPy text file containing
+            the atomic positions at the start of the dynamics block.
+
         auto_fix_minimise: bool
             Whether or not to automatically run minimisation if the
             trajectory exits with an error in the first few steps.
@@ -1768,6 +1824,14 @@ def run(
                 else:
                     save_velocities = False
 
+            if save_crash_report is None:
+                if self._d._map.specified("save_crash_report"):
+                    save_crash_report = (
+                        self._d._map["save_crash_report"].value().as_bool()
+                    )
+                else:
+                    save_crash_report = False
+
             self._d.run(
                 time=time,
                 save_frequency=save_frequency,
@@ -1778,6 +1842,7 @@ def run(
                 save_velocities=save_velocities,
                 save_frame_on_exit=save_frame_on_exit,
                 save_energy_on_exit=save_energy_on_exit,
+                save_crash_report=save_crash_report,
                 auto_fix_minimise=auto_fix_minimise,
                 num_energy_neighbours=num_energy_neighbours,
                 null_energy=null_energy,

tests/mol/test_dynamics.py

@@ -81,7 +81,7 @@ def test_cutoff_options(ala_mols):
     "openmm" not in sr.convert.supported_formats(),
     reason="openmm support is not available",
 )
-def test_sample_frequency(ala_mols):
+def test_sample_frequency(ala_mols, openmm_platform):
     """
     Test that energies and frames are saved at the correct frequency.
     """
@@ -92,7 +92,7 @@ def test_sample_frequency(ala_mols):
 
     mols = ala_mols
 
-    d = mols.dynamics(platform="Reference", timestep="1 fs")
+    d = mols.dynamics(platform=openmm_platform, timestep="1 fs")
 
     # Create a list of lambda windows.
     lambdas = [0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
@@ -145,3 +145,51 @@ def test_sample_frequency(ala_mols):
 
     # Check that the trajectory has 10 frames.
     assert new_mols.num_frames() == 10
+
+
+@pytest.mark.skipif(
+    "openmm" not in sr.convert.supported_formats(),
+    reason="openmm support is not available",
+)
+def test_crash_report(merged_ethane_methanol, openmm_platform):
+    """
+    Test that energies and frames are saved at the correct frequency.
+    """
+
+    import os
+    import glob
+    import tempfile
+    from sire.base import ProgressBar
+
+    ProgressBar.set_silent()
+
+    mols = merged_ethane_methanol.clone()
+    mols = sr.morph.link_to_reference(mols)
+
+    d = mols.dynamics(platform=openmm_platform)
+
+    # Run a short simulation within a temporary directory.
+    tmpdir = tempfile.TemporaryDirectory()
+
+    # Save the current directory.
+    old_dir = os.getcwd()
+
+    try:
+        # Change to the temporary directory.
+        os.chdir(tmpdir.name)
+
+        # Run a short simulation, forcing a crash.
+        d.run("1ps", save_crash_report=True)
+
+        # Glob for the crash report files.
+        crash_log = glob.glob("crash_*.log")
+        crash_system = glob.glob("system_*.xml")
+        crash_positions = glob.glob("positions_*.txt")
+
+        # Make sure we have one of each file.
+        assert len(crash_log) == 1
+        assert len(crash_system) == 1
+        assert len(crash_positions) == 1
+    finally:
+        # Change back to the old directory.
+        os.chdir(old_dir)