Backport fix from PR #223

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"
         )