From 7d02f7020e44c10278a31a94d807af91641422b1 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 29 Jun 2026 15:39:26 +0100 Subject: [PATCH 01/23] Update Sire development pin. --- pixi.toml | 16 ++++++++++------ recipes/biosimspace/recipe.yaml | 4 ++-- 2 files changed, 12 insertions(+), 8 deletions(-) diff --git a/pixi.toml b/pixi.toml index 25747b305..d95194edf 100644 --- a/pixi.toml +++ b/pixi.toml @@ -25,9 +25,9 @@ rdkit = "*" [target.linux-64.dependencies] # main -sire = ">=2026.1.0,<2026.2.0" +#sire = ">=2026.1.0,<2026.2.0" # devel -#sire = "==2026.2.0.dev" +sire = "==2026.2.0.dev" ambertools = ">=22" gromacs = "*" alchemlyb = "*" @@ -36,13 +36,17 @@ mdanalysis = "*" [target.linux-aarch64.dependencies] +# main +#sire = ">=2026.1.0,<2026.2.0" +# devel +sire = "==2026.2.0.dev" gromacs = "*" [target.osx-arm64.dependencies] # main -#sire = ">=2025.4.0,<2025.5.0" +#sire = ">=2026.1.0,<2026.2.0" # devel -sire = "==2026.1.0.dev" +sire = "==2026.2.0.dev" ambertools = ">=22" alchemlyb = "*" mdtraj = "*" @@ -50,9 +54,9 @@ mdanalysis = "*" [target.win-64.dependencies] # main -#sire = ">=2025.4.0,<2025.5.0" +#sire = ">=2026.1.0,<2026.2.0" # devel -sire = "==2026.1.0.dev" +sire = "==2026.2.0.dev" alchemlyb = "*" mdtraj = "*" mdanalysis = "*" diff --git a/recipes/biosimspace/recipe.yaml b/recipes/biosimspace/recipe.yaml index abdf4e06a..abdcd7712 100644 --- a/recipes/biosimspace/recipe.yaml +++ b/recipes/biosimspace/recipe.yaml @@ -43,9 +43,9 @@ requirements: - pyyaml - rdkit # main - - sire >=2026.1.0,<2026.2.0 + #- sire >=2026.1.0,<2026.2.0 # devel - #- sire ==2026.1.0.dev + - sire ==2026.2.0.dev - if: not aarch64 then: - alchemlyb From c7b1d50cb5485dfee5e36810f28cc747abfc0b02 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 29 Jun 2026 16:52:15 +0100 Subject: [PATCH 02/23] Add note regarding PyMBAR JAX issues. --- README.rst | 3 +++ 1 file changed, 3 insertions(+) diff --git a/README.rst b/README.rst index 7bf08bd71..47030af89 100644 --- a/README.rst +++ b/README.rst @@ -184,3 +184,6 @@ along with the BioSimSpace version number. This can be found by running: import BioSimSpace as BSS print(BSS.__version__) + +* If you experience ``JAX`` issues when using ``BioSimSpace.FreeEnergy.Relative.analyse``, try + setting the ``PYMBAR_DISABLE_JAX`` environment variable to ``1``. From 8b86bac65574942f20ba1fd1b7ecf2c28a579aaa Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 13 Jul 2026 16:33:51 +0100 Subject: [PATCH 03/23] Add script to autogenerate roadmap page for website. [ci skip] --- doc/Makefile | 8 +- doc/generate_roadmap.py | 347 +++++++++++++++++++++++++++++++++++++++ doc/source/changelog.rst | 68 ++++++++ doc/source/index.rst | 8 + 4 files changed, 429 insertions(+), 2 deletions(-) create mode 100644 doc/generate_roadmap.py diff --git a/doc/Makefile b/doc/Makefile index 175004fb0..8d67d3d53 100644 --- a/doc/Makefile +++ b/doc/Makefile @@ -16,7 +16,7 @@ PAPEROPT_a4 = -D latex_paper_size=a4 PAPEROPT_letter = -D latex_paper_size=letter ALLSPHINXOPTS = -d $(BUILD_DIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) source -.PHONY: help clean html +.PHONY: help clean html roadmap #------------------------------------------------------------------------------ all: html @@ -24,11 +24,15 @@ help: @echo "Please use \`make ' where is one of" @echo " html to make standalone HTML files" @echo " latex to make LaTeX files, you can set PAPER=a4 or PAPER=letter" + @echo " roadmap to regenerate source/roadmap.rst from changelog.rst" clean: -rm -rf $(BUILD_DIR)/* source/generated source/api/generated -html: +roadmap: + python3 generate_roadmap.py + +html: roadmap mkdir -p $(BUILD_DIR)/html $(BUILD_DIR)/doctrees $(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILD_DIR)/html $(FILES) @echo diff --git a/doc/generate_roadmap.py b/doc/generate_roadmap.py new file mode 100644 index 000000000..8c1d49d9f --- /dev/null +++ b/doc/generate_roadmap.py @@ -0,0 +1,347 @@ +#!/usr/bin/env python3 +"""Generate source/roadmap.rst from two sources: + +* Shipped: entries tagged ".. roadmap::" in source/changelog.rst. The + CHANGELOG is the single source of truth for what shipped and when - to + make a bullet appear here, add a ".. roadmap::" comment directly beneath + it (see changelog.rst for existing examples). +* Planned: open GitHub issues labelled "roadmap", fetched from the public + GitHub REST API (no auth needed, no ``gh`` CLI required). + +Nothing else needs to change - re-run this script and rebuild the docs. + +Usage: + python generate_roadmap.py +""" + +import json +import re +import urllib.error +import urllib.request +from pathlib import Path + +SOURCE = Path(__file__).parent / "source" / "changelog.rst" +OUTPUT = Path(__file__).parent / "source" / "roadmap.rst" + +GITHUB_REPO = "openbiosim/biosimspace" +GITHUB_LABEL = "roadmap" + +VERSION_LINKED_RE = re.compile(r"^`(?P[\w.]+) <(?P[^>]+)>`_ - (?P.+)$") +VERSION_PLAIN_RE = re.compile(r"^(?P\d[\w.]*) - (?P.+)$") +UNDERLINE_RE = re.compile(r"^-{5,}\s*$") + +ATTRIBUTION_RE = re.compile(r"\(`@[\w-]+[^)]*\)") +PR_RE = re.compile(r"\(`#(?P\d+) <(?P[^>]+)>`__\)\.?\s*$") +ROLE_RE = re.compile(r":\w+:`([^<`]+)(?: <[^>]+>)?`") +LITERAL_RE = re.compile(r"``([^`]+)``") +LINK_RE = re.compile(r"`([^`<]+) <([^>]+)>`_+") +BARE_INTERPRETED_RE = re.compile(r"`([^`<>]+)`") + + +def clean_text(text): + """Strip RST markup from a changelog bullet and pull out its PR link.""" + pr = None + match = PR_RE.search(text) + if match: + pr = (match.group("num"), match.group("url")) + text = text[: match.start()].rstrip() + + text = ATTRIBUTION_RE.sub("", text) + text = ROLE_RE.sub(r"\1", text) + text = LITERAL_RE.sub(r"\1", text) + text = LINK_RE.sub(lambda m: f'{m.group(1)}', text) + text = BARE_INTERPRETED_RE.sub(r"\1", text) + text = re.sub(r"\s+", " ", text).strip() + text = text.rstrip(" .") + return text, pr + + +def parse_changelog(path): + """Return a list of {version, url, date, items} dicts, newest first.""" + lines = path.read_text().splitlines() + releases = [] + current = None + i, n = 0, len(lines) + + while i < n: + line = lines[i] + m = VERSION_LINKED_RE.match(line) + m_plain = None if m else VERSION_PLAIN_RE.match(line) + + if (m or m_plain) and i + 1 < n and UNDERLINE_RE.match(lines[i + 1]): + group = m or m_plain + current = { + "version": group.group("version"), + "url": group.group("url") if m else None, + "date": group.group("date"), + "items": [], + } + releases.append(current) + i += 2 + continue + + if current is not None and line.startswith("* "): + bullet_lines = [line[2:]] + j = i + 1 + while j < n and lines[j].strip() and not lines[j].startswith(("* ", "..")) \ + and not (UNDERLINE_RE.match(lines[j])): + bullet_lines.append(lines[j].strip()) + j += 1 + bullet_text = " ".join(bullet_lines).strip() + + k = j + while k < n and lines[k].strip() == "": + k += 1 + tagged = k < n and lines[k].strip() == ".. roadmap::" + if tagged: + text, pr = clean_text(bullet_text) + current["items"].append({"text": text, "pr": pr}) + k += 1 + i = k + continue + + i += 1 + + return [r for r in releases if r["items"]] + + +def fetch_planned_issues(): + """Fetch open issues labelled GITHUB_LABEL via the public GitHub API. + + Returns an empty list (with a warning) rather than raising, so a + network hiccup or rate limit never blocks a docs build. + """ + url = ( + f"https://api.github.com/repos/{GITHUB_REPO}/issues" + f"?labels={GITHUB_LABEL}&state=open&per_page=100" + ) + request = urllib.request.Request( + url, + headers={ + "User-Agent": "biosimspace-roadmap-generator", + "Accept": "application/vnd.github+json", + }, + ) + try: + with urllib.request.urlopen(request, timeout=10) as response: + data = json.load(response) + except (urllib.error.URLError, urllib.error.HTTPError, TimeoutError, OSError) as exc: + print(f"Warning: could not fetch planned roadmap issues from GitHub ({exc}); leaving section empty.") + return [] + + # The issues endpoint also returns pull requests - exclude those. + return [ + {"number": item["number"], "title": item["title"], "url": item["html_url"]} + for item in data + if "pull_request" not in item + ] + + +def render_shipped_html(releases): + rows = [] + for release in releases: + if release["url"]: + heading = f'{release["version"]}' + else: + heading = release["version"] + + items_html = [] + for item in release["items"]: + pr_link = "" + if item["pr"]: + num, url = item["pr"] + pr_link = f' #{num}' + items_html.append(f'
  • {item["text"]}{pr_link}
    • ') + + rows.append(f""" +
    + +
    +
    + {heading} + {release["date"]} +
    +
      + {"".join(items_html)} +
    +
    +
    """) + + return f"""
    +{"".join(rows)} +
    +""" + + +def render_planned_html(issues): + issues_url = ( + f"https://github.com/{GITHUB_REPO}/issues" + f"?q=is%3Aissue+is%3Aopen+label%3A{GITHUB_LABEL}" + ) + if not issues: + return f"""
    +

    + Nothing tagged yet - see the + open issues labelled "{GITHUB_LABEL}". +

    +
    +""" + + items_html = "".join( + f'
  • {issue["title"]} ' + f'#{issue["number"]}
    • ' + for issue in issues + ) + return f"""
    +
      + {items_html} +
    +

    + Tracked via issues labelled + "{GITHUB_LABEL}". +

    +
    +""" + + +HEADER = """.. THIS FILE IS AUTO-GENERATED - DO NOT EDIT BY HAND. +.. Run ``python generate_roadmap.py`` from the ``doc`` directory to +.. regenerate it from changelog.rst and open GitHub issues. + +Roadmap +======= + +A sparse, curated timeline of the capabilities added to :mod:`BioSimSpace` +over time, plus what's currently planned. This isn't every change - see the +:doc:`changelog` for the full, linear history of every fix and tweak - just +the entries that mark a new piece of scientific or engineering capability. + +Planned +------- + +""" + +SHIPPED_HEADER = """ +Shipped +------- + +""" + +STYLE = """ +.. raw:: html + + +""" + + +def indent(text, prefix=" "): + return "\n".join(prefix + line if line else line for line in text.splitlines()) + + +def as_raw_block(html): + return ".. raw:: html\n\n" + indent(html) + "\n" + + +def main(): + releases = parse_changelog(SOURCE) + planned = fetch_planned_issues() + + parts = [ + HEADER, + as_raw_block(render_planned_html(planned)), + SHIPPED_HEADER, + as_raw_block(render_shipped_html(releases)), + STYLE, + ] + OUTPUT.write_text("\n".join(parts)) + print( + f"Wrote {len(planned)} planned issue(s), " + f"{len(releases)} release(s) ({sum(len(r['items']) for r in releases)} shipped entries) " + f"to {OUTPUT}" + ) + + +if __name__ == "__main__": + main() diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 502684d19..469e759e9 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -19,8 +19,12 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Update the ``prepareFEP`` helper script to generate ``SOMD1`` and ``SOMD2`` inputs simultaneously (`#516 `__). * Remove cross-bond angle and torsion terms for ring-making/breaking perturbations (`#517 `__). * Added functionaltiy for adding ions to an existing system (`#518 `__). + + .. roadmap:: * Expose the ``max_path`` and ``max_ring_size`` kwargs used for ring break and size change detection in the :func:`BioSimSpace.Align.generateNetwork ` function (`#520 `__). * Add support for user-defined "region-of-interest" merges (`#522 `__). + + .. roadmap:: * Add pins to handle OpenForceField Python version deprecations (`#524 `__). `2025.4.0 `_ - Feb 17 2026 @@ -28,6 +32,8 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Fixed centre of mass restraints for alchemical transfer method (ATM) simulations (`@mb2055 `__) (`#471 `__). * Add experimental :class:`ReplicaSystem` class to speed up handling of replica exchange simulations (`#473 `__). + + .. roadmap:: * Removed lazy imports from sub-modules that don't use Sire (`#475 `__). * Allow translation of a custom ``coordinates`` property for perturbable molecules (`#477 `__). * Added a kwarg to make zeroing of LJ sigma values for ghost atoms optional (`#482 `__). @@ -47,13 +53,20 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Added functionality for quickly getting and setting the coordinates array of a :class:`System ` (`#465 `__). * Reduce depdency import overheads by switching from module level imports to function and method level (`#466 `__). * Fully switch over to using the new Sire Python API naming convention, allowing BioSimSpace to be used within Sire (`#466 `__). + + .. roadmap:: * Fixed logic used for setting the GPU device index for :class:`Process.OpenMM ` (`#468 `__). `2025.2.0 `_ - Oct 08 2025 ------------------------------------------------------------------------------------------------- * Add support for ``SOMD2`` FEP analysis using data frames with different numbers of samples (`#415 `__). + + .. roadmap:: + * Add support for the ABCG2 charge method for GAFF parameterisation (`#421 `__). + + .. roadmap:: * Fixed f-string formatting error in FEP analysis exception message (`#423 `__). * Fixed FEP energy trajectory slicing when intitial sample time is non-zero (`#424 `__). * Add workaround for incompatibility between ``ParmEd`` and ``NumPy`` 2.3 (`#428 `__). @@ -73,9 +86,14 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Improved robustness of formal charge inference when reading molecules from PDB or SDF files (`#393 `__). * Make sure the system extracted from AMBER trajectory frames during free-energy perturbation simulations is in the original, unsquashed format (`#403 `__). * Add support for the ``ff19SB`` force field and OPC water (`#406 `__). + + .. roadmap:: + * Allow creation of ``SOMD`` perturbation files without modification to ghost atom bonded terms (`#407 `__). * Support analysis of ``SOMD2`` energy trajectories with time varying lambda sampling (`#408 `__). + .. roadmap:: + `2024.4.1 `_ - Feb 14 2025 ------------------------------------------------------------------------------------------------- @@ -101,6 +119,8 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Fixed alignment of monatomic molecules (`#313 `__ and (`#333 `__). * Expose missing ``extra_args`` keyword argument for the :class:`Process.Somd ` class (`#319 `__). * Add support for the Alchemical Transfer Method (ATM) (`@mb2055 `_) (`#327 `__). + + .. roadmap:: * Fixed :meth:`system.updateMolecules ` method when updating multiple molecules at once (`#336 `__). * Added a new :meth:`system.removeBox ` method to remove the box from a system (`#338 `__). * Fixed bug when using position restraints with OpenMM for perturbable systems (`#341 `__). @@ -118,6 +138,8 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Detect dummy atoms by checking ``element`` *and* ``ambertype`` properties when creating ``SOMD`` pert files (`#289 `__). * Add missing ``match_water`` kwarg to ``prepareFEP`` node (`#292 `__). * Add protein free-energy perturbation functionality (`@akalpokas `__). + + .. roadmap:: * Ensure that the LJ sigma parameter for perturbed atoms is non-zero (`#295 `__). * Fixed return type docstrings for functions in the :mod:`BioSimSpace.Parameters` module (`#298 `__). * Don't use ``sire.legacy.Base.wrap`` with the ``file_format`` property to avoid (incorrect) auto string to unit conversion of ``mol2`` to moles squared (`#300 `__). @@ -129,8 +151,13 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Switch to using Langevin integrator for GROMACS free energy simulations (`#264 `__). * Add support for clearing and disabling the IO file cache (`#266 `__). * Add support for using ``openff-nagl`` to generate partial charges (`#267 `__). + + .. roadmap:: + * Fixed non-reproducible search for backbone restraint atom indices (`#270 `__). * Add support for AMBER as an alchemical free-energy simulation engine (`#272 `__). + + .. roadmap:: * Switch to using ``os.path.join`` to generate directory file names (`#276 `__). * Make sure the ``fileformat`` property is preserved when creating single molecule systems (`#276 `__). * Add a ``getRestraintType`` method to the base protocol that returns ``None`` (`#276 `__). @@ -155,6 +182,8 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Add support for detecting nucleic acid backbones (`@fjclark `_) (`#189 `__). * Added SOMD and GROMACS support for multiple distance restraints for ABFE calculations (`#178 `__). + .. roadmap:: + `2023.4.1 `_ - Dec 14 2023 ------------------------------------------------------------------------------------------------- @@ -178,7 +207,12 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Add support for computing trajectory RMSDs using Sire backend (`#152 `__). * Add support for setting up systems containing crystal waters (`#154 `__). + + .. roadmap:: + * Add unified free-energy perturbation analysis using ``alchemlyb`` (`@annamherz `_) (`#155 `__). + + .. roadmap:: * Fix handling of connectivity changes during molecular perturbations (`#157 `__). * Fix issues related to new shared properties in Sire (`#160 `__). * Fix issues in SOMD perturbation files for absolute binding free-energy simulations (`@fjclark `_) (`#164 `__). @@ -214,6 +248,8 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Recenter molecules following vacuum simulation with GROMACS to avoid precision overflow with molecular coordinates on write (`#95 `__). * Fix expected angles used in unit test following updates to triclinic box code in Sire (`#99 `__). * Add absolute binding free-energy support for SOMD (`@fjclark `_) (`#104 `__). + + .. roadmap:: * Avoid streaming issues when reading binary AMBER restart files for a single frame (`#105 `__). * Improve overlap matrix plotting functionality (`@fjclark `_) (`#107 `__). * Handle updates to Sire parser format naming (`#108 `__). @@ -228,6 +264,8 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Catch exception when vacuum system has a cartesian space (`#120 `__). * Add support for Sire as a trajectory backend (`#121 `__). + .. roadmap:: + `2023.2.2 `_ - May 15 2023 ------------------------------------------------------------------------------------------------- @@ -274,6 +312,8 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Refactor code to use a unified :class:`WorkDir ` class to simplify the creation of working directories (`#2 `__). * Added :meth:`isSame ` method to compare systems using a sub-set of system and molecular properties. This improves our file caching support, allowing a user to exclude properties when comparing cached systems prior to write, e.g. ignoring coordinates and velocities, if those are the only things that differ between the systems `(#3 `__). * Added the initial version of :mod:`BioSimSpace.Convert `, which provides support for converting between native `BioSimSpace`, `Sire `__, and `RDKit `__ objects (`#9 `__). + + .. roadmap:: * Fixed several formatting issues with the website documentation. `2023.1.1 `_ - Feb 07 2023 @@ -289,6 +329,8 @@ within the biomolecular simulation community. Our software is hosted via the `Op * Wrapped the new `sire.load `__ function to allow loading of URLs. * Add basic file caching support to avoid re-writing files for the same molecular system. * Added :data:`BioSimSpace._Config` sub-package to simplify the generation of configuration files for molecular dynamics engines and improve flexiblity. (Adapted from code written by `@msuruzhon `_.) + + .. roadmap:: * Deprecated ``BioSimSpace.IO.glob`` since globbing is now performed automatically. * Autoformat entire codebase using `black `__. * Fix issues following Sire 2023 API updates. @@ -308,6 +350,8 @@ GitHub organisation. The following releases were made during that time. * Added wrapper for ``Sire.Units.GeneralUnit``. * Improved interoperability of ``BioSimSpace.Trajectory`` sub-package. * Added ``BioSimSpace.Sandpit`` for experimental features from external collaborators. + + .. roadmap:: * Added functionality to check for molecules in a ``BioSimSpace.System``. * Added functionality to extract atoms and residues by absolute index. * Allow continuation for GROMACS equilibration simulations. (`@kexul `_) @@ -338,10 +382,21 @@ GitHub organisation. The following releases were made during that time. * Added basic support for cleaning PDB files with `pdb4amber `_ prior to read. * Added basic support for exporting BioSimSpace Nodes as Common Workflow Language wrappers. * Added support for parameterising molecules using OpenForceField. + + .. roadmap:: + * Added support for using SMILES strings for input to parameterisation functions. * Added support for funnel metadynamics simulations (`@dlukauskis `_). + + .. roadmap:: + * Added support for steered molecular dynamics simulations (`@AdeleLip `_). + + .. roadmap:: + * Added support for generating perturbation networks using LOMAP (`@JenkeScheen `_). + + .. roadmap:: * Fixed bug affecting certain improper/dihedral terms in SOMD perturbation file writer. * Numerous performance improvements, particularly involving the manipulation and combination of molecular systems. @@ -376,6 +431,8 @@ GitHub organisation. The following releases were made during that time. * Switched to using `RDKit `_ for maximum common substructure (MCS) mappings. * Handle perturbable molecules for non free-energy protocols with SOMD and GROMACS. * Added basic metadynamics functionality with support for distance and torsion collective variables. + + .. roadmap:: * Added support for inferring formal charge of molecules. * Numerous MCS mapping fixes and improvements. Thanks to `@maxkuhn `_, `@dlukauskis `_, and `@ptosco `_ for help testing and debugging. * Added Dockerfile to build thirdparty packages required by the BioSimSpace notebook server. @@ -387,9 +444,20 @@ GitHub organisation. The following releases were made during that time. ------------------------------------------------------------------------------------------------- * Added support for parameterising proteins and ligands. + + .. roadmap:: + * Added support for solvating molecular systems. + + .. roadmap:: + * Molecular dynamics drivers updated to support SOMD and GROMACS. + + .. roadmap:: + * Support free energy perturbation simulations with SOMD and GROMACS. + + .. roadmap:: * Added Azure Pipeline to automatically build, test, document, and deploy BioSimSpace. * Created automatic Conda package pipeline. diff --git a/doc/source/index.rst b/doc/source/index.rst index ccc1cae59..f30f9db15 100644 --- a/doc/source/index.rst +++ b/doc/source/index.rst @@ -105,6 +105,14 @@ Contributing code_of_conduct +Roadmap +======= + +.. toctree:: + :maxdepth: 1 + + roadmap + Changelog ========= From abbfb7b8f3e893ca2a0c2a2b1f2b58362a336087 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 28 Jul 2026 11:35:30 +0100 Subject: [PATCH 04/23] Use pip_check: false to avoid NumPy pin errors with AmberTools. --- recipes/biosimspace/recipe.yaml | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/recipes/biosimspace/recipe.yaml b/recipes/biosimspace/recipe.yaml index abdcd7712..964f5097d 100644 --- a/recipes/biosimspace/recipe.yaml +++ b/recipes/biosimspace/recipe.yaml @@ -56,6 +56,10 @@ tests: - python: imports: - BioSimSpace + # AmberTools can be pulled in transitively, and installs tools into + # site-packages whose metadata still pins numpy <2, so 'pip check' + # fails even though nothing here uses them. + pip_check: false - script: - if: unix then: PYTHONPATH=. pytest -vvv --color=yes --import-mode=importlib ./tests From 98fe610533a5b1a7258c1ddec2b3a8d8f592ef75 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Jul 2026 11:24:02 +0100 Subject: [PATCH 05/23] Apply translation to correct object. [closes #539] --- .../Exscientia/_SireWrappers/_molecule.py | 4 +- src/BioSimSpace/_SireWrappers/_molecule.py | 4 +- .../Exscientia/_SireWrappers/test_molecule.py | 51 +++++++++++++++++++ tests/Sandpit/Exscientia/conftest.py | 11 ++++ tests/_SireWrappers/test_molecule.py | 51 +++++++++++++++++++ 5 files changed, 115 insertions(+), 6 deletions(-) diff --git a/src/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py b/src/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py index 654be76bf..e421d0bd6 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py +++ b/src/BioSimSpace/Sandpit/Exscientia/_SireWrappers/_molecule.py @@ -1472,9 +1472,7 @@ def translate(self, vector, property_map={}): if self._sire_object.has_property(coord_prop): _property_map["coordinates"] = coord_prop mol = ( - self._sire_object.move() - .translate(_SireMaths.Vector(vec), _property_map) - .commit() + mol.move().translate(_SireMaths.Vector(vec), _property_map).commit() ) else: diff --git a/src/BioSimSpace/_SireWrappers/_molecule.py b/src/BioSimSpace/_SireWrappers/_molecule.py index 75880eddb..662cb252a 100644 --- a/src/BioSimSpace/_SireWrappers/_molecule.py +++ b/src/BioSimSpace/_SireWrappers/_molecule.py @@ -1412,9 +1412,7 @@ def translate(self, vector, property_map={}): if self._sire_object.has_property(coord_prop): _property_map["coordinates"] = coord_prop mol = ( - self._sire_object.move() - .translate(_SireMaths.Vector(vec), _property_map) - .commit() + mol.move().translate(_SireMaths.Vector(vec), _property_map).commit() ) else: diff --git a/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py b/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py index 89ae620a4..eef91f38e 100644 --- a/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py +++ b/tests/Sandpit/Exscientia/_SireWrappers/test_molecule.py @@ -173,3 +173,54 @@ def test_makeCompatibleWith_regression(): # Make sure the energies are approximately equal. assert nrg_compatible == pytest.approx(nrg_leap, rel=1e-5) + + +def test_translate_custom_coordinates_property(perturbable_system): + """ + Regresson test to ensure that a molecule with a custom + coordinates property gets translated correctly. + """ + + from sire.maths import Vector + + # Extract a copy of the perturbable molecule. + mol = perturbable_system.getPerturbableMolecules()[0].copy() + + # Copy the "coordinates0" property to a custom property installed + # "test" + cursor = mol._sire_object.cursor() + cursor["test"] = cursor["coordinates0"] + mol._sire_object = cursor.commit() + + # Store the existing coordinates. + coords0 = mol._sire_object.property("coordinates0").to_vector() + coords1 = mol._sire_object.property("coordinates1").to_vector() + coords_test = mol._sire_object.property("test").to_vector() + + # Translate the molecule. + mol.translate(3 * [BSS.Units.Length.angstrom], property_map={"coordinates": "test"}) + + # Get the new coordinates. + new_coords0 = mol._sire_object.property("coordinates0").to_vector() + new_coords1 = mol._sire_object.property("coordinates1").to_vector() + new_coords_test = mol._sire_object.property("test").to_vector() + + # Create a vector of the displacement. + v_ref = Vector(1.0, 1.0, 1.0) + + # Check that the coordinates have been translated correctly. + for c_new, c_old in zip(new_coords0, coords0): + v = c_new - c_old + assert pytest.approx(v.x().value(), abs=1e-5) == v_ref.x().value() + assert pytest.approx(v.y().value(), abs=1e-5) == v_ref.y().value() + assert pytest.approx(v.z().value(), abs=1e-5) == v_ref.z().value() + for c_new, c_old in zip(new_coords1, coords1): + v = c_new - c_old + assert pytest.approx(v.x().value(), abs=1e-5) == v_ref.x().value() + assert pytest.approx(v.y().value(), abs=1e-5) == v_ref.y().value() + assert pytest.approx(v.z().value(), abs=1e-5) == v_ref.z().value() + for c_new, c_old in zip(new_coords_test, coords_test): + v = c_new - c_old + assert pytest.approx(v.x().value(), abs=1e-5) == v_ref.x().value() + assert pytest.approx(v.y().value(), abs=1e-5) == v_ref.y().value() + assert pytest.approx(v.z().value(), abs=1e-5) == v_ref.z().value() diff --git a/tests/Sandpit/Exscientia/conftest.py b/tests/Sandpit/Exscientia/conftest.py index cc4cebb38..82268b3b3 100644 --- a/tests/Sandpit/Exscientia/conftest.py +++ b/tests/Sandpit/Exscientia/conftest.py @@ -169,3 +169,14 @@ def merged_benzene_pyrrole(benzene, pyrrole, mapping_benzene_pyrrole): allow_ring_breaking=True, allow_ring_size_change=True, ) + + +@pytest.fixture(scope="module") +def perturbable_system(): + """A vacuum perturbable system.""" + return BSS.IO.readPerturbableSystem( + f"{url}/perturbable_system0.prm7", + f"{url}/perturbable_system0.rst7", + f"{url}/perturbable_system1.prm7", + f"{url}/perturbable_system1.rst7", + ) diff --git a/tests/_SireWrappers/test_molecule.py b/tests/_SireWrappers/test_molecule.py index bb447ec6e..2b4dbf16e 100644 --- a/tests/_SireWrappers/test_molecule.py +++ b/tests/_SireWrappers/test_molecule.py @@ -145,3 +145,54 @@ def test_makeCompatibleWith_regression(): # Make sure the energies are approximately equal. assert nrg_compatible == pytest.approx(nrg_leap, rel=1e-5) + + +def test_translate_custom_coordinates_property(perturbable_system): + """ + Regresson test to ensure that a molecule with a custom + coordinates property gets translated correctly. + """ + + from sire.maths import Vector + + # Extract a copy of the perturbable molecule. + mol = perturbable_system.getPerturbableMolecules()[0].copy() + + # Copy the "coordinates0" property to a custom property installed + # "test" + cursor = mol._sire_object.cursor() + cursor["test"] = cursor["coordinates0"] + mol._sire_object = cursor.commit() + + # Store the existing coordinates. + coords0 = mol._sire_object.property("coordinates0").to_vector() + coords1 = mol._sire_object.property("coordinates1").to_vector() + coords_test = mol._sire_object.property("test").to_vector() + + # Translate the molecule. + mol.translate(3 * [BSS.Units.Length.angstrom], property_map={"coordinates": "test"}) + + # Get the new coordinates. + new_coords0 = mol._sire_object.property("coordinates0").to_vector() + new_coords1 = mol._sire_object.property("coordinates1").to_vector() + new_coords_test = mol._sire_object.property("test").to_vector() + + # Create a vector of the displacement. + v_ref = Vector(1.0, 1.0, 1.0) + + # Check that the coordinates have been translated correctly. + for c_new, c_old in zip(new_coords0, coords0): + v = c_new - c_old + assert pytest.approx(v.x().value(), abs=1e-5) == v_ref.x().value() + assert pytest.approx(v.y().value(), abs=1e-5) == v_ref.y().value() + assert pytest.approx(v.z().value(), abs=1e-5) == v_ref.z().value() + for c_new, c_old in zip(new_coords1, coords1): + v = c_new - c_old + assert pytest.approx(v.x().value(), abs=1e-5) == v_ref.x().value() + assert pytest.approx(v.y().value(), abs=1e-5) == v_ref.y().value() + assert pytest.approx(v.z().value(), abs=1e-5) == v_ref.z().value() + for c_new, c_old in zip(new_coords_test, coords_test): + v = c_new - c_old + assert pytest.approx(v.x().value(), abs=1e-5) == v_ref.x().value() + assert pytest.approx(v.y().value(), abs=1e-5) == v_ref.y().value() + assert pytest.approx(v.z().value(), abs=1e-5) == v_ref.z().value() From b8760dd25ba37254df5b13f4410550a8f306ab19 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Jul 2026 17:01:40 +0100 Subject: [PATCH 06/23] Use two-atom in_ring so unrelated rings don't flag a ring break. --- src/BioSimSpace/Align/_merge.py | 8 ++++---- src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/BioSimSpace/Align/_merge.py b/src/BioSimSpace/Align/_merge.py index 24b50649e..cac08e0e2 100644 --- a/src/BioSimSpace/Align/_merge.py +++ b/src/BioSimSpace/Align/_merge.py @@ -1589,10 +1589,10 @@ def _check_ring(conn0, conn1, idx0, idy0, idx1, idy1, max_path=50, max_ring_size # Supplementary check for rings larger than max_path: find_paths may only # find the direct-bond path and miss the long way around the ring, giving # n=1 instead of n≥2. Sire's in_ring has no path-length limit and - # correctly identifies ring membership in macrocycles. - if (conn0.in_ring(idx0) and conn0.in_ring(idy0)) != ( - conn1.in_ring(idx1) and conn1.in_ring(idy1) - ): + # correctly identifies ring membership in macrocycles. The two-atom + # overload asks whether the atoms share a ring, so a ring built entirely + # from dummy atoms, which breaks no bond between mapped atoms, is ignored. + if conn0.in_ring(idx0, idy0) != conn1.in_ring(idx1, idy1): return True, False # A direct bond was replaced by a ring path (or vice versa), leaving the diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py index c4102611f..f14002dca 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_merge.py @@ -1455,10 +1455,10 @@ def _check_ring(conn0, conn1, idx0, idy0, idx1, idy1, max_path=50, max_ring_size # Supplementary check for rings larger than max_path: find_paths may only # find the direct-bond path and miss the long way around the ring, giving # n=1 instead of n≥2. Sire's in_ring has no path-length limit and - # correctly identifies ring membership in macrocycles. - if (conn0.in_ring(idx0) and conn0.in_ring(idy0)) != ( - conn1.in_ring(idx1) and conn1.in_ring(idy1) - ): + # correctly identifies ring membership in macrocycles. The two-atom + # overload asks whether the atoms share a ring, so a ring built entirely + # from dummy atoms, which breaks no bond between mapped atoms, is ignored. + if conn0.in_ring(idx0, idy0) != conn1.in_ring(idx1, idy1): return True, False # A direct bond was replaced by a ring path (or vice versa), leaving the From 646fbef39fc4072cb3f28313dc5427ee73590970 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Jul 2026 17:21:07 +0100 Subject: [PATCH 07/23] Add mcs_kwargs to allow the user to configure the RDKit MCS. --- src/BioSimSpace/Align/_align.py | 68 +++++++++++++++---- .../Sandpit/Exscientia/Align/_align.py | 63 +++++++++++++---- 2 files changed, 105 insertions(+), 26 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 2602bb5cc..666cfe44d 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -25,6 +25,7 @@ __email__ = "lester.hedges@gmail.com" __all__ = [ + "defaultMCSOptions", "generateNetwork", "matchAtoms", "viewMapping", @@ -702,6 +703,29 @@ def generateNetwork( return edges, scores +def defaultMCSOptions(): + """ + Return the default options used for the RDKit maximum common substructure + search. These can be overridden using the 'mcs_kwargs' argument of + :class:`matchAtoms `. + + Returns + ------- + + options : dict + The default RDKit MCS options. + """ + return { + "atomCompare": _rdFMCS.AtomCompare.CompareAny, + "bondCompare": _rdFMCS.BondCompare.CompareAny, + "completeRingsOnly": True, + "ringMatchesRingOnly": True, + "matchChiralTag": False, + "matchValences": False, + "maximizeBonds": False, + } + + def matchAtoms( molecule0, molecule1, @@ -719,6 +743,7 @@ def matchAtoms( prune_atom_types=False, property_map0={}, property_map1={}, + mcs_kwargs={}, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -775,6 +800,11 @@ def matchAtoms( option is only relevant to MCS performed using RDKit and will be ignored when falling back on Sire. + mcs_kwargs : dict + A dictionary of keyword arguments used to override the defaults + passed to the RDKit MCS search. This option is only relevant to MCS + performed using RDKit and will be ignored when falling back on Sire. + roi : list The region of interest to match. Consists of a list of ROI residue indices. @@ -881,6 +911,7 @@ def matchAtoms( prune_atom_types=prune_atom_types, property_map0=property_map0, property_map1=property_map1, + mcs_kwargs=mcs_kwargs, ) else: return _roiMatch( @@ -912,6 +943,7 @@ def _matchAtoms( prune_atom_types=False, property_map0={}, property_map1={}, + mcs_kwargs={}, ): import sys as _sys @@ -975,6 +1007,9 @@ def _matchAtoms( if not isinstance(complete_rings_only, bool): raise TypeError("'complete_rings_only' must be of type 'bool'") + if not isinstance(mcs_kwargs, dict): + raise TypeError("'mcs_kwargs' must be of type 'dict'") + if type(max_scoring_matches) is not int: raise TypeError("'max_scoring_matches' must be of type 'int'") @@ -1012,24 +1047,21 @@ def _matchAtoms( _Convert.toRDKit(mol1, property_map=property_map1), ] + # Default MCS options, overridden by anything in 'mcs_kwargs'. The + # timeout is applied last so that it can't be overridden. + mcs_options = defaultMCSOptions() + mcs_options["completeRingsOnly"] = complete_rings_only + mcs_options.update(mcs_kwargs) + mcs_options["timeout"] = timeout + # Generate the MCS match. - mcs = _rdFMCS.FindMCS( - mols, - atomCompare=_rdFMCS.AtomCompare.CompareAny, - bondCompare=_rdFMCS.BondCompare.CompareAny, - completeRingsOnly=complete_rings_only, - ringMatchesRingOnly=True, - matchChiralTag=False, - matchValences=False, - maximizeBonds=False, - timeout=timeout, - ) + mcs = _rdFMCS.FindMCS(mols, **mcs_options) # Get the common substructure as a SMARTS string. mcs_smarts = _Chem.MolFromSmarts(mcs.smartsString) - except: - raise RuntimeError("RDKit MCS mapping failed!") + except Exception as e: + raise RuntimeError(f"RDKit MCS mapping failed: {e}") # Score the mappings and return them in sorted order (best to worst). mappings, scores = _score_rdkit_mappings( @@ -1066,6 +1098,9 @@ def _matchAtoms( "Using Sire MCS. Ignoring unsupported 'complete_rings_only' option!" ) + if mcs_kwargs: + _warnings.warn("Using Sire MCS. Ignoring unsupported 'mcs_kwargs' options!") + # Convert timeout to a Sire Unit. timeout = timeout * _SireUnits.second @@ -2079,6 +2114,7 @@ def merge( roi=None, property_map0={}, property_map1={}, + mcs_kwargs={}, **kwargs, ): """ @@ -2130,6 +2166,11 @@ def merge( A dictionary that maps "properties" in molecule1 to their user defined values. + mcs_kwargs : dict + A dictionary of keyword arguments used to override the defaults + passed to the RDKit MCS search. This is only used when 'mapping' + is None, i.e. when a mapping is autogenerated. + Returns ------- @@ -2208,6 +2249,7 @@ def merge( molecule1, property_map0=property_map0, property_map1=property_map1, + mcs_kwargs=mcs_kwargs, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index 2f7170fcb..6ff60101f 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -25,6 +25,7 @@ __email__ = "lester.hedges@gmail.com" __all__ = [ + "defaultMCSOptions", "generateNetwork", "matchAtoms", "viewMapping", @@ -702,6 +703,29 @@ def generateNetwork( return edges, scores +def defaultMCSOptions(): + """ + Return the default options used for the RDKit maximum common substructure + search. These can be overridden using the 'mcs_kwargs' argument of + :class:`matchAtoms `. + + Returns + ------- + + options : dict + The default RDKit MCS options. + """ + return { + "atomCompare": _rdFMCS.AtomCompare.CompareAny, + "bondCompare": _rdFMCS.BondCompare.CompareAny, + "completeRingsOnly": True, + "ringMatchesRingOnly": True, + "matchChiralTag": False, + "matchValences": False, + "maximizeBonds": False, + } + + def matchAtoms( molecule0, molecule1, @@ -717,6 +741,7 @@ def matchAtoms( max_scoring_matches=1000, property_map0={}, property_map1={}, + mcs_kwargs={}, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -770,6 +795,11 @@ def matchAtoms( option is only relevant to MCS performed using RDKit and will be ignored when falling back on Sire. + mcs_kwargs : dict + A dictionary of keyword arguments used to override the defaults + passed to the RDKit MCS search. This option is only relevant to MCS + performed using RDKit and will be ignored when falling back on Sire. + prune_perturbed_constraints : bool Whether to remove hydrogen atoms that are perturbed to heavy atoms from the mapping. This is True for AMBER by default and False for @@ -930,24 +960,21 @@ def matchAtoms( _Convert.toRDKit(molecule1, property_map=property_map1), ] + # Default MCS options, overridden by anything in 'mcs_kwargs'. The + # timeout is applied last so that it can't be overridden. + mcs_options = defaultMCSOptions() + mcs_options["completeRingsOnly"] = complete_rings_only + mcs_options.update(mcs_kwargs) + mcs_options["timeout"] = timeout + # Generate the MCS match. - mcs = _rdFMCS.FindMCS( - mols, - atomCompare=_rdFMCS.AtomCompare.CompareAny, - bondCompare=_rdFMCS.BondCompare.CompareAny, - completeRingsOnly=complete_rings_only, - ringMatchesRingOnly=True, - matchChiralTag=False, - matchValences=False, - maximizeBonds=False, - timeout=timeout, - ) + mcs = _rdFMCS.FindMCS(mols, **mcs_options) # Get the common substructure as a SMARTS string. mcs_smarts = _Chem.MolFromSmarts(mcs.smartsString) - except: - raise RuntimeError("RDKit MCS mapping failed!") + except Exception as e: + raise RuntimeError(f"RDKit MCS mapping failed: {e}") # Score the mappings and return them in sorted order (best to worst). mappings, scores = _score_rdkit_mappings( @@ -984,6 +1011,9 @@ def matchAtoms( "Using Sire MCS. Ignoring unsupported 'complete_rings_only' option!" ) + if mcs_kwargs: + _warnings.warn("Using Sire MCS. Ignoring unsupported 'mcs_kwargs' options!") + # Convert timeout to a Sire Unit. timeout = timeout * _SireUnits.second @@ -1370,6 +1400,7 @@ def merge( roi=None, property_map0={}, property_map1={}, + mcs_kwargs={}, **kwargs, ): """ @@ -1417,6 +1448,11 @@ def merge( A dictionary that maps "properties" in molecule1 to their user defined values. + mcs_kwargs : dict + A dictionary of keyword arguments used to override the defaults + passed to the RDKit MCS search. This is only used when 'mapping' + is None, i.e. when a mapping is autogenerated. + Returns ------- @@ -1488,6 +1524,7 @@ def merge( molecule1, property_map0=property_map0, property_map1=property_map1, + mcs_kwargs=mcs_kwargs, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) From ab5eb3b2f95b99b78d023ab12df0219b4b53fed8 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 29 Jul 2026 17:22:42 +0100 Subject: [PATCH 08/23] Add unit tests for mcs_kwargs and updated ring-break check. --- tests/Align/test_align.py | 97 ++++++++++++++++---- tests/Sandpit/Exscientia/Align/test_align.py | 64 +++++++++++++ 2 files changed, 145 insertions(+), 16 deletions(-) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index be01565e5..ef7fc2804 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -674,10 +674,10 @@ def test_roi_flex_align(protein_inputs): def test_empty_custom_roi_mapping(): # mut contains a proline mutation at position 15 wt = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_wt_flare_processed.pdb") )[0] mut = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_mut_flare_processed.pdb") )[0] # use the custom_roi_map to specify that residue 15 in the WT protein should be @@ -691,15 +691,16 @@ def test_empty_custom_roi_mapping(): for atom_idx in roi_res_idx: assert atom_idx not in mapping.keys() + @pytest.mark.skipif(has_amber is False, reason="Requires AMBER to be installed.") def test_custom_roi_ring_break_merge(): # wt contains a leucine at position 15 # mut contains a proline at position 15 wt = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_wt_flare_processed.pdb") )[0] mut = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_mut_flare_processed.pdb") )[0] wt = BSS.Parameters.ff14SB(wt, ensure_compatible=False).getMolecule() @@ -743,13 +744,14 @@ def test_custom_roi_ring_break_merge(): assert n_bonds_created == 1 assert n_bonds_annihilated == 0 + @pytest.mark.skipif(has_amber is False, reason="Requires AMBER to be installed.") def test_custom_roi_map_invalid_outside_roi(): wt = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_wt_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_wt_flare_processed.pdb") )[0] mut = BSS.IO.readMolecules( - BSS.IO.expand(BSS.tutorialUrl(), f"1choFH_apo_mut_flare_processed.pdb") + BSS.IO.expand(BSS.tutorialUrl(), "1choFH_apo_mut_flare_processed.pdb") )[0] wt = BSS.Parameters.ff14SB(wt, ensure_compatible=False).getMolecule() @@ -761,7 +763,6 @@ def test_custom_roi_map_invalid_outside_roi(): molecule0=wt, molecule1=mut, roi=[15], - custom_roi_map={ 0: 0, 1: 1, @@ -1313,9 +1314,9 @@ def test_ring_breaking_cross_bond_cleanup(): mol_info.atom_idx(p.atom3()).value(), } for a, b in changing: - assert not ( - a in atoms and b in atoms - ), f"improper{suffix} spans absent bond ({a},{b})" + assert not (a in atoms and b in atoms), ( + f"improper{suffix} spans absent bond ({a},{b})" + ) # Check that the ring-breaking and ring-making bond properties are set. def _read_pairs(prop_name): @@ -1326,9 +1327,73 @@ def _read_pairs(prop_name): stored_breaking = _read_pairs("ring_breaking_bonds") stored_making = _read_pairs("ring_making_bonds") - assert ( - stored_breaking == ring_breaking - ), f"ring_breaking_bonds property mismatch: {stored_breaking} != {ring_breaking}" - assert ( - stored_making == ring_making - ), f"ring_making_bonds property mismatch: {stored_making} != {ring_making}" + assert stored_breaking == ring_breaking, ( + f"ring_breaking_bonds property mismatch: {stored_breaking} != {ring_breaking}" + ) + assert stored_making == ring_making, ( + f"ring_making_bonds property mismatch: {stored_making} != {ring_making}" + ) + + +@pytest.fixture(scope="session") +def ejm31(): + return BSS.IO.readMolecules( + [f"{url}/lig_ejm31.prm7.bz2", f"{url}/lig_ejm31.rst7.bz2"] + ).getMolecules()[0] + + +@pytest.fixture(scope="session") +def jmc28(): + return BSS.IO.readMolecules( + [f"{url}/lig_jmc28.prm7.bz2", f"{url}/lig_jmc28.rst7.bz2"] + ).getMolecules()[0] + + +def test_default_mcs_options(): + # The MCS defaults should be discoverable, and ring matching is on. + options = BSS.Align.defaultMCSOptions() + assert options["ringMatchesRingOnly"] is True + assert options["completeRingsOnly"] is True + + # The returned dictionary is a copy, so mutating it has no side effects. + options["ringMatchesRingOnly"] = False + assert BSS.Align.defaultMCSOptions()["ringMatchesRingOnly"] is True + + +def test_mcs_kwargs_ring_matches_ring_only(ejm31, jmc28): + # Perturbing a methyl to a 2-methylcyclopropyl. Atom 19 is the methyl + # carbon in ejm31 and the ring carbon bonded to the carbonyl in jmc28. + + # By default an acyclic atom can't map onto a ring atom, so the whole + # substituent is unmapped. + mapping = BSS.Align.matchAtoms(ejm31, jmc28) + assert 19 not in mapping + + # Allowing the match maps the two carbons onto each other, along with one + # of the methyl hydrogens. + mapping = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False} + ) + assert mapping[19] == 19 + assert len(mapping) == 30 + + # Only two hydrogens are removed and the ring is grown from dummy atoms. + assert sorted(set(range(32)) - set(mapping)) == [27, 28] + + +def test_mcs_kwargs_merge(ejm31, jmc28): + # The options are used when merge autogenerates a mapping. + merged = BSS.Align.merge(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + sire_mol = merged._sire_object + + # A ring grown entirely from dummy atoms breaks no bond between mapped + # atoms, so the merge doesn't require 'allow_ring_breaking' and the end + # states have the same number of bonds. + assert sire_mol.num_atoms() == 41 + assert len(sire_mol.property("bond0").potentials()) == len( + sire_mol.property("bond1").potentials() + ) + + # No ring is broken or made, so neither property is set. + assert not sire_mol.has_property("ring_breaking_bonds") + assert not sire_mol.has_property("ring_making_bonds") diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 18f6c86f0..6c85a1ee0 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -904,3 +904,67 @@ def test_ring_opening_and_size_change(ligands, mapping): BSS.Align.merge( m0, m1, mapping, allow_ring_breaking=True, allow_ring_size_change=True ) + + +@pytest.fixture(scope="session") +def ejm31(): + return BSS.IO.readMolecules( + [f"{url}/lig_ejm31.prm7.bz2", f"{url}/lig_ejm31.rst7.bz2"] + ).getMolecules()[0] + + +@pytest.fixture(scope="session") +def jmc28(): + return BSS.IO.readMolecules( + [f"{url}/lig_jmc28.prm7.bz2", f"{url}/lig_jmc28.rst7.bz2"] + ).getMolecules()[0] + + +def test_default_mcs_options(): + # The MCS defaults should be discoverable, and ring matching is on. + options = BSS.Align.defaultMCSOptions() + assert options["ringMatchesRingOnly"] is True + assert options["completeRingsOnly"] is True + + # The returned dictionary is a copy, so mutating it has no side effects. + options["ringMatchesRingOnly"] = False + assert BSS.Align.defaultMCSOptions()["ringMatchesRingOnly"] is True + + +def test_mcs_kwargs_ring_matches_ring_only(ejm31, jmc28): + # Perturbing a methyl to a 2-methylcyclopropyl. Atom 19 is the methyl + # carbon in ejm31 and the ring carbon bonded to the carbonyl in jmc28. + + # By default an acyclic atom can't map onto a ring atom, so the whole + # substituent is unmapped. + mapping = BSS.Align.matchAtoms(ejm31, jmc28) + assert 19 not in mapping + + # Allowing the match maps the two carbons onto each other, along with one + # of the methyl hydrogens. + mapping = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False} + ) + assert mapping[19] == 19 + assert len(mapping) == 30 + + # Only two hydrogens are removed and the ring is grown from dummy atoms. + assert sorted(set(range(32)) - set(mapping)) == [27, 28] + + +def test_mcs_kwargs_merge(ejm31, jmc28): + # The options are used when merge autogenerates a mapping. + merged = BSS.Align.merge(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + sire_mol = merged._sire_object + + # A ring grown entirely from dummy atoms breaks no bond between mapped + # atoms, so the merge doesn't require 'allow_ring_breaking' and the end + # states have the same number of bonds. + assert sire_mol.num_atoms() == 41 + assert len(sire_mol.property("bond0").potentials()) == len( + sire_mol.property("bond1").potentials() + ) + + # No ring is broken or made, so neither property is set. + assert not sire_mol.has_property("ring_breaking_bonds") + assert not sire_mol.has_property("ring_making_bonds") From b049629e0780dd2a72fd2e6a2b95cce7901c679c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Thu, 30 Jul 2026 19:22:30 +0100 Subject: [PATCH 09/23] Warn when a mapping stops short at a pairable attachment point. --- src/BioSimSpace/Align/_align.py | 157 +++++++++++++++++- .../Sandpit/Exscientia/Align/_align.py | 150 ++++++++++++++++- tests/Align/test_align.py | 13 ++ tests/Sandpit/Exscientia/Align/test_align.py | 13 ++ 4 files changed, 331 insertions(+), 2 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 666cfe44d..2706fcb8e 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -896,7 +896,7 @@ def matchAtoms( """ if roi is None: - return _matchAtoms( + result = _matchAtoms( molecule0=molecule0, molecule1=molecule1, scoring_function=scoring_function, @@ -913,6 +913,58 @@ def matchAtoms( property_map1=property_map1, mcs_kwargs=mcs_kwargs, ) + + # Only check a mapping generated from the defaults. If the user has + # configured the MCS then both the baseline and our idea of a sensible + # mapping may not match their intent. This also stops the retry below + # from recursing, since it passes 'mcs_kwargs'. Skip when a prematch is + # given, since the retry could then fall back on Sire MCS, which + # ignores 'mcs_kwargs', making the comparison meaningless. + if not mcs_kwargs and not prematch: + best = result[0] if return_scores else result + if isinstance(best, list): + best = best[0] if best else {} + + # Attachment points where the MCS stopped on both sides. + flagged = ( + _flag_unmapped_attachments(molecule0, molecule1, best) if best else [] + ) + + if flagged: + # Retry with ring matching relaxed to see if it does better. + retry = matchAtoms( + molecule0, + molecule1, + scoring_function=scoring_function, + prematch=prematch, + timeout=timeout, + complete_rings_only=complete_rings_only, + max_scoring_matches=max_scoring_matches, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + prune_atom_types=prune_atom_types, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, + ) + + # Only trust the retry if it extends the mapping, i.e. keeps + # every existing pair and adds sensible ones. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension(molecule0, molecule1, best, retry) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides of " + f"atom(s) {flagged} " + f"in molecule0. Relaxing 'ringMatchesRingOnly' gives a " + f"common core of {len(retry)} rather than {len(best)}. " + f"Consider passing " + f"mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + + return result else: return _roiMatch( molecule0, @@ -927,6 +979,109 @@ def matchAtoms( ) +def _flag_unmapped_attachments(molecule0, molecule1, mapping): + """ + Internal function to find mapped atoms that have an unmapped heavy atom + neighbour of a common element in both molecules. These are attachment + points where the MCS stopped on both sides, which usually means that a + pairable atom was missed. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The atom mapping between the two molecules. + + Returns + ------- + + flagged : [int] + The indices of the flagged atoms in molecule0. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + conn0 = mol0.property("connectivity") + conn1 = mol1.property("connectivity") + + def _heavy_elements(mol, conn, idx, mapped): + elements = set() + for i in conn.connections_to(_SireMol.AtomIdx(idx)): + if i.value() not in mapped: + protons = mol.atom(i).property("element").num_protons() + if protons > 1: + elements.add(protons) + return elements + + mapped1 = set(mapping.values()) + flagged = [] + + for idx0, idx1 in mapping.items(): + elements0 = _heavy_elements(mol0, conn0, idx0, mapping) + if not elements0: + continue + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1) + if elements0 & elements1: + flagged.append(idx0) + + return flagged + + +def _is_sensible_extension(molecule0, molecule1, mapping, extended): + """ + Internal function to test whether the atoms that 'extended' adds relative + to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to + pair a heavy atom with a hydrogen, which grows the common core without + improving the mapping. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The original atom mapping. + + extended : dict + The larger atom mapping to test. + + Returns + ------- + + is_sensible : bool + Whether the added atoms pair heavy with heavy and hydrogen with + hydrogen. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + for idx0, idx1 in extended.items(): + if idx0 not in mapping: + protons0 = ( + mol0.atom(_SireMol.AtomIdx(idx0)).property("element").num_protons() + ) + protons1 = ( + mol1.atom(_SireMol.AtomIdx(idx1)).property("element").num_protons() + ) + if (protons0 > 1) != (protons1 > 1): + return False + + return True + + def _matchAtoms( molecule0, molecule1, diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index 6ff60101f..ce8a05ab9 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -703,6 +703,109 @@ def generateNetwork( return edges, scores +def _flag_unmapped_attachments(molecule0, molecule1, mapping): + """ + Internal function to find mapped atoms that have an unmapped heavy atom + neighbour of a common element in both molecules. These are attachment + points where the MCS stopped on both sides, which usually means that a + pairable atom was missed. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The atom mapping between the two molecules. + + Returns + ------- + + flagged : [int] + The indices of the flagged atoms in molecule0. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + conn0 = mol0.property("connectivity") + conn1 = mol1.property("connectivity") + + def _heavy_elements(mol, conn, idx, mapped): + elements = set() + for i in conn.connections_to(_SireMol.AtomIdx(idx)): + if i.value() not in mapped: + protons = mol.atom(i).property("element").num_protons() + if protons > 1: + elements.add(protons) + return elements + + mapped1 = set(mapping.values()) + flagged = [] + + for idx0, idx1 in mapping.items(): + elements0 = _heavy_elements(mol0, conn0, idx0, mapping) + if not elements0: + continue + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1) + if elements0 & elements1: + flagged.append(idx0) + + return flagged + + +def _is_sensible_extension(molecule0, molecule1, mapping, extended): + """ + Internal function to test whether the atoms that 'extended' adds relative + to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to + pair a heavy atom with a hydrogen, which grows the common core without + improving the mapping. + + Parameters + ---------- + + molecule0 : :class:`Molecule ` + The first molecule. + + molecule1 : :class:`Molecule ` + The second molecule. + + mapping : dict + The original atom mapping. + + extended : dict + The larger atom mapping to test. + + Returns + ------- + + is_sensible : bool + Whether the added atoms pair heavy with heavy and hydrogen with + hydrogen. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + for idx0, idx1 in extended.items(): + if idx0 not in mapping: + protons0 = ( + mol0.atom(_SireMol.AtomIdx(idx0)).property("element").num_protons() + ) + protons1 = ( + mol1.atom(_SireMol.AtomIdx(idx1)).property("element").num_protons() + ) + if (protons0 > 1) != (protons1 > 1): + return False + + return True + + def defaultMCSOptions(): """ Return the default options used for the RDKit maximum common substructure @@ -948,7 +1051,9 @@ def matchAtoms( mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - # Convert the timeout to seconds and take the value as an integer. + # Convert the timeout to seconds and take the value as an integer, keeping + # the original for any onward call that expects a Time. + orig_timeout = timeout timeout = int(timeout.seconds().value()) # Use RDKkit to find the maximum common substructure. @@ -1076,6 +1181,49 @@ def matchAtoms( _prune_crossing_constraints(molecule0, molecule1, x) for x in mappings ] + # Warn if the mapping stopped short at an attachment point where a pairable + # atom exists. Only check a mapping generated from the defaults. If the + # user has configured the MCS then both the baseline and our idea of a + # sensible mapping may not match their intent. This also stops the retry + # below from recursing, since it passes 'mcs_kwargs'. Skip when a prematch + # is given, since the retry could then fall back on Sire MCS, which ignores + # 'mcs_kwargs', making the comparison meaningless. + if not mcs_kwargs and not prematch and mappings: + flagged = _flag_unmapped_attachments(molecule0, molecule1, mappings[0]) + + if flagged: + # Retry with ring matching relaxed to see if it does better. + retry = matchAtoms( + molecule0, + molecule1, + engine=engine, + scoring_function=scoring_function, + prematch=prematch, + timeout=orig_timeout, + complete_rings_only=complete_rings_only, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + max_scoring_matches=max_scoring_matches, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, + ) + + # Only trust the retry if it extends the mapping, i.e. keeps every + # existing pair and adds sensible ones. + if ( + len(retry) > len(mappings[0]) + and set(mappings[0].items()) <= set(retry.items()) + and _is_sensible_extension(molecule0, molecule1, mappings[0], retry) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides of " + f"atom(s) {flagged} in molecule0. Relaxing " + f"'ringMatchesRingOnly' gives a common core of " + f"{len(retry)} rather than {len(mappings[0])}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + if matches == 1: if return_scores: return (mappings[0], scores[0]) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index ef7fc2804..aa2ba2af8 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1,4 +1,5 @@ import sys +import warnings import pytest import sire as sr @@ -1397,3 +1398,15 @@ def test_mcs_kwargs_merge(ejm31, jmc28): # No ring is broken or made, so neither property is set. assert not sire_mol.has_property("ring_breaking_bonds") assert not sire_mol.has_property("ring_making_bonds") + + +def test_unmapped_attachment_warning(ejm31, jmc28): + # The default mapping stops at the carbonyl carbon (atom 17), leaving + # heavy atoms unmapped on both sides, so we should be told about it. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.matchAtoms(ejm31, jmc28) + + # No warning once the option has been set explicitly. + with warnings.catch_warnings(): + warnings.simplefilter("error") + BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 6c85a1ee0..b6597eba0 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -1,4 +1,5 @@ import sys +import warnings import pytest from sire.legacy.Maths import Vector @@ -968,3 +969,15 @@ def test_mcs_kwargs_merge(ejm31, jmc28): # No ring is broken or made, so neither property is set. assert not sire_mol.has_property("ring_breaking_bonds") assert not sire_mol.has_property("ring_making_bonds") + + +def test_unmapped_attachment_warning(ejm31, jmc28): + # The default mapping stops at the carbonyl carbon (atom 17), leaving + # heavy atoms unmapped on both sides, so we should be told about it. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.matchAtoms(ejm31, jmc28) + + # No warning once the option has been set explicitly. + with warnings.catch_warnings(): + warnings.simplefilter("error") + BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) From 8279dd577c35fe325c12daae21494e91db06ef4e Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 12:34:01 +0100 Subject: [PATCH 10/23] Harden the mapping diagnostic against crashes and custom properties. --- src/BioSimSpace/Align/_align.py | 187 +++++++++++++----- .../Sandpit/Exscientia/Align/_align.py | 166 ++++++++++++---- tests/Align/test_align.py | 16 +- tests/Sandpit/Exscientia/Align/test_align.py | 16 +- 4 files changed, 291 insertions(+), 94 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 2706fcb8e..fc04962c9 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -921,49 +921,76 @@ def matchAtoms( # given, since the retry could then fall back on Sire MCS, which # ignores 'mcs_kwargs', making the comparison meaningless. if not mcs_kwargs and not prematch: - best = result[0] if return_scores else result - if isinstance(best, list): - best = best[0] if best else {} - - # Attachment points where the MCS stopped on both sides. - flagged = ( - _flag_unmapped_attachments(molecule0, molecule1, best) if best else [] - ) - - if flagged: - # Retry with ring matching relaxed to see if it does better. - retry = matchAtoms( - molecule0, - molecule1, - scoring_function=scoring_function, - prematch=prematch, - timeout=timeout, - complete_rings_only=complete_rings_only, - max_scoring_matches=max_scoring_matches, - prune_perturbed_constraints=prune_perturbed_constraints, - prune_crossing_constraints=prune_crossing_constraints, - prune_atom_types=prune_atom_types, - property_map0=property_map0, - property_map1=property_map1, - mcs_kwargs={"ringMatchesRingOnly": False}, + # This is a diagnostic, so it must never be able to break a call + # that would otherwise have succeeded. + try: + best = result[0] if return_scores else result + if isinstance(best, list): + best = best[0] if best else {} + + # Attachment points where the MCS stopped on both sides. + flagged = ( + _flag_unmapped_attachments( + molecule0, molecule1, best, property_map0, property_map1 + ) + if best + else [] ) - # Only trust the retry if it extends the mapping, i.e. keeps - # every existing pair and adds sensible ones. - if ( - len(retry) > len(best) - and set(best.items()) <= set(retry.items()) - and _is_sensible_extension(molecule0, molecule1, best, retry) - ): - _warnings.warn( - f"Mapping leaves heavy atoms unmapped on both sides of " - f"atom(s) {flagged} " - f"in molecule0. Relaxing 'ringMatchesRingOnly' gives a " - f"common core of {len(retry)} rather than {len(best)}. " - f"Consider passing " - f"mcs_kwargs={{'ringMatchesRingOnly': False}}." + if flagged: + # Retry with ring matching relaxed to see if it does + # better. Note that the RDKit documentation implies that + # 'completeRingsOnly' forces 'ringMatchesRingOnly', which + # would make this a no-op. It doesn't: as of RDKit + # 2026.03.4 the relaxed search still returns a larger MCS + # with 'completeRingsOnly' enabled. If a future RDKit + # changes this, the feature will silently stop firing. + # + # 'matches' and 'return_scores' are passed explicitly so + # that 'retry' is always a plain dict. The comparison below + # relies on it. 'prematch' is omitted since the enclosing + # guard means it's always empty. + retry = matchAtoms( + molecule0, + molecule1, + scoring_function=scoring_function, + matches=1, + return_scores=False, + timeout=timeout, + complete_rings_only=complete_rings_only, + max_scoring_matches=max_scoring_matches, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + prune_atom_types=prune_atom_types, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, ) + # Only trust the retry if it extends the mapping, i.e. + # keeps every existing pair and adds sensible ones. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension( + molecule0, + molecule1, + best, + retry, + property_map0, + property_map1, + ) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides " + f"of atom(s) {_format_flagged(flagged)} in molecule0. " + f"Relaxing 'ringMatchesRingOnly' gives a common core " + f"of {len(retry)} rather than {len(best)}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + except Exception as e: + _warnings.warn(f"Unable to check the quality of the mapping: {e}") + return result else: return _roiMatch( @@ -979,7 +1006,37 @@ def matchAtoms( ) -def _flag_unmapped_attachments(molecule0, molecule1, mapping): +def _format_flagged(flagged, max_show=5): + """ + Internal helper to format a list of atom indices for use in a warning + message, truncating so that the message stays readable for large + molecules. + + Parameters + ---------- + + flagged : [int] + The indices to format. + + max_show : int + The maximum number of indices to show. + + Returns + ------- + + string : str + The formatted indices. + """ + if len(flagged) <= max_show: + return str(flagged) + else: + shown = ", ".join(str(x) for x in flagged[:max_show]) + return f"[{shown}, ... ({len(flagged)} in total)]" + + +def _flag_unmapped_attachments( + molecule0, molecule1, mapping, property_map0={}, property_map1={} +): """ Internal function to find mapped atoms that have an unmapped heavy atom neighbour of a common element in both molecules. These are attachment @@ -998,6 +1055,16 @@ def _flag_unmapped_attachments(molecule0, molecule1, mapping): mapping : dict The atom mapping between the two molecules. + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + Returns ------- @@ -1008,33 +1075,42 @@ def _flag_unmapped_attachments(molecule0, molecule1, mapping): mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - conn0 = mol0.property("connectivity") - conn1 = mol1.property("connectivity") - def _heavy_elements(mol, conn, idx, mapped): + # Build the connectivity explicitly, since the molecules aren't guaranteed + # to have a stored "connectivity" property. + conn0 = _SireMol.Connectivity(mol0, _SireMol.CovalentBondHunter()) + conn1 = _SireMol.Connectivity(mol1, _SireMol.CovalentBondHunter()) + + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + + def _heavy_elements(mol, conn, idx, mapped, element): elements = set() for i in conn.connections_to(_SireMol.AtomIdx(idx)): if i.value() not in mapped: - protons = mol.atom(i).property("element").num_protons() + protons = mol.atom(i).property(element).num_protons() if protons > 1: elements.add(protons) return elements + mapped0 = set(mapping) mapped1 = set(mapping.values()) flagged = [] for idx0, idx1 in mapping.items(): - elements0 = _heavy_elements(mol0, conn0, idx0, mapping) + elements0 = _heavy_elements(mol0, conn0, idx0, mapped0, element0) if not elements0: continue - elements1 = _heavy_elements(mol1, conn1, idx1, mapped1) + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1, element1) if elements0 & elements1: flagged.append(idx0) return flagged -def _is_sensible_extension(molecule0, molecule1, mapping, extended): +def _is_sensible_extension( + molecule0, molecule1, mapping, extended, property_map0={}, property_map1={} +): """ Internal function to test whether the atoms that 'extended' adds relative to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to @@ -1056,6 +1132,16 @@ def _is_sensible_extension(molecule0, molecule1, mapping, extended): extended : dict The larger atom mapping to test. + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + Returns ------- @@ -1068,13 +1154,16 @@ def _is_sensible_extension(molecule0, molecule1, mapping, extended): mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + for idx0, idx1 in extended.items(): if idx0 not in mapping: protons0 = ( - mol0.atom(_SireMol.AtomIdx(idx0)).property("element").num_protons() + mol0.atom(_SireMol.AtomIdx(idx0)).property(element0).num_protons() ) protons1 = ( - mol1.atom(_SireMol.AtomIdx(idx1)).property("element").num_protons() + mol1.atom(_SireMol.AtomIdx(idx1)).property(element1).num_protons() ) if (protons0 > 1) != (protons1 > 1): return False diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index ce8a05ab9..b60937331 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -703,7 +703,37 @@ def generateNetwork( return edges, scores -def _flag_unmapped_attachments(molecule0, molecule1, mapping): +def _format_flagged(flagged, max_show=5): + """ + Internal helper to format a list of atom indices for use in a warning + message, truncating so that the message stays readable for large + molecules. + + Parameters + ---------- + + flagged : [int] + The indices to format. + + max_show : int + The maximum number of indices to show. + + Returns + ------- + + string : str + The formatted indices. + """ + if len(flagged) <= max_show: + return str(flagged) + else: + shown = ", ".join(str(x) for x in flagged[:max_show]) + return f"[{shown}, ... ({len(flagged)} in total)]" + + +def _flag_unmapped_attachments( + molecule0, molecule1, mapping, property_map0={}, property_map1={} +): """ Internal function to find mapped atoms that have an unmapped heavy atom neighbour of a common element in both molecules. These are attachment @@ -722,6 +752,16 @@ def _flag_unmapped_attachments(molecule0, molecule1, mapping): mapping : dict The atom mapping between the two molecules. + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + Returns ------- @@ -732,33 +772,42 @@ def _flag_unmapped_attachments(molecule0, molecule1, mapping): mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - conn0 = mol0.property("connectivity") - conn1 = mol1.property("connectivity") - def _heavy_elements(mol, conn, idx, mapped): + # Build the connectivity explicitly, since the molecules aren't guaranteed + # to have a stored "connectivity" property. + conn0 = _SireMol.Connectivity(mol0, _SireMol.CovalentBondHunter()) + conn1 = _SireMol.Connectivity(mol1, _SireMol.CovalentBondHunter()) + + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + + def _heavy_elements(mol, conn, idx, mapped, element): elements = set() for i in conn.connections_to(_SireMol.AtomIdx(idx)): if i.value() not in mapped: - protons = mol.atom(i).property("element").num_protons() + protons = mol.atom(i).property(element).num_protons() if protons > 1: elements.add(protons) return elements + mapped0 = set(mapping) mapped1 = set(mapping.values()) flagged = [] for idx0, idx1 in mapping.items(): - elements0 = _heavy_elements(mol0, conn0, idx0, mapping) + elements0 = _heavy_elements(mol0, conn0, idx0, mapped0, element0) if not elements0: continue - elements1 = _heavy_elements(mol1, conn1, idx1, mapped1) + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1, element1) if elements0 & elements1: flagged.append(idx0) return flagged -def _is_sensible_extension(molecule0, molecule1, mapping, extended): +def _is_sensible_extension( + molecule0, molecule1, mapping, extended, property_map0={}, property_map1={} +): """ Internal function to test whether the atoms that 'extended' adds relative to 'mapping' pair like with like. The MCS uses CompareAny, so it is free to @@ -780,6 +829,16 @@ def _is_sensible_extension(molecule0, molecule1, mapping, extended): extended : dict The larger atom mapping to test. + property_map0 : dict + A dictionary that maps "properties" in molecule0 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + + property_map1 : dict + A dictionary that maps "properties" in molecule1 to their user + defined values. This allows the user to refer to properties with + their own naming scheme, e.g. { "charge" : "my-charge" } + Returns ------- @@ -792,13 +851,16 @@ def _is_sensible_extension(molecule0, molecule1, mapping, extended): mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() + element0 = property_map0.get("element", "element") + element1 = property_map1.get("element", "element") + for idx0, idx1 in extended.items(): if idx0 not in mapping: protons0 = ( - mol0.atom(_SireMol.AtomIdx(idx0)).property("element").num_protons() + mol0.atom(_SireMol.AtomIdx(idx0)).property(element0).num_protons() ) protons1 = ( - mol1.atom(_SireMol.AtomIdx(idx1)).property("element").num_protons() + mol1.atom(_SireMol.AtomIdx(idx1)).property(element1).num_protons() ) if (protons0 > 1) != (protons1 > 1): return False @@ -1189,41 +1251,63 @@ def matchAtoms( # is given, since the retry could then fall back on Sire MCS, which ignores # 'mcs_kwargs', making the comparison meaningless. if not mcs_kwargs and not prematch and mappings: - flagged = _flag_unmapped_attachments(molecule0, molecule1, mappings[0]) - - if flagged: - # Retry with ring matching relaxed to see if it does better. - retry = matchAtoms( - molecule0, - molecule1, - engine=engine, - scoring_function=scoring_function, - prematch=prematch, - timeout=orig_timeout, - complete_rings_only=complete_rings_only, - prune_perturbed_constraints=prune_perturbed_constraints, - prune_crossing_constraints=prune_crossing_constraints, - max_scoring_matches=max_scoring_matches, - property_map0=property_map0, - property_map1=property_map1, - mcs_kwargs={"ringMatchesRingOnly": False}, + # This is a diagnostic, so it must never be able to break a call that + # would otherwise have succeeded. + try: + best = mappings[0] + flagged = _flag_unmapped_attachments( + molecule0, molecule1, best, property_map0, property_map1 ) - # Only trust the retry if it extends the mapping, i.e. keeps every - # existing pair and adds sensible ones. - if ( - len(retry) > len(mappings[0]) - and set(mappings[0].items()) <= set(retry.items()) - and _is_sensible_extension(molecule0, molecule1, mappings[0], retry) - ): - _warnings.warn( - f"Mapping leaves heavy atoms unmapped on both sides of " - f"atom(s) {flagged} in molecule0. Relaxing " - f"'ringMatchesRingOnly' gives a common core of " - f"{len(retry)} rather than {len(mappings[0])}. Consider " - f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + if flagged: + # Retry with ring matching relaxed to see if it does better. + # Note that the RDKit documentation implies that + # 'completeRingsOnly' forces 'ringMatchesRingOnly', which would + # make this a no-op. It doesn't: as of RDKit 2026.03.4 the + # relaxed search still returns a larger MCS with + # 'completeRingsOnly' enabled. If a future RDKit changes this, + # the feature will silently stop firing. + # + # 'matches' and 'return_scores' are passed explicitly so that + # 'retry' is always a plain dict. The comparison below relies + # on it. 'prematch' is omitted since the enclosing guard means + # it's always empty. + retry = matchAtoms( + molecule0, + molecule1, + engine=engine, + scoring_function=scoring_function, + matches=1, + return_scores=False, + timeout=orig_timeout, + complete_rings_only=complete_rings_only, + prune_perturbed_constraints=prune_perturbed_constraints, + prune_crossing_constraints=prune_crossing_constraints, + max_scoring_matches=max_scoring_matches, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, ) + # Only trust the retry if it extends the mapping, i.e. keeps + # every existing pair and adds sensible ones. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension( + molecule0, molecule1, best, retry, property_map0, property_map1 + ) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides " + f"of atom(s) {_format_flagged(flagged)} in molecule0. " + f"Relaxing 'ringMatchesRingOnly' gives a common core " + f"of {len(retry)} rather than {len(best)}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + except Exception as e: + _warnings.warn(f"Unable to check the quality of the mapping: {e}") + if matches == 1: if return_scores: return (mappings[0], scores[0]) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index aa2ba2af8..bcc620917 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1406,7 +1406,19 @@ def test_unmapped_attachment_warning(ejm31, jmc28): with pytest.warns(UserWarning, match="ringMatchesRingOnly"): BSS.Align.matchAtoms(ejm31, jmc28) - # No warning once the option has been set explicitly. + # No warning once the option has been set explicitly. Only promote the + # warning we care about, so that unrelated warnings from RDKit or Sire + # don't fail the test. with warnings.catch_warnings(): - warnings.simplefilter("error") + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + + # The returned mapping must be unchanged by the check, since it is only + # meant to be an observation. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + checked = BSS.Align.matchAtoms(ejm31, jmc28) + unchecked = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + assert checked == unchecked diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index b6597eba0..9b4349315 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -977,7 +977,19 @@ def test_unmapped_attachment_warning(ejm31, jmc28): with pytest.warns(UserWarning, match="ringMatchesRingOnly"): BSS.Align.matchAtoms(ejm31, jmc28) - # No warning once the option has been set explicitly. + # No warning once the option has been set explicitly. Only promote the + # warning we care about, so that unrelated warnings from RDKit or Sire + # don't fail the test. with warnings.catch_warnings(): - warnings.simplefilter("error") + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") BSS.Align.matchAtoms(ejm31, jmc28, mcs_kwargs={"ringMatchesRingOnly": False}) + + # The returned mapping must be unchanged by the check, since it is only + # meant to be an observation. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + checked = BSS.Align.matchAtoms(ejm31, jmc28) + unchecked = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + assert checked == unchecked From 30895c10cacff6fde734751b29bb6ee52531fd1c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 12:53:14 +0100 Subject: [PATCH 11/23] Check the mapping before pruning, not after. --- src/BioSimSpace/Align/_align.py | 157 +++++++++--------- .../Sandpit/Exscientia/Align/_align.py | 50 +++--- 2 files changed, 104 insertions(+), 103 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index fc04962c9..40e20c867 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -896,7 +896,7 @@ def matchAtoms( """ if roi is None: - result = _matchAtoms( + return _matchAtoms( molecule0=molecule0, molecule1=molecule1, scoring_function=scoring_function, @@ -913,85 +913,6 @@ def matchAtoms( property_map1=property_map1, mcs_kwargs=mcs_kwargs, ) - - # Only check a mapping generated from the defaults. If the user has - # configured the MCS then both the baseline and our idea of a sensible - # mapping may not match their intent. This also stops the retry below - # from recursing, since it passes 'mcs_kwargs'. Skip when a prematch is - # given, since the retry could then fall back on Sire MCS, which - # ignores 'mcs_kwargs', making the comparison meaningless. - if not mcs_kwargs and not prematch: - # This is a diagnostic, so it must never be able to break a call - # that would otherwise have succeeded. - try: - best = result[0] if return_scores else result - if isinstance(best, list): - best = best[0] if best else {} - - # Attachment points where the MCS stopped on both sides. - flagged = ( - _flag_unmapped_attachments( - molecule0, molecule1, best, property_map0, property_map1 - ) - if best - else [] - ) - - if flagged: - # Retry with ring matching relaxed to see if it does - # better. Note that the RDKit documentation implies that - # 'completeRingsOnly' forces 'ringMatchesRingOnly', which - # would make this a no-op. It doesn't: as of RDKit - # 2026.03.4 the relaxed search still returns a larger MCS - # with 'completeRingsOnly' enabled. If a future RDKit - # changes this, the feature will silently stop firing. - # - # 'matches' and 'return_scores' are passed explicitly so - # that 'retry' is always a plain dict. The comparison below - # relies on it. 'prematch' is omitted since the enclosing - # guard means it's always empty. - retry = matchAtoms( - molecule0, - molecule1, - scoring_function=scoring_function, - matches=1, - return_scores=False, - timeout=timeout, - complete_rings_only=complete_rings_only, - max_scoring_matches=max_scoring_matches, - prune_perturbed_constraints=prune_perturbed_constraints, - prune_crossing_constraints=prune_crossing_constraints, - prune_atom_types=prune_atom_types, - property_map0=property_map0, - property_map1=property_map1, - mcs_kwargs={"ringMatchesRingOnly": False}, - ) - - # Only trust the retry if it extends the mapping, i.e. - # keeps every existing pair and adds sensible ones. - if ( - len(retry) > len(best) - and set(best.items()) <= set(retry.items()) - and _is_sensible_extension( - molecule0, - molecule1, - best, - retry, - property_map0, - property_map1, - ) - ): - _warnings.warn( - f"Mapping leaves heavy atoms unmapped on both sides " - f"of atom(s) {_format_flagged(flagged)} in molecule0. " - f"Relaxing 'ringMatchesRingOnly' gives a common core " - f"of {len(retry)} rather than {len(best)}. Consider " - f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." - ) - except Exception as e: - _warnings.warn(f"Unable to check the quality of the mapping: {e}") - - return result else: return _roiMatch( molecule0, @@ -1279,7 +1200,10 @@ def _matchAtoms( mol0 = molecule0._getSireObject() mol1 = molecule1._getSireObject() - # Convert the timeout to seconds and take the value as an integer. + # Convert the timeout to seconds and take the value as an integer. Keep + # the original, since the mapping check below re-enters this function, + # which expects a Time object. + orig_timeout = timeout timeout = int(timeout.seconds().value()) # Use RDKkit to find the maximum common substructure. @@ -1397,6 +1321,77 @@ def _matchAtoms( property_map1, ) + # Warn if the mapping stopped short at an attachment point where a pairable + # atom exists. This is done before the pruning below, since pruning deletes + # correctly mapped heavy atom pairs, which manufactures exactly the + # signature that the check looks for. Only check a mapping generated from + # the defaults. If the user has configured the MCS then both the baseline + # and our idea of a sensible mapping may not match their intent. This also + # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a + # prematch is given, since the retry could then fall back on Sire MCS, + # which ignores 'mcs_kwargs', making the comparison meaningless. + if not mcs_kwargs and not prematch and mappings: + # This is a diagnostic, so it must never be able to break a call that + # would otherwise have succeeded. + try: + best = mappings[0] + + # Attachment points where the MCS stopped on both sides. + flagged = _flag_unmapped_attachments( + molecule0, molecule1, best, property_map0, property_map1 + ) + + if flagged: + # Retry with ring matching relaxed to see if it does better. + # Note that the RDKit documentation implies that + # 'completeRingsOnly' forces 'ringMatchesRingOnly', which would + # make this a no-op. It doesn't: as of RDKit 2026.03.4 the + # relaxed search still returns a larger MCS with + # 'completeRingsOnly' enabled. If a future RDKit changes this, + # the feature will silently stop firing. + # + # Pruning is disabled so that the retry is compared like for + # like against the unpruned mapping above. 'matches' and + # 'return_scores' are passed explicitly so that 'retry' is + # always a plain dict, which the comparison below relies on. + # 'prematch' is omitted since the enclosing guard means it's + # always empty. + retry = _matchAtoms( + molecule0=molecule0, + molecule1=molecule1, + scoring_function=scoring_function, + matches=1, + return_scores=False, + timeout=orig_timeout, + complete_rings_only=complete_rings_only, + max_scoring_matches=max_scoring_matches, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, + prune_atom_types=False, + property_map0=property_map0, + property_map1=property_map1, + mcs_kwargs={"ringMatchesRingOnly": False}, + ) + + # Only trust the retry if it extends the mapping, i.e. keeps + # every existing pair and adds sensible ones. + if ( + len(retry) > len(best) + and set(best.items()) <= set(retry.items()) + and _is_sensible_extension( + molecule0, molecule1, best, retry, property_map0, property_map1 + ) + ): + _warnings.warn( + f"Mapping leaves heavy atoms unmapped on both sides " + f"of atom(s) {_format_flagged(flagged)} in molecule0. " + f"Relaxing 'ringMatchesRingOnly' gives a common core " + f"of {len(retry)} rather than {len(best)}. Consider " + f"passing mcs_kwargs={{'ringMatchesRingOnly': False}}." + ) + except Exception as e: + _warnings.warn(f"Unable to check the quality of the mapping: {e}") + # Optionally post-process the MCS for use with AMBER. if prune_perturbed_constraints: mappings = [ diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index b60937331..ec9450823 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -1233,28 +1233,22 @@ def matchAtoms( property_map1, ) - # Optionally post-process the MCS. - if prune_perturbed_constraints: - mappings = [ - _prune_perturbed_constraints(molecule0, molecule1, x) for x in mappings - ] - if prune_crossing_constraints: - mappings = [ - _prune_crossing_constraints(molecule0, molecule1, x) for x in mappings - ] - # Warn if the mapping stopped short at an attachment point where a pairable - # atom exists. Only check a mapping generated from the defaults. If the - # user has configured the MCS then both the baseline and our idea of a - # sensible mapping may not match their intent. This also stops the retry - # below from recursing, since it passes 'mcs_kwargs'. Skip when a prematch - # is given, since the retry could then fall back on Sire MCS, which ignores - # 'mcs_kwargs', making the comparison meaningless. + # atom exists. This is done before the pruning below, since pruning deletes + # correctly mapped heavy atom pairs, which manufactures exactly the + # signature that the check looks for. Only check a mapping generated from + # the defaults. If the user has configured the MCS then both the baseline + # and our idea of a sensible mapping may not match their intent. This also + # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a + # prematch is given, since the retry could then fall back on Sire MCS, + # which ignores 'mcs_kwargs', making the comparison meaningless. if not mcs_kwargs and not prematch and mappings: # This is a diagnostic, so it must never be able to break a call that # would otherwise have succeeded. try: best = mappings[0] + + # Attachment points where the MCS stopped on both sides. flagged = _flag_unmapped_attachments( molecule0, molecule1, best, property_map0, property_map1 ) @@ -1268,10 +1262,12 @@ def matchAtoms( # 'completeRingsOnly' enabled. If a future RDKit changes this, # the feature will silently stop firing. # - # 'matches' and 'return_scores' are passed explicitly so that - # 'retry' is always a plain dict. The comparison below relies - # on it. 'prematch' is omitted since the enclosing guard means - # it's always empty. + # Pruning is disabled so that the retry is compared like for + # like against the unpruned mapping above. 'matches' and + # 'return_scores' are passed explicitly so that 'retry' is + # always a plain dict, which the comparison below relies on. + # 'prematch' is omitted since the enclosing guard means it's + # always empty. retry = matchAtoms( molecule0, molecule1, @@ -1281,8 +1277,8 @@ def matchAtoms( return_scores=False, timeout=orig_timeout, complete_rings_only=complete_rings_only, - prune_perturbed_constraints=prune_perturbed_constraints, - prune_crossing_constraints=prune_crossing_constraints, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, max_scoring_matches=max_scoring_matches, property_map0=property_map0, property_map1=property_map1, @@ -1308,6 +1304,16 @@ def matchAtoms( except Exception as e: _warnings.warn(f"Unable to check the quality of the mapping: {e}") + # Optionally post-process the MCS. + if prune_perturbed_constraints: + mappings = [ + _prune_perturbed_constraints(molecule0, molecule1, x) for x in mappings + ] + if prune_crossing_constraints: + mappings = [ + _prune_crossing_constraints(molecule0, molecule1, x) for x in mappings + ] + if matches == 1: if return_scores: return (mappings[0], scores[0]) From 0dc6b77a5eb544d6585e66779280092b34df2833 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 13:09:53 +0100 Subject: [PATCH 12/23] Add negative, unit and branch coverage for the mapping check. --- tests/Align/test_align.py | 171 ++++++++++++++++++ tests/Sandpit/Exscientia/Align/test_align.py | 178 ++++++++++++++++++- 2 files changed, 348 insertions(+), 1 deletion(-) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index bcc620917..6c51a44b6 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1422,3 +1422,174 @@ def test_unmapped_attachment_warning(ejm31, jmc28): ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() ) assert checked == unchecked + + +def test_unmapped_attachment_no_warning(ejm31): + # A molecule mapped to itself leaves nothing unmapped, so there is + # nothing to flag. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, ejm31) + + +@pytest.mark.skipif( + not has_antechamber or not has_tleap, + reason="Requires antechamber and tLEaP to be installed.", +) +def test_unmapped_attachment_no_warning_r_group(monkeypatch): + """ + Regression test for the check running on the pruned mapping. Pruning + deletes correctly mapped heavy atom pairs, which looks identical to an + MCS that stopped short, so these ordinary R-group edits used to be + flagged and pay for a second MCS search for nothing. + + No warning was ever emitted for them, since the gate rejected the retry, + so the flagged atoms are spied on directly rather than the warning. + """ + from BioSimSpace.Align import _align + + pairs = [ + ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ] + + flagged = [] + original = _align._flag_unmapped_attachments + + def _spy(*args, **kwargs): + result = original(*args, **kwargs) + flagged.append(result) + return result + + monkeypatch.setattr(_align, "_flag_unmapped_attachments", _spy) + + for smiles0, smiles1 in pairs: + molecule0 = BSS.Parameters.gaff2(smiles0).getMolecule() + molecule1 = BSS.Parameters.gaff2(smiles1).getMolecule() + + del flagged[:] + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + warnings.filterwarnings("error", message=".*Unable to check.*") + BSS.Align.matchAtoms( + molecule0, + molecule1, + prune_perturbed_constraints=True, + prune_crossing_constraints=True, + ) + + # The check should run once, on the unpruned mapping, and find + # nothing. A second entry would mean the retry had run too. + assert flagged == [[]] + + +def test_unmapped_attachment_warning_other_branches(ejm31, jmc28): + # The check must cope with the alternative return shapes, and must not + # change what is returned in any of them. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + reference = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + + mappings = BSS.Align.matchAtoms(ejm31, jmc28, matches=5) + assert isinstance(mappings, list) + assert mappings[0] == reference + + mapping, score = BSS.Align.matchAtoms(ejm31, jmc28, return_scores=True) + assert mapping == reference + + mappings, scores = BSS.Align.matchAtoms( + ejm31, jmc28, matches=5, return_scores=True + ) + assert len(mappings) == len(scores) + assert mappings[0] == reference + + # A prematch skips the check, since the retry could fall back on the Sire + # MCS, which ignores 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28, prematch={0: 0}) + + +def test_flag_unmapped_attachments(ejm31): + """ + Unit test for the attachment point check, using hand-built mappings so + that no MCS search is involved. + """ + from sire.legacy import Mol as _SireMol + + from BioSimSpace.Align._align import _flag_unmapped_attachments + + sire_mol = ejm31._sire_object + connectivity = _SireMol.Connectivity(sire_mol, _SireMol.CovalentBondHunter()) + + # Map the molecule onto itself. Nothing is unmapped, so nothing is flagged. + identity = {x: x for x in range(sire_mol.num_atoms())} + assert _flag_unmapped_attachments(ejm31, ejm31, identity) == [] + + # Find a heavy atom with a heavy atom neighbour. + for atom in sire_mol.atoms(): + idx = atom.index().value() + if atom.property("element").num_protons() == 1: + continue + neighbours = [ + i.value() + for i in connectivity.connections_to(_SireMol.AtomIdx(idx)) + if sire_mol.atom(i).property("element").num_protons() > 1 + ] + if neighbours: + break + + # Drop the neighbour from the mapping. It is now an unmapped heavy atom + # of the same element on both sides of 'idx', so 'idx' is flagged. + truncated = dict(identity) + del truncated[neighbours[0]] + assert idx in _flag_unmapped_attachments(ejm31, ejm31, truncated) + + # Hydrogens are ignored, so dropping one flags nothing. + hydrogen = next( + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ) + truncated = dict(identity) + del truncated[hydrogen] + assert _flag_unmapped_attachments(ejm31, ejm31, truncated) == [] + + +def test_is_sensible_extension(ejm31): + """ + Unit test for the heavy/hydrogen check on the atoms that the retry adds. + """ + from BioSimSpace.Align._align import _is_sensible_extension + + sire_mol = ejm31._sire_object + + heavy = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() > 1 + ] + hydrogens = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ] + + mapping = {heavy[0]: heavy[0]} + + # Heavy to heavy and hydrogen to hydrogen are both sensible. + extended = dict(mapping) + extended[heavy[1]] = heavy[1] + extended[hydrogens[0]] = hydrogens[0] + assert _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Pairing a heavy atom with a hydrogen is not. + extended = dict(mapping) + extended[heavy[1]] = hydrogens[0] + assert not _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Adding nothing is trivially sensible. + assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 9b4349315..ea3c78ae5 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -7,7 +7,12 @@ from sire.legacy.Mol import AtomIdx, Element, PartialMolecule import BioSimSpace.Sandpit.Exscientia as BSS -from tests.Sandpit.Exscientia.conftest import has_antechamber, has_openff, url +from tests.Sandpit.Exscientia.conftest import ( + has_antechamber, + has_openff, + has_tleap, + url, +) @pytest.fixture(scope="session") @@ -993,3 +998,174 @@ def test_unmapped_attachment_warning(ejm31, jmc28): ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() ) assert checked == unchecked + + +def test_unmapped_attachment_no_warning(ejm31): + # A molecule mapped to itself leaves nothing unmapped, so there is + # nothing to flag. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, ejm31) + + +@pytest.mark.skipif( + not has_antechamber or not has_tleap, + reason="Requires antechamber and tLEaP to be installed.", +) +def test_unmapped_attachment_no_warning_r_group(monkeypatch): + """ + Regression test for the check running on the pruned mapping. Pruning + deletes correctly mapped heavy atom pairs, which looks identical to an + MCS that stopped short, so these ordinary R-group edits used to be + flagged and pay for a second MCS search for nothing. + + No warning was ever emitted for them, since the gate rejected the retry, + so the flagged atoms are spied on directly rather than the warning. + """ + from BioSimSpace.Sandpit.Exscientia.Align import _align + + pairs = [ + ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ] + + flagged = [] + original = _align._flag_unmapped_attachments + + def _spy(*args, **kwargs): + result = original(*args, **kwargs) + flagged.append(result) + return result + + monkeypatch.setattr(_align, "_flag_unmapped_attachments", _spy) + + for smiles0, smiles1 in pairs: + molecule0 = BSS.Parameters.gaff2(smiles0).getMolecule() + molecule1 = BSS.Parameters.gaff2(smiles1).getMolecule() + + del flagged[:] + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + warnings.filterwarnings("error", message=".*Unable to check.*") + BSS.Align.matchAtoms( + molecule0, + molecule1, + prune_perturbed_constraints=True, + prune_crossing_constraints=True, + ) + + # The check should run once, on the unpruned mapping, and find + # nothing. A second entry would mean the retry had run too. + assert flagged == [[]] + + +def test_unmapped_attachment_warning_other_branches(ejm31, jmc28): + # The check must cope with the alternative return shapes, and must not + # change what is returned in any of them. + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + reference = BSS.Align.matchAtoms( + ejm31, jmc28, mcs_kwargs=BSS.Align.defaultMCSOptions() + ) + + mappings = BSS.Align.matchAtoms(ejm31, jmc28, matches=5) + assert isinstance(mappings, list) + assert mappings[0] == reference + + mapping, score = BSS.Align.matchAtoms(ejm31, jmc28, return_scores=True) + assert mapping == reference + + mappings, scores = BSS.Align.matchAtoms( + ejm31, jmc28, matches=5, return_scores=True + ) + assert len(mappings) == len(scores) + assert mappings[0] == reference + + # A prematch skips the check, since the retry could fall back on the Sire + # MCS, which ignores 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28, prematch={0: 0}) + + +def test_flag_unmapped_attachments(ejm31): + """ + Unit test for the attachment point check, using hand-built mappings so + that no MCS search is involved. + """ + from sire.legacy import Mol as _SireMol + + from BioSimSpace.Sandpit.Exscientia.Align._align import _flag_unmapped_attachments + + sire_mol = ejm31._sire_object + connectivity = _SireMol.Connectivity(sire_mol, _SireMol.CovalentBondHunter()) + + # Map the molecule onto itself. Nothing is unmapped, so nothing is flagged. + identity = {x: x for x in range(sire_mol.num_atoms())} + assert _flag_unmapped_attachments(ejm31, ejm31, identity) == [] + + # Find a heavy atom with a heavy atom neighbour. + for atom in sire_mol.atoms(): + idx = atom.index().value() + if atom.property("element").num_protons() == 1: + continue + neighbours = [ + i.value() + for i in connectivity.connections_to(_SireMol.AtomIdx(idx)) + if sire_mol.atom(i).property("element").num_protons() > 1 + ] + if neighbours: + break + + # Drop the neighbour from the mapping. It is now an unmapped heavy atom + # of the same element on both sides of 'idx', so 'idx' is flagged. + truncated = dict(identity) + del truncated[neighbours[0]] + assert idx in _flag_unmapped_attachments(ejm31, ejm31, truncated) + + # Hydrogens are ignored, so dropping one flags nothing. + hydrogen = next( + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ) + truncated = dict(identity) + del truncated[hydrogen] + assert _flag_unmapped_attachments(ejm31, ejm31, truncated) == [] + + +def test_is_sensible_extension(ejm31): + """ + Unit test for the heavy/hydrogen check on the atoms that the retry adds. + """ + from BioSimSpace.Sandpit.Exscientia.Align._align import _is_sensible_extension + + sire_mol = ejm31._sire_object + + heavy = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() > 1 + ] + hydrogens = [ + a.index().value() + for a in sire_mol.atoms() + if a.property("element").num_protons() == 1 + ] + + mapping = {heavy[0]: heavy[0]} + + # Heavy to heavy and hydrogen to hydrogen are both sensible. + extended = dict(mapping) + extended[heavy[1]] = heavy[1] + extended[hydrogens[0]] = hydrogens[0] + assert _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Pairing a heavy atom with a hydrogen is not. + extended = dict(mapping) + extended[heavy[1]] = hydrogens[0] + assert not _is_sensible_extension(ejm31, ejm31, mapping, extended) + + # Adding nothing is trivially sensible. + assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) From 63c5704f60117b766af9d3b3d3545453443a7657 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 13:26:03 +0100 Subject: [PATCH 13/23] Only check the mapping where the advice can be acted on. --- src/BioSimSpace/Align/_align.py | 27 +++++++++-- .../Sandpit/Exscientia/Align/_align.py | 22 +++++++-- tests/Align/test_align.py | 46 +++++++++++++++++++ tests/Sandpit/Exscientia/Align/test_align.py | 28 +++++++++++ 4 files changed, 117 insertions(+), 6 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 40e20c867..9f1c0580b 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -744,6 +744,7 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + _check_mapping=True, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -912,6 +913,7 @@ def matchAtoms( property_map0=property_map0, property_map1=property_map1, mcs_kwargs=mcs_kwargs, + _check_mapping=_check_mapping, ) else: return _roiMatch( @@ -1109,6 +1111,7 @@ def _matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + _check_mapping=True, ): import sys as _sys @@ -1330,9 +1333,12 @@ def _matchAtoms( # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a # prematch is given, since the retry could then fall back on Sire MCS, # which ignores 'mcs_kwargs', making the comparison meaningless. - if not mcs_kwargs and not prematch and mappings: + if _check_mapping and not mcs_kwargs and not prematch and mappings: # This is a diagnostic, so it must never be able to break a call that - # would otherwise have succeeded. + # would otherwise have succeeded. The warning itself is emitted outside + # the guard, since it would otherwise be swallowed and re-reported as a + # failure whenever the user has promoted warnings to errors. + message = None try: best = mappings[0] @@ -1382,7 +1388,7 @@ def _matchAtoms( molecule0, molecule1, best, retry, property_map0, property_map1 ) ): - _warnings.warn( + message = ( f"Mapping leaves heavy atoms unmapped on both sides " f"of atom(s) {_format_flagged(flagged)} in molecule0. " f"Relaxing 'ringMatchesRingOnly' gives a common core " @@ -1392,6 +1398,9 @@ def _matchAtoms( except Exception as e: _warnings.warn(f"Unable to check the quality of the mapping: {e}") + if message is not None: + _warnings.warn(message) + # Optionally post-process the MCS for use with AMBER. if prune_perturbed_constraints: mappings = [ @@ -1732,6 +1741,9 @@ def _roiMatch( mapping = matchAtoms( res0_extracted, res1_extracted, + # The mapping check would report indices that are local to the + # extracted residue, not to molecule0 as its message claims. + _check_mapping=False, ) # Look up the absolute atom indices in the molecule if not using a custom ROI mapping. @@ -1955,6 +1967,9 @@ def _rmsdAlign(molecule0, molecule1, mapping=None, property_map0={}, property_ma molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Extract the Sire molecule from each BioSimSpace molecule. @@ -2157,6 +2172,9 @@ def _flexAlign( molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Convert the mapping to AtomIdx key:value pairs. @@ -2614,6 +2632,9 @@ def viewMapping( molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index ec9450823..554e002c2 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -907,6 +907,7 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + _check_mapping=True, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -1242,9 +1243,12 @@ def matchAtoms( # stops the retry from recursing, since it passes 'mcs_kwargs'. Skip when a # prematch is given, since the retry could then fall back on Sire MCS, # which ignores 'mcs_kwargs', making the comparison meaningless. - if not mcs_kwargs and not prematch and mappings: + if _check_mapping and not mcs_kwargs and not prematch and mappings: # This is a diagnostic, so it must never be able to break a call that - # would otherwise have succeeded. + # would otherwise have succeeded. The warning itself is emitted outside + # the guard, since it would otherwise be swallowed and re-reported as a + # failure whenever the user has promoted warnings to errors. + message = None try: best = mappings[0] @@ -1294,7 +1298,7 @@ def matchAtoms( molecule0, molecule1, best, retry, property_map0, property_map1 ) ): - _warnings.warn( + message = ( f"Mapping leaves heavy atoms unmapped on both sides " f"of atom(s) {_format_flagged(flagged)} in molecule0. " f"Relaxing 'ringMatchesRingOnly' gives a common core " @@ -1304,6 +1308,9 @@ def matchAtoms( except Exception as e: _warnings.warn(f"Unable to check the quality of the mapping: {e}") + if message is not None: + _warnings.warn(message) + # Optionally post-process the MCS. if prune_perturbed_constraints: mappings = [ @@ -1415,6 +1422,9 @@ def rmsdAlign(molecule0, molecule1, mapping=None, property_map0={}, property_map molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Extract the Sire molecule from each BioSimSpace molecule. @@ -1572,6 +1582,9 @@ def flexAlign( molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) # Convert the mapping to AtomIdx key:value pairs. @@ -1905,6 +1918,9 @@ def viewMapping( molecule1, property_map0=property_map0, property_map1=property_map1, + # This function doesn't take 'mcs_kwargs', so the advice from the + # mapping check can't be acted on here. + _check_mapping=False, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index 6c51a44b6..fa4295048 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1593,3 +1593,49 @@ def test_is_sensible_extension(ejm31): # Adding nothing is trivially sensible. assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) + + +def test_unmapped_attachment_check_suppressed(ejm31, jmc28): + """ + The check should only fire where the user can act on its advice, i.e. + where 'mcs_kwargs' can be passed through. It should also never fire on + the ROI path, where the flagged indices would be local to the extracted + residue rather than to molecule0 as the message claims. + """ + # These functions don't take 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.rmsdAlign(ejm31, jmc28) + BSS.Align.flexAlign(ejm31, jmc28) + + # 'merge' does, so the check stays on. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.merge(ejm31, jmc28, force=True) + + +def test_unmapped_attachment_check_suppressed_roi(protein_inputs): + # The ROI path maps each residue of interest separately, so any flagged + # indices would be local to that residue rather than to molecule0. + proteins, protein_mapping, roi = protein_inputs + p0 = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"{proteins}_mut_peptide.pdb") + )[0] + p1 = BSS.IO.readMolecules( + BSS.IO.expand(BSS.tutorialUrl(), f"{proteins}_wt_peptide.pdb") + )[0] + + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + assert BSS.Align.matchAtoms(p0, p1, roi=roi) == protein_mapping + + +def test_unmapped_attachment_warning_not_swallowed(ejm31, jmc28): + # Promoting the warning to an error must surface the warning itself, not + # a report that the check failed. The warning is emitted outside the + # try/except that guards the check for exactly this reason. + # Anchored, since the wrapped "Unable to check ..." message quotes the + # original and would otherwise match too. + with pytest.raises(UserWarning, match=r"^Mapping leaves heavy atoms"): + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28) diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index ea3c78ae5..0e413b826 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -1169,3 +1169,31 @@ def test_is_sensible_extension(ejm31): # Adding nothing is trivially sensible. assert _is_sensible_extension(ejm31, ejm31, mapping, dict(mapping)) + + +def test_unmapped_attachment_check_suppressed(ejm31, jmc28): + """ + The check should only fire where the user can act on its advice, i.e. + where 'mcs_kwargs' can be passed through. + """ + # These functions don't take 'mcs_kwargs'. + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.rmsdAlign(ejm31, jmc28) + BSS.Align.flexAlign(ejm31, jmc28) + + # 'merge' does, so the check stays on. + with pytest.warns(UserWarning, match="ringMatchesRingOnly"): + BSS.Align.merge(ejm31, jmc28, force=True) + + +def test_unmapped_attachment_warning_not_swallowed(ejm31, jmc28): + # Promoting the warning to an error must surface the warning itself, not + # a report that the check failed. The warning is emitted outside the + # try/except that guards the check for exactly this reason. + # Anchored, since the wrapped "Unable to check ..." message quotes the + # original and would otherwise match too. + with pytest.raises(UserWarning, match=r"^Mapping leaves heavy atoms"): + with warnings.catch_warnings(): + warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") + BSS.Align.matchAtoms(ejm31, jmc28) From ea14be897f42c2725d20a94b00d1a1872d0dacd3 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 14:33:59 +0100 Subject: [PATCH 14/23] Restore the viewMapping check and tidy the private parameter. --- src/BioSimSpace/Align/_align.py | 19 +++++++++++-------- .../Sandpit/Exscientia/Align/_align.py | 18 ++++++++++-------- tests/Align/test_align.py | 8 ++++++-- tests/Sandpit/Exscientia/Align/test_align.py | 8 ++++++-- 4 files changed, 33 insertions(+), 20 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 9f1c0580b..7b1391d70 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -744,6 +744,7 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + *, _check_mapping=True, ): """ @@ -1111,6 +1112,7 @@ def _matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + *, _check_mapping=True, ): import sys as _sys @@ -1334,10 +1336,11 @@ def _matchAtoms( # prematch is given, since the retry could then fall back on Sire MCS, # which ignores 'mcs_kwargs', making the comparison meaningless. if _check_mapping and not mcs_kwargs and not prematch and mappings: - # This is a diagnostic, so it must never be able to break a call that - # would otherwise have succeeded. The warning itself is emitted outside - # the guard, since it would otherwise be swallowed and re-reported as a - # failure whenever the user has promoted warnings to errors. + # This is a diagnostic, so a failure inside it is reported rather than + # raised. Note that this only holds while warnings are warnings: if the + # user has promoted them to errors then either notice below will raise + # out of here. The warning itself is emitted outside the guard, since + # it would otherwise be swallowed and re-reported as a failure. message = None try: best = mappings[0] @@ -1377,10 +1380,13 @@ def _matchAtoms( property_map0=property_map0, property_map1=property_map1, mcs_kwargs={"ringMatchesRingOnly": False}, + _check_mapping=False, ) # Only trust the retry if it extends the mapping, i.e. keeps - # every existing pair and adds sensible ones. + # every existing pair and adds sensible ones. The subset test is + # sensitive to symmetry: relabelling a symmetric ring the other + # way round discards the retry even though it is equivalent. if ( len(retry) > len(best) and set(best.items()) <= set(retry.items()) @@ -2632,9 +2638,6 @@ def viewMapping( molecule1, property_map0=property_map0, property_map1=property_map1, - # This function doesn't take 'mcs_kwargs', so the advice from the - # mapping check can't be acted on here. - _check_mapping=False, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index 554e002c2..d4be9b2e6 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -907,6 +907,7 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, + *, _check_mapping=True, ): """ @@ -1244,10 +1245,11 @@ def matchAtoms( # prematch is given, since the retry could then fall back on Sire MCS, # which ignores 'mcs_kwargs', making the comparison meaningless. if _check_mapping and not mcs_kwargs and not prematch and mappings: - # This is a diagnostic, so it must never be able to break a call that - # would otherwise have succeeded. The warning itself is emitted outside - # the guard, since it would otherwise be swallowed and re-reported as a - # failure whenever the user has promoted warnings to errors. + # This is a diagnostic, so a failure inside it is reported rather than + # raised. Note that this only holds while warnings are warnings: if the + # user has promoted them to errors then either notice below will raise + # out of here. The warning itself is emitted outside the guard, since + # it would otherwise be swallowed and re-reported as a failure. message = None try: best = mappings[0] @@ -1287,10 +1289,13 @@ def matchAtoms( property_map0=property_map0, property_map1=property_map1, mcs_kwargs={"ringMatchesRingOnly": False}, + _check_mapping=False, ) # Only trust the retry if it extends the mapping, i.e. keeps - # every existing pair and adds sensible ones. + # every existing pair and adds sensible ones. The subset test is + # sensitive to symmetry: relabelling a symmetric ring the other + # way round discards the retry even though it is equivalent. if ( len(retry) > len(best) and set(best.items()) <= set(retry.items()) @@ -1918,9 +1923,6 @@ def viewMapping( molecule1, property_map0=property_map0, property_map1=property_map1, - # This function doesn't take 'mcs_kwargs', so the advice from the - # mapping check can't be acted on here. - _check_mapping=False, ) molecule0 = rmsdAlign(molecule0, molecule1, mapping) diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index fa4295048..dca7b31be 100644 --- a/tests/Align/test_align.py +++ b/tests/Align/test_align.py @@ -1450,7 +1450,9 @@ def test_unmapped_attachment_no_warning_r_group(monkeypatch): pairs = [ ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("CCc1ccccc1", "CCCc1ccccc1"), # ethyl -> propyl ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ("O=C(N)c1ccccc1", "O=C(N)c1ccccc1Cl"), # hydrogen -> chlorine ] flagged = [] @@ -1602,11 +1604,13 @@ def test_unmapped_attachment_check_suppressed(ejm31, jmc28): the ROI path, where the flagged indices would be local to the extracted residue rather than to molecule0 as the message claims. """ - # These functions don't take 'mcs_kwargs'. + # 'rmsdAlign' doesn't take 'mcs_kwargs'. 'flexAlign' is suppressed for the + # same reason, but isn't exercised here since it needs fkcombu. Nor is + # 'viewMapping', which keeps the check but returns early outside a + # notebook, before it ever reaches 'matchAtoms'. with warnings.catch_warnings(): warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") BSS.Align.rmsdAlign(ejm31, jmc28) - BSS.Align.flexAlign(ejm31, jmc28) # 'merge' does, so the check stays on. with pytest.warns(UserWarning, match="ringMatchesRingOnly"): diff --git a/tests/Sandpit/Exscientia/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 0e413b826..6d265a166 100644 --- a/tests/Sandpit/Exscientia/Align/test_align.py +++ b/tests/Sandpit/Exscientia/Align/test_align.py @@ -1026,7 +1026,9 @@ def test_unmapped_attachment_no_warning_r_group(monkeypatch): pairs = [ ("Cc1ccccc1", "CCc1ccccc1"), # methyl -> ethyl + ("CCc1ccccc1", "CCCc1ccccc1"), # ethyl -> propyl ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ("O=C(N)c1ccccc1", "O=C(N)c1ccccc1Cl"), # hydrogen -> chlorine ] flagged = [] @@ -1176,11 +1178,13 @@ def test_unmapped_attachment_check_suppressed(ejm31, jmc28): The check should only fire where the user can act on its advice, i.e. where 'mcs_kwargs' can be passed through. """ - # These functions don't take 'mcs_kwargs'. + # 'rmsdAlign' doesn't take 'mcs_kwargs'. 'flexAlign' is suppressed for the + # same reason, but isn't exercised here since it needs fkcombu. Nor is + # 'viewMapping', which keeps the check but returns early outside a + # notebook, before it ever reaches 'matchAtoms'. with warnings.catch_warnings(): warnings.filterwarnings("error", message=".*ringMatchesRingOnly.*") BSS.Align.rmsdAlign(ejm31, jmc28) - BSS.Align.flexAlign(ejm31, jmc28) # 'merge' does, so the check stays on. with pytest.warns(UserWarning, match="ringMatchesRingOnly"): From 155728e8303b5a9a8f21e8260a532843bd018ebd Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 14:41:34 +0100 Subject: [PATCH 15/23] Note that the subset test is sensitive to equivalent relabellings. --- src/BioSimSpace/Align/_align.py | 4 ++-- src/BioSimSpace/Sandpit/Exscientia/Align/_align.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 7b1391d70..28060aef8 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -1385,8 +1385,8 @@ def _matchAtoms( # Only trust the retry if it extends the mapping, i.e. keeps # every existing pair and adds sensible ones. The subset test is - # sensitive to symmetry: relabelling a symmetric ring the other - # way round discards the retry even though it is equivalent. + # sensitive to relabelling: an equivalent mapping that traverses + # a ring the other way, or permutes hydrogens, is discarded. if ( len(retry) > len(best) and set(best.items()) <= set(retry.items()) diff --git a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index d4be9b2e6..db4675e34 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py @@ -1294,8 +1294,8 @@ def matchAtoms( # Only trust the retry if it extends the mapping, i.e. keeps # every existing pair and adds sensible ones. The subset test is - # sensitive to symmetry: relabelling a symmetric ring the other - # way round discards the retry even though it is equivalent. + # sensitive to relabelling: an equivalent mapping that traverses + # a ring the other way, or permutes hydrogens, is discarded. if ( len(retry) > len(best) and set(best.items()) <= set(retry.items()) From 934131b2e52cd921f24fff25b276944cd7fa5b8c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Fri, 31 Jul 2026 17:47:55 +0100 Subject: [PATCH 16/23] Take the mapping check flag off the public matchAtoms signature. --- src/BioSimSpace/Align/_align.py | 10 +++------- src/BioSimSpace/FreeEnergy/_atm.py | 9 +++++++-- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 28060aef8..72f0f50d7 100644 --- a/src/BioSimSpace/Align/_align.py +++ b/src/BioSimSpace/Align/_align.py @@ -744,8 +744,6 @@ def matchAtoms( property_map0={}, property_map1={}, mcs_kwargs={}, - *, - _check_mapping=True, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -914,7 +912,6 @@ def matchAtoms( property_map0=property_map0, property_map1=property_map1, mcs_kwargs=mcs_kwargs, - _check_mapping=_check_mapping, ) else: return _roiMatch( @@ -1105,7 +1102,6 @@ def _matchAtoms( timeout=5 * _Units.Time.second, complete_rings_only=True, max_scoring_matches=1000, - roi=None, prune_perturbed_constraints=False, prune_crossing_constraints=False, prune_atom_types=False, @@ -1744,7 +1740,7 @@ def _roiMatch( ) mapping = None else: - mapping = matchAtoms( + mapping = _matchAtoms( res0_extracted, res1_extracted, # The mapping check would report indices that are local to the @@ -1968,7 +1964,7 @@ def _rmsdAlign(molecule0, molecule1, mapping=None, property_map0={}, property_ma # Get the best match atom mapping. else: - mapping = matchAtoms( + mapping = _matchAtoms( molecule0, molecule1, property_map0=property_map0, @@ -2173,7 +2169,7 @@ def _flexAlign( # Get the best match atom mapping. else: - mapping = matchAtoms( + mapping = _matchAtoms( molecule0, molecule1, property_map0=property_map0, diff --git a/src/BioSimSpace/FreeEnergy/_atm.py b/src/BioSimSpace/FreeEnergy/_atm.py index b1f126f28..e7d237e51 100644 --- a/src/BioSimSpace/FreeEnergy/_atm.py +++ b/src/BioSimSpace/FreeEnergy/_atm.py @@ -593,7 +593,10 @@ def _makeSystemFromThree(protein, ligand_bound, ligand_free, displacement): BioSimSpace._SireWrappers.System The system for the ATM simulation. """ - from ..Align import matchAtoms as _matchAtoms + # The private form is used so that the mapping check can be switched + # off below. It is otherwise identical to 'matchAtoms', which is a + # dispatcher on 'roi' and is never given one here. + from ..Align._align import _matchAtoms from ..Align import rmsdAlign as _rmsdAlign from ..Types import Vector as _Vector @@ -683,7 +686,9 @@ def _findTranslationVector(system, displacement, protein, ligand): out_of_protein = displacement.value() * initial_normal_vector return out_of_protein - mapping = _matchAtoms(ligand_free, ligand_bound) + # ATM doesn't expose 'mcs_kwargs', nor a way to pass in a mapping, so + # the advice from the mapping check can't be acted on here. + mapping = _matchAtoms(ligand_free, ligand_bound, _check_mapping=False) ligand_free_aligned = _rmsdAlign(ligand_free, ligand_bound, mapping) prot_lig1 = (protein + ligand_bound).toSystem() From 4aaeefc74288e581e55fa30cffb88d08faa5b137 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 3 Aug 2026 16:20:53 +0100 Subject: [PATCH 17/23] Use native lazy import implementation. --- pixi.toml | 1 - recipes/biosimspace/recipe.yaml | 1 - .../Sandpit/Exscientia/__init__.py | 48 +++++++++---------- src/BioSimSpace/__init__.py | 44 ++++++++--------- 4 files changed, 44 insertions(+), 50 deletions(-) diff --git a/pixi.toml b/pixi.toml index d95194edf..34a4e009c 100644 --- a/pixi.toml +++ b/pixi.toml @@ -8,7 +8,6 @@ python = ">=3.10" configargparse = "*" ipywidgets = "*" kcombu_bss = "*" -lazy_import = "*" loguru = "*" lomap2 = "*" networkx = "*" diff --git a/recipes/biosimspace/recipe.yaml b/recipes/biosimspace/recipe.yaml index 964f5097d..a847c9bb8 100644 --- a/recipes/biosimspace/recipe.yaml +++ b/recipes/biosimspace/recipe.yaml @@ -22,7 +22,6 @@ requirements: - configargparse - ipywidgets - kcombu_bss - - lazy_import - loguru - lomap2 - networkx diff --git a/src/BioSimSpace/Sandpit/Exscientia/__init__.py b/src/BioSimSpace/Sandpit/Exscientia/__init__.py index eecc8a01f..554eb8856 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/__init__.py +++ b/src/BioSimSpace/Sandpit/Exscientia/__init__.py @@ -238,35 +238,31 @@ def _isVerbose(): # Lazy load submodules if possible. if _can_lazy_import: - import lazy_import as _lazy_import - - Align = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Align") - Box = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Box") - Convert = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Convert") - FreeEnergy = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.FreeEnergy") - Gateway = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Gateway") - IO = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.IO") - Metadynamics = _lazy_import.lazy_module( - "BioSimSpace.Sandpit.Exscientia.Metadynamics" - ) - MD = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.MD") - Node = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Node") - Notebook = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Notebook") - Parameters = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Parameters") - Process = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Process") - Protocol = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Protocol") - Solvent = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Solvent") - Stream = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Stream") - Trajectory = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Trajectory") - Types = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Types") - Units = _lazy_import.lazy_module("BioSimSpace.Sandpit.Exscientia.Units") - _SireWrappers = _lazy_import.lazy_module( - "BioSimSpace.Sandpit.Exscientia._SireWrappers" - ) + from sire._lazy_import import lazy_module as _lazy_module + + Align = _lazy_module("BioSimSpace.Sandpit.Exscientia.Align") + Box = _lazy_module("BioSimSpace.Sandpit.Exscientia.Box") + Convert = _lazy_module("BioSimSpace.Sandpit.Exscientia.Convert") + FreeEnergy = _lazy_module("BioSimSpace.Sandpit.Exscientia.FreeEnergy") + Gateway = _lazy_module("BioSimSpace.Sandpit.Exscientia.Gateway") + IO = _lazy_module("BioSimSpace.Sandpit.Exscientia.IO") + Metadynamics = _lazy_module("BioSimSpace.Sandpit.Exscientia.Metadynamics") + MD = _lazy_module("BioSimSpace.Sandpit.Exscientia.MD") + Node = _lazy_module("BioSimSpace.Sandpit.Exscientia.Node") + Notebook = _lazy_module("BioSimSpace.Sandpit.Exscientia.Notebook") + Parameters = _lazy_module("BioSimSpace.Sandpit.Exscientia.Parameters") + Process = _lazy_module("BioSimSpace.Sandpit.Exscientia.Process") + Protocol = _lazy_module("BioSimSpace.Sandpit.Exscientia.Protocol") + Solvent = _lazy_module("BioSimSpace.Sandpit.Exscientia.Solvent") + Stream = _lazy_module("BioSimSpace.Sandpit.Exscientia.Stream") + Trajectory = _lazy_module("BioSimSpace.Sandpit.Exscientia.Trajectory") + Types = _lazy_module("BioSimSpace.Sandpit.Exscientia.Types") + Units = _lazy_module("BioSimSpace.Sandpit.Exscientia.Units") + _SireWrappers = _lazy_module("BioSimSpace.Sandpit.Exscientia._SireWrappers") from . import _Exceptions, _Utils - del _lazy_import + del _lazy_module else: from . import ( IO, diff --git a/src/BioSimSpace/__init__.py b/src/BioSimSpace/__init__.py index e82999159..26e28ce4f 100644 --- a/src/BioSimSpace/__init__.py +++ b/src/BioSimSpace/__init__.py @@ -238,31 +238,31 @@ def _isVerbose(): # Lazy import submodules if possible. if _can_lazy_import: - import lazy_import as _lazy_import - - Align = _lazy_import.lazy_module("BioSimSpace.Align") - Box = _lazy_import.lazy_module("BioSimSpace.Box") - Convert = _lazy_import.lazy_module("BioSimSpace.Convert") - FreeEnergy = _lazy_import.lazy_module("BioSimSpace.FreeEnergy") - Gateway = _lazy_import.lazy_module("BioSimSpace.Gateway") - IO = _lazy_import.lazy_module("BioSimSpace.IO") - Metadynamics = _lazy_import.lazy_module("BioSimSpace.Metadynamics") - MD = _lazy_import.lazy_module("BioSimSpace.MD") - Node = _lazy_import.lazy_module("BioSimSpace.Node") - Notebook = _lazy_import.lazy_module("BioSimSpace.Notebook") - Parameters = _lazy_import.lazy_module("BioSimSpace.Parameters") - Process = _lazy_import.lazy_module("BioSimSpace.Process") - Protocol = _lazy_import.lazy_module("BioSimSpace.Protocol") - Solvent = _lazy_import.lazy_module("BioSimSpace.Solvent") - Stream = _lazy_import.lazy_module("BioSimSpace.Stream") - Trajectory = _lazy_import.lazy_module("BioSimSpace.Trajectory") - Types = _lazy_import.lazy_module("BioSimSpace.Types") - Units = _lazy_import.lazy_module("BioSimSpace.Units") - _SireWrappers = _lazy_import.lazy_module("BioSimSpace._SireWrappers") + from sire._lazy_import import lazy_module as _lazy_module + + Align = _lazy_module("BioSimSpace.Align") + Box = _lazy_module("BioSimSpace.Box") + Convert = _lazy_module("BioSimSpace.Convert") + FreeEnergy = _lazy_module("BioSimSpace.FreeEnergy") + Gateway = _lazy_module("BioSimSpace.Gateway") + IO = _lazy_module("BioSimSpace.IO") + Metadynamics = _lazy_module("BioSimSpace.Metadynamics") + MD = _lazy_module("BioSimSpace.MD") + Node = _lazy_module("BioSimSpace.Node") + Notebook = _lazy_module("BioSimSpace.Notebook") + Parameters = _lazy_module("BioSimSpace.Parameters") + Process = _lazy_module("BioSimSpace.Process") + Protocol = _lazy_module("BioSimSpace.Protocol") + Solvent = _lazy_module("BioSimSpace.Solvent") + Stream = _lazy_module("BioSimSpace.Stream") + Trajectory = _lazy_module("BioSimSpace.Trajectory") + Types = _lazy_module("BioSimSpace.Types") + Units = _lazy_module("BioSimSpace.Units") + _SireWrappers = _lazy_module("BioSimSpace._SireWrappers") from . import _Exceptions, _Utils - del _lazy_import + del _lazy_module else: from . import ( IO, From 6b494d3f2a3f0370554c68f4225d4eb69d16a06c Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 24 Aug 2026 09:12:28 +0100 Subject: [PATCH 18/23] Match tar files by extension. [closes #547] --- src/BioSimSpace/Gateway/_requirements.py | 8 ++++---- .../Sandpit/Exscientia/Gateway/_requirements.py | 8 ++++---- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/BioSimSpace/Gateway/_requirements.py b/src/BioSimSpace/Gateway/_requirements.py index 0ca478fb0..fe6c1b850 100644 --- a/src/BioSimSpace/Gateway/_requirements.py +++ b/src/BioSimSpace/Gateway/_requirements.py @@ -2037,13 +2037,13 @@ def _unarchive(name): else: dir = _os.path.splitext(name)[0] + "/" - # List of supported tar file formats. - tarfiles = ["tar.gz", "tar.bz2", "tar"] + # List of supported tar file extensions. + tar_exts = [".tar.gz", ".tar.bz2", ".tar"] # Check whether this is a tar compressed file. - for tar_name in tarfiles: + for ext in tar_exts: # Found a match. - if tar_name in name.lower(): + if name.lower().endswith(ext): # The list of decompressed files. files = [] diff --git a/src/BioSimSpace/Sandpit/Exscientia/Gateway/_requirements.py b/src/BioSimSpace/Sandpit/Exscientia/Gateway/_requirements.py index 0ca478fb0..fe6c1b850 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Gateway/_requirements.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Gateway/_requirements.py @@ -2037,13 +2037,13 @@ def _unarchive(name): else: dir = _os.path.splitext(name)[0] + "/" - # List of supported tar file formats. - tarfiles = ["tar.gz", "tar.bz2", "tar"] + # List of supported tar file extensions. + tar_exts = [".tar.gz", ".tar.bz2", ".tar"] # Check whether this is a tar compressed file. - for tar_name in tarfiles: + for ext in tar_exts: # Found a match. - if tar_name in name.lower(): + if name.lower().endswith(ext): # The list of decompressed files. files = [] From 01ef8127655110239a5b5fc545fd8bd75282596f Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Wed, 26 Aug 2026 14:48:31 +0100 Subject: [PATCH 19/23] Add determine_bond_orders kwarg to RDKit conversion function. --- src/BioSimSpace/Convert/_convert.py | 20 ++++++++++++++++++- .../Sandpit/Exscientia/Convert/_convert.py | 20 ++++++++++++++++++- 2 files changed, 38 insertions(+), 2 deletions(-) diff --git a/src/BioSimSpace/Convert/_convert.py b/src/BioSimSpace/Convert/_convert.py index c71ead665..817909703 100644 --- a/src/BioSimSpace/Convert/_convert.py +++ b/src/BioSimSpace/Convert/_convert.py @@ -194,6 +194,13 @@ def to(obj, format="biosimspace", property_map={}, **kwargs): raise TypeError("'force_stereo_inference' must be of type 'bool'.") property_map["force_stereo_inference"] = _SireBase.wrap(force_stereo_inference) + # Check for determine_bond_orders in kwargs. + if "determine_bond_orders" in kwargs: + determine_bond_orders = kwargs["determine_bond_orders"] + if not isinstance(determine_bond_orders, bool): + raise TypeError("'determine_bond_orders' must be of type 'bool'.") + property_map["determine_bond_orders"] = _SireBase.wrap(determine_bond_orders) + # Special handling for OpenMM conversion. Currently this is a one-way (toOpenMM) # conversion only and is only supported for specific Sire and BioSimSpace types. if format == "openmm": @@ -517,7 +524,9 @@ def toOpenMM(obj, property_map={}): ) -def toRDKit(obj, force_stereo_inference=False, property_map={}): +def toRDKit( + obj, force_stereo_inference=False, determine_bond_orders=True, property_map={} +): """ Convert an object to RDKit format. @@ -532,6 +541,11 @@ def toRDKit(obj, force_stereo_inference=False, property_map={}): stereochemistry present in the input object. This is useful when the object has been loaded from a file with invalid stereochemistry. + bool : determine_bond_orders + Whether to use RDKit's determineBondOrders function when bond orders + need to be inferred. This is more robust than the internal heuristic, + but can be slow for large molecules, e.g. proteins. + property_map : dict A dictionary that maps system "properties" to their user defined values. This allows the user to refer to properties with their @@ -548,10 +562,14 @@ def toRDKit(obj, force_stereo_inference=False, property_map={}): if not isinstance(force_stereo_inference, bool): raise TypeError("'force_stereo_inference' must be of type 'bool'.") + if not isinstance(determine_bond_orders, bool): + raise TypeError("'determine_bond_orders' must be of type 'bool'.") + if not isinstance(property_map, dict): raise TypeError("'property_map' must be of type 'dict'.") property_map["force_stereo_inference"] = _SireBase.wrap(force_stereo_inference) + property_map["determine_bond_orders"] = _SireBase.wrap(determine_bond_orders) return to(obj, format="rdkit", property_map=property_map) diff --git a/src/BioSimSpace/Sandpit/Exscientia/Convert/_convert.py b/src/BioSimSpace/Sandpit/Exscientia/Convert/_convert.py index c71ead665..817909703 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Convert/_convert.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Convert/_convert.py @@ -194,6 +194,13 @@ def to(obj, format="biosimspace", property_map={}, **kwargs): raise TypeError("'force_stereo_inference' must be of type 'bool'.") property_map["force_stereo_inference"] = _SireBase.wrap(force_stereo_inference) + # Check for determine_bond_orders in kwargs. + if "determine_bond_orders" in kwargs: + determine_bond_orders = kwargs["determine_bond_orders"] + if not isinstance(determine_bond_orders, bool): + raise TypeError("'determine_bond_orders' must be of type 'bool'.") + property_map["determine_bond_orders"] = _SireBase.wrap(determine_bond_orders) + # Special handling for OpenMM conversion. Currently this is a one-way (toOpenMM) # conversion only and is only supported for specific Sire and BioSimSpace types. if format == "openmm": @@ -517,7 +524,9 @@ def toOpenMM(obj, property_map={}): ) -def toRDKit(obj, force_stereo_inference=False, property_map={}): +def toRDKit( + obj, force_stereo_inference=False, determine_bond_orders=True, property_map={} +): """ Convert an object to RDKit format. @@ -532,6 +541,11 @@ def toRDKit(obj, force_stereo_inference=False, property_map={}): stereochemistry present in the input object. This is useful when the object has been loaded from a file with invalid stereochemistry. + bool : determine_bond_orders + Whether to use RDKit's determineBondOrders function when bond orders + need to be inferred. This is more robust than the internal heuristic, + but can be slow for large molecules, e.g. proteins. + property_map : dict A dictionary that maps system "properties" to their user defined values. This allows the user to refer to properties with their @@ -548,10 +562,14 @@ def toRDKit(obj, force_stereo_inference=False, property_map={}): if not isinstance(force_stereo_inference, bool): raise TypeError("'force_stereo_inference' must be of type 'bool'.") + if not isinstance(determine_bond_orders, bool): + raise TypeError("'determine_bond_orders' must be of type 'bool'.") + if not isinstance(property_map, dict): raise TypeError("'property_map' must be of type 'dict'.") property_map["force_stereo_inference"] = _SireBase.wrap(force_stereo_inference) + property_map["determine_bond_orders"] = _SireBase.wrap(determine_bond_orders) return to(obj, format="rdkit", property_map=property_map) From 6c8d5ffe8d99241202c996ae62349ccec7ac02bf Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Tue, 1 Sep 2026 15:28:02 +0100 Subject: [PATCH 20/23] Fix typos and inconsistencies. [ci skip] --- README.rst | 23 +++++++++-------------- doc/source/install.rst | 5 ++--- 2 files changed, 11 insertions(+), 17 deletions(-) diff --git a/README.rst b/README.rst index 47030af89..e36742b74 100644 --- a/README.rst +++ b/README.rst @@ -2,7 +2,7 @@ ==================================================== .. image:: https://github.com/openbiosim/biosimspace/actions/workflows/devel.yaml/badge.svg - :target: https://github.com/openbiosim/biosimspace/actions?query=workflow%3ARelease-Devel + :target: https://github.com/openbiosim/biosimspace/actions?query=workflow%3A%22Release+Devel%22 :alt: Build status .. image:: https://anaconda.org/openbiosim/biosimspace/badges/downloads.svg @@ -29,9 +29,9 @@ for biomolecular simulation. With it you can: * Start, stop, and monitor molecular simulation processes within interactive Python environments. Citation |DOI for Citing BioSimSpace| -===================================== +------------------------------------- -If you use BioSimSpace in any scientific software, please cite the following paper: :: +If you use BioSimSpace in any scientific work, please cite the following paper: :: @article{Hedges2019, doi = {10.21105/joss.01831}, @@ -79,19 +79,15 @@ To install the latest development version you can use: conda create -n openbiosim-dev -c conda-forge -c openbiosim/label/dev biosimspace conda activate openbiosim-dev -When updating the development version it is generally advised to update `Sire `_ +When updating the development version it is generally advised to update `Sire `__ at the same time: .. code-block:: bash conda update -c conda-forge -c openbiosim/label/dev biosimspace sire -Unless you add the required channels to your Conda configuration, then you'll -need to add them when updating, e.g., for the development package: - -.. code-block:: bash - - conda update -c conda-forge -c openbiosim/label/dev biosimspace +Unless you add the required channels to your Conda configuration, you'll need to +pass them on the command line when updating, as shown above. Installing from source (standalone) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ @@ -133,12 +129,11 @@ You may also want to install optional dependencies, such as ``ambertools`` and If you need OpenCL support (e.g. for OpenMM), note that pixi does not run conda post-link scripts, so the ``ocl-icd-system`` symlink won't be created -automatically. After creating the environment, run the following once to fix -this: +automatically. This applies to either of the routes above. From within the +activated environment, run the following once to fix this: .. code-block:: bash - pixi shell -e dev ln -s /etc/OpenCL/vendors "${CONDA_PREFIX}/etc/OpenCL/vendors/ocl-icd-system" Once finished, you can test the installation by running: @@ -176,7 +171,7 @@ Issues Please report bugs and other issues using the GitHub `issue tracker `__. When reporting issues please try to include a minimal code snippet that reproduces -the problem. Additional files can be also be uploaded as an archive, e.g. a zip +the problem. Additional files can also be uploaded as an archive, e.g. a zip file. Please also report the branch on which you are experiencing the issue, along with the BioSimSpace version number. This can be found by running: diff --git a/doc/source/install.rst b/doc/source/install.rst index 2d5297945..0e8b7242f 100644 --- a/doc/source/install.rst +++ b/doc/source/install.rst @@ -245,12 +245,11 @@ You may also want to install optional dependencies, such as ``ambertools`` and If you need OpenCL support (e.g. for OpenMM), note that pixi does not run conda post-link scripts, so the ``ocl-icd-system`` symlink won't be created -automatically. After creating the environment, run the following once to fix -this: +automatically. This applies to either of the routes above. From within the +activated environment, run the following once to fix this: .. code-block:: bash - pixi shell -e dev ln -s /etc/OpenCL/vendors "${CONDA_PREFIX}/etc/OpenCL/vendors/ocl-icd-system" Once finished, you can test the installation by running: From 6d008b1f553159d3741b357a79e9dade17ebe4dc Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 7 Sep 2026 10:31:08 +0100 Subject: [PATCH 21/23] Add functionality for parsing CMAP energy records. --- src/BioSimSpace/Process/_amber.py | 68 +++++++++++++++++++ src/BioSimSpace/Process/_gromacs.py | 41 +++++++++++ .../Sandpit/Exscientia/Process/_amber.py | 68 +++++++++++++++++++ .../Sandpit/Exscientia/Process/_gromacs.py | 41 +++++++++++ tests/Process/test_single_point_energy.py | 49 +++++++++++++ .../Process/test_single_point_energy.py | 49 +++++++++++++ 6 files changed, 316 insertions(+) diff --git a/src/BioSimSpace/Process/_amber.py b/src/BioSimSpace/Process/_amber.py index e7c94c254..f8dd37342 100644 --- a/src/BioSimSpace/Process/_amber.py +++ b/src/BioSimSpace/Process/_amber.py @@ -1413,6 +1413,74 @@ def getCurrentDihedralEnergy(self, time_series=False, region=0, soft_core=False) time_series=time_series, region=region, soft_core=soft_core, block=False ) + def getCMAPEnergy(self, time_series=False, region=0, soft_core=False, block="AUTO"): + """ + Get the CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + region : int + The region to which the record corresponds. There will only be more + than one region for FreeEnergy protocols, where 1 indicates the second + TI region. + + soft_core : bool + Whether to get the record for the soft-core part of the system for the + chosen region. + + block : bool + Whether to block until the process has finished running. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + from .. import Units as _Units + + return self.getRecord( + "CMAP", + time_series=time_series, + unit=_Units.Energy.kcal_per_mol, + region=region, + soft_core=soft_core, + block=block, + ) + + def getCurrentCMAPEnergy(self, time_series=False, region=0, soft_core=False): + """ + Get the current CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + region : int + The region to which the record corresponds. There will only be more + than one region for FreeEnergy protocols, where 1 indicates the second + TI region. + + soft_core : bool + Whether to get the record for the soft-core part of the system for the + chosen region. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + return self.getCMAPEnergy( + time_series=time_series, region=region, soft_core=soft_core, block=False + ) + def getElectrostaticEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): diff --git a/src/BioSimSpace/Process/_gromacs.py b/src/BioSimSpace/Process/_gromacs.py index bdcb49f0b..16c8e6f68 100644 --- a/src/BioSimSpace/Process/_gromacs.py +++ b/src/BioSimSpace/Process/_gromacs.py @@ -1382,6 +1382,47 @@ def getCurrentImproperEnergy(self, time_series=False): """ return self.getImproperEnergy(time_series, block=False) + def getCMAPEnergy(self, time_series=False, block="AUTO"): + """ + Get the CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + block : bool + Whether to block until the process has finished running. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + from .. import Units as _Units + + return self.getRecord("CMAPDIH", time_series, _Units.Energy.kj_per_mol, block) + + def getCurrentCMAPEnergy(self, time_series=False): + """ + Get the current CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + return self.getCMAPEnergy(time_series, block=False) + def getLennardJones14(self, time_series=False, block="AUTO"): """ Get the Lennard-Jones energy between atoms 1 and 4. diff --git a/src/BioSimSpace/Sandpit/Exscientia/Process/_amber.py b/src/BioSimSpace/Sandpit/Exscientia/Process/_amber.py index add46a4fb..36a52762b 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Process/_amber.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Process/_amber.py @@ -1497,6 +1497,74 @@ def getCurrentDihedralEnergy(self, time_series=False, region=0, soft_core=False) time_series=time_series, region=region, soft_core=soft_core, block=False ) + def getCMAPEnergy(self, time_series=False, region=0, soft_core=False, block="AUTO"): + """ + Get the CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + region : int + The region to which the record corresponds. There will only be more + than one region for FreeEnergy protocols, where 1 indicates the second + TI region. + + soft_core : bool + Whether to get the record for the soft-core part of the system for the + chosen region. + + block : bool + Whether to block until the process has finished running. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + from .. import Units as _Units + + return self.getRecord( + "CMAP", + time_series=time_series, + unit=_Units.Energy.kcal_per_mol, + region=region, + soft_core=soft_core, + block=block, + ) + + def getCurrentCMAPEnergy(self, time_series=False, region=0, soft_core=False): + """ + Get the current CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + region : int + The region to which the record corresponds. There will only be more + than one region for FreeEnergy protocols, where 1 indicates the second + TI region. + + soft_core : bool + Whether to get the record for the soft-core part of the system for the + chosen region. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + return self.getCMAPEnergy( + time_series=time_series, region=region, soft_core=soft_core, block=False + ) + def getElectrostaticEnergy( self, time_series=False, region=0, soft_core=False, block="AUTO" ): diff --git a/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py b/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py index c90c47914..b7125b955 100644 --- a/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py +++ b/src/BioSimSpace/Sandpit/Exscientia/Process/_gromacs.py @@ -1440,6 +1440,47 @@ def getCurrentImproperEnergy(self, time_series=False): """ return self.getImproperEnergy(time_series, block=False) + def getCMAPEnergy(self, time_series=False, block="AUTO"): + """ + Get the CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + block : bool + Whether to block until the process has finished running. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + from .. import Units as _Units + + return self.getRecord("CMAPDIH", time_series, _Units.Energy.kj_per_mol, block) + + def getCurrentCMAPEnergy(self, time_series=False): + """ + Get the current CMAP energy. + + Parameters + ---------- + + time_series : bool + Whether to return a list of time series records. + + Returns + ------- + + energy : :class:`Energy ` + The CMAP energy. + """ + return self.getCMAPEnergy(time_series, block=False) + def getLennardJones14(self, time_series=False, block="AUTO"): """ Get the Lennard-Jones energy between atoms 1 and 4. diff --git a/tests/Process/test_single_point_energy.py b/tests/Process/test_single_point_energy.py index 99c6e501f..b86dda785 100644 --- a/tests/Process/test_single_point_energy.py +++ b/tests/Process/test_single_point_energy.py @@ -12,6 +12,16 @@ def ubiquitin_system(): ) +@pytest.fixture(scope="module") +def cmap_system(): + """An ff19SB system, which carries CMAP backbone correction terms.""" + import sire as sr + + return BSS._SireWrappers.System( + sr.load_test_files("zero_k_torsions.gro", "zero_k_torsions.top")._system + ) + + @pytest.mark.skipif( has_amber is False or has_gromacs is False, reason="Requires that both AMBER and GROMACS are installed.", @@ -100,3 +110,42 @@ def test_amber_gromacs_triclinic(ubiquitin_system): nrg_amb = process_amb.getDihedralEnergy().kj_per_mol().value() nrg_gmx = process_gmx.getDihedralEnergy().kj_per_mol().value() assert nrg_amb == pytest.approx(nrg_gmx, rel=1e-2) + + +@pytest.mark.skipif( + has_amber is False or has_gromacs is False, + reason="Requires that both AMBER and GROMACS are installed.", +) +def test_amber_gromacs_cmap(cmap_system): + """Single point CMAP energy comparison between AMBER and GROMACS.""" + + # Create a single-step minimisation protocol. + protocol = BSS.Protocol.Minimisation(steps=1) + + # Create a process to run with AMBER. + process_amb = BSS.Process.Amber(cmap_system, protocol) + + # Create a process to run with GROMACS. + process_gmx = BSS.Process.Gromacs( + cmap_system, protocol, extra_options={"nsteps": 0} + ) + + # Run the AMBER process and wait for it to finish. + process_amb.start() + process_amb.wait() + + # Run the GROMACS process and wait for it to finish. + process_gmx.start() + process_gmx.wait() + + # Compare CMAP energies. (In kJ / mol) + nrg_amb = process_amb.getCMAPEnergy() + nrg_gmx = process_gmx.getCMAPEnergy() + + # The comparison is meaningless if either engine didn't report the term. + assert nrg_amb is not None + assert nrg_gmx is not None + + assert nrg_amb.kj_per_mol().value() == pytest.approx( + nrg_gmx.kj_per_mol().value(), rel=1e-2 + ) diff --git a/tests/Sandpit/Exscientia/Process/test_single_point_energy.py b/tests/Sandpit/Exscientia/Process/test_single_point_energy.py index 3fe986c0d..5fb6e089e 100644 --- a/tests/Sandpit/Exscientia/Process/test_single_point_energy.py +++ b/tests/Sandpit/Exscientia/Process/test_single_point_energy.py @@ -12,6 +12,16 @@ def system(): ) +@pytest.fixture(scope="session") +def cmap_system(): + """An ff19SB system, which carries CMAP backbone correction terms.""" + import sire as sr + + return BSS._SireWrappers.System( + sr.load_test_files("zero_k_torsions.gro", "zero_k_torsions.top")._system + ) + + @pytest.mark.skipif( has_amber is False or has_gromacs is False or has_pyarrow is False, reason="Requires that AMBER, GROMACS, and pyarrow are installed.", @@ -96,3 +106,42 @@ def test_amber_gromacs_triclinic(system): nrg_amb = process_amb.getDihedralEnergy().kj_per_mol().value() nrg_gmx = process_gmx.getDihedralEnergy().kj_per_mol().value() assert nrg_amb == pytest.approx(nrg_gmx, rel=1e-2) + + +@pytest.mark.skipif( + has_amber is False or has_gromacs is False or has_pyarrow is False, + reason="Requires that AMBER, GROMACS, and pyarrow are installed.", +) +def test_amber_gromacs_cmap(cmap_system): + """Single point CMAP energy comparison between AMBER and GROMACS.""" + + # Create a single-step minimisation protocol. + protocol = BSS.Protocol.Minimisation(steps=1) + + # Create a process to run with AMBER. + process_amb = BSS.Process.Amber(cmap_system, protocol) + + # Create a process to run with GROMACS. + process_gmx = BSS.Process.Gromacs( + cmap_system, protocol, extra_options={"nsteps": 0} + ) + + # Run the AMBER process and wait for it to finish. + process_amb.start() + process_amb.wait() + + # Run the GROMACS process and wait for it to finish. + process_gmx.start() + process_gmx.wait() + + # Compare CMAP energies. (In kJ / mol) + nrg_amb = process_amb.getCMAPEnergy() + nrg_gmx = process_gmx.getCMAPEnergy() + + # The comparison is meaningless if either engine didn't report the term. + assert nrg_amb is not None + assert nrg_gmx is not None + + assert nrg_amb.kj_per_mol().value() == pytest.approx( + nrg_gmx.kj_per_mol().value(), rel=1e-2 + ) From 4defa62198dbe24d5f1e03e347bc15d87102cb8b Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 14 Sep 2026 14:39:46 +0100 Subject: [PATCH 22/23] Update Sire version. --- pixi.toml | 4 ++-- recipes/biosimspace/recipe.yaml | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/pixi.toml b/pixi.toml index 34a4e009c..502aee409 100644 --- a/pixi.toml +++ b/pixi.toml @@ -24,9 +24,9 @@ rdkit = "*" [target.linux-64.dependencies] # main -#sire = ">=2026.1.0,<2026.2.0" +sire = ">=2026.2.0,<2026.3.0" # devel -sire = "==2026.2.0.dev" +#sire = "==2026.3.0.dev" ambertools = ">=22" gromacs = "*" alchemlyb = "*" diff --git a/recipes/biosimspace/recipe.yaml b/recipes/biosimspace/recipe.yaml index a847c9bb8..c3b0cc079 100644 --- a/recipes/biosimspace/recipe.yaml +++ b/recipes/biosimspace/recipe.yaml @@ -42,9 +42,9 @@ requirements: - pyyaml - rdkit # main - #- sire >=2026.1.0,<2026.2.0 + - sire >=2026.2.0,<2026.3.0 # devel - - sire ==2026.2.0.dev + #- sire ==2026.3.0.dev - if: not aarch64 then: - alchemlyb From 477dba056adc402a66a2960a3fb50300b37e6f16 Mon Sep 17 00:00:00 2001 From: Lester Hedges Date: Mon, 14 Sep 2026 14:40:16 +0100 Subject: [PATCH 23/23] Update CHANGELOG. --- doc/source/changelog.rst | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/doc/source/changelog.rst b/doc/source/changelog.rst index 469e759e9..ecc1630d7 100644 --- a/doc/source/changelog.rst +++ b/doc/source/changelog.rst @@ -9,6 +9,18 @@ company supporting open-source development of fostering academic/industrial coll within the biomolecular simulation community. Our software is hosted via the `OpenBioSim` `GitHub `__ organisation. +`2026.2.0 `_ - Sep 15 2026 +------------------------------------------------------------------------------------------------- + +* Fixed translation of molecules via a custom coordinates property (`#540 `__). +* Improve detection of ring-breaking perturbations during merge (`#542 `__). +* Exposed ``RDKit`` MCS options via a keyword argument (`#542 `__). +* Added a diagnostic to flag sub-optimal MCS mappings and suggest alternative ``mcs_kwargs`` (`#544 `__). +* Switched to a native lazy import system to avoid duplicate class objects across process boundaries and parallelisation issues (`#546 `__). +* Match tar files by common extensions, rather than using a partial fuzzy match (`#548 `__). +* Exposed the new ``determine_bond_orders`` kwarg in the :func:`toRDKit` function (`#550 `__). +* Added functionality to parse AMBER and GROMACS CMAP energy records (`#552 `__). + `2026.1.0 `_ - Jun 29 2026 -------------------------------------------------------------------------------------------------