diff --git a/CHANGELOG.md b/CHANGELOG.md index 34661f4..4e927a0 100644 --- a/CHANGELOG.md +++ b/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 -------------------------------------------------------------------------------------- diff --git a/src/somd2/config/_config.py b/src/somd2/config/_config.py index 1d45d7b..d712e40 100644 --- a/src/somd2/config/_config.py +++ b/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 diff --git a/src/somd2/runner/_base.py b/src/somd2/runner/_base.py index 16fec0d..e79fc72 100644 --- a/src/somd2/runner/_base.py +++ b/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" )