diff --git a/.github/renovate.json b/.github/renovate.json new file mode 100644 index 0000000..6e94515 --- /dev/null +++ b/.github/renovate.json @@ -0,0 +1,8 @@ +{ + "$schema": "https://docs.renovatebot.com/renovate-schema.json", + "extends": [ + "github>CCPBioSim/.github//renovate/default", + "github>CCPBioSim/.github//renovate/python" + ], + "reviewers": ["harryswift01", "jimboid"] +} diff --git a/.github/workflows/coveralls-main.yaml b/.github/workflows/coveralls-main.yaml new file mode 100644 index 0000000..c561366 --- /dev/null +++ b/.github/workflows/coveralls-main.yaml @@ -0,0 +1,35 @@ +name: Coveralls Main + +on: + push: + branches: [main] + + workflow_dispatch: + +permissions: + contents: read + +concurrency: + group: coveralls-main-${{ github.ref }} + cancel-in-progress: true + +jobs: + coveralls: + name: Build main coverage baseline + uses: CCPBioSim/.github/.github/workflows/ci-test-conda.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: read + checks: write + statuses: write + with: + package: prepmd + python-versions: '["3.13"]' + coverage-python-version: "3.13" + os: '["ubuntu-26.04", "macos-26"]' + environment-file: environment.yaml + activate-environment: prep + channels: salilab,conda-forge + test-path: test + pytest-args: -v + coverage: true + timeout-minutes: 360 diff --git a/.github/workflows/pr.yaml b/.github/workflows/pr.yaml new file mode 100644 index 0000000..e66d71f --- /dev/null +++ b/.github/workflows/pr.yaml @@ -0,0 +1,63 @@ +name: Pull Request CI + +on: + pull_request: + +permissions: + contents: read + +concurrency: + group: pr-ci-${{ github.ref }} + cancel-in-progress: true + +jobs: + unit: + name: Unit + uses: CCPBioSim/.github/.github/workflows/ci-test-conda.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: read + checks: write + statuses: write + with: + package: prepmd + python-versions: '["3.13"]' + coverage-python-version: "3.13" + os: '["ubuntu-26.04", "macos-26"]' + environment-file: environment.yaml + activate-environment: prep + channels: salilab,conda-forge + test-path: test + pytest-args: -v + coverage: true + timeout-minutes: 360 + + docs: + name: Docs + needs: unit + uses: CCPBioSim/.github/.github/workflows/ci-docs.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: read + with: + python-versions: '["3.13"]' + build-dir: build + timeout-minutes: 20 + + pre-commit: + name: Pre-commit + needs: unit + uses: CCPBioSim/.github/.github/workflows/ci-pre-commit.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: read + with: + python-version: "3.13" + timeout-minutes: 15 + + conda-build: + name: Conda build + needs: unit + uses: CCPBioSim/.github/.github/workflows/ci-conda-build.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: read + with: + python-versions: '["3.13"]' + timeout-minutes: 360 diff --git a/.github/workflows/python-app.yml b/.github/workflows/python-app.yml deleted file mode 100644 index 7a370d3..0000000 --- a/.github/workflows/python-app.yml +++ /dev/null @@ -1,87 +0,0 @@ -# This workflow will install Python dependencies, run tests and lint with a single version of Python -# For more information see: https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python - -name: prepmd CI - -on: - push: - branches: [main] - pull_request: - schedule: - - cron: '0 8 * * 1' - workflow_dispatch: - -jobs: - tests: - name: Run tests - runs-on: ${{ matrix.os }} - strategy: - fail-fast: false - matrix: - os: [ubuntu-24.04, macos-15] # no antechamber on windows-2025 :( - python-version: ["3.13"] - - #env: - #KEY_MODELLER: ${{ secrets.MODELLER_KEY }} - - steps: - - name: Checkout repo - uses: actions/checkout@v5.0.0 - - - name: Cache conda packages - uses: actions/cache@v4 - with: - path: ~/.conda/pkgs - key: ${{ runner.os }}-conda-${{ hashFiles('environment.yaml') }} - restore-keys: | - ${{ runner.os }}-conda- - - - name: Set up Mambaforge - uses: conda-incubator/setup-miniconda@v3 - with: - use-mamba: true - python-version: ${{ matrix.python-version }} - environment-file: environment.yaml - activate-environment: prep - auto-update-conda: true - auto-activate-base: false - - - name: Verify environment - shell: bash -l {0} - run: | - conda info - conda list - - - name: Install package and test dependencies - shell: bash -l {0} - run: | - which pip - pip install -e . - pip install pytest-cov - pip install pdbe-mmcif-validator - conda list - - - name: Run test suite with coverage - shell: bash -l {0} - run: | - pytest -v --cov=prepmd --cov-report=term-missing --cov-report=lcov - - - name: Coveralls GitHub Action - uses: coverallsapp/github-action@5cbfd81b66ca5d10c19b062c04de0199c215fb6e # v2.3.7 - with: - github-token: ${{ secrets.GITHUB_TOKEN }} - file: coverage.lcov - flag-name: ${{ matrix.os }}-py${{ matrix.python-version }} - parallel: true - fail-on-error: false - - coveralls: - name: Finalise Coveralls - needs: tests - runs-on: ubuntu-24.04 - steps: - - name: Coveralls Finished - uses: coverallsapp/github-action@5cbfd81b66ca5d10c19b062c04de0199c215fb6e # v2.3.7 - with: - github-token: ${{ secrets.GITHUB_TOKEN }} - parallel-finished: true diff --git a/.github/workflows/release.yaml b/.github/workflows/release.yaml new file mode 100644 index 0000000..e0b4d7b --- /dev/null +++ b/.github/workflows/release.yaml @@ -0,0 +1,63 @@ +name: Release prepmd + +on: + workflow_dispatch: + inputs: + version: + description: Version to release, in X.Y.Z form. + required: true + type: string + dry-run: + description: Check and build everything but do not open the version PR, tag, release or publish. + type: boolean + default: false + +permissions: + contents: read + +jobs: + prepare: + name: Prepare + uses: CCPBioSim/.github/.github/workflows/release-prepare.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: write + pull-requests: write + with: + version: ${{ inputs.version }} + dry-run: ${{ inputs.dry-run }} + version-files: >- + [ + {"path": "prepmd/__init__.py", "pattern": "^__version__ =.*", "replace": "__version__ = \"{version}\""}, + {"path": "CITATION.cff", "pattern": "^version:.*", "replace": "version: {version}"}, + {"path": "CITATION.cff", "pattern": "^date-released:.*", "replace": "date-released: '{date}'"} + ] + secrets: + AUTO_PR_MERGE: ${{ secrets.AUTO_PR_MERGE }} + ZULIP_BOT_EMAIL: ${{ secrets.ZULIP_BOT_EMAIL }} + ZULIP_BOT_API_KEY: ${{ secrets.ZULIP_BOT_API_KEY }} + ZULIP_SITE: ${{ secrets.ZULIP_SITE }} + + website: + name: Website + needs: prepare + uses: CCPBioSim/.github/.github/workflows/release-website.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: read + with: + version: ${{ inputs.version }} + tag: ${{ needs.prepare.outputs.tag }} + dry-run: ${{ inputs.dry-run }} + secrets: + BUILD_TOKEN: ${{ secrets.BUILD_TOKEN }} + + conda: + name: Conda + needs: prepare + uses: CCPBioSim/.github/.github/workflows/release-conda.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: read + with: + tag: ${{ needs.prepare.outputs.tag || inputs.version }} + dry-run: ${{ inputs.dry-run }} + secrets: + ANACONDA_TOKEN: ${{ secrets.ANACONDA_TOKEN }} diff --git a/.github/workflows/release.yml b/.github/workflows/release.yml deleted file mode 100644 index 2307f83..0000000 --- a/.github/workflows/release.yml +++ /dev/null @@ -1,285 +0,0 @@ -name: Release - -on: - workflow_dispatch: - inputs: - version: - required: true - default: 'x.y.z' - -permissions: - contents: write - pull-requests: write - -jobs: - checks: - name: Version check - runs-on: ubuntu-24.04 - steps: - - name: Checkout repository - id: repo - uses: actions/checkout@8e8c483db84b4bee98b60c0593521ed34d9990e8 # v6 - with: - fetch-depth: 0 - - - name: Set up Python - uses: actions/setup-python@a309ff8b426b58ec0e2a45f0f869d46889d02405 # v6.2.0 - with: - python-version: 3.13.0 - - - name: Get latest release tag - id: latestreleased - run: | - PREVIOUS_VERSION=$(git tag --sort=-version:refname | head -n 1) - PREVIOUS_VERSION=${PREVIOUS_VERSION#v} - - if [ -z "$PREVIOUS_VERSION" ]; then - PREVIOUS_VERSION="0.0.0" - fi - - echo "release_tag=$PREVIOUS_VERSION" >> "$GITHUB_OUTPUT" - echo "$PREVIOUS_VERSION" - - - name: version comparison - id: compare - run: | - pip3 install semver - output=$(pysemver compare ${{ steps.latestreleased.outputs.release_tag }} ${{ github.event.inputs.version }}) - if [ "$output" -ge 0 ]; then exit 1; fi - - version: - name: prepare ${{ github.event.inputs.version }} - needs: checks - runs-on: ubuntu-24.04 - steps: - - - name: checkout - uses: actions/checkout@8e8c483db84b4bee98b60c0593521ed34d9990e8 # v6 - - - name: Change version in repo and CITATION.cff - run: | - # Update Python package version - sed -i "s/__version__ =.*/__version__ = \"${{ github.event.inputs.version }}\"/g" prepmd/__init__.py - - # Update CITATION.cff version and date-released - if [ -f CITATION.cff ]; then - sed -i -E "s/^(version:\s*).*/\1${{ github.event.inputs.version }}/" CITATION.cff - sed -i -E "s/^(date-released:\s*).*/\1'$(date -u +%F)'/" CITATION.cff - fi - - - name: send PR - id: pr_id - uses: peter-evans/create-pull-request@c0f553fe549906ede9cf27b5156039d195d2ece0 # v8.1.0 - with: - commit-message: Update version to ${{ github.event.inputs.version }} - branch: version-update - title: "Update to version ${{ github.event.inputs.version }}" - body: | - Update version - - Update the __init__.py with new release - - Update CITATION.cff version & date-released - - Auto-generated by [CI] - committer: version-updater - author: version-updater - base: main - signoff: false - draft: false - - - name: auto approve review - uses: hmarr/auto-approve-action@f0939ea97e9205ef24d872e76833fa908a770363 # v4.0.0 - with: - pull-request-number: ${{ steps.pr_id.outputs.pull-request-number }} - review-message: "Auto approved version bump PR" - github-token: ${{ secrets.AUTO_PR_MERGE }} - - - name: merge PR - run: gh pr merge --merge --delete-branch --auto "${{ steps.pr_id.outputs.pull-request-number }}" - env: - GH_TOKEN: ${{ secrets.GITHUB_TOKEN }} - - tag: - name: tag release - needs: version - runs-on: ubuntu-24.04 - steps: - - name: Checkout repository - uses: actions/checkout@8e8c483db84b4bee98b60c0593521ed34d9990e8 # v6 - with: - ref: main - - - name: tag v${{ github.event.inputs.version }} - run: | - git config user.name github-actions - git config user.email github-actions@github.com - git tag ${{ github.event.inputs.version }} - git push origin tag ${{ github.event.inputs.version }} - - release: - name: make github release - needs: tag - runs-on: ubuntu-24.04 - steps: - - - name: create release - id: create_release - uses: softprops/action-gh-release@3d0d9888cb7fd7b750713d6e236d1fcb99157228 # v3.0.2 - with: - name: v${{ github.event.inputs.version }} - generate_release_notes: true - tag_name: ${{ github.event.inputs.version }} - - - name: fetch generated changelog - id: changelog - if: steps.create_release.outcome == 'success' - continue-on-error: true - env: - GH_TOKEN: ${{ secrets.GITHUB_TOKEN }} - run: | - # generate_release_notes builds the changelog server-side, so pull - # it back down here rather than trying to reconstruct it locally. - gh release view "${{ github.event.inputs.version }}" --repo "${{ github.repository }}" --json body -q .body > changelog.txt - - - name: notify Zulip of release - if: steps.create_release.outcome == 'success' - continue-on-error: true - env: - ZULIP_BOT_EMAIL: ${{ secrets.ZULIP_BOT_EMAIL }} - ZULIP_BOT_API_KEY: ${{ secrets.ZULIP_BOT_API_KEY }} - ZULIP_SITE: ${{ secrets.ZULIP_SITE }} - REPO: ${{ github.event.repository.name }} - TAG: v${{ github.event.inputs.version }} - RELEASE_URL: ${{ steps.create_release.outputs.url }} - run: | - CONTENT=$(cat < { "Full Name": "orcid url" } - authors=$(yq -o=json '.authors' "$CFF_FILE" | jq -c ' - map({ - key: ( - if .name then .name - else ([."given-names", ."family-names"] | map(select(. != null)) | join(" ")) - end - ), - value: (.orcid // "") - }) | from_entries') - - payload=$(jq -nc \ - --arg name "$NAME" \ - --arg github "$github_url" \ - --arg docs "$docs_url" \ - --arg version "$VERSION" \ - --argjson meta "$meta" \ - --argjson authors "$authors" \ - '{ - name: $name, - meta: ({github: $github, docs: $docs} + $meta + {version: $version, authors: $authors}) - }') - - echo "payload: $payload" - echo "json=$payload" >> "$GITHUB_OUTPUT" - - - name: Post version to website metadata repo - uses: peter-evans/repository-dispatch@28959ce8df70de7be546dd1250a005dd32156697 # v4.0.1 - with: - token: ${{ secrets.BUILD_TOKEN }} - repository: ccpbiosim/ccpbiosim.github.io - event-type: software - client-payload: ${{ steps.payload.outputs.json }} - - conda: - name: publish conda to anaconda.org - needs: [tag, release] - runs-on: ubuntu-24.04 - - steps: - - name: Checkout - uses: actions/checkout@v4 - with: - ref: main - - - name: Set version in conda recipe - run: | - sed -i "s/^{% set version = \".*\" %}$/{% set version = \"${{ github.event.inputs.version }}\" %}/" conda-recipe/meta.yaml - grep '{% set version' conda-recipe/meta.yaml - - - name: Setup Conda (conda-forge only) - uses: conda-incubator/setup-miniconda@v3 - with: - miniforge-variant: Miniforge3 - python-version: "3.13.2" - auto-activate: true - activate-environment: base - conda-remove-defaults: true - channels: conda-forge - channel-priority: strict - - - name: Install build tools - shell: bash -l {0} - run: | - conda install -y \ - --override-channels \ - --channel conda-forge \ - conda-build \ - anaconda-client - - conda config --set anaconda_upload no - conda list conda-build - conda build --help - - - name: Build package - shell: bash -l {0} - run: | - PKG_PATH=$(conda build conda-recipe --output --override-channels --channel conda-forge) - echo "PKG_PATH=$PKG_PATH" >> "$GITHUB_ENV" - - conda build conda-recipe --override-channels --channel conda-forge - - test -f "$PKG_PATH" - echo "Built: $PKG_PATH" - - - name: Upload to Anaconda - shell: bash -l {0} - env: - ANACONDA_API_TOKEN: ${{ secrets.ANACONDA_TOKEN }} - run: | - anaconda -t "$ANACONDA_API_TOKEN" upload "$PKG_PATH" --user CCPBioSim --force diff --git a/.github/workflows/update-copyright-years-in-license-file.yaml b/.github/workflows/update-copyright-years-in-license-file.yaml new file mode 100644 index 0000000..3ba8a34 --- /dev/null +++ b/.github/workflows/update-copyright-years-in-license-file.yaml @@ -0,0 +1,16 @@ +name: Update copyright year(s) in license file + +on: + schedule: + - cron: '0 3 1 1 *' # 03:00 AM on January 1 + workflow_dispatch: + +permissions: + contents: read + +jobs: + update-license-year: + uses: CCPBioSim/.github/.github/workflows/maintenance-license-year.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: write + pull-requests: write diff --git a/.github/workflows/weekly-docs.yaml b/.github/workflows/weekly-docs.yaml new file mode 100644 index 0000000..32ac63f --- /dev/null +++ b/.github/workflows/weekly-docs.yaml @@ -0,0 +1,24 @@ +name: Weekly Docs Build + +on: + schedule: + - cron: '0 8 * * 1' + workflow_dispatch: + +permissions: + contents: read + +concurrency: + group: weekly-docs-${{ github.ref }} + cancel-in-progress: true + +jobs: + docs: + name: Docs build + uses: CCPBioSim/.github/.github/workflows/ci-docs.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: read + with: + python-versions: '["3.13"]' + os: '["ubuntu-26.04", "macos-26"]' + build-dir: build diff --git a/.github/workflows/weekly-tests.yaml b/.github/workflows/weekly-tests.yaml new file mode 100644 index 0000000..24ea750 --- /dev/null +++ b/.github/workflows/weekly-tests.yaml @@ -0,0 +1,33 @@ +name: Weekly Unit Tests + +on: + schedule: + - cron: '0 8 * * 1' + workflow_dispatch: + +permissions: + contents: read + +concurrency: + group: weekly-tests-${{ github.ref }} + cancel-in-progress: true + +jobs: + unit: + name: Unit + uses: CCPBioSim/.github/.github/workflows/ci-test-conda.yaml@2bf63339f2515ec54777e952e8db3d8674d4e5ca # 0.14.0 + permissions: + contents: read + checks: write + statuses: write + with: + package: prepmd + python-versions: '["3.13"]' + coverage-python-version: "3.13" + os: '["ubuntu-26.04", "macos-26"]' + environment-file: environment.yaml + activate-environment: prep + channels: salilab,conda-forge + test-path: test + pytest-args: -v + timeout-minutes: 360 diff --git a/.gitignore b/.gitignore index e80b8a1..c5bccaf 100644 --- a/.gitignore +++ b/.gitignore @@ -182,9 +182,9 @@ cython_debug/ .abstra/ # Visual Studio Code -# Visual Studio Code specific template is maintained in a separate VisualStudioCode.gitignore +# Visual Studio Code specific template is maintained in a separate VisualStudioCode.gitignore # that can be found at https://github.com/github/gitignore/blob/main/Global/VisualStudioCode.gitignore -# and can be added to the global gitignore or merged into this file. However, if you prefer, +# and can be added to the global gitignore or merged into this file. However, if you prefer, # you could uncomment the following to ignore the entire vscode folder # .vscode/ diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml new file mode 100644 index 0000000..672e0d6 --- /dev/null +++ b/.pre-commit-config.yaml @@ -0,0 +1,39 @@ +exclude: ^test/test_data/ + +repos: + - repo: https://github.com/astral-sh/ruff-pre-commit + rev: v0.17.2 + hooks: + - id: ruff-check + args: [--fix] + - id: ruff-format + + - repo: https://github.com/pre-commit/pre-commit-hooks + rev: v6.0.0 + hooks: + - id: check-added-large-files + - id: check-merge-conflict + - id: check-yaml + exclude: ^conda-recipe/meta\.yaml$ + - id: check-toml + - id: check-case-conflict + - id: check-ast + - id: end-of-file-fixer + - id: trailing-whitespace + + - repo: https://github.com/rstcheck/rstcheck + rev: v6.3.0 + hooks: + - id: rstcheck + name: rstcheck (docstrings only) + files: \.py$ + args: ["--report-level", "warning"] + + - repo: local + hooks: + - id: sphinx-docs + name: sphinx-build (pre-push, warnings as errors) + entry: bash -lc 'make -C docs clean && make -C docs html SPHINXOPTS="-W --keep-going"' + language: system + pass_filenames: false + stages: [pre-push] diff --git a/.readthedocs.yaml b/.readthedocs.yaml index f99bb46..6bc548f 100644 --- a/.readthedocs.yaml +++ b/.readthedocs.yaml @@ -6,7 +6,7 @@ version: 2 build: os: ubuntu-24.04 tools: - python: "3.14" + python: "3.13" sphinx: configuration: docs/source/conf.py @@ -17,4 +17,3 @@ python: path: . extra_requirements: - docs - - myst-parser diff --git a/README.md b/README.md index ead6da3..5a5d24e 100644 --- a/README.md +++ b/README.md @@ -3,8 +3,8 @@ prepmd | Category | Badges | |----------------|--------| -| **Build** | [![prepmd CI](https://github.com/CCPBioSim/prepmd/actions/workflows/python-app.yml/badge.svg)](https://github.com/CCPBioSim/prepmd/actions/workflows/python-app.yml) | -| **Documentation** |[![Docs - Status](https://app.readthedocs.org/projects/prepmd/badge/?version=latest)](https://prepmd.readthedocs.io/en/latest/?badge=latest) | +| **Build** | [![PR Checks](https://github.com/CCPBioSim/prepmd/actions/workflows/pr.yaml/badge.svg)](https://github.com/CCPBioSim/prepmd/actions/workflows/pr.yaml) [![Weekly Tests](https://github.com/CCPBioSim/prepmd/actions/workflows/weekly-tests.yaml/badge.svg)](https://github.com/CCPBioSim/prepmd/actions/workflows/weekly-tests.yaml) | +| **Documentation** | [![Weekly Docs](https://github.com/CCPBioSim/prepmd/actions/workflows/weekly-docs.yaml/badge.svg)](https://github.com/CCPBioSim/prepmd/actions/workflows/weekly-docs.yaml) [![Docs - Status](https://app.readthedocs.org/projects/prepmd/badge/?version=latest)](https://prepmd.readthedocs.io/en/latest/?badge=latest) | | **Citation** | [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.21536438.svg)](https://doi.org/10.5281/zenodo.21536438) | | **Anaconda** | [![Anaconda.org](https://anaconda.org/CCPBioSim/prepmd/badges/version.svg)](https://anaconda.org/CCPBioSim/prepmd/) [![Last Updated](https://anaconda.org/CCPBioSim/prepmd/badges/latest_release_date.svg)](https://anaconda.org/CCPBioSim/prepmd) [![Platforms](https://anaconda.org/CCPBioSim/prepmd/badges/platforms.svg)](https://anaconda.org/CCPBioSim/prepmd) [![License](https://anaconda.org/CCPBioSim/prepmd/badges/license.svg)](https://anaconda.org/CCPBioSim/prepmd) [![Downloads](https://anaconda.org/CCPBioSim/prepmd/badges/downloads.svg)](https://anaconda.org/CCPBioSim/prepmd)| | **Quality** | [![Coverage Status](https://coveralls.io/repos/github/CCPBioSim/prepmd/badge.svg?branch=main)](https://coveralls.io/github/CCPBioSim/prepmd?branch=main) | diff --git a/conda-recipe/meta.yaml b/conda-recipe/meta.yaml index 59fc058..4e00b7c 100644 --- a/conda-recipe/meta.yaml +++ b/conda-recipe/meta.yaml @@ -59,4 +59,4 @@ about: and the OpenFF Toolkit. dev_url: https://github.com/CCPBioSim/prepmd doc_url: https://prepmd.readthedocs.io/en/latest/ - doc_source_url: https://github.com/CCPBioSim/prepmd/blob/main/README.md \ No newline at end of file + doc_source_url: https://github.com/CCPBioSim/prepmd/blob/main/README.md diff --git a/docs/source/Developer-Notes.md b/docs/source/Developer-Notes.md index 28358b1..771dc6d 100644 --- a/docs/source/Developer-Notes.md +++ b/docs/source/Developer-Notes.md @@ -18,7 +18,7 @@ Dependencies Manual install -------------- * Install [Conda](https://github.com/conda-forge/miniforge?tab=readme-ov-file#install) (if you don't already have it) -* Clone this repo and enter the folder: `git clone https://github.com/CCPBioSim/prepmd.git && cd prepmd` +* Clone this repo and enter the folder: `git clone https://github.com/CCPBioSim/prepmd.git && cd prepmd` * Run `conda env create --name prepmd --file environment.yaml && conda activate prepmd && pip install -e .` * For the MODELLER part of the workflow to work, you need a [modeller license key](https://salilab.org/modeller/registration.html) and add it to modeller's config.py file. If you use conda, the key will be in `envs/prep/lib/modeller-10.7/modlib/modeller/config.py` relative to the path where conda is installed. * After installing, run `pytest` to run tests. @@ -50,4 +50,4 @@ Style * All user-facing code (e.g. CLI) should have rigorous input validation and descriptive error messages * Pure functions, minimal coupling and imperative style preferred. Code is mostly WET to reduce coupling where possible, with abstractions only being made where it's obviously necessary to do so. Functions are allowed to be long to reduce indirection and avoid hiding complexity. In particular, longer functions aren't broken up in a way that would require shorter functions to accept objects with a specific state as arguments. * prep.py and run.py are both stateful in evil ways; MODELLER writes tons of junk to the working directory, so prep.py changes the working directory to a temporary folder and then copies output files to a path that the user has specified. This necessitates a bit of global state, and also juggling different paths and accounting for relative/absolute paths provided by the user. -* run.py is even more stateful - OpenMM objects are pretty stateful, and runmd can re-initialise openmm system objects, keep and restore old co-ordinates (if the simulation becomes numerically unstable) and also recursively call itself to try and fix metadynamics problems. This looks ugly as sin, but it's worth it for those features. \ No newline at end of file +* run.py is even more stateful - OpenMM objects are pretty stateful, and runmd can re-initialise openmm system objects, keep and restore old co-ordinates (if the simulation becomes numerically unstable) and also recursively call itself to try and fix metadynamics problems. This looks ugly as sin, but it's worth it for those features. diff --git a/docs/source/Installation.md b/docs/source/Installation.md index 919e6fe..bd4b8f2 100644 --- a/docs/source/Installation.md +++ b/docs/source/Installation.md @@ -4,4 +4,4 @@ Installation * Install [Conda](https://github.com/conda-forge/miniforge?tab=readme-ov-file#install) (if you don't already have it). If you have an existing conda install, make sure it can install packages from conda-forge. * Recommended: create a new virtual environment: `conda env create --name prepmd && conda activate prepmd` * Install prepmd from the CCPBioSim conda channel: `conda install -c ccpbiosim prepmd` -* Optional: add your [modeller license key](https://salilab.org/modeller/registration.html) by running `prep-license `. If you don't have a MODELLER license key, prepmd will instead use pdbfixer's loop builder. \ No newline at end of file +* Optional: add your [modeller license key](https://salilab.org/modeller/registration.html) by running `prep-license `. If you don't have a MODELLER license key, prepmd will instead use pdbfixer's loop builder. diff --git a/docs/source/Prep-Structure.md b/docs/source/Prep-Structure.md index 46a48f5..8b0b360 100644 --- a/docs/source/Prep-Structure.md +++ b/docs/source/Prep-Structure.md @@ -15,7 +15,7 @@ Steps in the `prepmd` workflow: * The final, mimimised structure file will be written out. Note: if ligands are present, the non-minimised structure will be written instead - this is to allow the user to choose which ligand files to include in their final structure, which can be minimised using `runmd`. ## prepmd command-line reference -* Use `prepmd --help` for a full list of parameters. +* Use `prepmd --help` for a full list of parameters. ## Worked examples diff --git a/docs/source/Python-API.md b/docs/source/Python-API.md index 4c18ba6..16dea2b 100644 --- a/docs/source/Python-API.md +++ b/docs/source/Python-API.md @@ -5,11 +5,14 @@ Python API ```python from prepmd.add_modeller_license import setup_license + setup_license("") # only needs to be done once from prepmd.prep import prep + prep("6xov", "6xov.cif", "working_dir") from prepmd.run import run + run("6xov.cif", traj_out="traj.xtc") ``` diff --git a/docs/source/Quickstart.md b/docs/source/Quickstart.md index 8f3fb18..4d76db4 100644 --- a/docs/source/Quickstart.md +++ b/docs/source/Quickstart.md @@ -3,4 +3,4 @@ Quickstart `prepmd 6xov 6xov_processed.cif` will download the structure for PDB entry `6xov`, process it and write it to `6xov_processed.cif`. -`runmd 6xov_processed.cif --traj_out traj.xtc --md_steps 5000` will minimise and run a simulation of 6xov_processed.cif writing a trajectory to `traj_out.xtc`, for 5000 steps. By default, `runmd` uses a minimal set of simulation parameters, which aren't likely to be accurate - check the `runmd` section of this documentation for more options. \ No newline at end of file +`runmd 6xov_processed.cif --traj_out traj.xtc --md_steps 5000` will minimise and run a simulation of 6xov_processed.cif writing a trajectory to `traj_out.xtc`, for 5000 steps. By default, `runmd` uses a minimal set of simulation parameters, which aren't likely to be accurate - check the `runmd` section of this documentation for more options. diff --git a/docs/source/Run-with-runmd.md b/docs/source/Run-with-runmd.md index c66adca..a542c79 100644 --- a/docs/source/Run-with-runmd.md +++ b/docs/source/Run-with-runmd.md @@ -10,7 +10,7 @@ Running MD simulations with runmd * If the run is a metadynamics run, and the metadynamics collective variables aren't minimised, the simulation will restart. ## runmd command-line reference -* Use `runmd --help` for a full list of parameters. +* Use `runmd --help` for a full list of parameters. ## Example usage diff --git a/docs/source/Whats-next.md b/docs/source/Whats-next.md index e3638e1..4569b4b 100644 --- a/docs/source/Whats-next.md +++ b/docs/source/Whats-next.md @@ -1,4 +1,4 @@ What next? ========== * Though you can run simple MD simulations, minimisations and validation with `prepmd`, for more in-depth MD we recommend using software such as GROMACS, AMBER, NAMD and OpenMM. -* If you're looking to generate an atomistic structure file that matches your EM map as closely as possible, you can use a flexible fitting tool such as [TEMPy-ReFF](https://gitlab.com/topf-lab/tempy-reff). \ No newline at end of file +* If you're looking to generate an atomistic structure file that matches your EM map as closely as possible, you can use a flexible fitting tool such as [TEMPy-ReFF](https://gitlab.com/topf-lab/tempy-reff). diff --git a/docs/source/_static/logos/prepmd-logo-black-text.svg b/docs/source/_static/logos/prepmd-logo-black-text.svg index e496d91..9ac883b 100644 --- a/docs/source/_static/logos/prepmd-logo-black-text.svg +++ b/docs/source/_static/logos/prepmd-logo-black-text.svg @@ -1 +1 @@ -prepmd \ No newline at end of file +prepmd diff --git a/docs/source/_static/logos/prepmd-logo-black.svg b/docs/source/_static/logos/prepmd-logo-black.svg index 269ab49..02cba08 100644 --- a/docs/source/_static/logos/prepmd-logo-black.svg +++ b/docs/source/_static/logos/prepmd-logo-black.svg @@ -1 +1 @@ - \ No newline at end of file + diff --git a/docs/source/_static/logos/prepmd-logo-white-text.svg b/docs/source/_static/logos/prepmd-logo-white-text.svg index c7e3d47..fe877ad 100644 --- a/docs/source/_static/logos/prepmd-logo-white-text.svg +++ b/docs/source/_static/logos/prepmd-logo-white-text.svg @@ -1 +1 @@ -prepmd \ No newline at end of file +prepmd diff --git a/docs/source/_static/logos/prepmd-logo-white.svg b/docs/source/_static/logos/prepmd-logo-white.svg index 59a559e..6b2af39 100644 --- a/docs/source/_static/logos/prepmd-logo-white.svg +++ b/docs/source/_static/logos/prepmd-logo-white.svg @@ -1 +1 @@ - \ No newline at end of file + diff --git a/docs/source/conf.py b/docs/source/conf.py index 3e74e26..b510a07 100644 --- a/docs/source/conf.py +++ b/docs/source/conf.py @@ -54,8 +54,8 @@ # -- Napoleon docstring settings -------------------------------------------- -napoleon_google_docstring = False -napoleon_numpy_docstring = True +napoleon_google_docstring = True +napoleon_numpy_docstring = False napoleon_use_param = False napoleon_use_ivar = True @@ -63,7 +63,6 @@ exclude_patterns = [] - # -- Options for HTML output ------------------------------------------------- # https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output diff --git a/prepmd/__init__.py b/prepmd/__init__.py index eee2baa..fb56e00 100644 --- a/prepmd/__init__.py +++ b/prepmd/__init__.py @@ -8,4 +8,4 @@ and the OpenFF Toolkit. """ -__version__ = "1.0.2" \ No newline at end of file +__version__ = "1.0.2" diff --git a/prepmd/add_modeller_license.py b/prepmd/add_modeller_license.py index d94f265..0f7aba1 100644 --- a/prepmd/add_modeller_license.py +++ b/prepmd/add_modeller_license.py @@ -1,21 +1,23 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Change the modeller license key """ import importlib -from pathlib import Path import sys +from pathlib import Path # note: don't do this # i have to do it like this to make the conda build work -HELP_MSG = "Usage: prep-license LICENSE-KEY (will replace current key)" +HELP_MSG = "Usage: prep-license LICENSE-KEY (will replace current key)" + def setup_license(key): modeller_init_path = Path(importlib.util.find_spec("modeller").origin) - modeller_lib = modeller_init_path.parent.parent.absolute() / "modeller" / "config.py" + modeller_lib = ( + modeller_init_path.parent.parent.absolute() / "modeller" / "config.py" + ) with open(modeller_lib) as file: contents = file.readlines() splitline = contents[1].split("'") @@ -24,6 +26,7 @@ def setup_license(key): with open(modeller_lib, "w") as file: file.writelines(contents) + def entry_point(): if len(sys.argv) != 2: print(HELP_MSG) @@ -35,6 +38,6 @@ def entry_point(): setup_license(key) print("Updated modeller license info.") -if __name__ == "__main__": - entry_point() \ No newline at end of file +if __name__ == "__main__": + entry_point() diff --git a/prepmd/align_together.py b/prepmd/align_together.py index e4afb99..98979d1 100644 --- a/prepmd/align_together.py +++ b/prepmd/align_together.py @@ -1,5 +1,4 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Trim the ends off two PDB files until their indices match. Note: this is a big pile of spaghetti, it turns out doing this in a robust way @@ -9,20 +8,23 @@ errors when you do for element in reversed(object). """ -from Bio.PDB.PDBParser import PDBParser -from Bio.PDB.MMCIFParser import MMCIFParser - -#from Bio.PDB import StructureAlignment +# from Bio.PDB import StructureAlignment import Bio.PDB from Bio.PDB.mmcifio import MMCIFIO +from Bio.PDB.MMCIFParser import MMCIFParser +from Bio.PDB.PDBParser import PDBParser + try: import modeller + modeller_found = True -except: +except Exception: modeller_found = False pass -from prepmd import get_residues, model import argparse + +from prepmd import get_residues, model + # from Bio import Align """ @@ -42,28 +44,30 @@ def align_biopython(s1, s2): return alignment """ + def align_modeller(s1, s2, code1="AAAA", code2="BBBB", outfile="alignment.out"): print("Getting missing residues from input file...") pdb1_res = get_residues.get_residues_pdb(s1, code1) pdb2_res = get_residues.get_residues_pdb(s2, code2) # pdb_fullseq = prepmd.get_residues.get_fullseq_pdb(inmodel, code) - with open(code1+".seq", "w") as file: + with open(code1 + ".seq", "w") as file: file.write("".join(pdb1_res)) - with open(code2+".seq", "w") as file: + with open(code2 + ".seq", "w") as file: file.write("".join(pdb2_res)) if not modeller_found: raise ImportError("Problem with MODELLER license key or installation") env = modeller.environ() print("Aligning sequences...") aln = modeller.Alignment(env) - aln.append(file=code1+".seq", align_codes=code1) - aln.append(file=code2+".seq", align_codes=code2) + aln.append(file=code1 + ".seq", align_codes=code1) + aln.append(file=code2 + ".seq", align_codes=code2) aln.salign() print("Succesfully aligned.", end="") aln.write(file=outfile) -# given an anlignment fasta, create an index of which sequences have gaps in which places +# given an anlignment fasta, create an index of which sequences have gaps in which +# places def get_gaps(alignment_fasta_path): alignment_fasta = model.fasta(alignment_fasta_path) sequence1 = "".join(alignment_fasta[0][1]) @@ -71,7 +75,7 @@ def get_gaps(alignment_fasta_path): sequence1 = sequence1.replace("*", "").replace("/", "") # no asterisk sequence2 = sequence2.replace("*", "").replace("/", "") # no asterisk gaps = "" - for res1, res2 in zip(sequence1, sequence2): + for res1, res2 in zip(sequence1, sequence2, strict=False): if res1 == "-" and res2 == "-": gaps += "G" if res1 == "-" and res2 != "-": @@ -113,11 +117,12 @@ def mark_missing_ends(gaps): # Find the innermost gap at each end at mark anything further than that gap for deletion def mark_innermost_gap(gaps): - gaps_half1 = gaps[:len(gaps)//2] - gaps_half2 = gaps[len(gaps)//2:] + gaps_half1 = gaps[: len(gaps) // 2] + gaps_half2 = gaps[len(gaps) // 2 :] delete_half1_seq1, delete_half1_seq2 = mark_missing_ends( - ''.join(reversed(gaps_half1))) + "".join(reversed(gaps_half1)) + ) delete_half1_seq1.reverse() delete_half1_seq2.reverse() @@ -129,13 +134,14 @@ def mark_innermost_gap(gaps): return missing_seq1, missing_seq2 + # get the exact number of residues at each end to be deleted def get_ends_deletion(missing_seq): delete_start = None delete_end = None prev_res = True - if missing_seq[0] != False: + if missing_seq[0]: for res in range(len(missing_seq)): if missing_seq[res] != prev_res: delete_start = res # indexing/ob1 error! @@ -144,13 +150,13 @@ def get_ends_deletion(missing_seq): delete_no_start = delete_start prev_res = True - if missing_seq[-1] != False: + if missing_seq[-1]: for res in reversed(range(len(missing_seq))): if missing_seq[res] != prev_res: - delete_end = res+1 # indexing/ob1 error! + delete_end = res + 1 # indexing/ob1 error! break prev_res = missing_seq[res] - if type(delete_end) == int: + if type(delete_end) is int: delete_no_end = len(missing_seq) - delete_end else: delete_no_end = None @@ -161,9 +167,9 @@ def remove_from_model(structure, number_start, number_end, verbose=True): # you need to build this index because the biopython objects don't let you # iterate through them backwards or do anything with the loop variables models = [] - for model in structure: + for structure_model in structure: chains = [] - for chain in model: + for chain in structure_model: residues = [] for residue in chain: residues.append(residue.id) @@ -173,8 +179,12 @@ def remove_from_model(structure, number_start, number_end, verbose=True): number_removed = 0 to_remove = [] if number_end: - for model_idx, model in zip(reversed(range(len(models))), reversed(models)): - for chain_idx, chain in zip(reversed(range(len(model))), reversed(model)): + for model_idx, model_chains in zip( + reversed(range(len(models))), reversed(models), strict=False + ): + for chain_idx, chain in zip( + reversed(range(len(model_chains))), reversed(model_chains), strict=False + ): for residue in reversed(chain): if number_removed < number_end: to_remove.append((model_idx, chain_idx, residue)) @@ -182,80 +192,85 @@ def remove_from_model(structure, number_start, number_end, verbose=True): if number_start: number_removed = 0 - for model_idx, model in zip(range(len(models)), models): - for chain_idx, chain in zip(range(len(model)), model): + for model_idx, model_chains in zip(range(len(models)), models, strict=False): + for chain_idx, chain in zip( + range(len(model_chains)), model_chains, strict=False + ): for residue in chain: if number_removed < number_start: to_remove.append((model_idx, chain_idx, residue)) number_removed += 1 - for model, chain, res in to_remove: - curr_chain = structure[model].child_list[chain] + for model_idx, chain_idx, res in to_remove: + curr_chain = structure[model_idx].child_list[chain_idx] curr_chain.detach_child(res) - print("Removed "+str(len(to_remove))+" residues.") + print("Removed " + str(len(to_remove)) + " residues.") return structure -def align_together(seq1path, seq2path, seq1out, seq2out, - code1="AAAA", code2="BBBB", verbose=True): + +def align_together( + seq1path, seq2path, seq1out, seq2out, code1="AAAA", code2="BBBB", verbose=True +): align_modeller(seq1path, seq2path, outfile="alignment.out") gaps = get_gaps("alignment.out") - + if verbose: print("GAPS:") print(gaps) print("") - + missing_seq1, missing_seq2 = mark_innermost_gap(gaps) - + if verbose: print("Missing in sequence 1:") print_mask(missing_seq1) print("") print("Missing in sequence 2:") print_mask(missing_seq2) - + seq2_delete_start, seq2_delete_end = get_ends_deletion(missing_seq1) seq1_delete_start, seq1_delete_end = get_ends_deletion(missing_seq2) - + if verbose: - print("Delete in sequence 1: "+str([seq1_delete_start,seq1_delete_end])) - print("Delete in sequence 2: "+str([seq2_delete_start,seq2_delete_end])) - + print("Delete in sequence 1: " + str([seq1_delete_start, seq1_delete_end])) + print("Delete in sequence 2: " + str([seq2_delete_start, seq2_delete_end])) + if ".pdb" in seq1path and ".pdb" in seq2path: parser = PDBParser(PERMISSIVE=1) - if ".mmcif" in seq1path+seq2path or ".mmCif" in seq1path+seq2path: + if ".mmcif" in seq1path + seq2path or ".mmCif" in seq1path + seq2path: parser = MMCIFParser(PERMISSIVE=1) seq1_structure = parser.get_structure(code1, seq1path) seq2_structure = parser.get_structure(code2, seq2path) - + if verbose: - print("SEQ1 residues: "+str(len(list(seq1_structure.get_residues())))) - print("SEQ2 residues: "+str(len(list(seq2_structure.get_residues())))) - + print("SEQ1 residues: " + str(len(list(seq1_structure.get_residues())))) + print("SEQ2 residues: " + str(len(list(seq2_structure.get_residues())))) + seq1_structure = remove_from_model( - seq1_structure, seq1_delete_start, seq1_delete_end, verbose=True) + seq1_structure, seq1_delete_start, seq1_delete_end, verbose=True + ) seq2_structure = remove_from_model( - seq2_structure, seq2_delete_start, seq2_delete_end, verbose=True) - + seq2_structure, seq2_delete_start, seq2_delete_end, verbose=True + ) + seq1_length = len(list(seq1_structure.get_residues())) seq2_length = len(list(seq2_structure.get_residues())) - + if verbose: - print("SEQ1 residues: "+str(seq1_length)) - print("SEQ2 residues: "+str(seq2_length)) - + print("SEQ1 residues: " + str(seq1_length)) + print("SEQ2 residues: " + str(seq2_length)) + if seq1_length != seq2_length: - raise ValueError("The modified structures still don't have the same" - " length.") - - if ".mmcif" in seq1path+seq2path or ".mmCif" in seq1path+seq2path: + raise ValueError("The modified structures still don't have the same length.") + + if ".mmcif" in seq1path + seq2path or ".mmCif" in seq1path + seq2path: io = MMCIFIO() if ".pdb" in seq1path and ".pdb" in seq2path: io = Bio.PDB.PDBIO() - + io.set_structure(seq1_structure) io.save(seq1out) io.set_structure(seq2_structure) @@ -265,28 +280,39 @@ def align_together(seq1path, seq2path, seq1out, seq2out, def entry_point(): "CLI entry point function. Uses sys.argv and argparse args object." - parser = argparse.ArgumentParser(prog="aligntogether", - description="Trim two structure files so " - "that they are aligned and are the same " - "length.") + parser = argparse.ArgumentParser( + prog="aligntogether", + description="Trim two structure files so " + "that they are aligned and are the same " + "length.", + ) parser.add_argument("struct1in", help="Path to first input structure") parser.add_argument("struct2in", help="Path to second input structure") parser.add_argument("struct1out", help="Path to first output structure") parser.add_argument("struct2out", help="Path to second output structure") parser.add_argument( - "-c1", "--code1", help="PDB code for file 1 (only used as a label)", - default="AAAA") + "-c1", + "--code1", + help="PDB code for file 1 (only used as a label)", + default="AAAA", + ) parser.add_argument( - "-c2", "--code2", help="PDB code for file 2 (only used as a label)", - default="BBBB") + "-c2", + "--code2", + help="PDB code for file 2 (only used as a label)", + default="BBBB", + ) args = parser.parse_args() - align_together(args.struct1in, - args.struct2in, - args.struct1out, - args.struct2out, - args.c1, - args.c2) + align_together( + args.struct1in, + args.struct2in, + args.struct1out, + args.struct2out, + args.c1, + args.c2, + ) + if __name__ == "__main__": - entry_point() \ No newline at end of file + entry_point() diff --git a/prepmd/analysis.py b/prepmd/analysis.py index 892a204..e689623 100644 --- a/prepmd/analysis.py +++ b/prepmd/analysis.py @@ -1,32 +1,33 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Created on Mon May 11 14:39:36 2026 @author: rob """ -import MDAnalysis as mda -from prepmd.ligand import load_universe import numpy as np +from prepmd.ligand import load_universe + + def get_selection_traj(universe, selection): selection_traj = [] - for ts in universe.trajectory: + for _ts in universe.trajectory: current_pos = selection.positions selection_traj.append(current_pos) return np.array(selection_traj) def dist(a, b): - return ( (b[0] - a[0])**2 + (b[1] - a[1])**2 + (b[2] - a[2])**2 )**0.5 + return ((b[0] - a[0]) ** 2 + (b[1] - a[1]) ** 2 + (b[2] - a[2]) ** 2) ** 0.5 -def get_ligand_centroid_traj(top, traj, ignore=["UNK"], output_trajs=True, - print_dists = True): +def get_ligand_centroid_traj( + top, traj, ignore=("UNK",), output_trajs=True, print_dists=True +): u = load_universe(top) u2 = u.load_new(traj) - ligand = u2.select_atoms('not protein and not water') + ligand = u2.select_atoms("not protein and not water") ligands = ligand.split("residue") ligand_centroids = {} for ligand in ligands: @@ -34,18 +35,18 @@ def get_ligand_centroid_traj(top, traj, ignore=["UNK"], output_trajs=True, if res in ignore: continue all_traj = get_selection_traj(u2, ligand) - if np.shape(all_traj)[1] > 1: # more than 1 atom + if np.shape(all_traj)[1] > 1: # more than 1 atom centroid = np.mean(all_traj, axis=2) else: - centroid = all_traj[0] # only 1 atom + centroid = all_traj[0] # only 1 atom res_no = 2 while True: if res not in ligand_centroids: ligand_centroids[res] = centroid break else: - if res+"_"+str(res_no) not in ligand_centroids: - ligand_centroids[res+"_"+str(res_no)] = centroid + if res + "_" + str(res_no) not in ligand_centroids: + ligand_centroids[res + "_" + str(res_no)] = centroid break else: res_no += 1 @@ -55,13 +56,11 @@ def get_ligand_centroid_traj(top, traj, ignore=["UNK"], output_trajs=True, if print_dists: for name, centroid in ligand_centroids.items(): ligand_dist = dist(centroid[0], centroid[-1]) - print(name+" moved "+str(ligand_dist)+"A") - - return ligand_centroids + print(name + " moved " + str(ligand_dist) + "A") + return ligand_centroids - -#top = "/home/rob/temp/testnewligand/min_3pdt.pdb" -#traj = "/home/rob/temp/testnewligand/out.xtc" -#centroids = get_ligand_centroid_traj(top, traj, output_trajs=False) +# top = "/home/rob/temp/testnewligand/min_3pdt.pdb" +# traj = "/home/rob/temp/testnewligand/out.xtc" +# centroids = get_ligand_centroid_traj(top, traj, output_trajs=False) diff --git a/prepmd/download.py b/prepmd/download.py index a8461df..88e6923 100644 --- a/prepmd/download.py +++ b/prepmd/download.py @@ -1,42 +1,52 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Download data from the PDB and UNIPROT """ import urllib.request from os.path import sep + import requests def get_em_map(emdb_id, directory): """ Download a structure from the EMDB. + Args: emdb_id: id of the em map to download, a string directory: directory to download the file into, a string + returns: path to the downloaded file. """ emdb_id = str(emdb_id).replace("EMD-", "").replace("emd-", "") - url = "https://ftp.ebi.ac.uk/pub/databases/emdb/structures/EMD-"+str(emdb_id)+"/map/emd_"+str(emdb_id)+".map.gz" - destination = directory+sep+str(emdb_id)+".map.gz" + url = ( + "https://ftp.ebi.ac.uk/pub/databases/emdb/structures/EMD-" + + str(emdb_id) + + "/map/emd_" + + str(emdb_id) + + ".map.gz" + ) + destination = directory + sep + str(emdb_id) + ".map.gz" try: urllib.request.urlretrieve(url, destination) except urllib.error.HTTPError as e: if e.code == 404: - msg = "EMDB entry with ID "+emdb_id+" not found." - raise IOError(msg) - return directory+sep+str(emdb_id)+".map.gz" + msg = "EMDB entry with ID " + emdb_id + " not found." + raise OSError(msg) from e + return directory + sep + str(emdb_id) + ".map.gz" def get_structure(pdb_id, directory, file_format="mmCif", redo=False): """ Download a structure from the PDB. + Args: pdb_id: id of the pdb to download, a string directory: directory to download the file into, a string file_format: mmCif or pdb, a string + returns: path to the downloaded file. """ @@ -47,43 +57,58 @@ def get_structure(pdb_id, directory, file_format="mmCif", redo=False): format_str = "pdb" try: if redo: - url = "https://pdb-redo.eu/db/"+pdb_id+"/"+pdb_id+"_final"+"."+format_str + url = ( + "https://pdb-redo.eu/db/" + + pdb_id + + "/" + + pdb_id + + "_final" + + "." + + format_str + ) print(url) else: - url = "https://files.rcsb.org/download/"+pdb_id+"."+format_str - destination = directory+sep+pdb_id+"."+format_str + url = "https://files.rcsb.org/download/" + pdb_id + "." + format_str + destination = directory + sep + pdb_id + "." + format_str urllib.request.urlretrieve(url, destination) except urllib.error.HTTPError as e: if e.code == 404 and redo: - msg = "PDB with ID "+pdb_id+" not found in PDB-REDO." - raise IOError(msg) + msg = "PDB with ID " + pdb_id + " not found in PDB-REDO." + raise OSError(msg) from e r = requests.get(url.replace(".pdb", ".cif")) if r.status_code == 200: - msg = "No PDB for "+pdb_id + \ - " exists (but an mmcif structure does). " - "Run with --fmt cif to use it." - raise IOError(msg) + msg = ( + "No PDB for " + pdb_id + " exists (but an mmcif structure does). " + "Run with --fmt cif to use it." + ) + raise OSError(msg) from e else: raise e - return directory+sep+pdb_id+"."+format_str + return directory + sep + pdb_id + "." + format_str -def get_uniprot_sequence(pdb_id, merge_sequence=True, write_to_file=None, - verbose=False): +def get_uniprot_sequence( + pdb_id, merge_sequence=True, write_to_file=None, verbose=False +): """ For a given pdb id, find the fasta sequence for all chains from UNIPROT. + Args: pdb_id: the id of the pdb, a string merge_sequence: whether to merge the sequences together into a single fasta file, a bool write_to_file: path of file to write to, a string verbose: a bool, whether to write debug info out + Returns: the fasta sequences as a dictionry keyed by pdb id, or a string with the fasta sequence. """ - query = '''{ - entries(entry_ids:["'''+pdb_id+'''"]){ + query = ( + '''{ + entries(entry_ids:["''' + + pdb_id + + """"]){ polymer_entities { rcsb_id rcsb_polymer_entity_container_identifiers { @@ -94,25 +119,34 @@ def get_uniprot_sequence(pdb_id, merge_sequence=True, write_to_file=None, } } } - }''' - data = requests.get('https://data.rcsb.org/graphql?query='+query) + }""" + ) + data = requests.get("https://data.rcsb.org/graphql?query=" + query) accessions = {} - for item in data.json()['data']['entries']: - entities = item['polymer_entities'] + for item in data.json()["data"]["entries"]: + entities = item["polymer_entities"] for entity in entities: - if entity['rcsb_polymer_entity_container_identifiers']['reference_sequence_identifiers']: - for identifier in entity['rcsb_polymer_entity_container_identifiers']['reference_sequence_identifiers']: - accessions[int(entity["rcsb_id"].split("_")[1]) - ] = identifier['database_accession'] + if entity["rcsb_polymer_entity_container_identifiers"][ + "reference_sequence_identifiers" + ]: + for identifier in entity["rcsb_polymer_entity_container_identifiers"][ + "reference_sequence_identifiers" + ]: + accessions[int(entity["rcsb_id"].split("_")[1])] = identifier[ + "database_accession" + ] else: accessions[int(entity["rcsb_id"].split("_")[1])] = None # rna? if verbose: - print("Accessions: "+str(accessions)) + print("Accessions: " + str(accessions)) sequences = {} for key, value in accessions.items(): if value: - url = "https://rest.uniprot.org/uniprotkb/stream?compressed=false&format=fasta&query=(accession:" + \ - value+")" + url = ( + "https://rest.uniprot.org/uniprotkb/stream?compressed=false&format=fasta&query=(accession:" + + value + + ")" + ) seq = requests.get(url) sequences[key] = seq.text else: @@ -122,16 +156,16 @@ def get_uniprot_sequence(pdb_id, merge_sequence=True, write_to_file=None, return sequences names = [] combined_sequences = [] - for key, sequence in sequences.items(): + for sequence in sequences.values(): if sequence is None: continue - names.append(sequence.split("\n")[0].replace( - "sp|", "").replace(">", "")) + names.append(sequence.split("\n")[0].replace("sp|", "").replace(">", "")) seq_noheader = sequence.split("\n")[1:] combined_sequences.append("".join(seq_noheader).replace("\n", "")) - fasta = ">P1;"+pdb_id+"\n" + \ - ";".join(names)+"\n"+"/".join(combined_sequences) + fasta = ( + ">P1;" + pdb_id + "\n" + ";".join(names) + "\n" + "/".join(combined_sequences) + ) if write_to_file: with open(write_to_file, "w") as file: file.write(fasta) diff --git a/prepmd/fix.py b/prepmd/fix.py index 2c3cc9f..6350262 100644 --- a/prepmd/fix.py +++ b/prepmd/fix.py @@ -1,11 +1,10 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Fix structures using PDBFixer """ -from pdbfixer import PDBFixer from openmm.app import PDBFile, PDBxFile +from pdbfixer import PDBFixer def fix_nonstandard_cif(cif): @@ -39,6 +38,7 @@ def fix_nonstandard_cif(cif): with open(cif, "w") as file: file.writelines(outtext) + # for some reason, i can't get pdbfixer to remove hetatm and unk residues? def remove_hetatms_unk(pdb, out): lines = [] @@ -58,6 +58,7 @@ def remove_hetatms_unk(pdb, out): def fix(pdb, out, remove_heterogens=True, add_hydrogens=True, ph=7.0): """ Fix a pdb file with pdbfixer. + Arguments: pdb: path to a pdb or mmcif file, a string out: output file path, a string @@ -78,18 +79,20 @@ def fix(pdb, out, remove_heterogens=True, add_hydrogens=True, ph=7.0): if add_hydrogens: fixer.addMissingHydrogens(ph) if ".cif" in pdb or ".mmcif" in pdb: - PDBxFile.writeFile(fixer.topology, fixer.positions, open(out, 'w')) + PDBxFile.writeFile(fixer.topology, fixer.positions, open(out, "w")) else: - PDBFile.writeFile(fixer.topology, fixer.positions, open(out, 'w')) + PDBFile.writeFile(fixer.topology, fixer.positions, open(out, "w")) def restore_metadata_pdb(pdb, fixed_pdb): """ Copy metadata from one pdb file to another. Useful as the output of pdbfixer and other tools often doesn't contain REMARKS and such. + Args: pdb: original pdb file path containing metadata, a string fixed_pdb: new pdb file path, containing no metadata, a string + Returns: nothing, but updates the contents of fixed_pdb """ diff --git a/prepmd/get_residues.py b/prepmd/get_residues.py index 7258805..395da90 100644 --- a/prepmd/get_residues.py +++ b/prepmd/get_residues.py @@ -1,24 +1,26 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Read residue information from structure files """ + +from prepmd import ligand, util + NO_MODELLER = False try: - from modeller import * -except: + from modeller import Alignment, Environ, Model, log +except Exception: NO_MODELLER = True -from prepmd import util -from prepmd import ligand def get_residues_pdb(pdb, code, get_hetatms=False): """ Get the fasta sequence of residues in the ATOM entries of a PDB or mmCif file. + Args: pdb: path to pdb file, a string code: PDB code + Returns: the fasta sequence as a string """ @@ -31,8 +33,8 @@ def get_residues_pdb(pdb, code, get_hetatms=False): m = Model(e, file=pdb) aln = Alignment(e) aln.append_model(m, align_codes=code) - aln.write(file=code+'.seq') - with open(code+".seq") as file: + aln.write(file=code + ".seq") + with open(code + ".seq") as file: original_fasta = file.readlines() return original_fasta @@ -41,9 +43,11 @@ def get_fullseq_pdb(pdb, code, get_hetatms=False): """ Get the fasta sequence of residues in the SEQRES records of a PDB/mmCif file. + Args: pdb: path to pdb/mmcif file, a string code: PDB code + Returns: the fasta sequence as a string """ @@ -58,7 +62,7 @@ def get_fullseq_pdb(pdb, code, get_hetatms=False): if chain not in seqres.keys(): seqres[chain] = [] sequence = split[4:] - seqres[chain] += (sequence) + seqres[chain] += sequence if line.startswith("HET") and not hetatms_found and get_hetatms: hetatms_found = True @@ -75,7 +79,7 @@ def get_fullseq_pdb(pdb, code, get_hetatms=False): if chain != prev_chain: chains.append(chain) prev_chain = chain - + if seqres == {}: reading_seq = False with open(pdb) as file: @@ -115,7 +119,7 @@ def get_fullseq_pdb(pdb, code, get_hetatms=False): reading_seq = False if line.startswith("HET") and not hetatms_found and get_hetatms: hetatms_found = True - + mmcif_cross_check = False # mmcif cross-check # first, get order of chains, which for some reason is only contained in @@ -133,13 +137,18 @@ def get_fullseq_pdb(pdb, code, get_hetatms=False): chains.append(chain) prev_chain = chain - # 5O80 is five instances of the same chain and just gives a contextless sequence and five strand IDs - # 10ED is all over the place and lists a chain ID somewhere after the sequence but it's not consistent + # 5O80 is five instances of the same chain and just gives a contextless sequence + # and five strand IDs + # 10ED is all over the place and lists a chain ID somewhere after the sequence but + # it's not consistent # it does label them by entity ID, which is in the ATOM data - could use that? - # no, i can't, the sequence data is labeled and even formatted differently depending on whether the file has multiple entities or different instances of the same entity - # the only way this is all going to shake out is if I get the entity ID AND chain id of each residue, - # then get the entity ID OR chain id of - + # no, i can't, the sequence data is labeled and even formatted differently + # depending on whether the file has multiple entities or different instances of + # the same entity + # the only way this is all going to shake out is if I get the entity ID AND chain + # id of each residue, + # then get the entity ID OR chain id of + # this is valid for 2O8A reading_seq = False fasta_lines = [] @@ -156,11 +165,11 @@ def get_fullseq_pdb(pdb, code, get_hetatms=False): fasta_lines = "".join(fasta_lines).split(";") fasta_lines = list(filter(None, fasta_lines)) fasta_lines.pop(0) - - - - - + + + + + seq = [] seq_chains = [] for line in fasta_lines: @@ -172,7 +181,7 @@ def get_fullseq_pdb(pdb, code, get_hetatms=False): chain_dict = {} for chain, sequence in zip(seq_chains, seq): chain_dict[chain.strip()] = sequence - + fasta_joined = "" for chain in chains: #print(chain+": "+chain_dict[chain]) @@ -191,29 +200,30 @@ def get_fullseq_pdb(pdb, code, get_hetatms=False): # count hetatm residues and add an equivalent number of "." entries in # the fasta sequence universe = ligand.load_universe(pdb) - ligands = universe.select_atoms('not protein and not water') + ligands = universe.select_atoms("not protein and not water") num_hetatm_residues = len(ligands.split("residue")) last_key = sorted(seqres.keys())[-1] - seqres[last_key] += [num_hetatm_residues*"."] + seqres[last_key] += [num_hetatm_residues * "."] - if get_hetatms and not hetatms_found: - raise ValueError("Was told to retrieve hetatms from "+pdb+" but none " - "were found.") + raise ValueError( + "Was told to retrieve hetatms from " + pdb + " but none were found." + ) # convert seqres to fasta fastas = [] - for chain, reses in seqres.items(): + for reses in seqres.values(): fasta = "" for res in reses: if util.is_residue(res): fasta += util.three_to_one(res) fastas.append(fasta) fasta_joined = "/".join(fastas) - if fastas == []: - raise IOError("Couldn't get full sequence from contents of "+pdb+". " - "Does it contain a sequence?") + raise OSError( + "Couldn't get full sequence from contents of " + pdb + ". " + "Does it contain a sequence?" + ) chains = ":::::::::" - return ">P1;"+code+"_fill\n"+chains+"\n"+fasta_joined+"*" + return ">P1;" + code + "_fill\n" + chains + "\n" + fasta_joined + "*" diff --git a/prepmd/lib/__init__.py b/prepmd/lib/__init__.py index 580526c..1021186 100644 --- a/prepmd/lib/__init__.py +++ b/prepmd/lib/__init__.py @@ -1,10 +1,10 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Created on Tue Feb 17 13:53:45 2026 @author: rob """ -from . import icp -from . import mdaCIF +from . import icp, mdaCIF + +__all__ = ["icp", "mdaCIF"] diff --git a/prepmd/lib/icp.py b/prepmd/lib/icp.py index 28e8045..c65d5ef 100644 --- a/prepmd/lib/icp.py +++ b/prepmd/lib/icp.py @@ -8,18 +8,21 @@ import numpy as np from sklearn.neighbors import NearestNeighbors + def best_fit_transform(A, B): - ''' + """ Calculates the least-squares best-fit transform between corresponding 3D points A->B - Input: + + Args: A: Nx3 numpy array of corresponding 3D points B: Nx3 numpy array of corresponding 3D points + Returns: T: 4x4 homogeneous transformation matrix R: 3x3 rotation matrix t: 3x1 column vector for translation s: 3x1 column vector for scaling - ''' + """ assert len(A) == len(B) @@ -36,8 +39,8 @@ def best_fit_transform(A, B): # special reflection case if np.linalg.det(R) < 0: - Vt[2,:] *= -1 - R = np.dot(Vt.T, U.T) + Vt[2, :] *= -1 + R = np.dot(Vt.T, U.T) # translation t = centroid_B.T - np.dot(R, centroid_A.T) @@ -49,41 +52,47 @@ def best_fit_transform(A, B): return T, R, t + def nearest_neighbor(src, dst): - ''' + """ Find the nearest (Euclidean) neighbor in dst for each point in src - Input: + + Args: src: Nx3 array of points dst: Nx3 array of points - Output: + + Returns: distances: Euclidean distances of the nearest neighbor indices: dst indices of the nearest neighbor - ''' + """ neigh = NearestNeighbors(n_neighbors=1) neigh.fit(dst) distances, indices = neigh.kneighbors(src, return_distance=True) return distances.ravel(), indices.ravel() + def icp(A, B, init_pose=None, max_iterations=100, tolerance=1e-10): - ''' + """ The Iterative Closest Point method - Input: + + Args: A: Nx3 numpy array of source 3D points B: Nx3 numpy array of destination 3D point init_pose: 4x4 homogeneous transformation max_iterations: exit algorithm after max_iterations tolerance: convergence criteria - Output: + + Returns: T: final homogeneous transformation distances: Euclidean distances (errors) of the nearest neighbor - ''' + """ # make points homogeneous, copy them so as to maintain the originals - src = np.ones((4,A.shape[0])) - dst = np.ones((4,B.shape[0])) - src[0:3,:] = np.copy(A.T) - dst[0:3,:] = np.copy(B.T) + src = np.ones((4, A.shape[0])) + dst = np.ones((4, B.shape[0])) + src[0:3, :] = np.copy(A.T) + dst[0:3, :] = np.copy(B.T) # apply the initial pose estimation if init_pose is not None: @@ -92,25 +101,26 @@ def icp(A, B, init_pose=None, max_iterations=100, tolerance=1e-10): prev_error = 0 try: + for _i in range(max_iterations): + # find the nearest neighbours between the current source and destination + # points + distances, indices = nearest_neighbor(src[0:3, :].T, dst[0:3, :].T) - for i in range(max_iterations): - # find the nearest neighbours between the current source and destination points - distances, indices = nearest_neighbor(src[0:3,:].T, dst[0:3,:].T) - - # compute the transformation between the current source and nearest destination points - T, _, _ = best_fit_transform(src[0:3,:].T, dst[0:3,indices].T) + # compute the transformation between the current source and nearest + # destination points + T, _, _ = best_fit_transform(src[0:3, :].T, dst[0:3, indices].T) # update the current source src = T.dot(src) # check error mean_error = np.sum(distances) / distances.size - if abs(prev_error-mean_error) < tolerance: + if abs(prev_error - mean_error) < tolerance: break prev_error = mean_error # calculate final transformation - T, _, _ = best_fit_transform(A, src[0:3,:].T) + T, _, _ = best_fit_transform(A, src[0:3, :].T) return T, distances except ValueError as e: @@ -118,4 +128,4 @@ def icp(A, B, init_pose=None, max_iterations=100, tolerance=1e-10): return np.eye(4), 1, np.array([np.inf]) except np.linalg.linalg.LinAlgError as e: print(e) - return np.eye(4), 1, np.array([np.inf]) \ No newline at end of file + return np.eye(4), 1, np.array([np.inf]) diff --git a/prepmd/ligand.py b/prepmd/ligand.py index d154c29..75fbc91 100644 --- a/prepmd/ligand.py +++ b/prepmd/ligand.py @@ -1,5 +1,4 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Created on Fri Feb 27 16:49:54 2026 @@ -7,22 +6,22 @@ """ # setup -import MDAnalysis as mda -from rdkit import Chem import copy +import os -# run -from openff.toolkit import Molecule, Topology, ForceField +import MDAnalysis as mda +import numpy as np +from Bio.PDB import MMCIFIO, PDBParser +from MDAnalysis.converters.OpenMMParser import OpenMMTopologyParser +from MDAnalysis.coordinates.memory import MemoryReader from openff.interchange import Interchange from openff.interchange.components._packmol import UNIT_CUBE, pack_box + +# run +from openff.toolkit import ForceField, Molecule, Topology from openff.units import unit from openmm.app import pdbxfile -from MDAnalysis.converters.OpenMMParser import OpenMMTopologyParser -from MDAnalysis.coordinates.memory import MemoryReader -from openmm.app.modeller import Modeller -from Bio.PDB import PDBParser, MMCIFIO -import numpy as np -import os +from rdkit import Chem """ Okay, here's my full laundry list of complaints about rdkit @@ -52,9 +51,10 @@ def load_universe(pdb_path): file into mdanalysis, as mdanalysis does not natively support mmCif files, despite the fact that the pdb format has been deprecated for two years (!) at time of writing. - + Args: pdb_path: path to a pdb or cif file + Returns: an mdanalysis Universe object """ @@ -68,17 +68,20 @@ def load_universe(pdb_path): u = mda.Universe(pdb_path) return u + def pdb_back_to_mmcif(inpdb, outcif): - p = PDBParser() - struc = p.get_structure("", inpdb) - io = MMCIFIO() - io.set_structure(struc) - io.save(outcif) + p = PDBParser() + struc = p.get_structure("", inpdb) + io = MMCIFIO() + io.set_structure(struc) + io.save(outcif) + # setup ligand system + def get_size(atoms): - if type(atoms) != np.ndarray: + if type(atoms) is not np.ndarray: atommap = [] conf = atoms.GetConformer() for atom in atoms.GetAtoms(): @@ -89,76 +92,80 @@ def get_size(atoms): maxima = [] minima = [] for i in range(2): - maxima.append(max(atoms[:,i])) - minima.append(min(atoms[:,i])) + maxima.append(max(atoms[:, i])) + minima.append(min(atoms[:, i])) size = np.array(maxima) - np.array(minima) return size + def split_pdb_ligand(pdb_path, new_pdb_output=None, neutralise_radicals=True): """ Extract ligands from a structure file (pdb or mmcif) and write out a new, ligand/hetatm-free structure file and sdf files containing the ligands. - + Args: pdb_path - path to input pdb or mmcif file, a string new_pdb_output - filename for the new pdb file, a string. Otherwise the existing pdb file will be overwritten. neutralise_radicals - whether to remove free radicals. A bool. + Returns: a list of strings, where each string is the filename of an sdf file that has been written, which are named according to their residue. """ u = load_universe(pdb_path) - protein = u.select_atoms('protein') - ligand = u.select_atoms('not protein and not water') - if len(ligand) == None: + protein = u.select_atoms("protein") + ligand = u.select_atoms("not protein and not water") + if len(ligand) is None: print("No hetatms/ligands found.") return None else: print("Extracting hetatms/ligands...") ligands = ligand.split("residue") - + # aaaaa # this is a horrible hack because mdanalysis has some kind of bug with # the global state whereby if you run .convert_to more than once it writes # out junk data for every structure after the first ligand_files = [] for ligand in range(len(ligands)): - name = "h"+str(ligand)+".pdb" + name = "h" + str(ligand) + ".pdb" ligands[ligand].write(name) ligand_files.append(name) - - #ligands = [ligands[4]] - + + # ligands = [ligands[4]] + if not new_pdb_output: output_path = pdb_path else: output_path = new_pdb_output if ".cif" in pdb_path or ".mmcif" in pdb_path: - print("Warning: mdanalysis inexplicably does not support writing " - "mmcif files, so your structure will be written to a temporary " - "pdb file and then converted to cif. This process may not work " - "with all structures, make sure to so double-check the output!") - protein.write(output_path+".pdb") + print( + "Warning: mdanalysis inexplicably does not support writing " + "mmcif files, so your structure will be written to a temporary " + "pdb file and then converted to cif. This process may not work " + "with all structures, make sure to so double-check the output!" + ) + protein.write(output_path + ".pdb") p = PDBParser() - struc = p.get_structure("", output_path+".pdb") + struc = p.get_structure("", output_path + ".pdb") io = MMCIFIO() io.set_structure(struc) - io.save(output_path) # todo: replace w call to pdb_back_to_mmcif + io.save(output_path) # todo: replace w call to pdb_back_to_mmcif else: protein.write(output_path) ligand_filenames = [] ligands_written = {} - for ligand, ligand_file in zip(ligands, ligand_files): + for ligand, ligand_file in zip(ligands, ligand_files, strict=False): resname = "".join(list(set(ligand.resnames))) if resname in ligands_written.keys(): ligands_written[resname] += 1 - resname += "_"+str(ligands_written[resname]) + resname += "_" + str(ligands_written[resname]) else: ligands_written[resname] = 1 fname = save_ligand(ligand_file, resname) @@ -172,14 +179,15 @@ def TEMP_write(mol, fname): writer.close() -def save_ligand(ligand, resname, neutralise_radicals = True): - - if type(ligand) == str: +def save_ligand(ligand, resname, neutralise_radicals=True): + + if type(ligand) is str: u = load_universe(ligand) ligand = u.atoms - mol = ligand.convert_to('RDKIT', NoImplicit=False) - #TEMP_write(mol, resname+"_1.sdf") + mol = ligand.convert_to("RDKIT", NoImplicit=False) + + # TEMP_write(mol, resname+"_1.sdf") # need to create a coordmap or rdkit rotates and translates the ligand # for absolutely no reason. luckily getting the positions of atoms in # rdkit is very easy. just make a conformer object and then iterate @@ -195,8 +203,9 @@ def get_coord_atommap(mol): coordmap[idx] = pos atommap.append([idx, idx]) return coordmap, atommap + coordmap, atommap = get_coord_atommap(mol) - + # do not touch this or it WILL break # the molecule does not contain any hydrogens but in order to add # hydrogens you have to first run the removehs function, if you don't @@ -206,11 +215,10 @@ def get_coord_atommap(mol): mol = Chem.RemoveHs(mol) mol = Chem.AddHs(mol, addCoords=True) except Exception as e: - print("WARNING: Couldn't extract "+resname+" ("+str(e)+").") + print("WARNING: Couldn't extract " + resname + " (" + str(e) + ").") return "" - - #TEMP_write(mol, resname+"_2.sdf") + # TEMP_write(mol, resname+"_2.sdf") original_mol = copy.copy(mol) # remove radicals @@ -223,55 +231,62 @@ def get_coord_atommap(mol): a.SetNumRadicalElectrons(0) radicals_neutralised += rad - #TEMP_write(mol, resname+"_3.sdf") + # TEMP_write(mol, resname+"_3.sdf") mol = Chem.AddHs(mol, addCoords=True) - #TEMP_write(mol, resname+"_4.sdf") + # TEMP_write(mol, resname+"_4.sdf") # if you don't call embedmolecule here then the molecule will have # exactly one atom in completely the wrong place that will cause the # simulation to immediately detonate. however, if you call # embedmolecule without a coordmap it will rotate and translate the # ligand. it does this out of spite and hatred for molecular modellers - #coordmap, atommap = get_coord_atommap(mol) - + # coordmap, atommap = get_coord_atommap(mol) + unembedded_mol = copy.copy(mol) Chem.AllChem.EmbedMolecule(mol, coordMap=coordmap) - if np.mean(get_size(mol)) < np.mean(get_size(unembedded_mol))*0.75: - # mol = unembedded_mol - print("Warning: embedding "+resname+" with rdkit broke it.") - #TEMP_write(mol, resname+"_5.sdf") + if np.mean(get_size(mol)) < np.mean(get_size(unembedded_mol)) * 0.75: + # mol = unembedded_mol + print("Warning: embedding " + resname + " with rdkit broke it.") + # TEMP_write(mol, resname+"_5.sdf") # i have to pass this stupid atommap rdMolAlign otherwise it won't # align the molecules because there's 'no substructure match' despite # them being literally the SAME MOLECULE - rmsd = Chem.rdMolAlign.AlignMol(mol, original_mol, atomMap=atommap) - #print("rmsd: "+str(rmsd)) - #TEMP_write(mol, resname+"_6.sdf") + Chem.rdMolAlign.AlignMol(mol, original_mol, atomMap=atommap) + # print("rmsd: "+str(rmsd)) + # TEMP_write(mol, resname+"_6.sdf") Chem.AllChem.MMFFOptimizeMolecule(mol) - #Chem.AllChem.UFFOptimizeMolecule(mol) - - #TEMP_write(mol, resname+"_7.sdf") + # Chem.AllChem.UFFOptimizeMolecule(mol) + + # TEMP_write(mol, resname+"_7.sdf") - writer = Chem.SDWriter(resname+".sdf") + writer = Chem.SDWriter(resname + ".sdf") writer.write(mol) writer.close() print("---") - return resname+".sdf" - - -waters = { - "tip3p": "tip3p.offxml", - "tip4pew": "tip4p_ew.offxml" -} - - -def solvate_openff(interchange, protein, ligands, n_water=None, gutter_nm=1.5, - positive_ion_smile="[Na+]", negative_ion_smile="[Cl-]", - desired_charge_elementary=0, water_ff="tip3p", - solvent_smiles="O", ff="openff-2.2.0.offxml", - solvent_name="HOH"): + return resname + ".sdf" + + +waters = {"tip3p": "tip3p.offxml", "tip4pew": "tip4p_ew.offxml"} + + +def solvate_openff( + interchange, + protein, + ligands, + n_water=None, + gutter_nm=1.5, + positive_ion_smile="[Na+]", + negative_ion_smile="[Cl-]", + desired_charge_elementary=0, + water_ff="tip3p", + solvent_smiles="O", + ff="openff-2.2.0.offxml", + solvent_name="HOH", +): """ Solvate an openff system. + Args: interchange: an openff interchange containing the protein AND ligand. protein: protein loaded as an openff topology @@ -288,20 +303,21 @@ def solvate_openff(interchange, protein, ligands, n_water=None, gutter_nm=1.5, solvent_smiles: smiles string for the solvent molecule. ff: force field to use, A string. Can be any openff force field. solvent name: name of solvent molecule (for topology!) - a string + Returns: openff topology and an interchange for just the waters. """ - all_molecules = [protein]+ligands + all_molecules = [protein] + ligands - total_charge = round( - sum(interchange["Electrostatics"].charges.values()), 3) + total_charge = round(sum(interchange["Electrostatics"].charges.values()), 3) ligand_charge = 0 for ligand in ligands: ligand_charge += ligand.total_charge assert total_charge == protein.total_charge + ligand_charge, ( - f"Total charge of the system is {total_charge}, not {protein.total_charge + ligand.total_charge}" + f"Total charge of the system is {total_charge}, not " + f"{protein.total_charge + ligand.total_charge}" ) water = Molecule.from_smiles(solvent_smiles) @@ -313,8 +329,9 @@ def solvate_openff(interchange, protein, ligands, n_water=None, gutter_nm=1.5, elif total_charge <= desired_charge_elementary * unit.elementary_charge: ion = Molecule.from_smiles(positive_ion_smile) - num_ions = abs(int((total_charge/unit.elementary_charge) - - (desired_charge_elementary))) + num_ions = abs( + int((total_charge / unit.elementary_charge) - (desired_charge_elementary)) + ) ion.generate_conformers(n_conformers=1) @@ -327,7 +344,7 @@ def solvate_openff(interchange, protein, ligands, n_water=None, gutter_nm=1.5, box_vol_angstroms = np.array(box_vectors.to("angstrom").m.diagonal()) # one water per 30A^3, old JAWS default n_water = int(box_vol_angstroms.prod() / 70) - print("Adding "+str(n_water)+" water molecules.") + print("Adding " + str(n_water) + " water molecules.") docked_topology = Topology.from_molecules(all_molecules) @@ -348,22 +365,24 @@ def solvate_openff(interchange, protein, ligands, n_water=None, gutter_nm=1.5, def create_ligand_system( - protein_structure, # pdb finame - ligand_structures, # list of ligand filenames - n_water=None, # if None, calculate water automatically - water_ff="tip3p", # tip3p or tip4pew only supported - solvent_smiles="O", - gutter_nm=1, - positive_ion_smile="[Na+]", - negative_ion_smile="[Cl-]", - desired_charge_elementary=0, - ff="ff14sb_off_impropers_0.0.4.offxml", - ligand_ff="openff-2.2.0.offxml", - output_topology=None, - solvent_name="HOH", - ligands_names=None): + protein_structure, # pdb finame + ligand_structures, # list of ligand filenames + n_water=None, # if None, calculate water automatically + water_ff="tip3p", # tip3p or tip4pew only supported + solvent_smiles="O", + gutter_nm=1, + positive_ion_smile="[Na+]", + negative_ion_smile="[Cl-]", + desired_charge_elementary=0, + ff="ff14sb_off_impropers_0.0.4.offxml", + ligand_ff="openff-2.2.0.offxml", + output_topology=None, + solvent_name="HOH", + ligands_names=None, +): """ Set up an openmm system containing ligands using openff. + Args: protein_structure: structure file path, pdb or mmcif, a string ligand_structure: a list of strings where each string is a path to @@ -382,13 +401,14 @@ def create_ligand_system( solvent_name - solvent name to appear in the topology, a string ligands_names: names of ligands to be written to the topology. a list of strings + Returns: an openmm system, topology and positions object. """ print("Setting up ligand simulation...") if water_ff not in waters.keys(): - raise ValueError("Water must be one of: "+",".join(waters.keys())) + raise ValueError("Water must be one of: " + ",".join(waters.keys())) protein_with_crystal_water = Topology.from_pdb(protein_structure) protein = protein_with_crystal_water.molecule(0) @@ -399,9 +419,9 @@ def create_ligand_system( force_field=ff14sb, topology=protein.to_topology(), ) - + if not ligands_names: - ligands_names = [None]*len(ligand_structures) + ligands_names = [None] * len(ligand_structures) for ligand_fname in range(len(ligand_structures)): try: basename = os.path.basename(ligand_structures[ligand_fname]) @@ -409,40 +429,43 @@ def create_ligand_system( if "_" in resname: resname = resname.split("_")[0] ligands_names[ligand_fname] = resname - except: + except Exception: pass if None in ligands_names: - print("Warning: missing ligand name(s)! Will appear in topology " - "as 'UNK'.") + print("Warning: missing ligand name(s)! Will appear in topology as 'UNK'.") ligands = [] - for ligand, name in zip(ligand_structures, ligands_names): + for ligand, name in zip(ligand_structures, ligands_names, strict=False): current_ligand = Molecule.from_file(ligand) if name: - current_ligand.name = name # ADD THIS + current_ligand.name = name # ADD THIS for atom in current_ligand.atoms: atom.metadata["residue_name"] = name ligands.append(current_ligand) ligand_topology = Topology.from_molecules(ligands) - ligand_force_field = ForceField( - ligand_ff, waters[water_ff]) # includes tip3p + ligand_force_field = ForceField(ligand_ff, waters[water_ff]) # includes tip3p - ligand_intrcg = Interchange.from_smirnoff( - ligand_force_field, ligand_topology) + ligand_intrcg = Interchange.from_smirnoff(ligand_force_field, ligand_topology) interchange = protein_intrcg.combine(ligand_intrcg) print("Solvating...") - packed_topology, water_interchange = solvate_openff(interchange, protein, - ligands, n_water=n_water, - gutter_nm=gutter_nm, - positive_ion_smile=positive_ion_smile, negative_ion_smile=negative_ion_smile, - desired_charge_elementary=desired_charge_elementary, water_ff=water_ff, - solvent_smiles=solvent_smiles, - ff=ligand_ff, - solvent_name=solvent_name) + packed_topology, water_interchange = solvate_openff( + interchange, + protein, + ligands, + n_water=n_water, + gutter_nm=gutter_nm, + positive_ion_smile=positive_ion_smile, + negative_ion_smile=negative_ion_smile, + desired_charge_elementary=desired_charge_elementary, + water_ff=water_ff, + solvent_smiles=solvent_smiles, + ff=ligand_ff, + solvent_name=solvent_name, + ) solvated_interchange = interchange.combine(water_interchange) solvated_interchange.positions = packed_topology.get_positions() @@ -452,10 +475,12 @@ def create_ligand_system( solvated_interchange.to_pdb(output_topology) if output_topology: if ".mmcif" in output_topology or ".mmCif" in output_topology: - print("WARNING: OpenMM can't write mmCif topologies. Writing a pdb" - " and converting it back to mmCif. Sorry!") - solvated_interchange.to_pdb(output_topology+".pdb") - pdb_back_to_mmcif(output_topology+".pdb", output_topology) + print( + "WARNING: OpenMM can't write mmCif topologies. Writing a pdb" + " and converting it back to mmCif. Sorry!" + ) + solvated_interchange.to_pdb(output_topology + ".pdb") + pdb_back_to_mmcif(output_topology + ".pdb", output_topology) openmm_system = solvated_interchange.to_openmm() openmm_topology = solvated_interchange.to_openmm_topology() diff --git a/prepmd/metadynamics.py b/prepmd/metadynamics.py index 2577265..bc3fb05 100644 --- a/prepmd/metadynamics.py +++ b/prepmd/metadynamics.py @@ -1,23 +1,33 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Run metadynamics simulations - used in run.py """ -from openmm.unit import * -import numpy as np import MDAnalysis.analysis.rms as rms +import numpy as np +from openmm.unit import nanometers + class StopSimulation(Exception): - """ Exception raised in order to stop simulation """ + """Exception raised in order to stop simulation""" + pass + class RMSDIncrease(Exception): pass -class RMSDReporter(object): - def __init__(self, file, reportInterval, ref_positions, threshold_nm = 0.17, consecutive_frames=2): - self._out = open(file, 'w') + +class RMSDReporter: + def __init__( + self, + file, + reportInterval, + ref_positions, + threshold_nm=0.17, + consecutive_frames=2, + ): + self._out = open(file, "w") self._reportInterval = reportInterval self._ref_positions = ref_positions.value_in_unit(nanometers) self._threshold_nm = threshold_nm @@ -30,17 +40,19 @@ def __del__(self): self._out.close() def describeNextReport(self, simulation): - steps = self._reportInterval - simulation.currentStep%self._reportInterval - return {'steps': steps, 'periodic': None, 'include':['positions']} + steps = self._reportInterval - simulation.currentStep % self._reportInterval + return {"steps": steps, "periodic": None, "include": ["positions"]} def report(self, simulation, state): positions = state.getPositions(asNumpy=True).value_in_unit(nanometers) # i think this measurement doesn't align whereas rmsdforce does - #rmsd = np.sqrt(((positions - self._ref_positions)**2).sum(-1).mean()) - report_rmsd = rms.rmsd(positions, self._ref_positions, center=True, superposition=True) - print("rmsd = "+str(report_rmsd), end=" ") + # rmsd = np.sqrt(((positions - self._ref_positions)**2).sum(-1).mean()) + report_rmsd = rms.rmsd( + positions, self._ref_positions, center=True, superposition=True + ) + print("rmsd = " + str(report_rmsd), end=" ") self.rmsd_log.append(report_rmsd) - self._out.write(str(report_rmsd)+"\n") + self._out.write(str(report_rmsd) + "\n") if report_rmsd < self._threshold_nm: self.current_consecutive_frames += 1 else: @@ -48,23 +60,24 @@ def report(self, simulation, state): if self.current_consecutive_frames > self.max_consecutive_frames: print(". RMSD below threshold value, stopping...") raise StopSimulation() - if self.reports%5 == 0 and self.reports > 0: + if self.reports % 5 == 0 and self.reports > 0: if np.mean(np.gradient(self.rmsd_log)) > 0: raise RMSDIncrease self.reports += 1 + def get_representative_rmsd_frames(rmsd_log, num_frames): num_frames -= 2 frames = [0] change = 0 - for frame in range(len(rmsd_log)-1): - change += abs(rmsd_log[frame+1] - rmsd_log[frame]) + for frame in range(len(rmsd_log) - 1): + change += abs(rmsd_log[frame + 1] - rmsd_log[frame]) change_per_frame = change / num_frames change = 0 - for frame in range(len(rmsd_log)-1): - change += abs(rmsd_log[frame+1] - rmsd_log[frame]) + for frame in range(len(rmsd_log) - 1): + change += abs(rmsd_log[frame + 1] - rmsd_log[frame]) if change > change_per_frame: frames.append(frame) change -= change_per_frame - frames.append(len(rmsd_log)-1) - return frames \ No newline at end of file + frames.append(len(rmsd_log) - 1) + return frames diff --git a/prepmd/model.py b/prepmd/model.py index 02d2ff9..8e4a084 100644 --- a/prepmd/model.py +++ b/prepmd/model.py @@ -1,32 +1,35 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Calls to the MODELLER API """ -NO_MODELLER = False -try: - import modeller - from modeller.automodel import * -except: - NO_MODELLER = True -from Bio import Align import os import pathlib import shutil -from prepmd import get_residues -from prepmd import point_cloud import sys +from Bio import Align + +from prepmd import get_residues, point_cloud + +NO_MODELLER = False +try: + import modeller + from modeller.automodel import automodel +except Exception: + NO_MODELLER = True + placeholder = "sequence:::::::::" -def retitle_alignment(alignmentfile, title1, title2, - metadata1=placeholder, metadata2=placeholder): +def retitle_alignment( + alignmentfile, title1, title2, metadata1=placeholder, metadata2=placeholder +): """ Alignment files generated by BioPython aren't formatted in a way that modeller can easily read, so this function reformats them with correct titles and header info. + Args: alignmentfile: path to alignment file, a string title1: a string, title of the first sequence @@ -34,27 +37,29 @@ def retitle_alignment(alignmentfile, title1, title2, metadata1: a string, metadata for the first sequence metadata2: a string, metadata for the second sequence. The metadata strings must contain at least the substring ':::::::::'. + Returns: none, but overwrites alignmentfile """ titles = [title1, title2] metadata = [metadata1, metadata2] - with open(alignmentfile, "r") as file: + with open(alignmentfile) as file: alignment = file.readlines() new_alignment = [] for line_no in range(len(alignment)): modified_line = False if ">" in alignment[line_no]: - new_alignment.append(alignment[line_no].replace( - ">", ">"+"P1;"+titles.pop(0))) + new_alignment.append( + alignment[line_no].replace(">", ">" + "P1;" + titles.pop(0)) + ) modified_line = True curr_metadata = metadata.pop(0) if curr_metadata is not None: - new_alignment.append(curr_metadata+os.linesep) + new_alignment.append(curr_metadata + os.linesep) line = alignment[line_no] if ">" not in line and ":" not in line and ";" not in line: if not line.strip().endswith("*"): - new_alignment.append(line.strip()+"*"+os.linesep) + new_alignment.append(line.strip() + "*" + os.linesep) modified_line = True if modified_line is False: new_alignment.append(alignment[line_no]) @@ -65,9 +70,11 @@ def retitle_alignment(alignmentfile, title1, title2, def fasta(fasta_name, include_metadata=False): """ Read data in fasta sequence. + Args: fasta_name: filename of fasta file, a string include_metadata: bool, reads metadata if true + Returns: a list of lists, the first list is the sequence in the file, the second are the different blocks in the file: first the name/id, then @@ -80,7 +87,7 @@ def fasta(fasta_name, include_metadata=False): for line in file: if line[0] == ">": header = line.strip() - if fastadata != [] and fastadata != ['']: + if fastadata != [] and fastadata != [""]: out.append([header, fastadata]) fastadata = [] if ":" in line and ">" not in line: @@ -90,16 +97,16 @@ def fasta(fasta_name, include_metadata=False): continue if ":" not in line and ">" not in line: fastadata.append(line.strip()) - if fastadata != [] and fastadata != ['']: + if fastadata != [] and fastadata != [""]: out.append([header, "".join(fastadata), metadata]) if out == []: - raise IOError("Couldn't read FASTA data in "+fasta_name) + raise OSError("Couldn't read FASTA data in " + fasta_name) return out def validate_alignment(fastafile): fix = False - with open(fastafile, "r") as file: + with open(fastafile) as file: lines = file.readlines() comment = [] header = [] @@ -114,18 +121,18 @@ def validate_alignment(fastafile): else: sequence.append(line) for item in range(len(comment)): - #print(comment[item]+", "+str(">" in comment[item])) + # print(comment[item]+", "+str(">" in comment[item])) if ">" in comment[item]: comment[item] = comment[item].replace(">", "") fix = True - #print(comment) + # print(comment) if len(comment) != 1: fix = True - comment = [ "".join(comment).replace("\n", "")+"\n" ] + comment = ["".join(comment).replace("\n", "") + "\n"] if fix: with open(fastafile, "w") as file: - file.writelines(header+comment+sequence) - print("Fixed incorrect formatting of "+fastafile+".") + file.writelines(header + comment + sequence) + print("Fixed incorrect formatting of " + fastafile + ".") return @@ -133,10 +140,10 @@ def get_alignment_info(alignmentout): """ For a FASTA-formatted alignment file, get the number of residues filled, number of gaps, and largest gap. - + Args: alignmentout: path to an alignment file, a string - + Returns: total_resididues_filled - total number of residues to be added total_gaps_filled - number of gaps to be filled @@ -144,6 +151,7 @@ def get_alignment_info(alignmentout): filled_gaps - gaps max_gap - largest gap """ + def get_info(aln): total_missing = 0 total_gaps = 0 @@ -166,27 +174,35 @@ def get_info(aln): with open(alignmentout) as file: alignment = "".join(file.readlines()) aln_missing, aln_filled = alignment.split(">")[1:] - + # fix for hetatms if ".h." in aln_missing: aln_missing = aln_missing.replace(".h.", "") if "---" in aln_filled[-5:]: aln_filled = "".join(aln_filled.rsplit("---", 1)) - + missing_residues, missing_gaps, missing_lengths = get_info(aln_missing) filled_residues, filled_gaps, filled_lengths = get_info(aln_filled) try: max_gap = max(missing_lengths) except ValueError: - max_gap = 0 # no gaps + max_gap = 0 # no gaps total_residues_filled = missing_residues - filled_residues total_gaps_filled = missing_gaps - filled_gaps - if total_residues_filled < 0: - print("WARNING: filled structure contains fewer residues than the " - "original. This is normally a sign that something has gone " - "wrong.") - return total_residues_filled, total_gaps_filled, filled_residues, filled_gaps, max_gap + if total_residues_filled < 0: + print( + "WARNING: filled structure contains fewer residues than the " + "original. This is normally a sign that something has gone " + "wrong." + ) + return ( + total_residues_filled, + total_gaps_filled, + filled_residues, + filled_gaps, + max_gap, + ) def get_objective_functions(pdbs): @@ -195,9 +211,10 @@ def get_objective_functions(pdbs): function (which measures the quality of the model) and similarity (of the model sequence and the sequence used to fill the missing loops) for each pdb. - + Args: pdbs: a list of strings, paths to pdb files + Returns: scores, similarities, two dictionaries keyed by file path containing the scores and similarities. @@ -205,7 +222,7 @@ def get_objective_functions(pdbs): scores = {} similarities = {} for pdb in pdbs: - with open(pdb, "r") as file: + with open(pdb) as file: for line in file: # pdb if "OBJECTIVE FUNCTION" in line: @@ -224,51 +241,62 @@ def get_objective_functions(pdbs): return scores, similarities -def get_best_pdb(directory, exts=["pdb", "cif", "mmcif", "mmCif"], - em_map=None, em_contour_level=None): +def get_best_pdb( + directory, exts=("pdb", "cif", "mmcif", "mmCif"), em_map=None, em_contour_level=None +): """ For a directory, find all of the pdbs generated by modeller and select the one with the highest objective function (arbitrary metric used by modeller). + Args: directory: a string, the directory to scan ext: file extensions to check for (a list of strings) em_map: path to an EM density map file (a string). If this is set, the best PDB will be picked based on similarity to the map. em_contour_level: contour level for the EM map, a float. + Returns: path to the file with the highest objective function, a string """ - + pdbs = [] for ext in exts: - pdbs += list(pathlib.Path(directory).glob('*.'+ext)) + pdbs += list(pathlib.Path(directory).glob("*." + ext)) scores, similarities = get_objective_functions(pdbs) if em_map: em_scores = {} for pdb in pdbs: - em_scores[pdb] = point_cloud.score_pdb_map(pdb, em_map, - em_contour_level) + em_scores[pdb] = point_cloud.score_pdb_map(pdb, em_map, em_contour_level) if len(similarities) == 0: - raise ValueError("MODELLER could not optimise missing " - "loops. Check the output folder for more info.") + raise ValueError( + "MODELLER could not optimise missing " + "loops. Check the output folder for more info." + ) max_sim = max(similarities.values()) if max_sim >= 97: - print("Similiarity: "+str(max_sim)+"%.") + print("Similiarity: " + str(max_sim) + "%.") if max_sim < 97 and max_sim >= 90: - print("WARNING: best sequence similarity is "+str(max_sim)+"%. " - "Modelled loops could be wrong. Please check them!") + print( + "WARNING: best sequence similarity is " + str(max_sim) + "%. " + "Modelled loops could be wrong. Please check them!" + ) if max_sim < 90 and max_sim >= 80: - print("SERIOUS WARNING: best sequence similarity is only "+ - str(max_sim)+"%. Modelled loops may have the wrong sequence." ) + print( + "SERIOUS WARNING: best sequence similarity is only " + + str(max_sim) + + "%. Modelled loops may have the wrong sequence." + ) if max_sim < 80: - raise ValueError("Sequences do not match. Cannot reconstruct missing " - "loops. Please double-check your sequence data (by " - "default these are the SEQRES records from the input " - "structure.") + raise ValueError( + "Sequences do not match. Cannot reconstruct missing " + "loops. Please double-check your sequence data (by " + "default these are the SEQRES records from the input " + "structure." + ) if em_map: em_err = min(em_scores.values()) - print("EM map alignemnt error: "+str(round(em_err, 2))+"A") + print("EM map alignemnt error: " + str(round(em_err, 2)) + "A") if em_err > 15: raise ValueError("EM map and PDB are not the same structure.") return str(min(em_scores, key=em_scores.get)) @@ -277,11 +305,21 @@ def get_best_pdb(directory, exts=["pdb", "cif", "mmcif", "mmCif"], # note: the pdb file isn't a parameter, it must be called code.pdb -def fix_missing_residues(code, fastafile, alignmentout, inmodel, outmodel, - wdir, num_models=1, em_map=None, em_contour=None, - keep_hetatms = False): +def fix_missing_residues( + code, + fastafile, + alignmentout, + inmodel, + outmodel, + wdir, + num_models=1, + em_map=None, + em_contour=None, + keep_hetatms=False, +): """ For a given structure, fill in missing loops using modeller. + Args: code: the pdb id of the structure, a string fastafile: path to the fasta sequence, a string @@ -291,6 +329,7 @@ def fix_missing_residues(code, fastafile, alignmentout, inmodel, outmodel, outmodel: output structure file, a string wdir: working directory, a string, will be created if it doesn't exist num_models: how many models to generate, an int. + Returns: nothing, but writes out outmodel and wdir. """ @@ -307,7 +346,7 @@ def fix_missing_residues(code, fastafile, alignmentout, inmodel, outmodel, except shutil.SameFileError: pass - pdb_dirs = ['.', '../atom_files'] + pdb_dirs = [".", "../atom_files"] modeller.log.none() @@ -316,92 +355,101 @@ def fix_missing_residues(code, fastafile, alignmentout, inmodel, outmodel, model = modeller.Model(env, file=inmodel) aln = modeller.Alignment(env) aln.append_model(model, align_codes=code) - aln.write(file=code+'.seq') + aln.write(file=code + ".seq") print("Using provided fasta sequence to generate alignment.") - pdb_sequence = fasta(code+'.seq') + pdb_sequence = fasta(code + ".seq") original_sequence = fasta(fastafile) - aln_metadata = fasta(code+'.seq', include_metadata=True)[0][2] + aln_metadata = fasta(code + ".seq", include_metadata=True)[0][2] aligner = Align.PairwiseAligner() try: print("Aligning sequences...") - alignments = aligner.align( - pdb_sequence[0][1], original_sequence[0][1]) + alignments = aligner.align(pdb_sequence[0][1], original_sequence[0][1]) print("Succesfully aligned.") Align.write(alignments[-1], alignmentout, "fasta") # the file doesn't have names/labels for whatever reason - retitle_alignment(alignmentout, code, code + - "_fill", metadata1=aln_metadata) + retitle_alignment( + alignmentout, code, code + "_fill", metadata1=aln_metadata + ) except OverflowError: - print("Could not align sequences (normally this is because the " - "sequences in the PDB and the original protein are too " - "different)") + print( + "Could not align sequences (normally this is because the " + "sequences in the PDB and the original protein are too " + "different)" + ) fastafile = None print("Will try to use the sequence data in the input file...") if not fastafile: print("Getting missing residues from input file...") - pdb_res = get_residues.get_residues_pdb(inmodel, code, - get_hetatms = keep_hetatms) - pdb_fullseq = get_residues.get_fullseq_pdb(inmodel, code, - get_hetatms = keep_hetatms) - with open(code+".seq", "w") as file: + pdb_res = get_residues.get_residues_pdb(inmodel, code, get_hetatms=keep_hetatms) + pdb_fullseq = get_residues.get_fullseq_pdb( + inmodel, code, get_hetatms=keep_hetatms + ) + with open(code + ".seq", "w") as file: file.write("".join(pdb_res)) - with open(code+"_fill.seq", "w") as file: + with open(code + "_fill.seq", "w") as file: file.write(pdb_fullseq) env = modeller.environ() env.io.hetatm = True print("Aligning sequences...") aln = modeller.Alignment(env) - validate_alignment(code+".seq") - aln.append(file=code+".seq", align_codes=code) + validate_alignment(code + ".seq") + aln.append(file=code + ".seq", align_codes=code) try: - aln.append(file=code+"_fill.seq", align_codes=code+"_fill") - except OSError: - raise OSError("MODELLER could not read the structure file. If " - "you're using cif, try switching to PDB, as that " - "sometimes helps.") + aln.append(file=code + "_fill.seq", align_codes=code + "_fill") + except OSError as e: + raise OSError( + "MODELLER could not read the structure file. If " + "you're using cif, try switching to PDB, as that " + "sometimes helps." + ) from e aln.salign() print("Succesfully aligned.", end="") aln.write(file=alignmentout) - print(" Wrote alignment to "+alignmentout) - - # TODO: here, open the alignment file and verify that the header has 10 : in it - # also check that each > begins a line, and if there are more than two then get rid of any + print(" Wrote alignment to " + alignmentout) + + # TODO: here, open the alignment file and verify that the header has 10 : in + # it also check that each > begins a line, and if there are more than two then + # get rid of any env.io.atom_files_directory = pdb_dirs print("Modelling missing loops...") if num_models > 4: - print("(Creating "+str(num_models)+" models - might be slow!)") + print("(Creating " + str(num_models) + " models - might be slow!)") old_stdout = sys.stdout - f = open(os.devnull, 'w') + f = open(os.devnull, "w") sys.stdout = f - a = automodel(env, alnfile=alignmentout, - knowns=code, sequence=code+"_fill") + a = automodel(env, alnfile=alignmentout, knowns=code, sequence=code + "_fill") a.auto_align() if ".mmCif" in inmodel or ".cif" in inmodel or "Cif" in inmodel: a.set_output_model_format("MMCIF") if num_models > 1: - a.starting_model= 1 - a.ending_model = num_models - + a.starting_model = 1 + a.ending_model = num_models + a.make() sys.stdout = old_stdout - + if em_map: print("Ranking PDBs by similarity to EM map...") best_pdb = get_best_pdb(wdir, em_map=em_map, em_contour_level=em_contour) - + print("Finished modelling missing loops.") residues, gaps, remain_res, remain_gaps, max_gap = get_alignment_info(alignmentout) - print("Filled "+str(residues)+" residues in "+str(gaps)+" gaps.") + print("Filled " + str(residues) + " residues in " + str(gaps) + " gaps.") if remain_res > 0: - print("(still missing: "+str(remain_res) + - " residues in "+str(remain_gaps)+" gaps)") + print( + "(still missing: " + + str(remain_res) + + " residues in " + + str(remain_gaps) + + " gaps)" + ) if max_gap > 20: - print("Warning: big missing loop ("+str(max_gap)+" residues)") + print("Warning: big missing loop (" + str(max_gap) + " residues)") shutil.copy2(best_pdb, os.path.basename(outmodel)) diff --git a/prepmd/pdb2pqr.py b/prepmd/pdb2pqr.py index 353d3e9..d52eb3f 100644 --- a/prepmd/pdb2pqr.py +++ b/prepmd/pdb2pqr.py @@ -1,14 +1,15 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Created on Mon Feb 9 14:04:25 2026 @author: rob """ -import pdb2pqr import subprocess +import pdb2pqr + + def run_pdb2pqr(infile, outfile, ff="AMBER", ph=7): """ Run PDB2PQR. @@ -17,18 +18,28 @@ def run_pdb2pqr(infile, outfile, ff="AMBER", ph=7): infile - path to input file, a string outfile - path to output PQR file, a string ff - force field to use for calculations and residue naming + Returns: nothing, but writes 'outfile' """ - - command = [infile, outfile, "--titration-state-method=propka", - "--with-ph="+str(ph), "--ff="+ff, "--ffout="+ff, "--drop-water"] + + command = [ + infile, + outfile, + "--titration-state-method=propka", + "--with-ph=" + str(ph), + "--ff=" + ff, + "--ffout=" + ff, + "--drop-water", + ] if pdb2pqr.__version__ < "3.7.0": - print("WARNIONG: Old version of pdb2pqr detected, please install a " - "version after 3.7.0. If you installed everything via conda, please " - "install pdb2pqr via pip (pip install pdb2pqr) as the conda package is " - "abandoned.") - subprocess.run( ["pdb2pqr"]+command, capture_output=True, text=True) + print( + "WARNIONG: Old version of pdb2pqr detected, please install a " + "version after 3.7.0. If you installed everything via conda, please " + "install pdb2pqr via pip (pip install pdb2pqr) as the conda package is " + "abandoned." + ) + subprocess.run(["pdb2pqr"] + command, capture_output=True, text=True) return - - pdb2pqr.run_pdb2pqr(command) \ No newline at end of file + + pdb2pqr.run_pdb2pqr(command) diff --git a/prepmd/point_cloud.py b/prepmd/point_cloud.py index 577755f..116d0d4 100644 --- a/prepmd/point_cloud.py +++ b/prepmd/point_cloud.py @@ -1,23 +1,21 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Created on Fri Feb 13 15:41:18 2026 @author: rob """ -import os -from sklearn.neighbors import NearestNeighbors -import numpy as np -from prepmd.lib.icp import icp +import MDAnalysis as mda import mrcfile import numpy as np -import MDAnalysis as mda + +from prepmd.lib.icp import icp def to_point_cloud(mrcdata, voxel, contour_level): """ Convert an EM density map to a point cloud. + Args: mrcdata - ndarray of size resolutionXresolutionXresolution containg MRC data. @@ -26,6 +24,7 @@ def to_point_cloud(mrcdata, voxel, contour_level): dimensions. contour_level: density above which to add a point to the point cloud, a float. + Returns: point cloud as an ndarray with three columns (x, y, z) and a row for each point. @@ -33,33 +32,39 @@ def to_point_cloud(mrcdata, voxel, contour_level): point_cloud = [] for x in range(len(mrcdata)): for y in range(len(mrcdata[x])): - for z in range(len(mrcdata[x,y])): - if mrcdata[x,y,z] > contour_level: - point_cloud.append([x*voxel.x, y*voxel.y, z*voxel.z]) + for z in range(len(mrcdata[x, y])): + if mrcdata[x, y, z] > contour_level: + point_cloud.append([x * voxel.x, y * voxel.y, z * voxel.z]) return np.array(point_cloud) + def score_pdb_map(pdb, em_map, contour_level): """ For a given pdb and em_map, convert them to point clouds and score their similarity based on the error in an alignment between two point clouds. + Args: pdb - path to a pdb file, a string em_map - path to an an em map file for the same structure as that pdb, a string. contour_level: density above which to add a point to the point cloud, a float. + Returns: the error, a float """ with mrcfile.open(em_map) as mrc: - vsize = mrc.voxel_size point_cloud = to_point_cloud(mrc.data, mrc.voxel_size, 0.01) u = mda.Universe(pdb) atoms_pointcloud = u.atoms.positions - transformation, err, = icp(point_cloud, atoms_pointcloud) + ( + transformation, + err, + ) = icp(point_cloud, atoms_pointcloud) return np.mean(err) - -#score = score_pdb_map('/home/rob/temp/testtemp/6XOV_PGNMLW/6XOV_fill.B99990001.pdb', + + +# score = score_pdb_map('/home/rob/temp/testtemp/6XOV_PGNMLW/6XOV_fill.B99990001.pdb', # "/home/rob/temp/testtemp/emd_22281.map.gz", # 0.01) -#print(score) \ No newline at end of file +# print(score) diff --git a/prepmd/prep.py b/prepmd/prep.py index df381ee..3f8e72f 100644 --- a/prepmd/prep.py +++ b/prepmd/prep.py @@ -1,107 +1,155 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Prepare structures from the PDB for molecular dynamics. """ import argparse -import pathlib -import os +import copy import glob +import json +import os +import pathlib import random -import string import shutil -import copy -import json +import string from pathlib import Path -from prepmd import download -from prepmd import run -from prepmd import fix -from prepmd import model +from prepmd import download, fix, ligand, model, run from prepmd import pdb2pqr as pqr -from prepmd import ligand -parser = argparse.ArgumentParser(prog="prepmd", - description="Get structures from the PDB ready for " - "molecular dynamics runs") +parser = argparse.ArgumentParser( + prog="prepmd", + description="Get structures from the PDB ready for molecular dynamics runs", +) parser.add_argument("code", help="4 or 12-character PDB code") parser.add_argument("out", help="Output filename") parser.add_argument("-w", "--wdir", help="Working directory") -parser.add_argument("-d", "--directory", - help="Input directory (will automatically check for fasta" - " sequences and structure files here)") parser.add_argument( - "-f", "--fasta", help="Fasta sequence to use for filling loops") -parser.add_argument("-s", "--structure", - help="Input structure file (pdb or mmCif)") -parser.add_argument("-a", "--alignmentout", - help="Alignment output file from sequences aligned for loop filling", - default="alignment_out.fasta") -parser.add_argument("-fmt", "--dlformat", - help="Structure format to download (only used when no structure file is " - "provided)", - default=None) -parser.add_argument("-ph", "--ph", type=float, help="pH to target when adding" - " hydrogens", default=7.0) -parser.add_argument("-q", "--quiet", - help="Do not print debug info", action="store_true") -parser.add_argument("-e", "--fixstart", - help="Fix pdb at the start of the process, not the end", - action="store_true") -parser.add_argument("-dl", "--download", - help="Download the sequence from UNIPROT instead of using the pdb or" - " an external fasta file", action="store_true") -parser.add_argument("-m", "--leavemissing", - help="Don't restore missing atoms", action="store_true") -parser.add_argument("-p", "--pqr", - help="PQR output filename. If set, will output a PQR " - "as well as a pdb.") -parser.add_argument("-ff", "--pqrfmt", - help="Force field to use for creating the PQR file. Can be" - "AMBER or LAMMPS. Will only output a PQR file if this is " - "set.", default="AMBER") -parser.add_argument("-r", "--redo", help="Get PDB from PDB-REDO", - action="store_true") -parser.add_argument("-n", "--num", - help="Number of models to create with MODELLER. The best " - "model will be selected based on the MODELLER objective " - "function.", type=int, default=1) -parser.add_argument("-em", "--em_map", - help="If multiple models are being generated, the best " - "model will be selected based on agreement with this " - "EM density map instead of the objective function.") -parser.add_argument("-c", "--contour", - help="Contour level for the EM map.", type=float) -parser.add_argument("-he", "--hetatms", - help="Split out hetatms into sdf files", - action="store_true") -parser.add_argument("-nm", "--nomodeller", - help="Don't use MODELLER to fill missing residues. Will " - "use pdbfixer instead, which is less accurate.", - action="store_true") + "-d", + "--directory", + help="Input directory (will automatically check for fasta" + " sequences and structure files here)", +) +parser.add_argument("-f", "--fasta", help="Fasta sequence to use for filling loops") +parser.add_argument("-s", "--structure", help="Input structure file (pdb or mmCif)") +parser.add_argument( + "-a", + "--alignmentout", + help="Alignment output file from sequences aligned for loop filling", + default="alignment_out.fasta", +) +parser.add_argument( + "-fmt", + "--dlformat", + help="Structure format to download (only used when no structure file is provided)", + default=None, +) +parser.add_argument( + "-ph", "--ph", type=float, help="pH to target when adding hydrogens", default=7.0 +) +parser.add_argument( + "-q", "--quiet", help="Do not print debug info", action="store_true" +) +parser.add_argument( + "-e", + "--fixstart", + help="Fix pdb at the start of the process, not the end", + action="store_true", +) +parser.add_argument( + "-dl", + "--download", + help="Download the sequence from UNIPROT instead of using the pdb or" + " an external fasta file", + action="store_true", +) +parser.add_argument( + "-m", "--leavemissing", help="Don't restore missing atoms", action="store_true" +) +parser.add_argument( + "-p", + "--pqr", + help="PQR output filename. If set, will output a PQR as well as a pdb.", +) +parser.add_argument( + "-ff", + "--pqrfmt", + help="Force field to use for creating the PQR file. Can be" + "AMBER or LAMMPS. Will only output a PQR file if this is " + "set.", + default="AMBER", +) +parser.add_argument("-r", "--redo", help="Get PDB from PDB-REDO", action="store_true") +parser.add_argument( + "-n", + "--num", + help="Number of models to create with MODELLER. The best " + "model will be selected based on the MODELLER objective " + "function.", + type=int, + default=1, +) +parser.add_argument( + "-em", + "--em_map", + help="If multiple models are being generated, the best " + "model will be selected based on agreement with this " + "EM density map instead of the objective function.", +) +parser.add_argument("-c", "--contour", help="Contour level for the EM map.", type=float) +parser.add_argument( + "-he", "--hetatms", help="Split out hetatms into sdf files", action="store_true" +) +parser.add_argument( + "-nm", + "--nomodeller", + help="Don't use MODELLER to fill missing residues. Will " + "use pdbfixer instead, which is less accurate.", + action="store_true", +) def id_generator(size=6, chars=string.ascii_uppercase + string.digits): """ Generate a random 6-character ID. Used for naming scratch directories. + Args: size: length of string, an int chars; chars to use, a string + Returns: the string """ - return ''.join(random.choice(chars) for _ in range(size)) + return "".join(random.choice(chars) for _ in range(size)) -def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None, - alignmentout="alignment_out.fasta", download_format=None, - quiet=False, fix_after=True, download_sequence=False, - fix_missing_atoms=True, write_metadata="prepmeta.json", pqrff="AMBER", - pqr_out=None, redo=False, num_models=1, em_map=None, em_contour=None, - split_hetatms=False, no_modeller=False, ph=7.0): +def prep( + code, + outmodel, + workingdir, + folder=None, + fastafile=None, + inmodel=None, + alignmentout="alignment_out.fasta", + download_format=None, + quiet=False, + fix_after=True, + download_sequence=False, + fix_missing_atoms=True, + write_metadata="prepmeta.json", + pqrff="AMBER", + pqr_out=None, + redo=False, + num_models=1, + em_map=None, + em_contour=None, + split_hetatms=False, + no_modeller=False, + ph=7.0, +): """ Prepare a PDB/MMCIF structure file for simulation. + Args: code: pdb code, a string outmodel: path to output file, a string @@ -137,6 +185,7 @@ def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None, split_hetatms: if this is set to True, hetams will be written to sdf files, otherwise they will be removed. ph - ph to target when adding hydrogens, a float + Returns: nothing, but writes out a file to outmodel. """ @@ -144,24 +193,23 @@ def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None, # save all input parameters - will be written to a json file later locals_copy = copy.copy(locals()) locals_json = json.dumps(locals_copy) - + # create output file in the working dir, copy it later, save working dir final_dir = None if os.path.dirname(outmodel) != "": final_dir = os.path.dirname(outmodel) outmodel = os.path.basename(outmodel) - + # infer download format from output format if not download_format and not inmodel: if (".pdb") in outmodel: download_format = "pdb" print("No download format specified, downloading PDB.") - elif (".cif") in outmodel or ".mmcif" in outmodel or ( - ".pdbx") in outmodel: + elif (".cif") in outmodel or ".mmcif" in outmodel or (".pdbx") in outmodel: download_format = "mmCif" print("No download format specified, downloading mmCif.") - - #infer download format from input format + + # infer download format from input format if not download_format and inmodel: if (".pdb") in inmodel: download_format = "pdb" @@ -170,42 +218,57 @@ def prep(code, outmodel, workingdir, folder=None, fastafile=None, inmodel=None, # check that input/output are specified in the same format # i'm not against converting the files but it shouldn't happen implicitly - def in_string(substr, text): return text == None or substr in text.lower() - - if in_string(".pdb", inmodel) and in_string( - ".pdb", outmodel) and in_string("pdb", download_format): + def in_string(substr, text): + return text is None or substr in text.lower() + + if ( + in_string(".pdb", inmodel) + and in_string(".pdb", outmodel) + and in_string("pdb", download_format) + ): pass - elif in_string(".cif", inmodel) and in_string( - ".cif", outmodel) and in_string("cif", download_format): + elif ( + in_string(".cif", inmodel) + and in_string(".cif", outmodel) + and in_string("cif", download_format) + ): pass else: - raise IOError("Inputs and outputs are in different formats! Please " - "use only one format!)") + raise OSError( + "Inputs and outputs are in different formats! Please use only one format!)" + ) if pqr_out and download_format == "mmCif": - raise IOError("PQR generation requires the PDB format to be used. " - "Please specify an input and output PDB!") - + raise OSError( + "PQR generation requires the PDB format to be used. " + "Please specify an input and output PDB!" + ) + if not fix_after and ligand: - raise ValueError("Cannot extract ligands because they will be " - "removed by pdbfixer!") - + raise ValueError( + "Cannot extract ligands because they will be removed by pdbfixer!" + ) + # if we don't have modeller, pdbfixer will fill the residues if model.NO_MODELLER: no_modeller = True if no_modeller: modeller_shim = True - print("MODELLER was not found or has been disabled, so PDBFixer will " - "be used to fill missing residues") + print( + "MODELLER was not found or has been disabled, so PDBFixer will " + "be used to fill missing residues" + ) else: modeller_shim = False - + if modeller_shim and num_models > 1: - print("WARNING: multiple model generation is enabled but MODELLER is " - "not. Multiple models will NOT be generated. If you need your" - "output to match an EM map, you can still use another tool" - "such as tempy-reff.") - + print( + "WARNING: multiple model generation is enabled but MODELLER is " + "not. Multiple models will NOT be generated. If you need your" + "output to match an EM map, you can still use another tool" + "such as tempy-reff." + ) + if pqr_out: hydrogens = False else: @@ -217,13 +280,15 @@ def in_string(substr, text): return text == None or substr in text.lower() # copy input file to working dir if inmodel: - shutil.copyfile(inmodel, - workingdir+os.path.sep+code+"." - ""+download_format.replace("mmCif", "cif")) + shutil.copyfile( + inmodel, + workingdir + os.path.sep + code + "." + "" + download_format.replace("mmCif", "cif"), + ) # note: modeller requires the filename to be the same as the pdb code # so here we change the filename - inmodel = code+"."+download_format.replace("mmCif", "cif") - + inmodel = code + "." + download_format.replace("mmCif", "cif") + # dl em map file if em_map: print(em_map) @@ -231,8 +296,7 @@ def in_string(substr, text): return text == None or substr in text.lower() print("Downloading EM map...") em_map = download.get_em_map(em_map, workingdir) else: - shutil.copyfile(em_map, - workingdir+os.path.sep+os.path.basename(em_map)) + shutil.copyfile(em_map, workingdir + os.path.sep + os.path.basename(em_map)) # change working dir to temporary folder - note: this is an important bit # of state - prep doesn't ever leave the working dir, instead we make a @@ -243,32 +307,33 @@ def in_string(substr, text): return text == None or substr in text.lower() # check folder for strucutre/sequence files if inmodel is None and folder: - pdbs = glob.glob(folder+'/*.pdb') - cifs = glob.glob(folder+'/*.cif') + pdbs = glob.glob(folder + "/*.pdb") + cifs = glob.glob(folder + "/*.cif") models = pdbs + cifs if len(cifs) == 1: inmodel = cifs[0] if len(pdbs) == 1: inmodel = pdbs[0] - print("Found input model:"+str(inmodel)) + print("Found input model:" + str(inmodel)) if len(pdbs) > 1 or len(cifs) > 1: - raise IOError( - "Mulitple structure files in folder, please specify a" - " structure file") + raise OSError( + "Mulitple structure files in folder, please specify a structure file" + ) if len(models) == 0: - raise IOError("Couldn't find input structure file in folder") + raise OSError("Couldn't find input structure file in folder") # download structure - if code and inmodel == None and folder == None: + if code and inmodel is None and folder is None: print("Downloading structure file...") inmodel = download.get_structure( - code, workingdir, file_format=download_format, redo=redo) - print("Downloaded to "+inmodel) + code, workingdir, file_format=download_format, redo=redo + ) + print("Downloaded to " + inmodel) # fix if not fix_after: print("Fixing structure file...") - #fix.fix(inmodel, inmodel, fix_missing_atoms=fix_missing_atoms, + # fix.fix(inmodel, inmodel, fix_missing_atoms=fix_missing_atoms, # fix_missing_residues=modeller_shim) fix.fix(inmodel, inmodel, ph=ph, add_hydrogens=hydrogens) if modeller_shim and not split_hetatms: @@ -277,113 +342,137 @@ def in_string(substr, text): return text == None or substr in text.lower() # get fasta sequence fastas = None if fastafile is None and folder: - fastas = glob.glob(folder+'/*.fasta') + fastas = glob.glob(folder + "/*.fasta") if len(fastas) == 1: fastafile = fastas[0] if len(fastas) > 1: - raise IOError( + raise OSError( "Multiple sequence files found in folder, please specify a" - " sequence file") + " sequence file" + ) fastas = None - print("Found fasta file "+str(fastafile)+" in "+str(folder)) + print("Found fasta file " + str(fastafile) + " in " + str(folder)) if (fastafile is None and download_sequence) or ( - fastas is None and download_sequence): + fastas is None and download_sequence + ): try: print("Downloading sequence file(s)") sequences = download.get_uniprot_sequence(code) - fastafile = workingdir+"/"+code+".fasta" + fastafile = workingdir + "/" + code + ".fasta" with open(fastafile, "w") as file: file.write(sequences) except Exception as e: print(e) - raise IOError( - "No sequence files found in folder or specified, could not" - "download sequence") + raise OSError( + "No sequence files found in folder or specified, could not " + "download sequence" + ) from e # fix missing residues (call MODELLER) if not no_modeller: - model.fix_missing_residues(code, fastafile, alignmentout, inmodel, - outmodel, workingdir, num_models=num_models, - em_map=em_map, em_contour=em_contour, - keep_hetatms = split_hetatms) - + model.fix_missing_residues( + code, + fastafile, + alignmentout, + inmodel, + outmodel, + workingdir, + num_models=num_models, + em_map=em_map, + em_contour=em_contour, + keep_hetatms=split_hetatms, + ) + # strip heterogens, or split them out if the user wants no_sim_output = False if split_hetatms: ligands = ligand.split_pdb_ligand(inmodel) - print("Wrote hetatms/ligands to "+", ".join(ligands)) - no_sim_output = True # don't minimise the ligandless structure without + print("Wrote hetatms/ligands to " + ", ".join(ligands)) + no_sim_output = True # don't minimise the ligandless structure without # the ligands! if not split_hetatms: ligands = None - print("Note: all hetatms/ligands will be removed! Enable the 'split " - "hetatms/ligands' option to keep them.") - + print( + "Note: all hetatms/ligands will be removed! Enable the 'split " + "hetatms/ligands' option to keep them." + ) + if split_hetatms and not ligands: print("Warning: 'split hetatms' is enabled but no hetatms were found.") - + if fix_after: print("Fixing PDB...") - fix.fix(os.path.basename(outmodel), os.path.basename(outmodel), ph=ph, - add_hydrogens=hydrogens) + fix.fix( + os.path.basename(outmodel), + os.path.basename(outmodel), + ph=ph, + add_hydrogens=hydrogens, + ) if modeller_shim and not split_hetatms: - fix.remove_hetatms_unk(os.path.basename(outmodel), - os.path.basename(outmodel)) - + fix.remove_hetatms_unk( + os.path.basename(outmodel), os.path.basename(outmodel) + ) + if pqr_out: print("Creating PQR...") - pqr.run_pdb2pqr(outmodel, os.path.basename(pqr_out), ff=pqrff, ph=ph) + pqr.run_pdb2pqr(outmodel, os.path.basename(pqr_out), ff=pqrff, ph=ph) # if the PQR is written then don't output a minimised structure file - print("Simulating "+code) + print("Simulating " + code) if pqr_out: - run.run(os.path.basename(pqr_out), - minimised_structure_out=None, md_steps=None, - integrator_str="LangevinMiddleIntegrator", - pressure=None, test_sim_steps=250) - print("Note: minimised structure file has not been written out, as " - "this would replace the PQR.") + run.run( + os.path.basename(pqr_out), + minimised_structure_out=None, + md_steps=None, + integrator_str="LangevinMiddleIntegrator", + pressure=None, + test_sim_steps=250, + ) + print( + "Note: minimised structure file has not been written out, as " + "this would replace the PQR." + ) else: - run.test_sim(os.path.basename(outmodel), no_output = no_sim_output) + run.test_sim(os.path.basename(outmodel), no_output=no_sim_output) print("Done.") if ligand and no_sim_output: - print("Note: "+outmodel+" is not minimsed. This is because there" - " are possible ligands in the input file, and you may want to " - "add these back to the structure before minimising!") - + print( + "Note: " + outmodel + " is not minimsed. This is because there" + " are possible ligands in the input file, and you may want to " + "add these back to the structure before minimising!" + ) + with open(write_metadata, "w") as file: file.write(locals_json) - + # restore user's desired destination folder and handle abs/relative paths if final_dir: outmodel = os.path.join(final_dir, outmodel) if not os.path.isabs(final_dir): - outdir = run_dir+os.path.sep+final_dir+os.path.sep + outdir = run_dir + os.path.sep + final_dir + os.path.sep else: - outdir = final_dir+os.path.sep + outdir = final_dir + os.path.sep else: - outdir = run_dir+os.path.sep + outdir = run_dir + os.path.sep # copy all output files (structure, ligands, pqr, etc) to output directory if not pqr_out: - shutil.copyfile(os.path.basename(outmodel), - outdir+os.path.basename(outmodel)) - print("Wrote structure file to "+outdir+os.path.basename(outmodel)+".") - + shutil.copyfile(os.path.basename(outmodel), outdir + os.path.basename(outmodel)) + print("Wrote structure file to " + outdir + os.path.basename(outmodel) + ".") + # note: this copies to outdir which is found from outmodel so the pqr # currently can't go in a different folder - if pqr_out: - shutil.copyfile(os.path.basename(pqr_out), - outdir+os.path.basename(pqr_out)) - print("Wrote final PQR file to "+pqr_out+".") - - shutil.copyfile(write_metadata, outdir+write_metadata) - print("Wrote metadata to "+str(outdir+write_metadata)) - + if pqr_out: + shutil.copyfile(os.path.basename(pqr_out), outdir + os.path.basename(pqr_out)) + print("Wrote final PQR file to " + pqr_out + ".") + + shutil.copyfile(write_metadata, outdir + write_metadata) + print("Wrote metadata to " + str(outdir + write_metadata)) + if ligands and split_hetatms: for lig in ligands: - shutil.copyfile(lig, outdir+lig) - print("Wrote hetatms to "+outdir+lig) + shutil.copyfile(lig, outdir + lig) + print("Wrote hetatms to " + outdir + lig) def entry_point(): @@ -391,16 +480,30 @@ def entry_point(): args = parser.parse_args() fix_after = not args.fixstart if args.wdir is None: - args.wdir = args.code+"_"+id_generator(6) - prep(args.code, args.out, os.getcwd()+os.path.sep+args.wdir, - folder=args.directory, fastafile=args.fasta, inmodel=args.structure, - alignmentout=args.alignmentout, download_format=args.dlformat, - quiet=args.quiet, fix_after=fix_after, - download_sequence=args.download, fix_missing_atoms=args.leavemissing, - pqrff=args.pqrfmt, pqr_out=args.pqr, - redo=args.redo, num_models=args.num, em_map=args.em_map, - em_contour=args.contour, split_hetatms=args.hetatms, - no_modeller=args.nomodeller, ph=args.ph) + args.wdir = args.code + "_" + id_generator(6) + prep( + args.code, + args.out, + os.getcwd() + os.path.sep + args.wdir, + folder=args.directory, + fastafile=args.fasta, + inmodel=args.structure, + alignmentout=args.alignmentout, + download_format=args.dlformat, + quiet=args.quiet, + fix_after=fix_after, + download_sequence=args.download, + fix_missing_atoms=args.leavemissing, + pqrff=args.pqrfmt, + pqr_out=args.pqr, + redo=args.redo, + num_models=args.num, + em_map=args.em_map, + em_contour=args.contour, + split_hetatms=args.hetatms, + no_modeller=args.nomodeller, + ph=args.ph, + ) if __name__ == "__main__": diff --git a/prepmd/run.py b/prepmd/run.py index 8068083..c5b95c6 100644 --- a/prepmd/run.py +++ b/prepmd/run.py @@ -1,60 +1,110 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Use OpenMM to run simulations """ -from openmm.app import * -from openmm import * -from openmm.unit import * -from operator import methodcaller import argparse -import sys +import copy import json import os -import pathlib -import prepmd.metadynamics -import MDAnalysis as mda +import sys +from operator import methodcaller from sys import stdout -from prepmd import ligand -from prepmd import analysis + +import MDAnalysis as mda +from openmm import ( + CustomCVForce, + LangevinMiddleIntegrator, + MonteCarloBarostat, + OpenMMException, + Platform, + RMSDForce, + VariableLangevinIntegrator, +) +from openmm.app import ( + PME, + BiasVariable, + CutoffNonPeriodic, + CutoffPeriodic, + DCDReporter, + ForceField, + HBonds, + Metadynamics, + Modeller, + NoCutoff, + PDBFile, + PDBxFile, + Simulation, + StateDataReporter, + XTCReporter, +) +from openmm.unit import ( + Quantity, + amu, + angstrom, + bar, + kelvin, + kilocalories_per_mole, + kilojoule, + molar, + mole, + nanometer, + nanometers, + picosecond, + picoseconds, + quantity, +) + +import prepmd.metadynamics +from prepmd import analysis, ligand ff_lookup = { - "amber14,tip3p": ['amber14-all.xml', 'amber14/tip3p.xml'], - "amber14,spce": ['amber14-all.xml', 'amber14/spce.xml'], - "amber14,tip4pew": ['amber14-all.xml', 'amber14/tip4pew.xml'], - "charmm36,tip3p": ['charmm36.xml', 'charmm36/water.xml'], - "charmm36,tip4pew": ['charmm36.xml', 'charmm36/tip4pew.xml'], - "charmm36,spce": ['charmm36.xml', 'charmm36/spce.xml'], - "amoeba,tip3p": ['amoeba2018.xml'], - "amoeba,tip4pew": ['amoeba2018.xml'], - "amoeba,spce": ['amoeba2018.xml'], - "amber19,tip3p": ['amber19-all.xml' 'amber19/tip3p.xml'], - "amber19,tip4pew": ['amber19-all.xml' 'amber19/tip4pew.xml'], - "amber19,spce": ['amber19-all.xml' 'amber19/spce.xml'], - "amber14,None": ['amber14-all.xml'], - "charmm36,None": ['charmm36.xml'], - "amoeba,None": ['amoeba2018.xml'], - "amber19,None": ['amber19-all.xml'], - "amber19,implicit": ['amber19-all.xml', 'implicit/gbn2.xml'], - "amber14,implicit": ['amber14-all.xml', 'implicit/gbn2.xml',], - "charmm36,implicit": ['charmm36.xml', 'implicit/gbn2.xml'], - "amoeba,implicit": ['amoeba2018.xml', 'amoeba2009_gk.xml'], + "amber14,tip3p": ["amber14-all.xml", "amber14/tip3p.xml"], + "amber14,spce": ["amber14-all.xml", "amber14/spce.xml"], + "amber14,tip4pew": ["amber14-all.xml", "amber14/tip4pew.xml"], + "charmm36,tip3p": ["charmm36.xml", "charmm36/water.xml"], + "charmm36,tip4pew": ["charmm36.xml", "charmm36/tip4pew.xml"], + "charmm36,spce": ["charmm36.xml", "charmm36/spce.xml"], + "amoeba,tip3p": ["amoeba2018.xml"], + "amoeba,tip4pew": ["amoeba2018.xml"], + "amoeba,spce": ["amoeba2018.xml"], + "amber19,tip3p": ["amber19-all.xmlamber19/tip3p.xml"], + "amber19,tip4pew": ["amber19-all.xmlamber19/tip4pew.xml"], + "amber19,spce": ["amber19-all.xmlamber19/spce.xml"], + "amber14,None": ["amber14-all.xml"], + "charmm36,None": ["charmm36.xml"], + "amoeba,None": ["amoeba2018.xml"], + "amber19,None": ["amber19-all.xml"], + "amber19,implicit": ["amber19-all.xml", "implicit/gbn2.xml"], + "amber14,implicit": [ + "amber14-all.xml", + "implicit/gbn2.xml", + ], + "charmm36,implicit": ["charmm36.xml", "implicit/gbn2.xml"], + "amoeba,implicit": ["amoeba2018.xml", "amoeba2009_gk.xml"], } def make_forcefield(ff, solvent=None): - if ff+","+str(solvent) not in ff_lookup.keys(): - err = "Solvent "+str(solvent)+" not in list of available solvents." - validff = set([item[0] for item in list( - map(methodcaller("split", ","), ff_lookup.keys()))]) - validsolvent = set([item[1] for item in list( - map(methodcaller("split", ","), ff_lookup.keys()))]) + if ff + "," + str(solvent) not in ff_lookup.keys(): + err = "Solvent " + str(solvent) + " not in list of available solvents." + validff = set( + [ + item[0] + for item in list(map(methodcaller("split", ","), ff_lookup.keys())) + ] + ) + validsolvent = set( + [ + item[1] + for item in list(map(methodcaller("split", ","), ff_lookup.keys())) + ] + ) err += "Please pick a valid force field and solvent.\n" - err += "Valid force fields are "+str(validff)+"\n" - err += "Valid solvents are "+str(validsolvent)+"\n" + err += "Valid force fields are " + str(validff) + "\n" + err += "Valid solvents are " + str(validsolvent) + "\n" raise ValueError(err) - ffargs = ff_lookup[ff+","+str(solvent)] + ffargs = ff_lookup[ff + "," + str(solvent)] return ForceField(*ffargs) @@ -71,50 +121,58 @@ def check_platforms(): print("WARNING: no GPU platform found. Simulation will be slow!!!") -def test_sim(pdb, no_output = False): +def test_sim(pdb, no_output=False): if no_output: min_out = None else: min_out = pdb - run(pdb, minimised_structure_out=min_out, md_steps=None, - integrator_str="LangevinMiddleIntegrator", - pressure=None, test_sim_steps=250, no_output = no_output) - - -def run(pdb, - minimised_structure_out=None, - traj_out=None, - max_minimise_iterations=500, - minimise_error = None, #default = 0.001, - test_sim_steps=500, + run( + pdb, + minimised_structure_out=min_out, md_steps=None, - md_timestep=None, # default = 0.002*picoseconds, - forcefield_str=None, # or amber 14, amber19, amoeba - integrator_str="LangevinMiddleIntegrator", # or LangevinMiddleIntegrator - friction_coeff=1/picosecond, - temperature=300*kelvin, - minimise=True, - test_run=True, - fix_backbone=False, - constraints_str="Default", # Default, None, HBonds - solvent=None, - strip_solvent=False, - ionic_strength=0.0*molar, - pressure=None,# formerly 1*bar, - step=1000, - thermo_out_file="thermo.txt", - non_bonded_method_str="Default", - checkpoint_output="checkpoint.dat", - verbose=True, - write_params="params.json", - metadynamics_morph = None, # filename of structure to morph to - meta_rmsd_threshold_nm = 0.17, - no_output = False, - ligands = [], - ligands_names = None - ): + integrator_str="LangevinMiddleIntegrator", + pressure=None, + test_sim_steps=250, + no_output=no_output, + ) + + +def run( + pdb, + minimised_structure_out=None, + traj_out=None, + max_minimise_iterations=500, + minimise_error=None, # default = 0.001, + test_sim_steps=500, + md_steps=None, + md_timestep=None, # default = 0.002*picoseconds, + forcefield_str=None, # or amber 14, amber19, amoeba + integrator_str="LangevinMiddleIntegrator", # or LangevinMiddleIntegrator + friction_coeff=1 / picosecond, + temperature=300 * kelvin, + minimise=True, + test_run=True, + fix_backbone=False, + constraints_str="Default", # Default, None, HBonds + solvent=None, + strip_solvent=False, + ionic_strength=0.0 * molar, + pressure=None, # formerly 1*bar, + step=1000, + thermo_out_file="thermo.txt", + non_bonded_method_str="Default", + checkpoint_output="checkpoint.dat", + verbose=True, + write_params="params.json", + metadynamics_morph=None, # filename of structure to morph to + meta_rmsd_threshold_nm=0.17, + no_output=False, + ligands=None, + ligands_names=None, +): """ Run an MD simulation from a pdb/mmcif structure created with prepmd. + Args: minimised_structure_out - filename to write the final minimised structure to (string) @@ -150,7 +208,7 @@ def run(pdb, constraints_str - whether to constrain Hbonds or not. Possible values: None, HBonds. By default, HBonds will be constrained if the backbone is not fixed. A string. - solvent - solvent to use. Possible values: None (no solvent), tip3p, + solvent - solvent to use. Possible values: None (no solvent), tip3p, tip4pew, spce, implicit (gbn model, equivalent to AMBER igb=8). A string. write_solvent: whether to write solvent atoms to the minimised @@ -175,10 +233,13 @@ def run(pdb, ligand topology. Can run without this but ligands will show up as UNK """ + if ligands is None: + ligands = [] + # don't look at this locals_copy = copy.copy(locals()) for key, value in locals_copy.items(): - if type(value) == Quantity: + if type(value) is Quantity: locals_copy[key] = value.value_in_unit(value.unit) locals_json = json.dumps(locals_copy) @@ -189,77 +250,90 @@ def run(pdb, raise ValueError("Format of output trajectory must be DCD or XTC") if md_steps and traj_out and strip_solvent: - print("Warning: you are removing solvent molecules from the structure" - " file, but writing a trajectory, which will contain solvent" - "molecules. Ideally the trajectory and structure file should " - "have the same number of atoms.") + print( + "Warning: you are removing solvent molecules from the structure" + " file, but writing a trajectory, which will contain solvent" + "molecules. Ideally the trajectory and structure file should " + "have the same number of atoms." + ) if traj_out and not md_steps: - raise ValueError("Output trajectory specified, but not a number of " - "steps!") + raise ValueError("Output trajectory specified, but not a number of steps!") if md_steps and not traj_out: raise ValueError("No output trajectory file specified!") if (constraints_str == "HBonds" or constraints_str == HBonds) and fix_backbone: - raise ValueError("Can't fix the backbone with HBonds constraints! " - "Please disable constraints or unfix the backbone!") - + raise ValueError( + "Can't fix the backbone with HBonds constraints! " + "Please disable constraints or unfix the backbone!" + ) + # set some defaults - if forcefield_str == None: + if forcefield_str is None: if not ligands: forcefield_str = "charmm36" else: forcefield_str = "amber14" - - if solvent == None: + + if solvent is None: if not ligands: solvent = "implicit" else: solvent = "tip3p" - + # default timestep for ligand simulations is smaller if not md_timestep and not ligands: - md_timestep = 0.002*picoseconds + md_timestep = 0.002 * picoseconds if not md_timestep and ligands: - md_timestep = 0.00025*picoseconds + md_timestep = 0.00025 * picoseconds if not minimise_error and not ligands: minimise_error = 0.001 if not minimise_error and ligands: - minimise_error = 0.001*0.25 - #print("Timestep: "+str(md_timestep)) - #print("Min_err: "+str(minimise_error)) + minimise_error = 0.001 * 0.25 + # print("Timestep: "+str(md_timestep)) + # print("Min_err: "+str(minimise_error)) if forcefield_str == "amoeba" and not solvent: if non_bonded_method_str == "Default": non_bonded_method = NoCutoff - print("WARNING: AMOEBA does not support non-bonded cutoffs without" - " defining a periodic simulation box. For optimal " - "performance, use a solvated box or a different force" - " field.") + print( + "WARNING: AMOEBA does not support non-bonded cutoffs without" + " defining a periodic simulation box. For optimal " + "performance, use a solvated box or a different force" + " field." + ) else: - raise ValueError("AMOEBA can't be run with the specified non " - "bonded method without a simulation box. Try " - "running a solvated simulation instead.") + raise ValueError( + "AMOEBA can't be run with the specified non " + "bonded method without a simulation box. Try " + "running a solvated simulation instead." + ) if not solvent: print("WARNING: No solvent. Minimised structures may be wrong!!!") - + if fix_backbone and traj_out: - print("Warning: you are writing a trajectory but have constrained " - "the backbone. Your system won't move!") - + print( + "Warning: you are writing a trajectory but have constrained " + "the backbone. Your system won't move!" + ) + if minimise and not minimised_structure_out and not no_output: - minimised_structure_out = "min_"+os.path.basename(pdb) - print("No minimised structure output file specified. Writing " - "minimised structure to "+"min_"+pdb) + minimised_structure_out = "min_" + os.path.basename(pdb) + print( + "No minimised structure output file specified. Writing " + "minimised structure to " + "min_" + pdb + ) if ligands: if forcefield_str != "amber14": - raise ValueError("Error: Ligand simulations are only supported " - "with amber14") + raise ValueError( + "Error: Ligand simulations are only supported with amber14" + ) if solvent != "tip3p" and solvent != "tip4pew": - raise ValueError("Ligand simulations are only supported with" - " tip3p and tip4pew solvent") + raise ValueError( + "Ligand simulations are only supported with tip3p and tip4pew solvent" + ) if strip_solvent: raise ValueError("Can't strip solvent on a ligand system!") if constraints_str == "Default" or constraints_str == "HBonds": @@ -268,27 +342,33 @@ def run(pdb, if constraints_str == "None": ff = "ff14sb_off_impropers_0.0.4.offxml" ligand_ff = "openff_unconstrained-2.2.0.offxml" - if md_timestep >= 0.001*picoseconds or minimise_error >= 0.001: - print("WARNING: Timestep/error tolerance are high for a ligand " - "simulation. If the simulation crashes (e.g. due to " - "'Particle co-ordinate is NaN'), try reducing them.") + if md_timestep >= 0.001 * picoseconds or minimise_error >= 0.001: + print( + "WARNING: Timestep/error tolerance are high for a ligand " + "simulation. If the simulation crashes (e.g. due to " + "'Particle co-ordinate is NaN'), try reducing them." + ) if ".cif" in pdb or ".mmcif" in pdb: - writer = PDBxFile + writer = PDBxFile else: - writer = PDBFile - system, topology, positions = ligand.create_ligand_system(pdb, - ligands, - water_ff = solvent, - desired_charge_elementary = ionic_strength._value, - ff = ff, - ligand_ff = ligand_ff, - n_water = 5000, - ligands_names=ligands_names) - print("Note: topology written to "+str(minimised_structure_out)+" is " - "not minimised. This is because OpenMM cannot write a correct " - "topology for OpenFF ligands and solvents.") + writer = PDBFile + system, topology, positions = ligand.create_ligand_system( + pdb, + ligands, + water_ff=solvent, + desired_charge_elementary=ionic_strength._value, + ff=ff, + ligand_ff=ligand_ff, + n_water=5000, + ligands_names=ligands_names, + ) + print( + "Note: topology written to " + str(minimised_structure_out) + " is " + "not minimised. This is because OpenMM cannot write a correct " + "topology for OpenFF ligands and solvents." + ) # TODO: maybe write a separate file? - + else: # read in mmcif or pdb if ".cif" in pdb or ".mmcif" in pdb: @@ -297,35 +377,36 @@ def run(pdb, else: structure = PDBFile(pdb) writer = PDBFile - + modeller = Modeller(structure.topology, structure.positions) - + forcefield = make_forcefield(forcefield_str, solvent) - + if solvent and solvent != "implicit": print("Solvating system...") - modeller.addSolvent(forcefield, model=solvent, - ionicStrength=ionic_strength, - padding=1.0*nanometers) - + modeller.addSolvent( + forcefield, + model=solvent, + ionicStrength=ionic_strength, + padding=1.0 * nanometers, + ) + # can't set hbond constraints if also restraining whole backbone if constraints_str not in [None, "Default", HBonds, "HBonds"]: - raise ValueError( - "Constraints must be one of: Default, None, HBonds") + raise ValueError("Constraints must be one of: Default, None, HBonds") if constraints_str == "Default": if fix_backbone: constraints = None - print("As the backbone is fixed, no constraints will be " - "applied.") + print("As the backbone is fixed, no constraints will be applied.") else: constraints = HBonds print("No constraints specified. Constraining HBonds.") if constraints_str == "HBonds": constraints = HBonds - + # setup non bonded method - if none is chosen, select based on box if non_bonded_method_str == "Default": - if solvent == None or solvent == "implicit": + if solvent is None or solvent == "implicit": non_bonded_method = CutoffNonPeriodic else: non_bonded_method = PME @@ -338,58 +419,76 @@ def run(pdb, elif non_bonded_method == NoCutoff: pass else: - raise ValueError("Supply a valid nonbondedMethod! Could be " - "Default, PME, CutoffPeriodic or " - "CutoffNonPeriodic") + raise ValueError( + "Supply a valid nonbondedMethod! Could be " + "Default, PME, CutoffPeriodic or " + "CutoffNonPeriodic" + ) if non_bonded_method == PME and not solvent: - raise ValueError("Cannot run a PME simulation on a system with no " - "simulation box. Either change the " - "nonbondedMethod to 'NoCutoff' or add a solvent," - " which will initialise a box automatically.") - - #modeller.addHydrogens() - system = forcefield.createSystem(modeller.topology, - nonbondedMethod=non_bonded_method, - nonbondedCutoff=1*nanometer, - constraints=constraints) + raise ValueError( + "Cannot run a PME simulation on a system with no " + "simulation box. Either change the " + "nonbondedMethod to 'NoCutoff' or add a solvent," + " which will initialise a box automatically." + ) + + # modeller.addHydrogens() + system = forcefield.createSystem( + modeller.topology, + nonbondedMethod=non_bonded_method, + nonbondedCutoff=1 * nanometer, + constraints=constraints, + ) # constrain backbone by setting CA mass to zero if fix_backbone: for atom in modeller.topology.atoms(): - if atom.name == 'CA': - system.setParticleMass(atom.index, 0*amu) + if atom.name == "CA": + system.setParticleMass(atom.index, 0 * amu) if pressure: - print("Adding pressure coupling: "+str(pressure)) + print("Adding pressure coupling: " + str(pressure)) system.addForce(MonteCarloBarostat(pressure, temperature)) - + if metadynamics_morph: final = PDBFile(metadynamics_morph) ref_positions = final.getPositions(asNumpy=False) ref_positions_np = final.getPositions(asNumpy=True) r = RMSDForce(ref_positions) - rforce = CustomCVForce('1*r') - rforce.addCollectiveVariable('r', r) - rmsd = BiasVariable(force=rforce, minValue=0 * angstrom, - maxValue=500 * angstrom, biasWidth=1 * angstrom, - periodic=False) - meta = Metadynamics(system=system, variables=[rmsd], - temperature=300, biasFactor=2500.0, - height=10 * kilocalories_per_mole, - frequency=5, saveFrequency=50, biasDir='./') + rforce = CustomCVForce("1*r") + rforce.addCollectiveVariable("r", r) + rmsd = BiasVariable( + force=rforce, + minValue=0 * angstrom, + maxValue=500 * angstrom, + biasWidth=1 * angstrom, + periodic=False, + ) + meta = Metadynamics( + system=system, + variables=[rmsd], + temperature=300, + biasFactor=2500.0, + height=10 * kilocalories_per_mole, + frequency=5, + saveFrequency=50, + biasDir="./", + ) # set up integrator if integrator_str == "VariableLangevinIntegrator": - integrator = VariableLangevinIntegrator(temperature, - friction_coeff, minimise_error) + integrator = VariableLangevinIntegrator( + temperature, friction_coeff, minimise_error + ) elif integrator_str == "LangevinMiddleIntegrator": - integrator = LangevinMiddleIntegrator(temperature, - friction_coeff, md_timestep) + integrator = LangevinMiddleIntegrator(temperature, friction_coeff, md_timestep) else: - raise ValueError("Integrator must be one of: " - "VariableLangevinIntegrator," - " LangevinMiddleIntegrator") - + raise ValueError( + "Integrator must be one of: " + "VariableLangevinIntegrator," + " LangevinMiddleIntegrator" + ) + if ligands: simulation = Simulation(topology, system, integrator) simulation.context.setPositions(positions) @@ -402,34 +501,35 @@ def run(pdb, # run simulation # try and resolve steric clashes with variable langevin integrator # todo: this section could be tidier i think - strict_min = quantity.Quantity(value=2.5,unit=kilojoule/(nanometer*mole)) + strict_min = quantity.Quantity(value=2.5, unit=kilojoule / (nanometer * mole)) if minimise: try: print("Minimising...") - simulation.minimizeEnergy(maxIterations=max_minimise_iterations, - tolerance = strict_min) - curr_state = simulation.context.getState( - getPositions=True).getPositions(asNumpy=True) + simulation.minimizeEnergy( + maxIterations=max_minimise_iterations, tolerance=strict_min + ) + curr_state = simulation.context.getState(getPositions=True).getPositions( + asNumpy=True + ) except OpenMMException as e: if integrator_str != "VariableLangevinIntegrator": - print("Minimisation blew up. " - "Trying variable langevin integrator...") - integrator = VariableLangevinIntegrator(temperature, - friction_coeff, - minimise_error) + print("Minimisation blew up. Trying variable langevin integrator...") + integrator = VariableLangevinIntegrator( + temperature, friction_coeff, minimise_error + ) simulation = Simulation(topology, system, integrator) simulation.context.setPositions(positions) integrator_str = "VariableLangevinIntegrator" simulation.minimizeEnergy( - maxIterations=max_minimise_iterations, - tolerance = strict_min) + maxIterations=max_minimise_iterations, tolerance=strict_min + ) curr_state = simulation.context.getState( - getPositions=True).getPositions(asNumpy=True) + getPositions=True + ).getPositions(asNumpy=True) else: - print("Could not minimise even with variable langevin " - "integrator.") + print("Could not minimise even with variable langevin integrator.") raise e - + if test_run: try: print("Running test simulation...") @@ -437,25 +537,31 @@ def run(pdb, print("Test simulation ran successfully.") except OpenMMException as e: if integrator_str != "VariableLangevinIntegrator": - print("Simulation blew up, running with variable langevin " - "integrator...") - integrator = VariableLangevinIntegrator(temperature, - friction_coeff, - minimise_error) + print( + "Simulation blew up, running with variable langevin integrator..." + ) + integrator = VariableLangevinIntegrator( + temperature, friction_coeff, minimise_error + ) simulation = Simulation(topology, system, integrator) simulation.context.setPositions(positions) simulation.minimizeEnergy( - maxIterations = max_minimise_iterations, - tolerance = strict_min) + maxIterations=max_minimise_iterations, tolerance=strict_min + ) print("Fixed.") curr_state = simulation.context.getState( - getPositions=True).getPositions(asNumpy=True) + getPositions=True + ).getPositions(asNumpy=True) if md_steps: - print("Note: production MD will run with the variable " - "langevin integrator.") + print( + "Note: production MD will run with the variable " + "langevin integrator." + ) else: - print("Simulation blew up even while using variable langevin" - " integrator. This doesn't normally happen.") + print( + "Simulation blew up even while using variable langevin" + " integrator. This doesn't normally happen." + ) raise e with open(write_params, "w") as file: @@ -465,71 +571,108 @@ def run(pdb, modeller_out = Modeller(simulation.topology, curr_state) if strip_solvent: modeller_out.deleteWater() - writer.writeFile(modeller_out.topology, modeller_out.positions, - file=open(minimised_structure_out, "w"), keepIds=True) - print("Wrote minimised structure to "+str(minimised_structure_out)+".") + writer.writeFile( + modeller_out.topology, + modeller_out.positions, + file=open(minimised_structure_out, "w"), + keepIds=True, + ) + print("Wrote minimised structure to " + str(minimised_structure_out) + ".") else: if minimised_structure_out and not (minimise or test_run): - raise ValueError("Minimised structure output was requested, but" - " minimisation and test run were both skipped.") - + raise ValueError( + "Minimised structure output was requested, but" + " minimisation and test run were both skipped." + ) + if md_steps and traj_out and not metadynamics_morph: if ".dcd" in traj_out.lower(): simulation.reporters.append(DCDReporter(traj_out, step)) elif ".xtc" in traj_out.lower(): simulation.reporters.append(XTCReporter(traj_out, step)) if thermo_out_file: - simulation.reporters.append(StateDataReporter(thermo_out_file, - step, step=True, - potentialEnergy=True, - kineticEnergy=True, - elapsedTime=True, - temperature=True)) + simulation.reporters.append( + StateDataReporter( + thermo_out_file, + step, + step=True, + potentialEnergy=True, + kineticEnergy=True, + elapsedTime=True, + temperature=True, + ) + ) if verbose: - simulation.reporters.append(StateDataReporter(sys.stdout, - step, step=True, - speed=True, - remainingTime=True, - totalSteps=md_steps, - potentialEnergy=True, - temperature=True)) - - if md_steps and traj_out and not metadynamics_morph: + simulation.reporters.append( + StateDataReporter( + sys.stdout, + step, + step=True, + speed=True, + remainingTime=True, + totalSteps=md_steps, + potentialEnergy=True, + temperature=True, + ) + ) + + if md_steps and traj_out and not metadynamics_morph: print("Running production MD...") simulation.step(md_steps) simulation.saveCheckpoint(checkpoint_output) print("Done!") - print("Wrote trajectory to "+traj_out) + print("Wrote trajectory to " + traj_out) if thermo_out_file: - print("Wrote thermo info to "+thermo_out_file) - print("Wrote checkpoint to "+checkpoint_output) - + print("Wrote thermo info to " + thermo_out_file) + print("Wrote checkpoint to " + checkpoint_output) + if ligand and traj_out and minimised_structure_out: print("Running ligand analysis...") analysis.get_ligand_centroid_traj(minimised_structure_out, traj_out) print("Done.") - + # metadynamics run - + if metadynamics_morph and not traj_out: - print("Metadynamics run requested but no output trajectory specified. " - "Writing to metadynamics.xtc.") + print( + "Metadynamics run requested but no output trajectory specified. " + "Writing to metadynamics.xtc." + ) traj_out = "metadynamics.xtc" if metadynamics_morph and not md_steps: - print("Metadynamics run requested but no steps specified. Will run a " - "maximum of 5000000 steps.") + print( + "Metadynamics run requested but no steps specified. Will run a " + "maximum of 5000000 steps." + ) md_steps = 5000000 - + write_morph = False redo_run = False if metadynamics_morph: try: # Set up the integrator and simulation reporters simulation.reporters.append(XTCReporter(traj_out, 500)) - simulation.reporters.append(StateDataReporter(stdout, 500, step=True, - potentialEnergy=True, temperature=True, progress=True, remainingTime=True, totalSteps=5000000, separator='\t')) - simulation.reporters.append(prepmd.metadynamics.RMSDReporter( - "rmsd.txt", 500, ref_positions_np, threshold_nm=meta_rmsd_threshold_nm)) + simulation.reporters.append( + StateDataReporter( + stdout, + 500, + step=True, + potentialEnergy=True, + temperature=True, + progress=True, + remainingTime=True, + totalSteps=5000000, + separator="\t", + ) + ) + simulation.reporters.append( + prepmd.metadynamics.RMSDReporter( + "rmsd.txt", + 500, + ref_positions_np, + threshold_nm=meta_rmsd_threshold_nm, + ) + ) meta.step(simulation, md_steps) print("WARNING: Couldn't reach target structure.") except prepmd.metadynamics.StopSimulation: @@ -542,116 +685,281 @@ def run(pdb, pass if redo_run: - run(pdb, minimised_structure_out, traj_out, - max_minimise_iterations, minimise_error, - test_sim_steps, md_steps, md_timestep, forcefield_str, - integrator_str, friction_coeff, temperature, minimise, - test_run, fix_backbone, constraints_str, solvent, - strip_solvent, ionic_strength, pressure, step, - thermo_out_file, non_bonded_method_str, checkpoint_output, - verbose, write_params, metadynamics_morph) - + run( + pdb, + minimised_structure_out, + traj_out, + max_minimise_iterations, + minimise_error, + test_sim_steps, + md_steps, + md_timestep, + forcefield_str, + integrator_str, + friction_coeff, + temperature, + minimise, + test_run, + fix_backbone, + constraints_str, + solvent, + strip_solvent, + ionic_strength, + pressure, + step, + thermo_out_file, + non_bonded_method_str, + checkpoint_output, + verbose, + write_params, + metadynamics_morph, + ) + if write_morph: # write out pdbs and cifs as well final_traj = mda.Universe(pdb, traj_out) frames = prepmd.metadynamics.get_representative_rmsd_frames( - simulation.reporters[2].rmsd_log, 10) + simulation.reporters[2].rmsd_log, 10 + ) filtered = final_traj.trajectory[frames] with mda.Writer("morph.xtc", final_traj.atoms.n_atoms) as w: - for frame in filtered: + for _frame in filtered: w.write(final_traj) for frame in frames: - with mda.Writer("snapshot-"+str(frame)+".pdb") as w: + with mda.Writer("snapshot-" + str(frame) + ".pdb") as w: final_traj.trajectory[frame] w.write(final_traj) print("Done simulating.") - def entry_point(): "CLI entry point function. Uses sys.argv and argparse args object." - parser = argparse.ArgumentParser(prog="runmd", - description="Run MD simulations") + parser = argparse.ArgumentParser(prog="runmd", description="Run MD simulations") parser.add_argument("pdb", help="Input filename (pdb or mmcif)") parser.add_argument( - "-o", "--min_out", help="filename to write the final minimised structure to (pdb or mmCif)", default=None) + "-o", + "--min_out", + help="filename to write the final minimised structure to (pdb or mmCif)", + default=None, + ) + parser.add_argument( + "-tr", + "--traj_out", + help="Production MD trajectory output file. Can be in DCD or XTC format.", + default=None, + ) + parser.add_argument( + "-it", + "--minimise_iterations", + help="Maximum minimise iterations", + default=100, + type=int, + ) + parser.add_argument( + "-e", + "--minimise_err", + help=( + "error tolerance for variable langevin integrator. The value is arbitrary, " + "0.001 is a good starting point, increasing this will make the simulation " + "run faster at the expense of accuracy" + ), + default=None, + type=float, + ) + parser.add_argument( + "-test", + "--test_steps", + help=( + "how many steps to run of the test simulation. This isn't production MD, " + "this is just the simulation that checks that your structure doesn't have " + "any steric clashes" + ), + default=500, + type=int, + ) + parser.add_argument( + "-mdsteps", + "--md_steps", + help="how many steps to run production MD for", + default=None, + type=int, + ) + parser.add_argument( + "-ts", "--timestep", help="MD timestep in picoseconds", default=None, type=float + ) + parser.add_argument( + "-ff", + "--forcefield", + help=( + "Force field to use (can be charmm36, amber14, amoeba, or amber19 on newer " + "openmm versions" + ), + default=None, + ) + parser.add_argument( + "-i", + "--integrator", + help=( + "Numerical integration scheme. Can be one of: VariableLangevinIntegrator, " + "LangevinMiddleIntegrator" + ), + default="LangevinMiddleIntegrator", + ) + parser.add_argument( + "-f", + "--friction", + help=( + "Friction coefficient which couples the system to the heat bath. Units " + "1/picosecond." + ), + default=1, + type=float, + ) parser.add_argument( - "-tr", "--traj_out", help="Production MD trajectory output file. Can be in DCD or XTC format.", default=None) - parser.add_argument("-it", "--minimise_iterations", - help="Maximum minimise iterations", default=100, type=int) - parser.add_argument("-e", "--minimise_err", help="error tolerance for variable langevin integrator. The value is arbitrary, 0.001 is a good starting point, increasing this will make the simulation run faster at the expense of accuracy", default=None, type=float) - parser.add_argument("-test", "--test_steps", help="how many steps to run of the test simulation. This isn't production MD, this is just the simulation that checks that your structure doesn't have any steric clashes", default=500, type=int) - parser.add_argument("-mdsteps", "--md_steps", - help="how many steps to run production MD for", default=None, type=int) + "-t", + "--temperature", + help="Simulation temperature in Kelvin.", + default=300, + type=float, + ) parser.add_argument( - "-ts", "--timestep", help="MD timestep in picoseconds", default=None, type=float) - parser.add_argument("-ff", "--forcefield", - help="Force field to use (can be charmm36, amber14, amoeba, or amber19 on newer openmm versions", default=None) - parser.add_argument("-i", "--integrator", help="Numerical integration scheme. Can be one of: VariableLangevinIntegrator, LangevinMiddleIntegrator", - default="LangevinMiddleIntegrator") + "-nomin", "--no_minimise", help="Skip minimisation", action="store_true" + ) parser.add_argument( - "-f", "--friction", help="Friction coefficient which couples the system to the heat bath. Units 1/picosecond.", default=1, type=float) - parser.add_argument("-t", "--temperature", - help="Simulation temperature in Kelvin.", default=300, type=float) - parser.add_argument("-nomin", "--no_minimise", - help="Skip minimisation", action="store_true") - parser.add_argument("-notest", "--no_testrun", - help="Skip test run", action="store_true") - parser.add_argument("-fb", "--fix_backbone", - help="Fix the backbone", action="store_true") - parser.add_argument("-c", "--constraints", help="Constraints to apply. Possible values: None, HBonds. Constraining HBonds can make simulations run faster, at the cost of some accuracy. You can't constrain HBonds if you're also fixing the backbone, for obvious reasons.", default="Default") + "-notest", "--no_testrun", help="Skip test run", action="store_true" + ) parser.add_argument( - "-solv", "--solvent", help="Solvent to use. Possible values: tip3p, tip4pew, spce, implicit (GBn model, equivalent to AMBER igb=8). If no solvent is specified, implicit will be used.", default=None) - parser.add_argument("-ns", "--strip_solvent", - help="Don't write solvent molecules to the minimised output file", action="store_true") - parser.add_argument("-ion", "--ionic_strength", - help="Ionic strength of the solvent in molar.", default=0.0, type=float) + "-fb", "--fix_backbone", help="Fix the backbone", action="store_true" + ) parser.add_argument( - "-p", "--pressure", help="Pressure in bar, coupling via monte carlo barostat If pressure is not specified, no pressure coupling will be used.", default=None) + "-c", + "--constraints", + help=( + "Constraints to apply. Possible values: None, HBonds. Constraining HBonds " + "can make simulations run faster, at the cost of some accuracy. You can't " + "constrain HBonds if you're also fixing the backbone, for obvious reasons." + ), + default="Default", + ) parser.add_argument( - "-step", "--step", help="In production MD, write to output files (traj, thermo) every x steps.", default=1000, type=int) + "-solv", + "--solvent", + help=( + "Solvent to use. Possible values: tip3p, tip4pew, spce, implicit (GBn " + "model, equivalent to AMBER igb=8). If no solvent is specified, implicit " + "will be used." + ), + default=None, + ) parser.add_argument( - "-th", "--thermo", help="Thermo information output file.", default="thermo.txt") + "-ns", + "--strip_solvent", + help="Don't write solvent molecules to the minimised output file", + action="store_true", + ) parser.add_argument( - "-nb", "--nonbonded", help="Non-bonded interaction method. Can be: PME, CutoffPeriodic, CutoffNonPeriodic.", default="Default") - parser.add_argument("-ckpt", "--checkpoint", - help="Checkpoint output filename", default="checkpoint.dat") + "-ion", + "--ionic_strength", + help="Ionic strength of the solvent in molar.", + default=0.0, + type=float, + ) parser.add_argument( - "-q", "--quiet", help="Don't print info to the stdout", action="store_true") - parser.add_argument("-parm", "--write_params", - help="Write simulation params to a (json) file", default="params.json") - parser.add_argument("-m", "--metamorph", - help="Instead of running full MD, create a morph trajectory", default=None) - parser.add_argument("-l", "--ligand", action="append", - help="Add a ligand from an sdf file. Call multiple times to add multiple ligands.", default=None) - parser.add_argument("-ln", "--ligname", action="append", - help="Name the ligand. This is how ligands will be named in the topology file. Call once for each ligand.", default=None) + "-p", + "--pressure", + help=( + "Pressure in bar, coupling via monte carlo barostat If pressure is not " + "specified, no pressure coupling will be used." + ), + default=None, + ) + parser.add_argument( + "-step", + "--step", + help="In production MD, write to output files (traj, thermo) every x steps.", + default=1000, + type=int, + ) + parser.add_argument( + "-th", "--thermo", help="Thermo information output file.", default="thermo.txt" + ) + parser.add_argument( + "-nb", + "--nonbonded", + help=( + "Non-bonded interaction method. Can be: PME, CutoffPeriodic, " + "CutoffNonPeriodic." + ), + default="Default", + ) + parser.add_argument( + "-ckpt", + "--checkpoint", + help="Checkpoint output filename", + default="checkpoint.dat", + ) + parser.add_argument( + "-q", "--quiet", help="Don't print info to the stdout", action="store_true" + ) + parser.add_argument( + "-parm", + "--write_params", + help="Write simulation params to a (json) file", + default="params.json", + ) + parser.add_argument( + "-m", + "--metamorph", + help="Instead of running full MD, create a morph trajectory", + default=None, + ) + parser.add_argument( + "-l", + "--ligand", + action="append", + help=( + "Add a ligand from an sdf file. Call multiple times to add multiple " + "ligands." + ), + default=None, + ) + parser.add_argument( + "-ln", + "--ligname", + action="append", + help=( + "Name the ligand. This is how ligands will be named in the topology file. " + "Call once for each ligand." + ), + default=None, + ) args = parser.parse_args() minimise = not args.no_minimise testrun = not args.no_testrun verbose = not args.quiet - run(args.pdb, + run( + args.pdb, minimised_structure_out=args.min_out, traj_out=args.traj_out, max_minimise_iterations=args.minimise_iterations, minimise_error=args.minimise_err, test_sim_steps=args.test_steps, md_steps=args.md_steps, - md_timestep=args.timestep*picoseconds, + md_timestep=args.timestep * picoseconds, forcefield_str=args.forcefield, integrator_str=args.integrator, - friction_coeff=args.friction/picosecond, - temperature=args.temperature*kelvin, + friction_coeff=args.friction / picosecond, + temperature=args.temperature * kelvin, minimise=minimise, test_run=testrun, fix_backbone=args.fix_backbone, constraints_str=args.constraints, solvent=args.solvent, strip_solvent=args.strip_solvent, - ionic_strength=args.ionic_strength*molar, - pressure=args.pressure*bar, + ionic_strength=args.ionic_strength * molar, + pressure=args.pressure * bar, step=args.step, thermo_out_file=args.thermo, non_bonded_method_str=args.nonbonded, @@ -660,8 +968,8 @@ def entry_point(): write_params=args.write_params, metadynamics_morph=args.metamorph, ligands=args.ligand, - ligands_names=args.ligname - ) + ligands_names=args.ligname, + ) if __name__ == "__main__": diff --git a/prepmd/util.py b/prepmd/util.py index 2795e42..82b99c5 100644 --- a/prepmd/util.py +++ b/prepmd/util.py @@ -1,5 +1,4 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Utility functions """ @@ -34,8 +33,10 @@ def is_residue(resid): """ Check if a string is a valid residue. + Args: resid: the string + Returns: a boolean that is true if the string is a valid residue """ @@ -45,8 +46,10 @@ def is_residue(resid): def is_residue_sequence(sequence): """ Check if a sequence only contains valid residues + Args: sequence: an iterable + Returns: true if every item in the iterable is a valid residue, otherwise false """ @@ -60,14 +63,16 @@ def three_to_one(resid, ignore_non_standard=False): """ Convert the three-character residue representation found in SEQRES records to the one-character FASTA representation. + Args: resid: the residue, a string ignore_non_standard: ignore non-standard residues instead of throwing an error, a bool + Returns: the fasta residue, a string """ - if list(set(resid)) == ["."]: # hetatms + if list(set(resid)) == ["."]: # hetatms return resid if is_residue(resid): return codes[resid] @@ -81,8 +86,10 @@ def three_to_one_sequence(resids): """ Convert the three-character residue representation found in SEQRES records to the one-character FASTA representation. + Args: resid: an iterable of strings containing three-character residue codes + Returns: the residue sequence in FASTA format """ diff --git a/pyproject.toml b/pyproject.toml index 3cf32ab..a102bf8 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -33,6 +33,16 @@ Repository = "https://github.com/CCPBioSim/prepmd" Documentation = "https://prepmd.readthedocs.io/en/latest/" [project.optional-dependencies] +testing = [ + "pytest>=9.0,<10.0", + "pytest-cov>=7.0,<8.0", +] +pre-commit = [ + "pre-commit>=4.5,<5.0", + "ruff>=0.15.2,<1.0", + "pylint>=4.0,<5.0", + "rstcheck>=6.2,<7.0" +] docs = [ "sphinx>=9,<10", "nbsphinx>=0.9,<1.0", @@ -49,3 +59,14 @@ docs = [ prepmd = "prepmd.prep:entry_point" runmd = "prepmd.run:entry_point" prep-license = "prepmd.add_modeller_license:entry_point" + +[tool.ruff] +line-length = 88 +target-version = "py313" + +[tool.ruff.lint] +select = ["E", "F", "I", "B", "UP"] + +[tool.ruff.format] +quote-style = "double" +indent-style = "space" diff --git a/test/test_all.py b/test/test_all.py index f25d2cb..d64334d 100644 --- a/test/test_all.py +++ b/test/test_all.py @@ -1,25 +1,31 @@ #!/usr/bin/env python3 -# -*- coding: utf-8 -*- """ Integration tests """ +import os + import openmm.unit as unit -from prepmd.run import run + from prepmd.prep import prep -from prepmd.align_together import align_together -import os +from prepmd.run import run -def run_sim(code, path, fmt="pdb", ): - prep(code, - str(path)+os.path.sep+code+"."+fmt, - str(path)+os.path.sep+"testout"+os.path.sep+code+"_test", - download_format=fmt, no_modeller=True) +def run_sim( + code, + path, + fmt="pdb", +): + prep( + code, + str(path) + os.path.sep + code + "." + fmt, + str(path) + os.path.sep + "testout" + os.path.sep + code + "_test", + download_format=fmt, + no_modeller=True, + ) class TestPrep: - def test_9CS5(self, tmp_path): run_sim("9CS5", tmp_path, fmt="pdb") @@ -34,109 +40,156 @@ def test_7IB8(self, tmp_path): def test_9A9G(self, tmp_path): run_sim("9A9G", tmp_path, fmt="cif") - + def test_bestpdb(self, tmp_path): path = str(tmp_path) code = "6XOV" - prep(code, - str(path)+os.path.sep+code+"."+"pdb", - str(path)+os.path.sep+"testout"+os.path.sep+code+"_test", - download_format="pdb", - num_models=2, em_map="22281", em_contour=0.01, - no_modeller=True) - + prep( + code, + str(path) + os.path.sep + code + "." + "pdb", + str(path) + os.path.sep + "testout" + os.path.sep + code + "_test", + download_format="pdb", + num_models=2, + em_map="22281", + em_contour=0.01, + no_modeller=True, + ) + def test_pqr(self, tmp_path): path = str(tmp_path) code = "1UBQ" - prep(code, - str(path)+os.path.sep+code+"."+"pdb", - str(path)+os.path.sep+"testout"+os.path.sep+code+"_test", - download_format="pdb", - pqr_out=str(path)+os.path.sep+code+"."+"pqr", - no_modeller=True) + prep( + code, + str(path) + os.path.sep + code + "." + "pdb", + str(path) + os.path.sep + "testout" + os.path.sep + code + "_test", + download_format="pdb", + pqr_out=str(path) + os.path.sep + code + "." + "pqr", + no_modeller=True, + ) + # removed: 6TY4, 6XOV, 9I3U, 8RTO (too slow!) file_path = os.path.realpath(__file__) sep = os.path.sep -test_file = os.path.dirname(file_path)+sep+"test_data"+sep+"101M_proc.cif" +test_file = os.path.dirname(file_path) + sep + "test_data" + sep + "101M_proc.cif" class TestRun: - def test_basic(self, tmp_path): - run(test_file, - traj_out=str(tmp_path)+"101M_proc.xtc", md_steps=100, step=25, - solvent=None, pressure=None, minimise=False, test_run=False, - write_params=str(tmp_path)+sep+"params.json", - thermo_out_file=str(tmp_path)+sep+"thermo.txt", - checkpoint_output=str(tmp_path)+sep+"checkpoint.dat") + run( + test_file, + traj_out=str(tmp_path) + "101M_proc.xtc", + md_steps=100, + step=25, + solvent=None, + pressure=None, + minimise=False, + test_run=False, + write_params=str(tmp_path) + sep + "params.json", + thermo_out_file=str(tmp_path) + sep + "thermo.txt", + checkpoint_output=str(tmp_path) + sep + "checkpoint.dat", + ) def test_minimise_run(self, tmp_path): - run(test_file, minimised_structure_out=str(tmp_path)+"101m_prof_min.cif", - traj_out=str(tmp_path)+"101M_proc.xtc", md_steps=100, step=25, - solvent=None, pressure=None, - write_params=str(tmp_path)+sep+"params.json", - thermo_out_file=str(tmp_path)+sep+"thermo.txt", - checkpoint_output=str(tmp_path)+sep+"checkpoint.dat") + run( + test_file, + minimised_structure_out=str(tmp_path) + "101m_prof_min.cif", + traj_out=str(tmp_path) + "101M_proc.xtc", + md_steps=100, + step=25, + solvent=None, + pressure=None, + write_params=str(tmp_path) + sep + "params.json", + thermo_out_file=str(tmp_path) + sep + "thermo.txt", + checkpoint_output=str(tmp_path) + sep + "checkpoint.dat", + ) def test_variable_langevin(self, tmp_path): - run(test_file, - traj_out=str(tmp_path)+"101M_proc.xtc", md_steps=10, step=2, + run( + test_file, + traj_out=str(tmp_path) + "101M_proc.xtc", + md_steps=10, + step=2, integrator_str="VariableLangevinIntegrator", - solvent="tip4pew", pressure=1.0*unit.bar, minimise=False, + solvent="tip4pew", + pressure=1.0 * unit.bar, + minimise=False, test_run=False, - write_params=str(tmp_path)+sep+"params.json", - thermo_out_file=str(tmp_path)+sep+"thermo.txt", - checkpoint_output=str(tmp_path)+sep+"checkpoint.dat") + write_params=str(tmp_path) + sep + "params.json", + thermo_out_file=str(tmp_path) + sep + "thermo.txt", + checkpoint_output=str(tmp_path) + sep + "checkpoint.dat", + ) def test_amber14(self, tmp_path): - run(test_file, - traj_out=str(tmp_path)+"101M_proc.xtc", md_steps=5, step=1, - solvent="tip3p", forcefield_str="amber14", minimise=False, + run( + test_file, + traj_out=str(tmp_path) + "101M_proc.xtc", + md_steps=5, + step=1, + solvent="tip3p", + forcefield_str="amber14", + minimise=False, + test_run=False, + md_timestep=0.001 * unit.picoseconds, + write_params=str(tmp_path) + sep + "params.json", + thermo_out_file=str(tmp_path) + sep + "thermo.txt", + checkpoint_output=str(tmp_path) + sep + "checkpoint.dat", + ) + + def test_fix_backbone(self, tmp_path): # huge memory leak???? + run( + test_file, + traj_out=str(tmp_path) + "101M_proc.xtc", + md_steps=5, + step=1, + solvent="tip4pew", + pressure=1.0 * unit.bar, + forcefield_str="amber14", + fix_backbone=True, + minimise=True, test_run=False, - md_timestep=0.001*unit.picoseconds, - write_params=str(tmp_path)+sep+"params.json", - thermo_out_file=str(tmp_path)+sep+"thermo.txt", - checkpoint_output=str(tmp_path)+sep+"checkpoint.dat") - - def test_fix_backbone(self, tmp_path): # huge memory leak???? - run(test_file, - traj_out=str(tmp_path)+"101M_proc.xtc", md_steps=5, step=1, - solvent="tip4pew", pressure=1.0*unit.bar, forcefield_str="amber14", - fix_backbone=True, minimise=True, test_run=False, - md_timestep=0.001*unit.picoseconds, - write_params=str(tmp_path)+sep+"params.json", - thermo_out_file=str(tmp_path)+sep+"thermo.txt", - checkpoint_output=str(tmp_path)+sep+"checkpoint.dat") - + md_timestep=0.001 * unit.picoseconds, + write_params=str(tmp_path) + sep + "params.json", + thermo_out_file=str(tmp_path) + sep + "thermo.txt", + checkpoint_output=str(tmp_path) + sep + "checkpoint.dat", + ) + def test_metamorph(self, tmp_path): - testpath = os.path.dirname(file_path)+sep+"test_data"+sep - run(testpath+"6xou_cropped.pdb", - metadynamics_morph = testpath+"6xov_cropped.pdb", - minimised_structure_out = str(tmp_path)+"_min.pdb", - meta_rmsd_threshold_nm = 0.90) + testpath = os.path.dirname(file_path) + sep + "test_data" + sep + run( + testpath + "6xou_cropped.pdb", + metadynamics_morph=testpath + "6xov_cropped.pdb", + minimised_structure_out=str(tmp_path) + "_min.pdb", + meta_rmsd_threshold_nm=0.90, + ) # note: the threshold being high means this isn't a very rigorous # test - this test always passes on my local machine but it stalls # on the github CI for unknown reasons - possibly to do with the RNG? def test_ligand(self, tmp_path): - testpath = os.path.dirname(file_path)+sep+"test_data"+sep - run(testpath+"3PDT_out.pdb", - traj_out=str(tmp_path)+"3PDT.xtc", md_steps=3, step=1, - ligands = [testpath+"ADP.sdf"], - write_params=str(tmp_path)+sep+"params.json", - thermo_out_file=str(tmp_path)+sep+"thermo.txt", - checkpoint_output=str(tmp_path)+sep+"checkpoint.dat", - minimise=True, test_run=True) - -test_aln1 = os.path.dirname(file_path)+sep+"test_data"+sep+"6xov_prep.pdb" -test_aln2 = os.path.dirname(file_path)+sep+"test_data"+sep+"6xou_prep.pdb" + testpath = os.path.dirname(file_path) + sep + "test_data" + sep + run( + testpath + "3PDT_out.pdb", + traj_out=str(tmp_path) + "3PDT.xtc", + md_steps=3, + step=1, + ligands=[testpath + "ADP.sdf"], + write_params=str(tmp_path) + sep + "params.json", + thermo_out_file=str(tmp_path) + sep + "thermo.txt", + checkpoint_output=str(tmp_path) + sep + "checkpoint.dat", + minimise=True, + test_run=True, + ) + + +test_aln1 = os.path.dirname(file_path) + sep + "test_data" + sep + "6xov_prep.pdb" +test_aln2 = os.path.dirname(file_path) + sep + "test_data" + sep + "6xou_prep.pdb" # can't be tested without a modeller key because it uses modeller for alignment -#class TestAlignTogether: -# +# class TestAlignTogether: +# # def test_align_together(self, tmp_path): -# align_together(test_aln1, test_aln2, +# align_together(test_aln1, test_aln2, # str(tmp_path)+sep+"6xov_cropped.pdb", -# str(tmp_path)+sep+"6xou_cropped.pdb", "6xov", "6xou") \ No newline at end of file +# str(tmp_path)+sep+"6xou_cropped.pdb", "6xov", "6xou") diff --git a/test/test_packaging.py b/test/test_packaging.py index 6ea6178..c4299a9 100644 --- a/test/test_packaging.py +++ b/test/test_packaging.py @@ -2,9 +2,10 @@ import importlib -import prepmd import pytest +import prepmd + def test_package_version_is_defined(): """The package should expose a valid version string.""" @@ -24,4 +25,4 @@ def test_command_entry_points_are_importable(module_name, function_name): """Every command declared in pyproject.toml should be importable.""" module = importlib.import_module(module_name) - assert callable(getattr(module, function_name)) \ No newline at end of file + assert callable(getattr(module, function_name))