diff --git a/.github/workflows/devel.yaml b/.github/workflows/devel.yaml index 3ab8dab..ba36088 100644 --- a/.github/workflows/devel.yaml +++ b/.github/workflows/devel.yaml @@ -3,7 +3,7 @@ name: Release Devel on: workflow_dispatch: push: - branches: [ devel ] + branches: [devel] jobs: build: @@ -18,11 +18,9 @@ jobs: - { name: "linux", os: "ubuntu-latest", shell: "bash -l {0}" } - { name: "macos", os: "macos-latest", shell: "bash -l {0}" } exclude: - # Exclude all but the latest Python from all - # but Linux - platform: { name: "macos", os: "macos-latest", shell: "bash -l {0}" } - python-version: "3.12" # MacOS can't run 3.12 yet... We want 3.10 and 3.11 + python-version: "3.12" # MacOS can't run 3.12 yet... We want 3.10 and 3.11 environment: name: ghostly-build defaults: @@ -32,30 +30,44 @@ jobs: SIRE_DONT_PHONEHOME: 1 SIRE_SILENT_PHONEHOME: 1 steps: - - uses: conda-incubator/setup-miniconda@v3 + # + - uses: actions/checkout@v4 with: - auto-update-conda: true - python-version: ${{ matrix.python-version }} - activate-environment: ghostly_build - miniforge-version: latest -# - - name: Clone the devel branch - run: git clone -b devel https://github.com/openbiosim/ghostly -# - - name: Setup Conda - run: conda install -y -c conda-forge boa anaconda-client packaging -# - - name: Update Conda recipe - run: python ${{ github.workspace }}/ghostly/actions/update_recipe.py -# - - name: Prepare build location - run: mkdir ${{ github.workspace }}/build -# - - name: Build Conda package using conda build - run: conda build -c conda-forge -c openbiosim/label/dev ${{ github.workspace }}/ghostly/recipes/ghostly -# - - name: Upload Conda package - run: python ${{ github.workspace }}/ghostly/actions/upload_package.py + fetch-depth: 0 + # + - name: Compute version info + shell: bash + run: python actions/update_recipe.py + # + - name: Create sdist + run: pip install build && python -m build --sdist && mv dist/*.tar.gz ghostly-source.tar.gz + working-directory: ${{ github.workspace }} + # + - name: Install pixi + uses: prefix-dev/setup-pixi@v0.9.4 + with: + run-install: false + # + - name: Install rattler-build + shell: bash + run: pixi global install rattler-build + # + - name: Write Python variant config + shell: bash + run: printf 'python:\n - "${{ matrix.python-version }}"\n' > "${{ github.workspace }}/python_variant.yaml" + # + - name: Build package using rattler-build + shell: bash + run: rattler-build build --recipe "${{ github.workspace }}/recipes/ghostly" -c conda-forge -c openbiosim/label/dev --variant-config "${{ github.workspace }}/python_variant.yaml" + # + - name: Install anaconda-client + shell: bash + run: python -m pip install anaconda-client + if: ${{ matrix.platform.name == 'linux' && matrix.python-version == '3.11' }} + # + - name: Upload package + shell: bash + run: python actions/upload_package.py env: ANACONDA_TOKEN: ${{ secrets.ANACONDA_TOKEN }} ANACONDA_LABEL: dev diff --git a/.github/workflows/main.yaml b/.github/workflows/main.yaml index 7ef9819..66e71d0 100644 --- a/.github/workflows/main.yaml +++ b/.github/workflows/main.yaml @@ -23,7 +23,7 @@ jobs: exclude: - platform: { name: "macos", os: "macos-latest", shell: "bash -l {0}" } - python-version: "3.12" # MacOS can't run 3.12 yet... + python-version: "3.12" # MacOS can't run 3.12 yet... environment: name: ghostly-build defaults: @@ -33,30 +33,45 @@ jobs: SIRE_DONT_PHONEHOME: 1 SIRE_SILENT_PHONEHOME: 1 steps: - - uses: conda-incubator/setup-miniconda@v3 + # + - uses: actions/checkout@v4 with: - auto-update-conda: true - python-version: ${{ matrix.python-version }} - activate-environment: ghostly_build - miniforge-version: latest -# - - name: Clone the main branch - run: git clone -b main https://github.com/openbiosim/ghostly -# - - name: Setup Conda - run: conda install -y -c conda-forge boa anaconda-client packaging -# - - name: Update Conda recipe - run: python ${{ github.workspace }}/ghostly/actions/update_recipe.py -# - - name: Prepare build location - run: mkdir ${{ github.workspace }}/build -# - - name: Build Conda package using conda build - run: conda build -c conda-forge -c openbiosim/label/main ${{ github.workspace }}/ghostly/recipes/ghostly -# - - name: Upload Conda package - run: python ${{ github.workspace }}/ghostly/actions/upload_package.py + ref: main + fetch-depth: 0 + # + - name: Compute version info + shell: bash + run: python actions/update_recipe.py + # + - name: Create sdist + run: pip install build && python -m build --sdist && mv dist/*.tar.gz ghostly-source.tar.gz + working-directory: ${{ github.workspace }} + # + - name: Install pixi + uses: prefix-dev/setup-pixi@v0.9.4 + with: + run-install: false + # + - name: Install rattler-build + shell: bash + run: pixi global install rattler-build + # + - name: Write Python variant config + shell: bash + run: printf 'python:\n - "${{ matrix.python-version }}"\n' > "${{ github.workspace }}/python_variant.yaml" + # + - name: Build package using rattler-build + shell: bash + run: rattler-build build --recipe "${{ github.workspace }}/recipes/ghostly" -c conda-forge -c openbiosim/label/main --variant-config "${{ github.workspace }}/python_variant.yaml" + # + - name: Install anaconda-client + shell: bash + run: python -m pip install anaconda-client + if: github.event.inputs.upload_packages == 'true' && matrix.platform.name == 'linux' && matrix.python-version == '3.11' + # + - name: Upload package + shell: bash + run: python actions/upload_package.py env: ANACONDA_TOKEN: ${{ secrets.ANACONDA_TOKEN }} ANACONDA_LABEL: main diff --git a/.github/workflows/pr.yaml b/.github/workflows/pr.yaml index 6b56abc..ad6c7e0 100644 --- a/.github/workflows/pr.yaml +++ b/.github/workflows/pr.yaml @@ -17,14 +17,12 @@ jobs: - { name: "linux", os: "ubuntu-latest", shell: "bash -l {0}" } - { name: "macos", os: "macos-latest", shell: "bash -l {0}" } exclude: - # Exclude all but the latest Python from all - # but Linux - platform: { name: "macos", os: "macos-latest", shell: "bash -l {0}" } python-version: "3.10" - platform: { name: "macos", os: "macos-latest", shell: "bash -l {0}" } - python-version: "3.12" # MacOS can't run 3.12 yet... + python-version: "3.12" # MacOS can't run 3.12 yet... environment: name: ghostly-build defaults: @@ -35,29 +33,38 @@ jobs: SIRE_SILENT_PHONEHOME: 1 REPO: "${{ github.event.pull_request.head.repo.full_name || github.repository }}" steps: - - uses: conda-incubator/setup-miniconda@v3 + # + - uses: actions/checkout@v4 with: - auto-update-conda: true - python-version: ${{ matrix.python-version }} - activate-environment: ghostly_build - miniforge-version: latest -# - - name: Clone the feature branch - run: git clone -b ${{ github.head_ref }} --single-branch https://github.com/${{ env.REPO }} ghostly -# - - name: Setup Conda - run: conda install -y -c conda-forge boa anaconda-client packaging -# - - name: Update Conda recipe - run: python ${{ github.workspace }}/ghostly/actions/update_recipe.py -# - - name: Prepare build location - run: mkdir ${{ github.workspace }}/build -# - - name: Build Conda package using conda build using main channel + fetch-depth: 0 + # + - name: Compute version info + shell: bash + run: python actions/update_recipe.py + # + - name: Create sdist + run: pip install build && python -m build --sdist && mv dist/*.tar.gz ghostly-source.tar.gz + working-directory: ${{ github.workspace }} + # + - name: Install pixi + uses: prefix-dev/setup-pixi@v0.9.4 + with: + run-install: false + # + - name: Install rattler-build + shell: bash + run: pixi global install rattler-build + # + - name: Write Python variant config + shell: bash + run: printf 'python:\n - "${{ matrix.python-version }}"\n' > "${{ github.workspace }}/python_variant.yaml" + # + - name: Build package using rattler-build (main channel) if: ${{ github.base_ref == 'main' }} - run: conda build -c conda-forge -c openbiosim/label/main ${{ github.workspace }}/ghostly/recipes/ghostly -# - - name: Build Conda package using conda build using dev channel + shell: bash + run: rattler-build build --recipe "${{ github.workspace }}/recipes/ghostly" -c conda-forge -c openbiosim/label/main --variant-config "${{ github.workspace }}/python_variant.yaml" + # + - name: Build package using rattler-build (dev channel) if: ${{ github.base_ref != 'main' }} - run: conda build -c conda-forge -c openbiosim/label/dev ${{ github.workspace }}/ghostly/recipes/ghostly + shell: bash + run: rattler-build build --recipe "${{ github.workspace }}/recipes/ghostly" -c conda-forge -c openbiosim/label/dev --variant-config "${{ github.workspace }}/python_variant.yaml" diff --git a/.gitignore b/.gitignore index e232e0e..5789c59 100644 --- a/.gitignore +++ b/.gitignore @@ -20,7 +20,6 @@ setup.err dist/ build/ ghostly.egg-info -src/ghostly/_version.py # Test output. output.yaml @@ -35,8 +34,8 @@ output.yaml # VSCode config .vscode/ -# Conda recipe (it is auto-generated) -recipes/ghostly/meta.yaml - # Sire cache files. cache/ + +# Auto-generated version file. +src/ghostly/_version.py diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..0043ca4 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,23 @@ +files: ^(src|tests)/ +exclude: ^tests/(input|output)/ + +repos: + # General file quality checks + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v5.0.0 + hooks: + - id: trailing-whitespace + - id: end-of-file-fixer + - id: check-added-large-files + args: [--maxkb=1000] # Prevent files larger than 1MB + - id: check-merge-conflict + + # Python formatting and linting + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.8.4 + hooks: + # Run the formatter + - id: ruff-format + # Run the linter (optional - remove if too strict) + - id: ruff + args: [--fix, --exit-zero] # Auto-fix but don't block commits diff --git a/MANIFEST.in b/MANIFEST.in new file mode 100644 index 0000000..96737d3 --- /dev/null +++ b/MANIFEST.in @@ -0,0 +1 @@ +graft tests diff --git a/README.md b/README.md index 3598949..e769f4e 100644 --- a/README.md +++ b/README.md @@ -30,41 +30,60 @@ See the [examples](examples) directory for more details. ## Installation -First create a conda environment using the provided environment file: +### Conda package + +Install `ghostly` directly from the `openbiosim` channel: ``` -conda env create -f environment.yaml +conda install -c conda-forge -c openbiosim ghostly ``` -(We recommend using [Miniforge](https://github.com/conda-forge/miniforge).) - -Now install `ghostly` into the environment: +Or, for the development version: ``` -conda activate ghostly -pip install . +conda install -c conda-forge -c openbiosim/label/dev ghostly ``` -Or, for an editable install (useful for development): +### Installing from source (standalone) + +To install from source using [pixi](https://pixi.sh), which will +automatically create an environment with all required dependencies +(including pre-built [Sire](https://github.com/OpenBioSim/sire) and +[BioSimSpace](https://github.com/OpenBioSim/biosimspace)): ``` -conda activate ghostly +git clone https://github.com/openbiosim/ghostly +cd ghostly +pixi install +pixi shell pip install -e . ``` -For an existing conda environment, you can also install `ghostly` directly from -the `openbiosim` channel: +### Installing from source (full OpenBioSim development) + +If you are developing across the full OpenBioSim stack, first install +[Sire](https://github.com/OpenBioSim/sire) from source by following the +instructions [here](https://github.com/OpenBioSim/sire#installation), then +activate its pixi environment: ``` -conda install -c conda-forge -c openbiosim ghostly +pixi shell --manifest-path /path/to/sire/pixi.toml -e dev ``` -Or, for the development version: +You may also need to install other packages from source, e.g. +[BioSimSpace](https://github.com/OpenBioSim/biosimspace): ``` -conda install -c conda-forge -c openbiosim/label/dev ghostly +pip install -e /path/to/biosimspace +``` + +Then install `ghostly` into the environment: + +``` +pip install -e . ``` +### Testing You should now have a `ghostly` executable in your path. To test, run: @@ -72,6 +91,24 @@ You should now have a `ghostly` executable in your path. To test, run: ghostly --help ``` +## Development + +Pre-commit hooks are used to ensure consistent code formatting and linting. +To set up pre-commit in your development environment: + +``` +pixi shell -e dev +pre-commit install +``` + +This will run [ruff](https://docs.astral.sh/ruff/) formatting and linting +checks automatically on each commit. To run the checks manually against all +files: + +``` +pre-commit run --all-files +``` + ## Usage Ghostly requires topology and coordinate files for the reference and perturbed molecules diff --git a/actions/update_recipe.py b/actions/update_recipe.py index 16276ed..cd9a0fc 100644 --- a/actions/update_recipe.py +++ b/actions/update_recipe.py @@ -1,50 +1,58 @@ -import sys +"""Compute git version info for rattler-build. + +This script computes GIT_DESCRIBE_TAG and GIT_DESCRIBE_NUMBER from the +git history and outputs them in GitHub Actions format for setting +environment variables. + +It also writes a _version.py file so that versioningit has a fallback +when .git is not available (e.g., when rattler-build excludes it). +""" + import os import subprocess +import sys -# Get the name of the script. script = os.path.abspath(sys.argv[0]) - -# we want to import the 'get_requirements' package from this directory -sys.path.insert(0, os.path.dirname(script)) - -# go up one directories to get the source directory -# (this script is in BioSimSpace/actions/) srcdir = os.path.dirname(os.path.dirname(script)) - -condadir = os.path.join(srcdir, "recipes", "ghostly") - -print(f"conda recipe in {condadir}") - -# Store the name of the recipe and template YAML files. -recipe = os.path.join(condadir, "meta.yaml") -template = os.path.join(condadir, "template.yaml") - gitdir = os.path.join(srcdir, ".git") def run_cmd(cmd): - p = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE) - return str(p.stdout.read().decode("utf-8")).lstrip().rstrip() - - -# Get the remote. -remote = run_cmd( - f"git --git-dir={gitdir} --work-tree={srcdir} config --get remote.origin.url" -) -print(remote) - -# Get the branch. -branch = run_cmd( - f"git --git-dir={gitdir} --work-tree={srcdir} rev-parse --abbrev-ref HEAD" -) -print(branch) - -lines = open(template, "r").readlines() - -with open(recipe, "w") as FILE: - for line in lines: - line = line.replace("GHOSTLY_REMOTE", remote) - line = line.replace("GHOSTLY_BRANCH", branch) - - FILE.write(line) + p = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE) + stdout, _ = p.communicate() + return stdout.decode("utf-8").strip() + + +# Get the full git describe output (e.g., "2024.1.0-5-gabcdef" or "2024.1.0") +describe = run_cmd(f"git --git-dir={gitdir} --work-tree={srcdir} describe --tags") + +if "-" in describe: + # Format: tag-number-hash (e.g., "2024.1.0-5-gabcdef") + parts = describe.rsplit("-", 2) + tag = parts[0] + number = parts[1] + rev = parts[2] # e.g., "gabcdef" + version = f"{tag}+{number}.{rev}" +else: + # Exactly on a tag + tag = describe + number = "0" + version = tag + +print(f"GIT_DESCRIBE_TAG={tag}") +print(f"GIT_DESCRIBE_NUMBER={number}") +print(f"Version={version}") + +# Write to GITHUB_ENV if running in GitHub Actions +github_env = os.environ.get("GITHUB_ENV") +if github_env: + with open(github_env, "a") as f: + f.write(f"GIT_DESCRIBE_TAG={tag}\n") + f.write(f"GIT_DESCRIBE_NUMBER={number}\n") + print("Exported to GITHUB_ENV") + +# Write _version.py for versioningit fallback +version_file = os.path.join(srcdir, "src", "ghostly", "_version.py") +with open(version_file, "w") as f: + f.write(f'__version__ = "{version}"\n') +print(f"Wrote {version_file}") diff --git a/actions/upload_package.py b/actions/upload_package.py index 504bea7..fc3435f 100644 --- a/actions/upload_package.py +++ b/actions/upload_package.py @@ -1,16 +1,18 @@ +"""Upload built packages to the openbiosim Anaconda Cloud channel.""" + import os import sys import glob +import subprocess script = os.path.abspath(sys.argv[0]) -# go up one directories to get the source directory -# (this script is in ghostly/actions/) +# Go up one directory to get the source directory. srcdir = os.path.dirname(os.path.dirname(script)) print(f"Ghostly source is in {srcdir}\n") -# Get the anaconda token to authorise uploads +# Get the anaconda token to authorise uploads. if "ANACONDA_TOKEN" in os.environ: conda_token = os.environ["ANACONDA_TOKEN"] else: @@ -22,42 +24,30 @@ else: conda_label = "dev" -# get the root conda directory -conda = os.environ["CONDA"] - -# Set the path to the conda-bld directory. -conda_bld = os.path.join(conda, "envs", "ghostly_build", "conda-bld") - -print(f"conda_bld = {conda_bld}") +# Search for rattler-build output first. +packages = glob.glob(os.path.join("output", "**", "*.conda"), recursive=True) -# Find the packages to upload -ghostly_pkg = glob.glob(os.path.join(conda_bld, "noarch", "ghostly-*.tar.bz2")) +# Fall back to conda-bld output. +if not packages: + if "CONDA" in os.environ: + conda = os.environ["CONDA"] + conda_bld = os.path.join(conda, "envs", "ghostly_build", "conda-bld") + packages = glob.glob( + os.path.join(conda_bld, "**", "ghostly-*.tar.bz2"), recursive=True + ) -if len(ghostly_pkg) == 0: +if not packages: print("No ghostly packages to upload?") sys.exit(-1) -packages = ghostly_pkg - -print(f"Uploading packages:") -print(" * ", "\n * ".join(packages)) - -packages = " ".join(packages) - - -def run_cmd(cmd): - import subprocess - - p = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE) - return str(p.stdout.read().decode("utf-8")).lstrip().rstrip() +print("Uploading packages:") +for pkg in packages: + print(f" * {pkg}") - -gitdir = os.path.join(srcdir, ".git") - -tag = run_cmd(f"git --git-dir={gitdir} --work-tree={srcdir} tag --contains") +packages_str = " ".join(packages) # Upload the packages to the openbiosim channel on Anaconda Cloud. -cmd = f"anaconda --token {conda_token} upload --user openbiosim --label {conda_label} --force {packages}" +cmd = f"anaconda --token {conda_token} upload --user openbiosim --label {conda_label} --force {packages_str}" print(f"\nUpload command:\n\n{cmd}\n") @@ -65,8 +55,12 @@ def run_cmd(cmd): print("Not uploading as the ANACONDA_TOKEN is not set!") sys.exit(-1) -output = run_cmd(cmd) -print(output) +def run_cmd(cmd): + p = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE) + return str(p.stdout.read().decode("utf-8")).lstrip().rstrip() + +output = run_cmd(cmd) +print(output) print("Package uploaded!") diff --git a/environment.yaml b/environment.yaml deleted file mode 100644 index 3fee8fc..0000000 --- a/environment.yaml +++ /dev/null @@ -1,9 +0,0 @@ -name: ghostly - -channels: - - conda-forge - - openbiosim/label/dev - -dependencies: - - biosimspace - - loguru diff --git a/pixi.toml b/pixi.toml new file mode 100644 index 0000000..1d59545 --- /dev/null +++ b/pixi.toml @@ -0,0 +1,34 @@ +[workspace] +name = "ghostly" +channels = ["conda-forge", "openbiosim/label/dev"] +platforms = ["linux-64", "linux-aarch64", "osx-arm64", "win-64"] + +[dependencies] +python = ">=3.10" +loguru = "*" + +[target.linux-64.dependencies] +biosimspace = "*" + +[target.linux-aarch64.dependencies] +# biosimspace/sire not available as conda packages on linux-aarch64; +# build from source first + +[target.osx-arm64.dependencies] +biosimspace = "*" + +[target.win-64.dependencies] +biosimspace = "*" + +[feature.test.dependencies] +pytest = "*" + +[feature.lint.dependencies] +pre-commit = "*" +rattler-build = "*" +ruff = "*" + +[environments] +default = [] +test = ["test"] +dev = ["test", "lint"] diff --git a/pyproject.toml b/pyproject.toml index de0298d..908968a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -32,3 +32,9 @@ distance-dirty = "{base_version}+{distance}.{vcs}{rev}.dirty" [tool.versioningit.write] file = "src/ghostly/_version.py" + +[tool.ruff.lint] +ignore = ["E402"] + +[tool.ruff.lint.per-file-ignores] +"tests/**" = ["F841"] diff --git a/recipes/ghostly/conda_build_config.yaml b/recipes/ghostly/conda_build_config.yaml deleted file mode 100644 index 3e8e203..0000000 --- a/recipes/ghostly/conda_build_config.yaml +++ /dev/null @@ -1,3 +0,0 @@ -pin_run_as_build: - sire: - max_pin: x.x diff --git a/recipes/ghostly/recipe.yaml b/recipes/ghostly/recipe.yaml new file mode 100644 index 0000000..dacf549 --- /dev/null +++ b/recipes/ghostly/recipe.yaml @@ -0,0 +1,50 @@ +context: + name: ghostly + +package: + name: ${{ name }} + version: ${{ env.get('GIT_DESCRIBE_TAG', default='PR') }} + +source: + path: ../../ghostly-source.tar.gz + +build: + number: ${{ env.get('GIT_DESCRIBE_NUMBER', default='0') }} + noarch: python + script: python -m pip install . --no-deps --ignore-installed -vv + +requirements: + host: + - pip + - python + - setuptools + - versioningit + run: + - biosimspace + - loguru + - python + +tests: + - python: + imports: + - ghostly + - script: + - pytest -vvv --color=yes --import-mode=importlib ./tests + files: + source: + - tests/ + requirements: + run: + - pytest + +about: + homepage: https://github.com/openbiosim/ghostly + license: GPL-3.0-or-later + license_file: LICENSE + summary: "Ghost atom bonded term modifications for alchemical free-energy simulations." + repository: https://github.com/openbiosim/ghostly + documentation: https://github.com/openbiosim/ghostly + +extra: + recipe-maintainers: + - lohedges diff --git a/recipes/ghostly/template.yaml b/recipes/ghostly/template.yaml deleted file mode 100644 index 475b629..0000000 --- a/recipes/ghostly/template.yaml +++ /dev/null @@ -1,83 +0,0 @@ -{% set name = "ghostly" %} - -package: - name: {{ name }} - version: {{ environ.get('GIT_DESCRIBE_TAG', 'PR') }} - -source: - git_url: GHOSTLY_REMOTE - git_tag: GHOSTLY_BRANCH - -build: - number: {{ environ.get('GIT_DESCRIBE_NUMBER', 0) }} - noarch: python - script: {{ PYTHON }} -m pip install . --no-deps --ignore-installed -vv - -requirements: - host: - - biosimspace - - loguru - - pip - - python - - setuptools - - sire - - versioningit - run: - - biosimspace - - loguru - - python - - sire - -test: - script_env: - - SIRE_DONT_PHONEHOME - - SIRE_SILENT_PHONEHOME - requires: - - black == 25 # [linux and x86_64 and py==311] - - pytest - - pytest-black # [linux and x86_64 and py==311] - imports: - - ghostly - source_files: - - src/ghostly - - tests - commands: - - pytest -vvv --color=yes --black src/ghostly # [linux and x86_64 and py==311] - - pytest -vvv --color=yes --import-mode=importlib tests - -about: - home: https://github.com/openbiosim/ghostly - license: GPL-3.0-or-later - license_file: '{{ environ["RECIPE_DIR"] }}/LICENSE' - summary: "Ghost atom bonded term modifications for alchemical free-energy simulations." - dev_url: https://github.com/openbiosim/ghostly - doc_url: https://github.com/openbiosim/ghostly - description: | - Ghostly is a package to perform modification of ghost (dummy) atom bonded - terms for alchemical free-energy calculations, using the approach described - in the paper "Dummy Atoms in Alchemical Free Energy Calculations" by - Fleck et al., JCTC, 2021, 17, 7, 4403-4419. These modifications were - designed to solve two key issues: - - * To ensure that ghost atoms only give a multiplicative contribution to the - partition function, which will cancel when computing double free-energy - differences. - * To avoid spurious coupling between the physical and ghost systems, which - can affect the equilibrium geometry of the physical system. - - To install: - - `conda install -c conda-forge -c openbiosim ghostly` - - To install the development version: - - `conda install -c conda-forge -c openbiosim/label/dev ghostly` - - When updating the development version it is generally advised to - update Sire at the same time: - - `conda install -c conda-forge -c openbiosim/label/dev ghostly sire` - -extra: - recipe-maintainers: - - lohedges diff --git a/src/ghostly/_cli.py b/src/ghostly/_cli.py index 39ebf38..59794fc 100644 --- a/src/ghostly/_cli.py +++ b/src/ghostly/_cli.py @@ -36,7 +36,6 @@ def run(): import argparse import json - import os import sys import BioSimSpace as BSS @@ -100,6 +99,19 @@ def run(): required=False, ) + parser.add_argument( + "--k-hard-ring", + type=str, + help=""" + The force constant to use when setting angle terms involving ghost + atoms to 90 degrees for bridge atoms that are in a ring. A lower + value reduces strain where the 90 degree target fights the natural + ring geometry. + """, + default="75 kcal/mol/rad**2", + required=False, + ) + parser.add_argument( "--k-soft", type=str, @@ -133,6 +145,48 @@ def run(): required=False, ) + parser.add_argument( + "--soften-anchors", + type=float, + help=""" + Scale factor for surviving mixed ghost/physical dihedral force + constants. 1.0 keeps the original force constants (Boresch + approach). 0.0 removes all mixed dihedrals entirely (SOMD1 + scheme). Intermediate values (e.g. 0.5) scale the force + constants, reducing the constraint on ghost group orientation. + Softening can prevent dynamics crashes at small lambda for + complex perturbations, particularly with multiple ghost groups. + """, + default=1.0, + required=False, + ) + + parser.add_argument( + "--stiffen-rotamers", + action=argparse.BooleanOptionalAction, + help=""" + Whether to replace rotamer anchor dihedrals with a stiff + single-well cosine potential. When a bridge-physical bond is + rotatable (not in a ring, sp3 bridge), surviving anchor + dihedrals can allow rotameric transitions of ghost atoms at + intermediate lambda. + """, + default=False, + required=False, + ) + + parser.add_argument( + "--k-rotamer", + type=str, + help=""" + The force constant for the replacement cosine well when + stiffening rotamer anchor dihedrals. The resulting barrier + height is 2 * k_rotamer. + """, + default="50 kcal/mol", + required=False, + ) + parser.add_argument( "--output-prefix", type=str, @@ -238,6 +292,12 @@ def run(): logger.error(f"An error occurred while parsing the k-hard value: {e}") sys.exit(1) + try: + k_hard_ring = sr.u(args.k_hard_ring) + except Exception as e: + logger.error(f"An error occurred while parsing the k-hard-ring value: {e}") + sys.exit(1) + try: k_soft = sr.u(args.k_soft) except Exception as e: @@ -248,10 +308,30 @@ def run(): if not k_hard.has_same_units(u): logger.error("k-hard must have units of kcal/mol/rad**2") sys.exit(1) + if not k_hard_ring.has_same_units(u): + logger.error("k-hard-ring must have units of kcal/mol/rad**2") + sys.exit(1) if not k_soft.has_same_units(u): logger.error("k-soft must have units of kcal/mol/rad**2") sys.exit(1) + # Validate soften-anchors. + if not 0.0 <= args.soften_anchors <= 1.0: + logger.error("soften-anchors must be between 0.0 and 1.0") + sys.exit(1) + + # Try to parse the k-rotamer value. + try: + k_rotamer = sr.u(args.k_rotamer) + except Exception as e: + logger.error(f"An error occurred while parsing the k-rotamer value: {e}") + sys.exit(1) + + u_energy = sr.u("kcal/mol") + if not k_rotamer.has_same_units(u_energy): + logger.error("k-rotamer must have units of kcal/mol") + sys.exit(1) + # Try to merge the reference and perturbed molecules. if args.system is None: try: @@ -272,9 +352,13 @@ def run(): system, modifications = modify( system, k_hard.value(), - k_soft.value(), - args.optimise_angles, - args.num_optimise, + k_hard_ring=k_hard_ring.value(), + k_soft=k_soft.value(), + optimise_angles=args.optimise_angles, + num_optimise=args.num_optimise, + soften_anchors=args.soften_anchors, + stiffen_rotamers=args.stiffen_rotamers, + k_rotamer=k_rotamer.value(), ) except Exception as e: logger.error( diff --git a/src/ghostly/_ghostly.py b/src/ghostly/_ghostly.py index e64641f..1d6d1c7 100644 --- a/src/ghostly/_ghostly.py +++ b/src/ghostly/_ghostly.py @@ -31,7 +31,7 @@ try: from somd2 import _logger -except: +except Exception: from loguru import logger as _logger import platform as _platform @@ -44,7 +44,17 @@ del _platform -def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): +def modify( + system, + k_hard=100, + k_hard_ring=75, + k_soft=5, + optimise_angles=True, + num_optimise=10, + soften_anchors=1.0, + stiffen_rotamers=False, + k_rotamer=50, +): """ Apply modifications to ghost atom bonded terms to avoid non-physical coupling between the ghost atoms and the physical region. @@ -59,6 +69,13 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): The force constant to use to when setting angle terms involving ghost atoms to 90 degrees to avoid flapping. (In kcal/mol/rad^2) + k_hard_ring : float, optional + The force constant to use when setting angle terms involving ghost + atoms to 90 degrees for bridge atoms that are in a ring. Ring bridges + suffer strain from the 90 degree target fighting the natural ring + geometry (~110 degrees). A lower force constant reduces this strain + while still preventing flapping. (In kcal/mol/rad^2) + k_soft : float, optional The force constant to use when setting angle terms involving ghost atoms for non-planar triple junctions. (In kcal/mol/rad^2) @@ -72,6 +89,32 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): value of the angle terms involving ghost atoms for non-planar triple junctions. + soften_anchors : float, optional + Scale factor for surviving mixed ghost/physical dihedral force + constants. Applied as a post-processing step after all per-bridge + junction handlers and residual removal passes. A value of 1.0 + (default) keeps the original force constants (Boresch approach). + A value of 0.0 removes all mixed dihedrals entirely (SOMD1 scheme). + Intermediate values (e.g. 0.5) scale the force constants, reducing + the constraint on ghost group orientation while preserving the + thermodynamic correction. Softening can prevent dynamics crashes + at small lambda for complex perturbations where ghost groups are + constrained too tightly, particularly when multiple ghost groups + share hub atoms. + + stiffen_rotamers : bool, optional + Whether to replace rotamer anchor dihedrals with a stiff single-well + cosine potential. When a bridge--physical bond is rotatable (not in a + ring, sp3 bridge), surviving anchor dihedrals can allow rotameric + transitions of ghost atoms at intermediate lambda. Enabling this + replaces those dihedrals with a single cosine well of depth + ``2 * k_rotamer``. Default is False (only log warnings). + + k_rotamer : float, optional + The force constant for the replacement cosine well when stiffening + rotamer anchor dihedrals (in kcal/mol). The resulting barrier height + is ``2 * k_rotamer``. Only used when ``stiffen_rotamers`` is True. + Returns ------- @@ -114,20 +157,21 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): modifications["lambda_0"] = { "removed_angles": [], "removed_dihedrals": [], - "stiffened_angles": [], + "stiffened_angles": {}, "softened_angles": {}, + "softened_dihedrals": [], + "stiffened_dihedrals": [], } modifications["lambda_1"] = { "removed_angles": [], "removed_dihedrals": [], - "stiffened_angles": [], + "stiffened_angles": {}, "softened_angles": {}, + "softened_dihedrals": [], + "stiffened_dihedrals": [], } for mol in pert_mols: - # Store the molecule info. - info = mol.info() - # Generate the end state connectivity objects. connectivity0 = _create_connectivity(_morph.link_to_reference(mol)) connectivity1 = _create_connectivity(_morph.link_to_perturbed(mol)) @@ -163,12 +207,16 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): else: bridges0[c].append(ghost) # Work out the indices of the other physical atoms that are connected to - # the bridge atoms, sorted by the atom index. + # the bridge atoms, sorted by the atom index. These are "core" physical + # atoms, i.e. they are physical in both end states. physical0 = {} for b in bridges0: physical0[b] = [] for c in connectivity0.connections_to(b): - if not _is_ghost(mol, [c])[0]: + if ( + not _is_ghost(mol, [c])[0] + and not _is_ghost(mol, [c], is_lambda1=True)[0] + ): physical0[b].append(c) for b in physical0: physical0[b].sort(key=lambda x: x.value()) @@ -186,7 +234,10 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): for b in bridges1: physical1[b] = [] for c in connectivity1.connections_to(b): - if not _is_ghost(mol, [c], is_lambda1=True)[0]: + if ( + not _is_ghost(mol, [c])[0] + and not _is_ghost(mol, [c], is_lambda1=True)[0] + ): physical1[b].append(c) for b in physical1: physical1[b].sort(key=lambda x: x.value()) @@ -219,6 +270,10 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): # Now process the bridges. + # Compute the set of bridge atom indices for anchor selection. + bridge_indices0 = set(b.value() for b in bridges0) + bridge_indices1 = set(b.value() for b in bridges1) + # First lambda = 0. for b in bridges0: # Determine the type of junction. @@ -227,7 +282,13 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): # Terminal junction. if junction == 1: mol = _terminal( - mol, b, bridges0[b], physical0[b], connectivity0, modifications + mol, + b, + bridges0[b], + physical0[b], + connectivity0, + modifications, + bridge_indices=bridge_indices0, ) # Dual junction. @@ -240,6 +301,8 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): connectivity0, modifications, k_hard=k_hard, + k_hard_ring=k_hard_ring, + bridge_indices=bridge_indices0, ) # Triple junction. @@ -249,11 +312,14 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): b, bridges0[b], physical0[b], + connectivity0, modifications, k_hard=k_hard, + k_hard_ring=k_hard_ring, k_soft=k_soft, optimise_angles=optimise_angles, num_optimise=num_optimise, + bridge_indices=bridge_indices0, ) # Higher order junction. @@ -271,6 +337,40 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): num_optimise=num_optimise, ) + # Remove any improper dihedrals connecting ghosts to the physical region. + mol = _remove_impropers(mol, ghosts0, modifications, is_lambda1=False) + + # Remove any residual ghost dihedrals not caught by the per-bridge + # junction handlers (cross-bridge and ghost-middle patterns). + mol = _remove_residual_ghost_dihedrals( + mol, ghosts0, modifications, is_lambda1=False + ) + + # Remove any angles where the central atom is ghost and both terminal + # atoms are physical (e.g. B1-G-B2 in ring-breaking topologies). + mol = _remove_ghost_centre_angles(mol, ghosts0, modifications, is_lambda1=False) + + # Soften any surviving mixed ghost/physical dihedrals. + mol = _soften_mixed_dihedrals( + mol, + ghosts0, + modifications, + soften_anchors=soften_anchors, + is_lambda1=False, + ) + + # Check for potential rotamer anchor dihedrals. + mol = _check_rotamer_anchors( + mol, + bridges0, + physical0, + ghosts0, + modifications, + is_lambda1=False, + stiffen=stiffen_rotamers, + k_rotamer=k_rotamer, + ) + # Now lambda = 1. for b in bridges1: junction = len(physical1[b]) @@ -284,6 +384,7 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): connectivity1, modifications, is_lambda1=True, + bridge_indices=bridge_indices1, ) elif junction == 2: @@ -295,7 +396,9 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): connectivity1, modifications, k_hard=k_hard, + k_hard_ring=k_hard_ring, is_lambda1=True, + bridge_indices=bridge_indices1, ) elif junction == 3: @@ -304,12 +407,15 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): b, bridges1[b], physical1[b], + connectivity1, modifications, k_hard=k_hard, + k_hard_ring=k_hard_ring, k_soft=k_soft, optimise_angles=optimise_angles, num_optimise=num_optimise, is_lambda1=True, + bridge_indices=bridge_indices1, ) # Higher order junction. @@ -328,6 +434,40 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): is_lambda1=True, ) + # Remove any improper dihedrals connecting ghosts to the physical region. + mol = _remove_impropers(mol, ghosts1, modifications, is_lambda1=True) + + # Remove any residual ghost dihedrals not caught by the per-bridge + # junction handlers (cross-bridge and ghost-middle patterns). + mol = _remove_residual_ghost_dihedrals( + mol, ghosts1, modifications, is_lambda1=True + ) + + # Remove any angles where the central atom is ghost and both terminal + # atoms are physical (e.g. B1-G-B2 in ring-breaking topologies). + mol = _remove_ghost_centre_angles(mol, ghosts1, modifications, is_lambda1=True) + + # Soften any surviving mixed ghost/physical dihedrals. + mol = _soften_mixed_dihedrals( + mol, + ghosts1, + modifications, + soften_anchors=soften_anchors, + is_lambda1=True, + ) + + # Check for potential rotamer anchor dihedrals. + mol = _check_rotamer_anchors( + mol, + bridges1, + physical1, + ghosts1, + modifications, + is_lambda1=True, + stiffen=stiffen_rotamers, + k_rotamer=k_rotamer, + ) + # Update the molecule in the system. system.update(mol) @@ -335,8 +475,113 @@ def modify(system, k_hard=100, k_soft=5, optimise_angles=True, num_optimise=10): return system, modifications +def _score_atom(mol, atom, bridge_indices): + """ + Score an atom based on its bridge and transmuting status. + + Lower scores are better: + + 0 - not a bridge, not transmuting (ideal) + 1 - transmuting but not a bridge + 2 - bridge but not transmuting + 3 - both bridge and transmuting (worst) + + Parameters + ---------- + + mol : sire.mol.Molecule + The perturbable molecule. + + atom : sire.legacy.Mol.AtomIdx + The atom to score. + + bridge_indices : set of int + The atom index values of all bridge atoms at the current end state. + + Returns + ------- + + score : int + The score (0-3). + """ + + score = 0 + + # Penalise bridge atoms (couples ghost groups). + if atom.value() in bridge_indices: + score += 2 + + # Penalise transmuting atoms (unstable reference frame). + a = mol.atom(atom) + if a.property("element0").symbol() != a.property("element1").symbol(): + score += 1 + + return score + + +def _select_anchor(mol, candidates, bridge_indices): + """ + Select the best anchor atom from physical2 candidates for a terminal + junction. + + The anchor atom is retained to constrain the ghost group's rotation via + surviving anchor dihedrals. A poor choice (e.g. a transmuting bridge atom) + can couple independent ghost groups or tie ghost constraints to element + transmutations. + + Candidates are scored using ``_score_atom()`` (lower is better). + Ties are broken by atom index (lowest first), preserving the previous + ``physical2[0]`` behaviour when all candidates score equally. + + Parameters + ---------- + + mol : sire.mol.Molecule + The perturbable molecule. + + candidates : List[sire.legacy.Mol.AtomIdx] + The physical2 candidates, sorted by atom index. + + bridge_indices : set of int + The atom index values of all bridge atoms at the current end state. + + Returns + ------- + + anchor : sire.legacy.Mol.AtomIdx + The best candidate to use as anchor. + """ + + best = candidates[0] + best_score = 3 # worst possible + + for cand in candidates: + score = _score_atom(mol, cand, bridge_indices) + + if score < best_score: + best_score = score + best = cand + if score == 0: + break # Can't do better. + + if best != candidates[0]: + _logger.debug( + f" Anchor selection: chose atom {best.value()} over " + f"{candidates[0].value()} (avoiding bridge/transmuting atom)" + ) + + return best + + def _terminal( - mol, bridge, ghosts, physical, connectivity, modifications, is_lambda1=False + mol, + bridge, + ghosts, + physical, + connectivity, + modifications, + is_lambda1=False, + bridge_indices=None, ): r""" Apply modifications to a terminal junction. @@ -376,6 +621,12 @@ def _terminal( is_lambda1 : bool, optional Whether the junction is at lambda = 1. + bridge_indices : set of int, optional + The atom index values of all bridge atoms at the current end state. + Used by ``_select_anchor`` to avoid choosing bridge or transmuting + atoms as anchors. When ``None``, falls back to first-by-index + selection. + Returns ------- @@ -388,6 +639,10 @@ def _terminal( f"{_lam_sym} = {int(is_lambda1)}:" ) + # Nothing to do if there are no ghost atoms. + if not ghosts: + return mol + # Store the molecular info. info = mol.info() @@ -420,17 +675,27 @@ def _terminal( # Initialise a container to store the updated dihedrals. new_dihedrals = _SireMM.FourAtomFunctions(mol.info()) - # Remove all dihedral terms for all but one of the physical atoms two atoms - # from the physical bridge atom. - physical2.pop(0) + # Select the best anchor atom and remove it from the list. Dihedrals + # through the remaining physical2 atoms will be removed; dihedrals through + # the anchor are kept to constrain the ghost group's orientation. + if bridge_indices is not None and len(physical2) > 1: + anchor = _select_anchor(mol, physical2, bridge_indices) + else: + anchor = physical2[0] + physical2.remove(anchor) + _logger.debug(f" Anchor atom: {anchor.value()}") for p in dihedrals.potentials(): idx0 = info.atom_idx(p.atom0()) idx1 = info.atom_idx(p.atom1()) idx2 = info.atom_idx(p.atom2()) idx3 = info.atom_idx(p.atom3()) - if (idx0 in physical2 and idx3 in ghosts) or ( + + # Cross-bridge dihedral: non-anchor physical2 at one end, ghost at other. + is_cross_bridge = (idx0 in physical2 and idx3 in ghosts) or ( idx3 in physical2 and idx0 in ghosts - ): + ) + + if is_cross_bridge: _logger.debug( f" Removing dihedral: [{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], {p.function()}" ) @@ -455,7 +720,9 @@ def _dual( connectivity, modifications, k_hard=100, + k_hard_ring=75, is_lambda1=False, + bridge_indices=None, ): r""" Apply modifications to a dual junction. @@ -496,9 +763,17 @@ def _dual( The force constant to use when setting angle terms involving ghost atoms to 90 degrees to avoid flapping. (In kcal/mol/rad^2) + k_hard_ring : float, optional + The force constant to use when the bridge atom is in a ring, where + the 90 degree target fights the natural ring geometry. (In kcal/mol/rad^2) + is_lambda1 : bool, optional Whether the junction is at lambda = 1. + bridge_indices : set of int, optional + The atom index values of all bridge atoms at the current end state. + Used to score physical neighbours. + Returns ------- @@ -511,6 +786,10 @@ def _dual( f"{_lam_sym} = {int(is_lambda1)}:" ) + # Nothing to do if there are no ghost atoms. + if not ghosts: + return mol + # Store the molecular info. info = mol.info() @@ -522,6 +801,34 @@ def _dual( mod_key = "lambda_0" suffix = "0" + # Check if the bridge atom is in a ring using RDKit. + if is_lambda1: + end_state_mol = _morph.link_to_perturbed(mol) + else: + end_state_mol = _morph.link_to_reference(mol) + + from sire.convert import to_rdkit + + try: + rdmol = to_rdkit(end_state_mol) + except Exception: + rdmol = None + + bridge_in_ring = ( + rdmol is not None and rdmol.GetAtomWithIdx(bridge.value()).IsInRing() + ) + + # Score the physical neighbours. + if bridge_indices is not None: + phys_scores = {p: _score_atom(mol, p, bridge_indices) for p in physical} + best_phys_score = min(phys_scores.values()) + else: + phys_scores = None + best_phys_score = 0 + + # Choose the force constant for angle stiffening. + k = k_hard_ring if bridge_in_ring else k_hard + # Single branch. if len(ghosts) == 1: _logger.debug(" Single branch:") @@ -583,13 +890,27 @@ def _dual( or idx0 in physical and idx2 in ghosts ): + # Identify the physical atom in this angle. + phys_atom = idx2 if idx0 in ghosts else idx0 + + # Skip stiffening through poorly-scoring atoms if better ones exist. + if phys_scores is not None and phys_scores[phys_atom] > best_phys_score: + new_angles.set(idx0, idx1, idx2, p.function()) + _logger.debug( + f" Skipping stiffening for angle " + f"[{idx0.value()}-{idx1.value()}-{idx2.value()}]: " + f"physical atom {phys_atom.value()} scores " + f"{phys_scores[phys_atom]} (transmuting/bridge)" + ) + continue + from math import pi from sire.legacy.CAS import Symbol theta0 = pi / 2.0 # Create the new angle function. - amber_angle = _SireMM.AmberAngle(k_hard, theta0) + amber_angle = _SireMM.AmberAngle(k, theta0) # Generate the new angle expression. expression = amber_angle.to_expression(Symbol("theta")) @@ -602,8 +923,10 @@ def _dual( f"{p.function()} --> {expression}" ) - ang_idx = (idx0.value(), idx1.value(), idx2.value()) - modifications[mod_key]["stiffened_angles"].append(ang_idx) + ang_idx = ",".join( + [str(i) for i in (idx0.value(), idx1.value(), idx2.value())] + ) + modifications[mod_key]["stiffened_angles"][ang_idx] = {"k": k} else: new_angles.set(idx0, idx1, idx2, p.function()) @@ -621,7 +944,7 @@ def _dual( else: _logger.debug(" Dual branch:") - # First, delete all bonded terms between atoms in two ghost branches. + # First, delete all bonded terms between atoms in the two ghost branches. # Get the end state bond functions. angles = mol.property("angle" + suffix) @@ -679,16 +1002,18 @@ def _dual( ) # Now treat the ghost branches individually. - for d in ghosts: + for ghost in ghosts: mol = _dual( mol, bridge, - [d], + [ghost], physical, connectivity, modifications, k_hard=k_hard, + k_hard_ring=k_hard_ring, is_lambda1=is_lambda1, + bridge_indices=bridge_indices, ) # Return the updated molecule. @@ -700,12 +1025,15 @@ def _triple( bridge, ghosts, physical, + connectivity, modifications, k_hard=100, + k_hard_ring=75, k_soft=5, optimise_angles=True, num_optimise=10, is_lambda1=False, + bridge_indices=None, ): r""" Apply modifications to a triple junction. @@ -735,6 +1063,9 @@ def _triple( physical : List[sire.legacy.Mol.AtomIdx] The list of physical atoms connected to the bridge atom. + connectivity : sire.legacy.MM.Connectivity + The connectivity of the molecule at the relevant end state. + modifications : dict A dictionary to store details of the modifications made. @@ -742,6 +1073,10 @@ def _triple( The force constant to use when setting angle terms involving ghost atoms to 90 degrees to avoid flapping. (In kcal/mol/rad^2) + k_hard_ring : float, optional + The force constant to use when the bridge atom is in a ring, where + the 90 degree target fights the natural ring geometry. (In kcal/mol/rad^2) + k_soft : float, optional The force constant to use when setting angle terms involving ghost atoms for non-planar triple junctions. (In kcal/mol/rad^2) @@ -757,6 +1092,10 @@ def _triple( is_lambda1 : bool, optional Whether the junction is at lambda = 1. + bridge_indices : set of int, optional + The atom index values of all bridge atoms at the current end state. + Used to score physical neighbours. + Returns ------- @@ -769,6 +1108,10 @@ def _triple( f"{_lam_sym} = {int(is_lambda1)}:" ) + # Nothing to do if there are no ghost atoms. + if not ghosts: + return mol + # Store the molecular info. info = mol.info() @@ -793,13 +1136,31 @@ def _triple( _logger.error(msg) raise RuntimeError(msg) - # Get the hybridisation of the bridge atom. - hybridisation = rdmol.GetAtomWithIdx(bridge.value()).GetHybridization() + # Get the hybridisation and ring membership of the bridge atom. + bridge_rdatom = rdmol.GetAtomWithIdx(bridge.value()) + hybridisation = bridge_rdatom.GetHybridization() + bridge_in_ring = bridge_rdatom.IsInRing() + + # Warn if the hybridisation is unspecified. + if hybridisation in (HybridizationType.UNSPECIFIED, HybridizationType.OTHER): + _logger.warning( + f"Unspecified hybridisation for bridge atom {bridge.value()} " + f"at {_lam_sym} = {int(is_lambda1)}. Defaulting to planar junction." + ) + + # Score the physical neighbours. + if bridge_indices is not None: + phys_scores = {p: _score_atom(mol, p, bridge_indices) for p in physical} + best_phys_score = min(phys_scores.values()) + else: + phys_scores = None + best_phys_score = 0 # Non-planar junction. - if ( - hybridisation == HybridizationType.SP3 - or hybridisation == HybridizationType.SP3D + if hybridisation in ( + HybridizationType.SP3, + HybridizationType.SP3D, + HybridizationType.SP3D2, ): _logger.debug(" Non-planar junction.") @@ -826,6 +1187,20 @@ def _triple( or idx2 in ghosts and idx0 in physical ): + # Identify the physical atom in this angle. + phys_atom = idx2 if idx0 in ghosts else idx0 + + # Skip softening through poorly-scoring atoms if better ones exist. + if phys_scores is not None and phys_scores[phys_atom] > best_phys_score: + new_angles.set(idx0, idx1, idx2, p.function()) + _logger.debug( + f" Skipping softening for angle " + f"[{idx0.value()}-{idx1.value()}-{idx2.value()}]: " + f"physical atom {phys_atom.value()} scores " + f"{phys_scores[phys_atom]} (transmuting/bridge)" + ) + continue + from sire.legacy.CAS import Symbol theta = Symbol("theta") @@ -922,6 +1297,7 @@ def _triple( theta0s[idx] = [] # Perform multiple minimisations to get an average for the theta0 values. + is_error = False for _ in range(num_optimise): # Minimise the molecule. min_mol = _morph.link_to_reference(mol) @@ -930,7 +1306,10 @@ def _triple( constraint="none", platform="cpu", ) - minimiser.run() + try: + minimiser.run() + except Exception: + is_error = True # Commit the changes. min_mol = minimiser.commit() @@ -939,9 +1318,15 @@ def _triple( for idx in angle_idxs: try: theta0s[idx].append(min_mol.angles(*idx).sizes()[0].to(_radian)) - except: + except Exception: raise ValueError(f"Could not find optimised angle term: {idx}") + if is_error: + _logger.warning( + "Minimisation failed to converge during angle optimisation." + ) + break + # Compute the mean and standard error. import numpy as _np @@ -1008,8 +1393,17 @@ def _triple( # First remove all bonded terms between one of the physical atoms # and the ghost group. - # Store the index of the first physical atom. - idx = physical[0] + # Choose the worst-scoring physical atom as the sacrificial one. + if phys_scores is not None: + idx = max(physical, key=lambda p: (phys_scores[p], p.value())) + if idx != physical[0]: + _logger.debug( + f" Sacrificial atom selection: chose atom {idx.value()} " + f"over {physical[0].value()} (worst score " + f"{phys_scores[idx]})" + ) + else: + idx = physical[0] # Get the end state bond functions. angles = mol.property("angle" + suffix) @@ -1056,52 +1450,8 @@ def _triple( else: new_dihedrals.set(idx0, idx1, idx2, idx3, p.function()) - # Next we modify the angle terms between the remaining physical and - # ghost atoms so that the equilibrium angle is 90 degrees. - new_new_angles = _SireMM.ThreeAtomFunctions(mol.info()) - for p in new_angles.potentials(): - idx0 = info.atom_idx(p.atom0()) - idx1 = info.atom_idx(p.atom1()) - idx2 = info.atom_idx(p.atom2()) - - if ( - idx0 in ghosts - and idx2 in physical[1:] - or idx0 in physical[1:] - and idx2 in ghosts - ): - from math import pi - from sire.legacy.CAS import Symbol - - theta0 = pi / 2.0 - - # Create the new angle function. - amber_angle = _SireMM.AmberAngle(k_hard, theta0) - - # Generate the new angle expression. - expression = amber_angle.to_expression(Symbol("theta")) - - # Set the equilibrium angle to 90 degrees. - new_new_angles.set(idx0, idx1, idx2, expression) - - _logger.debug( - f" Stiffening angle: [{idx0.value()}-{idx1.value()}-{idx2.value()}], " - f"{p.function()} --> {expression}" - ) - - ang_idx = (idx0.value(), idx1.value(), idx2.value()) - modifications[mod_key]["stiffened_angles"].append(ang_idx) - - else: - new_new_angles.set(idx0, idx1, idx2, p.function()) - # Update the molecule. - mol = ( - mol.edit() - .set_property("angle" + suffix, new_new_angles) - .molecule() - .commit() - ) + mol = mol.edit().set_property("angle" + suffix, new_angles).molecule().commit() mol = ( mol.edit() .set_property("dihedral" + suffix, new_dihedrals) @@ -1109,6 +1459,21 @@ def _triple( .commit() ) + # Next we treat the remaining terms as a dual junction. + remaining_physical = [p for p in physical if p != idx] + mol = _dual( + mol, + bridge, + ghosts, + remaining_physical, + connectivity, + modifications, + k_hard=k_hard, + k_hard_ring=k_hard_ring, + is_lambda1=is_lambda1, + bridge_indices=bridge_indices, + ) + # Return the updated molecule. return mol @@ -1181,6 +1546,10 @@ def _higher( f"{_lam_sym} = {int(is_lambda1)}:" ) + # Nothing to do if there are no ghost atoms. + if not ghosts: + return mol + # Store the molecular info. info = mol.info() @@ -1296,6 +1665,7 @@ def _higher( bridge, ghosts, physical, + connectivity, modifications, k_hard=k_hard, k_soft=k_soft, @@ -1305,6 +1675,617 @@ def _higher( ) +def _remove_impropers(mol, ghosts, modifications, is_lambda1=False): + """ + Remove improper dihedral terms that bridge the ghost and physical systems. + + Parameters + ---------- + + mol : sire.mol.Molecule + The perturbable molecule. + + ghosts : List[sire.legacy.Mol.AtomIdx] + The list of ghost atoms. + + modifications : dict + A dictionary to store details of the modifications made. + + is_lambda1 : bool, optional + Whether to remove the improper dihedrals at lambda = 1. + + Returns + ------- + + mol : sire.mol.Molecule + The updated molecule. + """ + + # Nothing to do if there are no ghost atoms. + if not ghosts: + return mol + + # Store the molecular info. + info = mol.info() + + # Get the end state bond functions. + suffix = "1" if is_lambda1 else "0" + impropers = mol.property("improper" + suffix) + + # Initialise a container to store the updated bonded terms. + new_impropers = _SireMM.FourAtomFunctions(mol.info()) + + # Loop over the improper dihedrals. + for p in impropers.potentials(): + idx0 = info.atom_idx(p.atom0()) + idx1 = info.atom_idx(p.atom1()) + idx2 = info.atom_idx(p.atom2()) + idx3 = info.atom_idx(p.atom3()) + + # Remove any improper dihedrals that bridge the ghost and physical systems. + if not all(idx in ghosts for idx in (idx0, idx1, idx2, idx3)) and any( + idx in ghosts for idx in (idx0, idx1, idx2, idx3) + ): + _logger.debug( + f" Removing improper dihedral: [{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], {p.function()}" + ) + dih_idx = (idx0.value(), idx1.value(), idx2.value(), idx3.value()) + dih_idx = ",".join([str(i) for i in dih_idx]) + key = "lambda_1" if is_lambda1 else "lambda_0" + modifications[key]["removed_dihedrals"].append(dih_idx) + else: + new_impropers.set(idx0, idx1, idx2, idx3, p.function()) + + # Set the updated impropers. + mol = ( + mol.edit().set_property("improper" + suffix, new_impropers).molecule().commit() + ) + + # Return the updated molecule. + return mol + + +def _remove_residual_ghost_dihedrals(mol, ghosts, modifications, is_lambda1=False): + r""" + Remove dihedral terms that couple ghost and physical regions but were + not caught by the per-bridge junction handlers. This covers two cases: + + 1. Cross-bridge: both terminal atoms are ghost and both middle atoms + are physical. This arises when two ghost groups have adjacent bridge + atoms. The dihedral DR1-X1-X2-DR2 escapes both junction handlers + because each handler only sees its own ghost group. + + DR1 DR2 + \ / + X1------X2 + / \ + R1 R2 + + Removed dihedral: DR1-X1-X2-DR2 + + 2. Ghost middle: both terminal atoms are physical but at least one + middle atom is ghost. This arises in ring-breaking topologies where + a ghost atom is bonded to two bridge atoms. + + R1 R2 + \ / + X1-DR-X2 + / \ + R3 R4 + + Removed dihedrals: e.g. R1-X1-DR-X2, X1-DR-X2-R2 + + Parameters + ---------- + + mol : sire.mol.Molecule + The perturbable molecule. + + ghosts : List[sire.legacy.Mol.AtomIdx] + The list of ghost atoms at the current end state. + + modifications : dict + A dictionary to store details of the modifications made. + + is_lambda1 : bool, optional + Whether to modify dihedrals at lambda = 1. + + Returns + ------- + + mol : sire.mol.Molecule + The updated molecule. + """ + + # Nothing to do if there are no ghost atoms. + if not ghosts: + return mol + + # Store the molecular info. + info = mol.info() + + # Get the end state property. + if is_lambda1: + mod_key = "lambda_1" + suffix = "1" + else: + mod_key = "lambda_0" + suffix = "0" + + # Get the end state dihedral functions. + dihedrals = mol.property("dihedral" + suffix) + + # Initialise a container to store the updated dihedral functions. + new_dihedrals = _SireMM.FourAtomFunctions(mol.info()) + + # Track whether any modifications were made. + modified = False + + # Loop over the dihedral potentials. + for p in dihedrals.potentials(): + idx0 = info.atom_idx(p.atom0()) + idx1 = info.atom_idx(p.atom1()) + idx2 = info.atom_idx(p.atom2()) + idx3 = info.atom_idx(p.atom3()) + + # Case 1: Both terminals ghost, both middles physical (cross-bridge). + cross_bridge = ( + idx0 in ghosts + and idx3 in ghosts + and idx1 not in ghosts + and idx2 not in ghosts + ) + + # Case 2: Both terminals physical, at least one middle ghost + # (ring-breaking). + ghost_middle = ( + idx0 not in ghosts + and idx3 not in ghosts + and (idx1 in ghosts or idx2 in ghosts) + ) + + if cross_bridge or ghost_middle: + _logger.debug( + f" Removing residual ghost dihedral: " + f"[{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], " + f"{p.function()}" + ) + dih_idx = (idx0.value(), idx1.value(), idx2.value(), idx3.value()) + dih_idx = ",".join([str(i) for i in dih_idx]) + modifications[mod_key]["removed_dihedrals"].append(dih_idx) + modified = True + else: + new_dihedrals.set(idx0, idx1, idx2, idx3, p.function()) + + # Set the updated dihedrals. + if modified: + mol = ( + mol.edit() + .set_property("dihedral" + suffix, new_dihedrals) + .molecule() + .commit() + ) + + # Return the updated molecule. + return mol + + +def _remove_ghost_centre_angles(mol, ghosts, modifications, is_lambda1=False): + r""" + Remove angle terms where the central atom is ghost and both terminal + atoms are physical. These can arise in ring-breaking topologies where + a ghost atom is bonded to two bridge atoms. The per-bridge junction + handlers only catch angles with ghost terminal atoms, not ghost central + atoms. + + R1 R2 + \ / + X1-DR-X2 + / \ + R3 R4 + + Removed angle: X1-DR-X2 + + Parameters + ---------- + + mol : sire.mol.Molecule + The perturbable molecule. + + ghosts : List[sire.legacy.Mol.AtomIdx] + The list of ghost atoms at the current end state. + + modifications : dict + A dictionary to store details of the modifications made. + + is_lambda1 : bool, optional + Whether to modify angles at lambda = 1. + + Returns + ------- + + mol : sire.mol.Molecule + The updated molecule. + """ + + # Nothing to do if there are no ghost atoms. + if not ghosts: + return mol + + # Store the molecular info. + info = mol.info() + + # Get the end state property. + if is_lambda1: + mod_key = "lambda_1" + suffix = "1" + else: + mod_key = "lambda_0" + suffix = "0" + + # Get the end state angle functions. + angles = mol.property("angle" + suffix) + + # Initialise a container to store the updated angle functions. + new_angles = _SireMM.ThreeAtomFunctions(mol.info()) + + # Track whether any modifications were made. + modified = False + + # Loop over the angle potentials. + for p in angles.potentials(): + idx0 = info.atom_idx(p.atom0()) + idx1 = info.atom_idx(p.atom1()) + idx2 = info.atom_idx(p.atom2()) + + # Remove any angle where the central atom is ghost and both + # terminal atoms are physical. + if idx1 in ghosts and idx0 not in ghosts and idx2 not in ghosts: + _logger.debug( + f" Removing ghost centre angle: " + f"[{idx0.value()}-{idx1.value()}-{idx2.value()}], " + f"{p.function()}" + ) + ang_idx = (idx0.value(), idx1.value(), idx2.value()) + ang_idx = ",".join([str(i) for i in ang_idx]) + modifications[mod_key]["removed_angles"].append(ang_idx) + modified = True + else: + new_angles.set(idx0, idx1, idx2, p.function()) + + # Set the updated angles. + if modified: + mol = mol.edit().set_property("angle" + suffix, new_angles).molecule().commit() + + # Return the updated molecule. + return mol + + +def _soften_mixed_dihedrals( + mol, ghosts, modifications, soften_anchors=1.0, is_lambda1=False +): + r""" + Soften surviving mixed ghost/physical dihedral terms by scaling their + force constants. This is a post-processing step that runs after all + per-bridge junction handlers and residual removal passes. + + A "mixed" dihedral is one that involves at least one ghost atom and at + least one physical atom. These dihedrals couple the ghost and physical + regions and can cause dynamics crashes at small lambda when the ghost + atoms start gaining softcore nonbonded interactions but are constrained + too tightly by bonded terms. + + When ``soften_anchors`` is 1.0 (default), no modifications are made. + When 0.0, all mixed dihedrals are removed. Intermediate values scale + the force constants. + + Parameters + ---------- + + mol : sire.mol.Molecule + The perturbable molecule. + + ghosts : List[sire.legacy.Mol.AtomIdx] + The list of ghost atoms at the current end state. + + modifications : dict + A dictionary to store details of the modifications made. + + soften_anchors : float, optional + Scale factor for mixed dihedral force constants (0.0 to 1.0). + + is_lambda1 : bool, optional + Whether to modify dihedrals at lambda = 1. + + Returns + ------- + + mol : sire.mol.Molecule + The updated molecule. + """ + + # Nothing to do if there are no ghost atoms or no softening is requested. + if not ghosts or soften_anchors >= 1.0: + return mol + + # Store the molecular info. + info = mol.info() + + # Get the end state property. + if is_lambda1: + mod_key = "lambda_1" + suffix = "1" + else: + mod_key = "lambda_0" + suffix = "0" + + # Get the end state dihedral functions. + dihedrals = mol.property("dihedral" + suffix) + + # Initialise a container to store the updated dihedral functions. + new_dihedrals = _SireMM.FourAtomFunctions(mol.info()) + + # Track whether any modifications were made. + modified = False + + # Loop over the dihedral potentials. + for p in dihedrals.potentials(): + idx0 = info.atom_idx(p.atom0()) + idx1 = info.atom_idx(p.atom1()) + idx2 = info.atom_idx(p.atom2()) + idx3 = info.atom_idx(p.atom3()) + + atoms = (idx0, idx1, idx2, idx3) + has_ghost = any(a in ghosts for a in atoms) + has_physical = any(a not in ghosts for a in atoms) + + if has_ghost and has_physical: + dih_idx_str = ",".join(str(a.value()) for a in atoms) + if soften_anchors > 0.0: + scaled = p.function() * soften_anchors + new_dihedrals.set(idx0, idx1, idx2, idx3, scaled) + _logger.debug( + f" Softening mixed dihedral: [{dih_idx_str}], " + f"scale={soften_anchors}" + ) + modifications[mod_key]["softened_dihedrals"].append(dih_idx_str) + else: + _logger.debug( + f" Removing mixed dihedral: [{dih_idx_str}], {p.function()}" + ) + modifications[mod_key]["removed_dihedrals"].append(dih_idx_str) + modified = True + else: + new_dihedrals.set(idx0, idx1, idx2, idx3, p.function()) + + # Set the updated dihedrals. + if modified: + mol = ( + mol.edit() + .set_property("dihedral" + suffix, new_dihedrals) + .molecule() + .commit() + ) + n_softened = len(modifications[mod_key]["softened_dihedrals"]) + lam = f"{_lam_sym}={int(is_lambda1)}" + if soften_anchors > 0.0: + _logger.debug( + f"Softened {n_softened} mixed ghost/physical dihedrals at " + f"{lam} (scale={soften_anchors})" + ) + else: + _logger.debug(f"Removed all mixed ghost/physical dihedrals at {lam}") + + # Return the updated molecule. + return mol + + +def _check_rotamer_anchors( + mol, + bridges, + physical, + ghosts, + modifications, + is_lambda1=False, + stiffen=False, + k_rotamer=50, +): + r""" + Detect and optionally stiffen rotamer anchor dihedrals. + + After ghost modifications, each bridge atom retains anchor dihedral(s) + that constrain the ghost group's rotation. If the bridge--physical bond + is not in a ring and the bridge atom is sp3-hybridised, the anchor + dihedral is likely a rotamer. This can allow the ghost group to flip + between rotameric states at intermediate lambda, degrading convergence + (see Boresch et al., JCTC 2021). + + R1 DR1 + \ / + \ / + X ---DR2 X is sp3, P-X is rotatable + / \ + / \ + R2 DR3 + + Rotamer anchors are always logged as warnings. When ``stiffen=True``, + affected dihedrals (those spanning the rotatable bond with at least one + ghost terminal) are replaced with a single n=1 cosine well: + + V(phi) = k [1 + cos(phi - pi)] + + which has a single minimum at phi = 0 (trans) and a barrier of 2k. + + .. note:: + + Stiffening is not currently enabled. When wiring in, add + ``stiffen_rotamers`` and ``k_rotamer`` parameters to ``modify()`` + and expose them through the CLI. The ``modifications`` dict will + also need a ``"stiffened_dihedrals"`` key initialised to an empty + list for each end state. + + Parameters + ---------- + + mol : sire.mol.Molecule + The perturbable molecule. + + bridges : dict + A dictionary mapping bridge atoms to their ghost neighbours. + + physical : dict + A dictionary mapping bridge atoms to their physical neighbours. + + ghosts : List[sire.legacy.Mol.AtomIdx] + The list of ghost atoms at the current end state. + + modifications : dict + A dictionary to store details of the modifications made. + + is_lambda1 : bool, optional + Whether to check/modify the lambda = 1 end state. + + stiffen : bool, optional + Whether to replace rotamer anchor dihedrals with a stiff cosine + well. If False (default), only log warnings. + + k_rotamer : float, optional + The force constant for the replacement cosine term. The resulting + barrier height is 2 * k_rotamer. (In kcal/mol) Only used when + ``stiffen=True``. The default of 50 is a placeholder and has not + been calibrated against simulation data. It should be benchmarked + before use. + + Returns + ------- + + mol : sire.mol.Molecule + The updated molecule (unchanged if ``stiffen=False``). + """ + + # Nothing to do if there are no bridges. + if not bridges: + return mol + + from rdkit.Chem import HybridizationType + from sire.convert import to_rdkit + + lam = int(is_lambda1) + + # Link the molecule to the desired end state and convert to RDKit. + if is_lambda1: + end_state_mol = _morph.link_to_perturbed(mol) + else: + end_state_mol = _morph.link_to_reference(mol) + + try: + rdmol = to_rdkit(end_state_mol) + except Exception as e: + _logger.warning(f"Failed to convert molecule to RDKit for rotamer check: {e}") + return mol + + # Identify rotatable bridge--physical bond pairs. + rotatable_bonds = set() + for bridge in bridges: + for p in physical[bridge]: + b_idx = bridge.value() + p_idx = p.value() + + rd_bond = rdmol.GetBondBetweenAtoms(p_idx, b_idx) + if rd_bond is None or rd_bond.IsInRing(): + continue + + rd_bridge = rdmol.GetAtomWithIdx(b_idx) + hybridisation = rd_bridge.GetHybridization() + + if hybridisation in ( + HybridizationType.SP3, + HybridizationType.SP3D, + HybridizationType.SP3D2, + ): + rotatable_bonds.add((p, bridge)) + if not stiffen: + _logger.warning( + f"Potential rotamer anchor at {_lam_sym} = {lam}: " + f"bond {p_idx}-{b_idx} is a rotatable sp3 bond. " + f"Surviving anchor dihedrals may allow rotameric " + f"transitions of ghost atoms." + ) + + # If not stiffening, or no rotatable bonds found, return early. + if not stiffen or not rotatable_bonds: + return mol + + # Stiffen the anchor dihedrals spanning the rotatable bonds. + + from math import pi + + from sire.legacy.CAS import Symbol + + # Get the end state property suffix. + if is_lambda1: + mod_key = "lambda_1" + suffix = "1" + else: + mod_key = "lambda_0" + suffix = "0" + + # Store the molecular info. + info = mol.info() + + # Get the end state dihedral functions. + dihedrals = mol.property("dihedral" + suffix) + + # Create the stiff single-well replacement: k[1 + cos(phi - pi)]. + # Minimum at phi = 0 (trans), barrier height = 2k. + replacement = _SireMM.AmberDihedral( + _SireMM.AmberDihPart(k_rotamer, 1, pi) + ).to_expression(Symbol("phi")) + + # Initialise a container to store the updated dihedral functions. + new_dihedrals = _SireMM.FourAtomFunctions(mol.info()) + + modified = False + + for p in dihedrals.potentials(): + idx0 = info.atom_idx(p.atom0()) + idx1 = info.atom_idx(p.atom1()) + idx2 = info.atom_idx(p.atom2()) + idx3 = info.atom_idx(p.atom3()) + + # Check if the central bond (idx1-idx2) is a rotatable bridge bond + # and at least one terminal atom is ghost. + bond_pair = (idx1, idx2) + bond_pair_rev = (idx2, idx1) + has_ghost_terminal = idx0 in ghosts or idx3 in ghosts + + if has_ghost_terminal and ( + bond_pair in rotatable_bonds or bond_pair_rev in rotatable_bonds + ): + _logger.debug( + f" Stiffening rotamer anchor dihedral: " + f"[{idx0.value()}-{idx1.value()}-{idx2.value()}-{idx3.value()}], " + f"{p.function()} --> {replacement}" + ) + new_dihedrals.set(idx0, idx1, idx2, idx3, replacement) + dih_idx = (idx0.value(), idx1.value(), idx2.value(), idx3.value()) + dih_idx = ",".join([str(i) for i in dih_idx]) + modifications[mod_key]["stiffened_dihedrals"].append(dih_idx) + modified = True + else: + new_dihedrals.set(idx0, idx1, idx2, idx3, p.function()) + + if modified: + mol = ( + mol.edit() + .set_property("dihedral" + suffix, new_dihedrals) + .molecule() + .commit() + ) + + return mol + + def _create_connectivity(mol): """ Create a connectivity object for an end state molecule. diff --git a/tests/test_ghostly.py b/tests/test_ghostly.py index a1c52ee..bf7709b 100644 --- a/tests/test_ghostly.py +++ b/tests/test_ghostly.py @@ -55,12 +55,9 @@ def test_hexane_to_propane(): ) # Check that the missing dihedrals are not in the new dihedrals. - assert ( - all( - check_dihedral(info, new_dihedrals.potentials(), *dihedral) - for dihedral in missing_dihedrals - ) - == False + assert not all( + check_dihedral(info, new_dihedrals.potentials(), *dihedral) + for dihedral in missing_dihedrals ) @@ -108,21 +105,15 @@ def test_toluene_to_pyridine(): info = mols[0].info() # Check that the missing dihedrals are in the original dihedrals. - assert ( - all( - check_dihedral(info, dihedrals.potentials(), *dihedral) - for dihedral in missing_dihedrals - ) - == True + assert all( + check_dihedral(info, dihedrals.potentials(), *dihedral) + for dihedral in missing_dihedrals ) # Check that the missing dihedrals are not in the new dihedrals. - assert ( - all( - check_dihedral(info, new_dihedrals.potentials(), *dihedral) - for dihedral in missing_dihedrals - ) - == False + assert not all( + check_dihedral(info, new_dihedrals.potentials(), *dihedral) + for dihedral in missing_dihedrals ) # Create a list of angle IDs for the modified angles. @@ -132,7 +123,8 @@ def test_toluene_to_pyridine(): ] # Functional form of the modified angles. - expression = "100 [theta - 1.5708]^2" + # Bridge atom 1 is in a ring, so k_hard_ring (75) is used. + expression = "75 [theta - 1.5708]^2" # Check that the original angles don't have the modified functional form. for p in angles.potentials(): @@ -202,21 +194,15 @@ def test_acetone_to_propenol(): info = mols[0].info() # Check that the missing dihedrals are in the original dihedrals at lambda = 0. - assert ( - all( - check_dihedral(info, dihedrals0.potentials(), *dihedral) - for dihedral in missing_dihedrals0 - ) - == True + assert all( + check_dihedral(info, dihedrals0.potentials(), *dihedral) + for dihedral in missing_dihedrals0 ) # Check that the missing dihedrals are not in the new dihedrals at lambda = 0. - assert ( - all( - check_dihedral(info, new_dihedrals0.potentials(), *dihedral) - for dihedral in missing_dihedrals0 - ) - == False + assert not all( + check_dihedral(info, new_dihedrals0.potentials(), *dihedral) + for dihedral in missing_dihedrals0 ) # Create dihedral IDs for the missing dihedrals at lambda = 1. @@ -226,47 +212,37 @@ def test_acetone_to_propenol(): ] # Check that the missing dihedrals are in the original dihedrals at lambda = 1. - assert ( - all( - check_dihedral(info, dihedrals1.potentials(), *dihedral) - for dihedral in missing_dihedrals1 - ) - == True + assert all( + check_dihedral(info, dihedrals1.potentials(), *dihedral) + for dihedral in missing_dihedrals1 ) # Check that the missing dihedrals are not in the new dihedrals at lambda = 1. - assert ( - all( - check_dihedral(info, new_dihedrals1.potentials(), *dihedral) - for dihedral in missing_dihedrals1 - ) - == False + assert not all( + check_dihedral(info, new_dihedrals1.potentials(), *dihedral) + for dihedral in missing_dihedrals1 ) # Create angle IDs for the removed angles at lambda = 1. + # Atom 9 is transmuting (worst score), so it is the sacrificial atom. removed_angles = [ - (AtomIdx(1), AtomIdx(3), AtomIdx(7)), + (AtomIdx(7), AtomIdx(3), AtomIdx(9)), ] # Check that the removed angles are in the original angles at lambda = 1. - assert ( - all(check_angle(info, angles1.potentials(), *angle) for angle in removed_angles) - == True + assert all( + check_angle(info, angles1.potentials(), *angle) for angle in removed_angles ) # Check that the removed angles are not in the new angles at lambda = 1. - assert ( - all( - check_angle(info, new_angles1.potentials(), *angle) - for angle in removed_angles - ) - == False + assert not all( + check_angle(info, new_angles1.potentials(), *angle) for angle in removed_angles ) # Create angle IDs for the modified angles at lambda = 1. modified_angles = [ + (AtomIdx(1), AtomIdx(3), AtomIdx(7)), (AtomIdx(7), AtomIdx(3), AtomIdx(8)), - (AtomIdx(7), AtomIdx(3), AtomIdx(9)), ] # Functional form of the modified angles. @@ -291,6 +267,195 @@ def test_acetone_to_propenol(): assert str(p.function()) == expression +def test_ejm49_to_ejm31(): + """ + Test ghost atom modifications for the TYK ligands EJM 49 to 31. This is assert + more complex perturbation. + """ + + # Load the system. Here pruned means that the atom mapping has pruned + # atoms where the constraint changes between the end states, which is + # what is used by OpenFE. + mols = sr.load_test_files("ejm49_ejm31_pruned.bss") + + # Store the orginal angles and dihedrals at lambda = 0 and lambda = 1. + angles0 = mols[0].property("angle0") + angles1 = mols[0].property("angle1") + dihedrals0 = mols[0].property("dihedral0") + dihedrals1 = mols[0].property("dihedral1") + improper0 = mols[0].property("improper0") + improper1 = mols[0].property("improper1") + + # Apply the ghost atom modifications. + new_mols, _ = modify(mols) + + # Get the new angles and dihedrals. + new_angles0 = new_mols[0].property("angle0") + new_angles1 = new_mols[0].property("angle1") + new_dihedrals0 = new_mols[0].property("dihedral0") + new_dihedrals1 = new_mols[0].property("dihedral1") + new_improper0 = new_mols[0].property("improper0") + new_improper1 = new_mols[0].property("improper1") + + # The number of angles should remain the same at lambda = 0. + assert angles0.num_functions() == new_angles0.num_functions() + + # The number of dihedrals should be five fewer at lambda = 0. + assert dihedrals0.num_functions() - 5 == new_dihedrals0.num_functions() + + # The number of impropers shoudld be three fewer at lambda = 0. + assert improper0.num_functions() - 3 == new_improper0.num_functions() + + # The number of angles should remain the same at lambda = 1. + assert angles1.num_functions() == new_angles1.num_functions() + + # The number of dihedrals should be four fewer at lambda = 1. + assert dihedrals1.num_functions() - 4 == new_dihedrals1.num_functions() + + # The number of impropers should be six fewer at lambda = 1. + assert improper1.num_functions() - 6 == new_improper1.num_functions() + + # Create dihedral IDs for the missing dihedrals at lambda = 0. + + from sire.legacy.Mol import AtomIdx + + missing_dihedrals0 = [ + (AtomIdx(18), AtomIdx(17), AtomIdx(39), AtomIdx(40)), + (AtomIdx(18), AtomIdx(17), AtomIdx(39), AtomIdx(41)), + (AtomIdx(18), AtomIdx(17), AtomIdx(39), AtomIdx(42)), + (AtomIdx(33), AtomIdx(16), AtomIdx(17), AtomIdx(39)), + (AtomIdx(14), AtomIdx(16), AtomIdx(17), AtomIdx(39)), + ] + + # Store the molecular info. + info = mols[0].info() + + # Check that the missing dihedrals are in the original dihedrals at lambda = 0. + assert all( + check_dihedral(info, dihedrals0.potentials(), *dihedral) + for dihedral in missing_dihedrals0 + ) + + # Check that the missing dihedrals are not in the new dihedrals at lambda = 0. + assert not all( + check_dihedral(info, new_dihedrals0.potentials(), *dihedral) + for dihedral in missing_dihedrals0 + ) + + # Create dihedral IDs for the missing dihedrals at lambda = 1. + missing_dihedrals1 = [ + (AtomIdx(18), AtomIdx(17), AtomIdx(20), AtomIdx(21)), + (AtomIdx(18), AtomIdx(17), AtomIdx(20), AtomIdx(25)), + (AtomIdx(20), AtomIdx(17), AtomIdx(16), AtomIdx(33)), + (AtomIdx(14), AtomIdx(16), AtomIdx(17), AtomIdx(20)), + ] + + # Check that the missing dihedrals are in the original dihedrals at lambda = 1. + assert all( + check_dihedral(info, dihedrals1.potentials(), *dihedral) + for dihedral in missing_dihedrals1 + ) + + # Check that the missing dihedrals are not in the new dihedrals at lambda = 1. + assert not all( + check_dihedral(info, new_dihedrals1.potentials(), *dihedral) + for dihedral in missing_dihedrals1 + ) + + # Create angle IDs for the modified angles at lambda = 0. + modified_angles0 = [ + (AtomIdx(16), AtomIdx(17), AtomIdx(39)), + (AtomIdx(18), AtomIdx(17), AtomIdx(39)), + ] + + # Functional form of the modified angles. + expression = "100 [theta - 1.5708]^2" + + # Check that the original angles don't have the modified functional form. + for p in angles0.potentials(): + idx0 = info.atom_idx(p.atom0()) + idx1 = info.atom_idx(p.atom1()) + idx2 = info.atom_idx(p.atom2()) + + if (idx0, idx1, idx2) in modified_angles0: + assert str(p.function()) != expression + + # Check that the modified angles have the correct functional form. + for p in new_angles0.potentials(): + idx0 = info.atom_idx(p.atom0()) + idx1 = info.atom_idx(p.atom1()) + idx2 = info.atom_idx(p.atom2()) + + if (idx0, idx1, idx2) in modified_angles0: + assert str(p.function()) == expression + + # Create angle IDs for the modified angles at lambda = 1. + modified_angles1 = [ + (AtomIdx(16), AtomIdx(17), AtomIdx(20)), + (AtomIdx(18), AtomIdx(17), AtomIdx(20)), + ] + + # Check that the original angles don't have the modified functional form. + for p in angles1.potentials(): + idx0 = info.atom_idx(p.atom0()) + idx1 = info.atom_idx(p.atom1()) + idx2 = info.atom_idx(p.atom2()) + + if (idx0, idx1, idx2) in modified_angles1: + assert str(p.function()) != expression + + # Check that the modified angles have the correct functional form. + for p in new_angles1.potentials(): + idx0 = info.atom_idx(p.atom0()) + idx1 = info.atom_idx(p.atom1()) + idx2 = info.atom_idx(p.atom2()) + + if (idx0, idx1, idx2) in modified_angles1: + assert str(p.function()) == expression + + # Create improper IDs for the missing impropers at lambda = 0. + missing_impropers0 = [ + (AtomIdx(17), AtomIdx(16), AtomIdx(18), AtomIdx(39)), + (AtomIdx(16), AtomIdx(39), AtomIdx(18), AtomIdx(17)), + (AtomIdx(17), AtomIdx(39), AtomIdx(16), AtomIdx(18)), + ] + + # Check that the missing impropers are in the original impropers at lambda = 0. + assert all( + check_improper(info, improper0.potentials(), *improper) + for improper in missing_impropers0 + ) + + # Check that the missing impropers are not in the new impropers at lambda = 0. + assert not all( + check_improper(info, new_improper0.potentials(), *improper) + for improper in missing_impropers0 + ) + + # Create improper IDs for the missing impropers at lambda = 1. + missing_impropers1 = [ + (AtomIdx(17), AtomIdx(25), AtomIdx(21), AtomIdx(20)), + (AtomIdx(17), AtomIdx(20), AtomIdx(16), AtomIdx(18)), + (AtomIdx(16), AtomIdx(20), AtomIdx(18), AtomIdx(17)), + (AtomIdx(17), AtomIdx(20), AtomIdx(16), AtomIdx(18)), + (AtomIdx(20), AtomIdx(25), AtomIdx(17), AtomIdx(21)), + (AtomIdx(16), AtomIdx(20), AtomIdx(18), AtomIdx(17)), + (AtomIdx(20), AtomIdx(17), AtomIdx(21), AtomIdx(25)), + ] + + # Check that the missing impropers are in the original impropers at lambda = 1. + assert all( + check_improper(info, improper1.potentials(), *improper) + for improper in missing_impropers1 + ) + + # Check that the missing impropers are not in the new impropers at lambda = 1. + assert not all( + check_improper(info, new_improper1.potentials(), *improper) + for improper in missing_impropers1 + ) + + def check_angle(info, potentials, idx0, idx1, idx2): """ Check if an angle potential is in a list of potentials. @@ -322,3 +487,20 @@ def check_dihedral(info, potentials, idx0, idx1, idx2, idx3): return True return False + + +def check_improper(info, potentials, idx0, idx1, idx2, idx3): + """ + Check if an improper potential is in a list of potentials. + """ + + for p in potentials: + if ( + idx0 == info.atom_idx(p.atom0()) + and idx1 == info.atom_idx(p.atom1()) + and idx2 == info.atom_idx(p.atom2()) + and idx3 == info.atom_idx(p.atom3()) + ): + return True + + return False