diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index f30a236c3..3fea1f740 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -31,6 +31,8 @@ organisation on `GitHub `__. * Add support for ``openmm.MonteCarloMembraneBarostat`` in Sire-to-OpenMM conversion. +* Added support for saving crash reports during Sire dynamics runs. + `2025.2.0 `__ - October 2025 -------------------------------------------------------------------------------------------- diff --git a/src/sire/mol/_dynamics.py b/src/sire/mol/_dynamics.py index b18e186c3..7dd72f130 100644 --- a/src/sire/mol/_dynamics.py +++ b/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, diff --git a/tests/mol/test_dynamics.py b/tests/mol/test_dynamics.py index 7b86fad4a..bca606c38 100644 --- a/tests/mol/test_dynamics.py +++ b/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)