Fix issue #222

CHANGELOG.md

@@ -9,6 +9,7 @@ Changelog
 * Buffer energy components and write them at checkpoint time, rather than rewriting the parquet file on every energy save, which cost a few milliseconds per replica per cycle and grew with the length of the run [#212](https://github.com/OpenBioSim/somd2/pull/212).
 * Silence Sire's progress bars when a runner is constructed rather than when `somd2` is imported, so that importing `somd2` as a library no longer changes how Sire reports progress [#215](https://github.com/OpenBioSim/somd2/pull/215).
 * Save the replica exchange state once at the end of a run rather than twice when the last cycle is a checkpoint cycle, and include the GCMC statistics in the final save [#218](https://github.com/OpenBioSim/somd2/pull/218).
+* Fall back to the Aldeghi Boresch restraint search protocol when the default RXRX protocol can't be used, e.g. for ligands with no N/O atoms to act as hydrogen-bond partners. Topology failures are now detected before the restraint search trajectory is run [#223](https://github.com/OpenBioSim/somd2/pull/223).
 
 [2026.2.0](https://github.com/openbiosim/somd2/compare/2026.1.0...2026.2.0) - Sep 2026
 --------------------------------------------------------------------------------------

src/somd2/config/_config.py

@@ -634,6 +634,8 @@ def __init__(
             Sire selection string for receptor anchor atom candidates used
             during automatic Boresch restraint generation. If None, the default
             backbone selection is used (CA, C, N atoms in non-water molecules).
+            Only used by the Aldeghi protocol, which is the fallback when the
+            default RXRX protocol fails.
 
         morse_hard_well_depth: str
             The well depth of the "hard" Morse potential that replaces the

src/somd2/runner/_base.py

@@ -1155,7 +1155,7 @@ def _generate_boresch_restraint(self, device=None):
         crash on every restart, whereas a fresh search may pick a different frame
         or anchor, and re-seeds ``self._system`` naturally.
         """
-        from sire.restraints import boresch_search
+        from sire.restraints import boresch_search, check_boresch_search
 
         restraint_file = str(self._config.output_directory / "abfe_restraint.s3")
 
@@ -1191,6 +1191,16 @@ def _generate_boresch_restraint(self, device=None):
             "No restraint supplied for ABFE. Running Boresch restraint search."
         )
 
+        protocol = "rxrx"
+        try:
+            check_boresch_search(self._system, protocol=protocol)
+        except ValueError as e:
+            _logger.warning(
+                f"RXRX Boresch restraint search cannot be used for this system: {e} "
+                "Falling back to the Aldeghi protocol."
+            )
+            protocol = "aldeghi"
+
         search_system = self._system
 
         if self._config.minimise:
@@ -1235,22 +1245,39 @@ def _generate_boresch_restraint(self, device=None):
         )
         search_system = dynamics.commit()
 
-        search_kwargs = {"temperature": self._config.temperature}
+        # The restraint lever must match the one used by the ABFE schedules.
+        search_kwargs = {
+            "temperature": self._config.temperature,
+            "restraint_lever": "split",
+        }
         if self._config.restraint_search_receptor_selection is not None:
             search_kwargs["receptor_selection"] = (
                 self._config.restraint_search_receptor_selection
             )
 
-        restraints, correction, starting_structure = boresch_search(
-            search_system, **search_kwargs
-        )
+        try:
+            restraints, correction, starting_structure = boresch_search(
+                search_system, protocol=protocol, **search_kwargs
+            )
+        except ValueError as e:
+            if protocol != "rxrx":
+                raise
+            _logger.warning(
+                f"RXRX Boresch restraint search failed: {e} "
+                "Falling back to the Aldeghi protocol."
+            )
+            protocol = "aldeghi"
+            restraints, correction, starting_structure = boresch_search(
+                search_system, protocol=protocol, **search_kwargs
+            )
 
         # Cache so it can be written into the energy trajectory parquet
         # metadata (see _checkpoint), letting analysis code automatically
         # apply the correction without needing to scan the logs.
         self._standard_state_correction = float(correction.to(_sr.units.kcal_per_mol))
         _logger.info(
-            f"Boresch restraint generated. Standard state correction: "
+            f"Boresch restraint generated using the {protocol} protocol. "
+            "Standard state correction: "
             f"{self._standard_state_correction:.4f} kcal mol-1"
         )