diff --git a/README.rst b/README.rst index 7bf08bd71..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: @@ -184,3 +179,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``. 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..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 ------------------------------------------------------------------------------------------------- @@ -19,8 +31,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 +44,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 +65,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 +98,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 +131,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 +150,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 +163,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 +194,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 +219,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 +260,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 +276,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 +324,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 +341,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 +362,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 +394,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 +443,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 +456,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 ========= 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: diff --git a/pixi.toml b/pixi.toml index 25747b305..502aee409 100644 --- a/pixi.toml +++ b/pixi.toml @@ -8,7 +8,6 @@ python = ">=3.10" configargparse = "*" ipywidgets = "*" kcombu_bss = "*" -lazy_import = "*" loguru = "*" lomap2 = "*" networkx = "*" @@ -25,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 = "*" @@ -36,13 +35,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 +53,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 5ddb5b831..c3b0cc079 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 @@ -43,9 +42,9 @@ requirements: - pyyaml - rdkit # main - - sire >=2026.1.0,<2026.2.0 + - sire >=2026.2.0,<2026.3.0 # devel - #- sire ==2026.1.0.dev + #- sire ==2026.3.0.dev - if: not aarch64 then: - alchemlyb diff --git a/src/BioSimSpace/Align/_align.py b/src/BioSimSpace/Align/_align.py index 2602bb5cc..72f0f50d7 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( @@ -896,6 +927,171 @@ def matchAtoms( ) +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 + 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. + + 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 + ------- + + flagged : [int] + The indices of the flagged atoms in molecule0. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + # 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() + 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, mapped0, element0) + if not elements0: + continue + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1, element1) + if elements0 & elements1: + flagged.append(idx0) + + return flagged + + +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 + 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. + + 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 + ------- + + 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() + + 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(element0).num_protons() + ) + protons1 = ( + mol1.atom(_SireMol.AtomIdx(idx1)).property(element1).num_protons() + ) + if (protons0 > 1) != (protons1 > 1): + return False + + return True + + def _matchAtoms( molecule0, molecule1, @@ -906,12 +1102,14 @@ 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, property_map0={}, property_map1={}, + mcs_kwargs={}, + *, + _check_mapping=True, ): import sys as _sys @@ -975,6 +1173,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'") @@ -1000,7 +1201,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. @@ -1012,24 +1216,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 +1267,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 @@ -1118,6 +1322,87 @@ 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 _check_mapping and not mcs_kwargs and not prematch and mappings: + # 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] + + # 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}, + _check_mapping=False, + ) + + # 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 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()) + and _is_sensible_extension( + molecule0, molecule1, best, retry, property_map0, property_map1 + ) + ): + 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 " + 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 message is not None: + _warnings.warn(message) + # Optionally post-process the MCS for use with AMBER. if prune_perturbed_constraints: mappings = [ @@ -1455,9 +1740,12 @@ def _roiMatch( ) mapping = None else: - mapping = matchAtoms( + 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. @@ -1676,11 +1964,14 @@ 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, 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. @@ -1878,11 +2169,14 @@ def _flexAlign( # Get the best match atom mapping. else: - mapping = matchAtoms( + mapping = _matchAtoms( molecule0, 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. @@ -2079,6 +2373,7 @@ def merge( roi=None, property_map0={}, property_map1={}, + mcs_kwargs={}, **kwargs, ): """ @@ -2130,6 +2425,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 +2508,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/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() 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/Align/_align.py b/src/BioSimSpace/Sandpit/Exscientia/Align/_align.py index 2f7170fcb..db4675e34 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,194 @@ def generateNetwork( return edges, scores +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 + 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. + + 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 + ------- + + flagged : [int] + The indices of the flagged atoms in molecule0. + """ + from sire.legacy import Mol as _SireMol + + mol0 = molecule0._getSireObject() + mol1 = molecule1._getSireObject() + + # 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() + 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, mapped0, element0) + if not elements0: + continue + elements1 = _heavy_elements(mol1, conn1, idx1, mapped1, element1) + if elements0 & elements1: + flagged.append(idx0) + + return flagged + + +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 + 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. + + 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 + ------- + + 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() + + 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(element0).num_protons() + ) + protons1 = ( + mol1.atom(_SireMol.AtomIdx(idx1)).property(element1).num_protons() + ) + if (protons0 > 1) != (protons1 > 1): + return False + + return True + + +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 +906,9 @@ def matchAtoms( max_scoring_matches=1000, property_map0={}, property_map1={}, + mcs_kwargs={}, + *, + _check_mapping=True, ): """ Find mappings between atom indices in molecule0 to those in molecule1. @@ -770,6 +962,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 @@ -918,7 +1115,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. @@ -930,24 +1129,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 +1180,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 @@ -1036,6 +1235,87 @@ 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 _check_mapping and not mcs_kwargs and not prematch and mappings: + # 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] + + # 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, + molecule1, + engine=engine, + scoring_function=scoring_function, + matches=1, + return_scores=False, + timeout=orig_timeout, + complete_rings_only=complete_rings_only, + prune_perturbed_constraints=False, + prune_crossing_constraints=False, + max_scoring_matches=max_scoring_matches, + 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. The subset test is + # 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()) + and _is_sensible_extension( + molecule0, molecule1, best, retry, property_map0, property_map1 + ) + ): + 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 " + 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 message is not None: + _warnings.warn(message) + # Optionally post-process the MCS. if prune_perturbed_constraints: mappings = [ @@ -1147,6 +1427,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. @@ -1304,6 +1587,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. @@ -1370,6 +1656,7 @@ def merge( roi=None, property_map0={}, property_map1={}, + mcs_kwargs={}, **kwargs, ): """ @@ -1417,6 +1704,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 +1780,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/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/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, diff --git a/tests/Align/test_align.py b/tests/Align/test_align.py index be01565e5..dca7b31be 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 @@ -674,10 +675,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 +692,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 +745,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 +764,6 @@ def test_custom_roi_map_invalid_outside_roi(): molecule0=wt, molecule1=mut, roi=[15], - custom_roi_map={ 0: 0, 1: 1, @@ -1313,9 +1315,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 +1328,318 @@ 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") + + +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. 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.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 + + +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 + ("CCc1ccccc1", "CCCc1ccccc1"), # ethyl -> propyl + ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ("O=C(N)c1ccccc1", "O=C(N)c1ccccc1Cl"), # hydrogen -> chlorine + ] + + 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)) + + +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. + """ + # '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) + + # '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/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/Align/test_align.py b/tests/Sandpit/Exscientia/Align/test_align.py index 18f6c86f0..6d265a166 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 @@ -6,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") @@ -904,3 +910,294 @@ 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") + + +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. 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.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 + + +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 + ("CCc1ccccc1", "CCCc1ccccc1"), # ethyl -> propyl + ("COc1ccccc1", "CCOc1ccccc1"), # methoxy -> ethoxy + ("O=C(N)c1ccccc1", "O=C(N)c1ccccc1Cl"), # hydrogen -> chlorine + ] + + 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)) + + +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. + """ + # '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) + + # '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) 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 + )