diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 2656b1d..753ddb0 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -17,7 +17,7 @@ jobs: - uses: actions/checkout@v7 - uses: astral-sh/setup-uv@v7 - run: uv python install ${{ matrix.python-version }} - - run: uv sync --extra test + - run: uv sync --extra test --extra pymol - run: uv run pytest lint: diff --git a/CHANGELOG.md b/CHANGELOG.md index 4dd43e9..7a9fd10 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,18 @@ All notable changes to DBSTEP. Versions follow [semantic versioning](https://semver.org/). -## Unreleased +## 2.2.0 — 2026-09-26 + +DBSTEP now covers everything the earlier `patonlab/sterimol` and `patonlab/wsterimol` codes did: CPK radii, cone angles of half-sandwich and phosphine ligands, Boltzmann weighting over separate QM outputs and a PyMOL plugin. + +### Added +- PyMOL plugin (`import dbstep.pymol_plugin` inside PyMOL): `dbstep_sterimol`, `dbstep_vbur` (with residue contributions coloured on the protein), `dbstep_cone`, `dbstep_vdw`, `dbstep_ensemble`, `dbstep_conformers` and `dbstep_style`, all taking PyMOL selections and drawing the results on the open structure. Tested headless against `pymol-open-source` (optional extra `pymol`). +- `from_coords(atoms, coords, ...)` and `parse_data.Structure`: measure structures held in memory, with optional residue metadata for `--decompose`; `aligned_coords`, `sterimol_vectors()` and `calculator.rigid_transform` expose DBSTEP's frame for drawing. +- `--radii cpk`: the CPK radii of Verloop's original Sterimol program, with atom types (sp3/aromatic carbon, single/double-bonded oxygen, ...) assigned from DFT-D3 style coordination numbers, so the earlier `patonlab/sterimol` code is no longer needed. Classic Sterimol values with this set match the original Fortran program to 0.01 Å on the benchmark set. +- `--cone`: Tolman cone angle and metal-to-centroid distance of a ligand (ring or donor atom, auto-detected or given with `--atom1`/`--atom2`), with the ligand's Sterimol parameters measured from the metal; `cone_angle` and `metal_centroid` columns in `results` and the CSV. Reproduces the half-sandwich analysis of `patonlab/sterimol` (173.97° for its [RhCp*Cl2(PMe3)] example with `--radii cpk`). +- `--boltzmann` over separate QM output files: energies are read through cclib (Gibbs free energy when available, else the SCF energy, in hartree), gzipped outputs are accepted, and several single-structure input files form one `ensemble` row. Reproduces the wSterimol example (wL 6.33, wB1 1.79, wB5 3.75), so `patonlab/wsterimol` is superseded by AQME/CREST + DBSTEP. +- `--energy-window` (kcal/mol): conformers above the window are left out of the Boltzmann average but stay listed with population 0; `boltzmann_average(runs, window=..., label=...)` in Python, with `energy_rel`, `in_window` and `energy_key` set on each run. +- `path` column (the input path as given) in `results`, the CSV and the contributions CSV, so same-named files from different folders stay distinguishable. ### Changed - `run.energy` set by `boltzmann_average` is now always in kcal/mol (it used to be in the input units). @@ -15,13 +26,6 @@ All notable changes to DBSTEP. Versions follow [semantic versioning](https://sem - `all_frames(frames=...)` no longer changes the caller's options object; `frames=0` from Python selects frame 0 instead of all frames. - `--residue all --atom1 N` uses the named atom to pick residues (previously only `--atom` was consulted); a default atom2 can no longer coincide with atom1 (e.g. `--atom CB`), and an explicit coincidence is an error. -### Added -- `--radii cpk`: the CPK radii of Verloop's original Sterimol program, with atom types (sp3/aromatic carbon, single/double-bonded oxygen, ...) assigned from DFT-D3 style coordination numbers, so the earlier `patonlab/sterimol` code is no longer needed. Classic Sterimol values with this set match the original Fortran program to 0.01 Å on the benchmark set. -- `--cone`: Tolman cone angle and metal-to-centroid distance of a ligand (ring or donor atom, auto-detected or given with `--atom1`/`--atom2`), with the ligand's Sterimol parameters measured from the metal; `cone_angle` and `metal_centroid` columns in `results` and the CSV. Reproduces the half-sandwich analysis of `patonlab/sterimol` (173.97° for its [RhCp*Cl2(PMe3)] example with `--radii cpk`). -- `--boltzmann` over separate QM output files: energies are read through cclib (Gibbs free energy when available, else the SCF energy, in hartree), gzipped outputs are accepted, and several single-structure input files form one `ensemble` row. Reproduces the wSterimol example (wL 6.33, wB1 1.79, wB5 3.75), so `patonlab/wsterimol` is superseded by AQME/CREST + DBSTEP. -- `--energy-window` (kcal/mol): conformers above the window are left out of the Boltzmann average but stay listed with population 0; `boltzmann_average(runs, window=..., label=...)` in Python, with `energy_rel`, `in_window` and `energy_key` set on each run. -- `path` column (the input path as given) in `results`, the CSV and the contributions CSV, so same-named files from different folders stay distinguishable. - ## 2.1.0 — 2026-09-25 ### Added diff --git a/CLAUDE.md b/CLAUDE.md index e011825..6b2fb66 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -16,11 +16,12 @@ DBSTEP (DFT-Based Steric Parameters) is a Python package for computing steric pa - `sterics.py` also holds `buried_vol_by_group` (per-residue %V_bur decomposition for `--decompose`) - `trajectory.py` — Frames of multi-structure files (`--frames` selection, frame counting) - `ensemble.py` — Boltzmann weighting over conformer ensembles (`--boltzmann`, `--energy-window`; energies from SDF data fields, xyz comments or QM outputs via cclib, pooled across single-structure files) - - `parse_data.py` — Input file parsing (xyz, sdf/mol, pdb with residue metadata, cube, cclib-supported formats) + - `parse_data.py` — Input file parsing (xyz, sdf/mol, pdb with residue metadata, cube, cclib-supported formats) and `Structure` for in-memory input (`from_coords`) - `constants.py` — Chemical constants (periodic table, Bondi/Charry-Tkatchenko/CPK radii, covalent radii, metals) - `cone.py` — Tolman cone angle, metal-centroid distance and ligand detection (`--cone`) - `radii.py` — Radius assignment per set; CPK atom typing from DFT-D3 coordination numbers (`--radii cpk`) - `graph.py` — 2D graph-based steric contribution calculations + - `pymol_plugin.py` — `dbstep_*` commands for PyMOL (selections in, CGO drawings out); importable without PyMOL for its pure helpers - `writer.py` — Output formatting and file writing (PyMOL scripts, xyz, CSV results) - `__init__.py` — Package init, `__version__`, `__all__` - `__main__.py` — Module entry point for `python -m dbstep` @@ -40,6 +41,8 @@ DBSTEP (DFT-Based Steric Parameters) is a Python package for computing steric pa - `test_trajectory.py` — `--frames` parsing, per-frame runs on `ala5_traj.pdb`, multi-frame xyz, CSV time series - `test_decompose.py` — `--decompose` per-residue contributions: sums to %V_bur, overlap sharing, CLI/CSV - `test_qm_ensemble.py` — Boltzmann weighting over separate (gzipped) Gaussian outputs vs the wSterimol example, `--energy-window`, pooling rule + - `test_from_coords.py` — In-memory `Structure` input, `--decompose` from metadata, `rigid_transform` and `sterimol_vectors` + - `test_pymol_plugin.py` — The PyMOL plugin against headless `pymol-open-source` (skipped when PyMOL is missing; `uv sync --extra pymol`) - `test_examples.py` — Executes the code cells of `examples/proteins_and_conformers.ipynb` - `test_ensemble.py` — SDF data fields, Boltzmann weights/averages, CLI `--boltzmann` on `sdf_files/ether_conformers.sdf` (AQME output) - `cube_files/` — Test cube file fixtures (keep these small; use `benzene_coarse.cube` / `Ne_medium.cube` for new tests) @@ -67,6 +70,11 @@ uv sync uv sync --extra dev ``` +### Install open-source PyMOL for the plugin tests (Python < 3.14) +``` +uv sync --extra pymol +``` + ### Install extras for the analysis scripts (matplotlib, pandas, rdkit, tqdm) ``` uv sync --group analysis diff --git a/README.md b/README.md index dd7512f..8a2c091 100644 --- a/README.md +++ b/README.md @@ -206,6 +206,26 @@ If no atoms are specified, the first two atoms in the file will be used as refer `-b` adds the buried volume of the ligand around the metal, `--boltzmann` averages cone angles over conformers, and `--csv` writes `cone_angle` and `metal_centroid` columns. In Python the object carries `cone_angle`, `metal_centroid`, `cone_sectors` (half angle per sector) and `ligand_atoms`. For dimers or several metals, choose the centre with `--atom1`. +### PyMOL plugin + +The same measurements are available inside PyMOL, on the structures you have open, with atoms picked by PyMOL selections and the results drawn on the molecule (this replaces the visual side of [wSterimol](https://github.com/patonlab/wsterimol)). Install DBSTEP into the Python that PyMOL uses (`pip install dbstep`; the open-source `pymol-open-source` wheel on PyPI works too), then in PyMOL or in your `.pymolrc`: + +``` +import dbstep.pymol_plugin +``` + +| Command | What it does | +|---|---| +| `dbstep_sterimol atom1, atom2 [, atom3, radii=bondi, measure=classic, selection=, name=sterimol]` | L, Bmin and Bmax from `atom1` along the `atom1`-`atom2` axis, drawn as a blue L axis and green/red Bmin/Bmax circles in the object's own frame | +| `dbstep_vbur atom1 [, radius=3.5, radii=bondi, selection=, decompose=1]` | %V_bur in a sphere around `atom1`; with `decompose` every residue's B-factor is set to its contribution and the contributing residues are coloured white to red | +| `dbstep_cone metal [, ligand, radii=cpk]` | Tolman cone angle, metal-centroid distance and the ligand's Sterimol parameters, drawn as a translucent cone plus the distance | +| `dbstep_vdw object [, radii=bondi, scale=1.0]` | Translucent van der Waals copy `object_vdw` with DBSTEP's radii (Bondi, Charry-Tkatchenko or CPK atom types) | +| `dbstep_ensemble files, atom1, atom2 [, radii, temperature, window, vbur=1]` | Boltzmann-weighted parameters over conformers (QM outputs or a multi-structure file), loaded into PyMOL with their populations as titles and a transparency that fades out minor conformers | +| `dbstep_conformers folder [, pattern=*.pdb]` | Load every structure of a folder into one group | +| `dbstep_style` | White background and light settings for figures | + +`selection` restricts the atoms that are measured; atoms named as `atom1`/`atom2` outside it are kept as zero-radius ghosts, so `dbstep_vbur /1a8o//A/186/CA, 3.5, polymer and not resi 186` measures the pocket around Thr186 without waters and without the residue itself. Every command returns the DBSTEP run object, so from the PyMOL Python prompt the full results (`results`, `contributions`, `cone_sectors`, ...) are available. From Python outside PyMOL the same in-memory route is `dbstep.Dbstep.from_coords(atoms, coords, ...)`. + ### Examples A notebook covering the protein, trajectory and conformer-ensemble workflows is at `examples/proteins_and_conformers.ipynb`. diff --git a/dbstep/Dbstep.py b/dbstep/Dbstep.py index c7e1c4e..9255388 100644 --- a/dbstep/Dbstep.py +++ b/dbstep/Dbstep.py @@ -47,6 +47,8 @@ def __init__(self, *args, **kwargs): self.tensor, self.tensor_grid = False, False # Cone angle parameters (--cone) self.cone_angle, self.metal_centroid, self.ligand_atoms, self.cone_sectors = None, None, None, None + # PyMOL CGO strings of the last Bmax/Bmin measurement (see sterimol_vectors) + self.cylinders = [] if "options" in kwargs: # Work on a copy: the calculation renumbers spec atoms (--noH/--exclude) and @@ -65,6 +67,8 @@ def __init__(self, *args, **kwargs): if isinstance(file, str): name, ext = os.path.splitext(file) + elif isinstance(file, parse_data.Structure): + name, ext = file.name, "array" else: name = file ext = "rdkit" @@ -168,6 +172,9 @@ def __init__(self, *args, **kwargs): # Rotate molecule to align atom1-atom2 bond along Z-axis self._orient_molecule(mol, options) + # the measured atoms in DBSTEP's frame (atom1 at the origin, axis along z); with self.coords this + # gives the rigid transform back to the input frame (calculator.rigid_transform), used by the PyMOL plugin + self.aligned_coords = np.array(mol.CARTESIANS) if options.surface == "vdw" else None # Recompute grid bounds after rotation so the grid covers the full rotated molecule if options.sterimol and options.surface == "vdw": @@ -226,7 +233,7 @@ def __init__(self, *args, **kwargs): print("\n L parameter is {:5.2f} Ang".format(L)) # Write PyMOL visualization files - if options.pymol and ext != "rdkit": + if options.pymol and ext not in ("rdkit", "array"): if options.sterimol: cylinders.append(" CYLINDER, 0., 0., 0., 0., 0., {:5.3f}, 0.1, 1.0, 1.0, 1.0, 0., 0.0, 1.0,".format(self.L)) writer.xyz_export(file, mol) @@ -393,6 +400,8 @@ def _compute(self, mol, file, options, origin, occ_grid, occ_vol, point_tree, gr if isinstance(file, str): fname = os.path.basename(file) + elif isinstance(file, parse_data.Structure): + fname = file.name else: try: fname = file.GetProp("_Name") @@ -448,6 +457,8 @@ def _compute(self, mol, file, options, origin, occ_grid, occ_vol, point_tree, gr sys.exit(" Can't use classic Sterimol with the isodensity surface. Use --measure grid or --surface vdw") Bmin_list.append(Bmin) Bmax_list.append(Bmax) + # CGO strings for the Bmax/Bmin directions in the aligned frame (see sterimol_vectors) + self.cylinders = list(cyl) if options.pymol: cylinders.extend(cyl) @@ -519,6 +530,20 @@ def _compute(self, mol, file, options, origin, occ_grid, occ_vol, point_tree, gr return spheres, cylinders + def sterimol_vectors(self): + """Directions of the last Bmax and Bmin measurements in the aligned frame, as (x, y, z) points on the + molecular surface, or None when no Sterimol parameters were computed.""" + if not self.cylinders: + return None + vectors = {} + for cylinder in self.cylinders: + values = [v.strip() for v in cylinder.split(",")] + if values[0].endswith("CYLINDER") and len(values) >= 14: + point = tuple(float(v) for v in values[4:7]) + red, green = float(values[11]), float(values[12]) + vectors["bmax" if red == 1.0 and green == 0.0 else "bmin"] = point + return vectors or None + def _decompose(self, mol, options, grid_axes, origin, R): """Per-residue contributions to %V_bur (--decompose), as a dict "A:45 LEU" -> percent.""" if options.surface != "vdw" or grid_axes is None: @@ -702,6 +727,22 @@ def format_result_row(row, options, fw): return fmt.format(label, row["atom1"], row["atom2"], row["bmin"], row["bmax"], row["L"]) + cone_cols +def from_coords(atoms, coords, name="structure", metadata=None, **kwargs): + """Measure a structure held in memory. Accepts the same keyword arguments as dbstep(). + + Args: + atoms: element symbols ("Bq" marks a zero-radius ghost atom) + coords: Cartesian coordinates in Angstrom, shape (n, 3) + name: label used in the results table + metadata: optional per-atom arrays as the PDB parser provides them (name, resname, chain, resseq, + icode, resid, element, record); with them --decompose works on in-memory structures + + Example: + result = from_coords(["H", "C", "H", "H", "H"], xyz, atom1=1, atom2=2, sterimol=True) + """ + return dbstep(parse_data.Structure(atoms, coords, name, metadata), **kwargs) + + def from_rdkit(mol, **kwargs): """Measure an RDKit molecule that carries a 3D conformer. Accepts the same keyword arguments as dbstep(). diff --git a/dbstep/__init__.py b/dbstep/__init__.py index 4218c6f..ab755e5 100644 --- a/dbstep/__init__.py +++ b/dbstep/__init__.py @@ -1,6 +1,6 @@ """DBSTEP: DFT-Based Steric Parameters.""" -__version__ = "2.1.0" +__version__ = "2.2.0" __all__ = [ "Dbstep", @@ -9,6 +9,9 @@ "parse_data", "writer", "constants", + "radii", + "cone", + "ensemble", ] from dbstep import Dbstep @@ -17,3 +20,6 @@ from dbstep import parse_data from dbstep import writer from dbstep import constants +from dbstep import radii +from dbstep import cone +from dbstep import ensemble diff --git a/dbstep/calculator.py b/dbstep/calculator.py index ff1bbe0..0ee0b03 100644 --- a/dbstep/calculator.py +++ b/dbstep/calculator.py @@ -189,3 +189,20 @@ def translate_dens(mol, options, xmin, xmax, ymin, ymax, zmin, zmax, xyz_max, or except Exception: sys.exit(" WARNING! Unable to find an atom (e.g. metal) to set at the origin") return [coords, cube_origin, xmin, xmax, ymin, ymax, zmin, zmax, xyz_max] + + +def rigid_transform(source, target): + """Rotation R and translation t with target = source @ R.T + t (Kabsch fit; exact for a rigid motion). + + Used to map DBSTEP's aligned frame (atom1 at the origin, atom1-atom2 along z) back onto the input + coordinates: original = (aligned - t) @ R. + """ + source = np.asarray(source, dtype=float).reshape(-1, 3) + target = np.asarray(target, dtype=float).reshape(-1, 3) + source_centre, target_centre = source.mean(axis=0), target.mean(axis=0) + h = (source - source_centre).T @ (target - target_centre) + u, _, vt = np.linalg.svd(h) + d = np.sign(np.linalg.det(vt.T @ u.T)) or 1.0 + rotation = vt.T @ np.diag([1.0, 1.0, d]) @ u.T + translation = target_centre - source_centre @ rotation.T + return rotation, translation diff --git a/dbstep/parse_data.py b/dbstep/parse_data.py index f907e26..c943e9d 100644 --- a/dbstep/parse_data.py +++ b/dbstep/parse_data.py @@ -185,6 +185,8 @@ def read_input(molecule, ext, options): mol = PDBParser(molecule, ext[1:], False, False, options.spec_atom_1, options.spec_atom_2, structure=structure) elif ext == "rdkit": mol = RDKitParser(molecule, False, False, options.spec_atom_1, options.spec_atom_2) + elif ext == "array": + mol = ArrayParser(molecule, False, False, options.spec_atom_1, options.spec_atom_2) else: mol = cclibParser(molecule, ext[1:], False, False, options.spec_atom_1, options.spec_atom_2) if getattr(options, "radii", "bondi") == "cpk": @@ -587,6 +589,43 @@ def parse_input(self): self.PROPERTIES.update(energies) +class Structure: + """A structure held in memory, accepted by dbstep() in place of a file name (see Dbstep.from_coords). + + Args: + atoms: element symbols + coords: Cartesian coordinates in Angstrom, shape (n, 3) + name: label used in the results table + metadata: optional per-atom arrays (record, name, resname, chain, resseq, icode, element, resid), + as the PDB parser provides them, which enable --decompose + properties: optional per-structure properties (e.g. an "energy" for Boltzmann weighting) + """ + + def __init__(self, atoms, coords, name="structure", metadata=None, properties=None): + self.atoms = [str(atom) for atom in atoms] + self.coords = np.asarray(coords, dtype=float).reshape(-1, 3) + if len(self.atoms) != len(self.coords): + raise ValueError("Structure: {} atoms but {} coordinates".format(len(self.atoms), len(self.coords))) + self.name = str(name) + self.metadata = dict(metadata or {}) + self.properties = dict(properties or {}) + + +class ArrayParser(DataParser): + """Wrap an in-memory Structure in the parser interface.""" + + def __init__(self, structure, noH, exclude, spec_atom_1, spec_atom_2): + super().__init__(structure, "array", noH, exclude, spec_atom_1, spec_atom_2) + + def parse_input(self): + structure = self._input + self.ATOMTYPES = list(structure.atoms) + self.CARTESIANS = [list(row) for row in structure.coords] + if structure.metadata: + self.METADATA = {key: np.array(list(values)) for key, values in structure.metadata.items()} + self.PROPERTIES.update(structure.properties) + + class RDKitParser(DataParser): """Extract coordinates and atom types from rdkit mol object diff --git a/dbstep/pymol_plugin.py b/dbstep/pymol_plugin.py new file mode 100644 index 0000000..200ecd1 --- /dev/null +++ b/dbstep/pymol_plugin.py @@ -0,0 +1,420 @@ +# -*- coding: UTF-8 -*- + +""" +pymol_plugin + +DBSTEP inside PyMOL: measure the structures you have open, pick atoms with PyMOL selections and see +the Sterimol axes, buried-volume spheres, residue contributions and cone angles drawn on the molecule. + +Setup (once, in the Python that PyMOL uses): ``pip install dbstep``. Then in PyMOL, or in your +.pymolrc: + + import dbstep.pymol_plugin + +(or ``run /path/to/dbstep/pymol_plugin.py``; ``python -m dbstep.pymol_plugin`` prints that path). + +Commands (all arguments are PyMOL selections or plain values, ``help dbstep_sterimol`` in PyMOL): + + dbstep_sterimol atom1, atom2 [, atom3, radii=bondi, measure=classic, grid=0.05, scale=1.0, noH=0, selection=, name=sterimol] + dbstep_vbur atom1 [, radius=3.5, radii=bondi, scale=1.0, noH=0, selection=, decompose=1, grid=0.05, name=vbur] + dbstep_cone metal [, ligand, radii=cpk, scale=1.0, selection=, name=cone] + dbstep_vdw object [, radii=bondi, scale=1.0, transparency=0.6] + dbstep_ensemble files, atom1, atom2 [, radii=bondi, temperature=298.15, window=, tag=, vbur=0, radius=3.5, group=ensemble] + dbstep_conformers folder [, pattern=*.pdb, group=conformers] + dbstep_style +""" + +import glob +import math +import os +import sys + +import numpy as np + +try: + from pymol import cmd + from pymol.cgo import ALPHA, BEGIN, COLOR, CONE, CYLINDER, END, LINEWIDTH, LINE_LOOP, SPHERE, VERTEX +except ImportError: # imported outside PyMOL (tests, --help): the helpers still work, nothing is registered + cmd = None + +from dbstep import Dbstep, calculator, ensemble, trajectory +from dbstep import radii as radii_sets +from dbstep.parse_data import Structure + +COLORS = {"L": (0.25, 0.45, 1.0), "bmin": (0.15, 0.75, 0.25), "bmax": (0.95, 0.2, 0.2), "sphere": (0.95, 0.45, 0.8), "cone": (1.0, 0.7, 0.15), "axis": (0.4, 0.4, 0.4)} +_ICODE_LETTERS = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz" + + +class PluginError(Exception): + pass + + +def _flag(value): + return str(value).strip().lower() in ("1", "true", "yes", "on") + + +def _object_of(selection): + """The single PyMOL object a selection belongs to.""" + objects = sorted({obj for obj, _ in cmd.index(selection)}) + if not objects: + raise PluginError("selection '{}' matches no atoms".format(selection)) + if len(objects) > 1: + raise PluginError("selection '{}' spans several objects ({}); pick one".format(selection, ", ".join(objects))) + return objects[0] + + +def _element(atom): + symbol = (atom.symbol or "").strip() + if symbol: + return symbol[0].upper() + symbol[1:].lower() + name = atom.name.strip("0123456789 ") + return name[:1].upper() if name else "X" + + +def _structure(object_name, selection="", extra=()): + """In-memory Structure of `object_name` restricted to `selection` (whole object when empty). + + Atoms of the `extra` selections (atom1, atom2, ...) are always part of the structure; when they + fall outside `selection` they are kept as zero-radius ghosts, so `selection="not resi 45"` with + atom1 in residue 45 measures the environment of the residue, like --exclude-self. + + Returns: + (Structure, list of PyMOL atom indices in structure order) + """ + kept = "({})".format(object_name) if not selection else "(({}) and ({}))".format(object_name, selection) + pieces = [kept] + ["({})".format(e) for e in extra if e] + model = cmd.get_model(" or ".join(pieces)) + if not model.atom: + raise PluginError("no atoms in '{}'".format(object_name)) + inside = {a.index for a in cmd.get_model(kept).atom} + atoms, coords, ids = [], [], [] + meta = {key: [] for key in ("record", "name", "resname", "chain", "resseq", "icode", "element", "resid")} + for a in model.atom: + element = _element(a) + atoms.append(element if a.index in inside else "Bq") + coords.append(a.coord) + ids.append(a.index) + resi = str(a.resi) + number = resi.rstrip(_ICODE_LETTERS) + meta["record"].append("HETATM" if a.hetatm else "ATOM") + meta["name"].append(a.name) + meta["resname"].append(a.resn) + meta["chain"].append(a.chain) + meta["resseq"].append(int(number) if number.lstrip("-").isdigit() else 0) + meta["icode"].append(resi[len(number):]) + meta["element"].append(element) + meta["resid"].append("{}:{}".format(a.chain, resi)) + return Structure(atoms, coords, object_name, meta), ids + + +def _positions(ids, selection): + """1-based positions (in structure order) of the atoms of a selection.""" + if not selection: + return [] + positions = [] + for _, index in cmd.index(selection): + if index not in ids: + raise PluginError("atom {} of '{}' is not part of the measured structure".format(index, selection)) + positions.append(ids.index(index) + 1) + if not positions: + raise PluginError("selection '{}' matches no atoms".format(selection)) + return positions + + +def _run(structure, **kwargs): + """Run DBSTEP on a Structure, printing its table to the PyMOL console.""" + Dbstep.dbstep._column_header_printed = False + Dbstep.dbstep._file_col_width = max(20, len(structure.name) + 2) + try: + return Dbstep.dbstep(structure, **kwargs) + except SystemExit as error: # DBSTEP reports problems through sys.exit; keep PyMOL alive + raise PluginError(str(error).strip()) from None + + +def to_original(run): + """Function mapping points of DBSTEP's aligned frame (atom1 at the origin, axis along z) back to the + coordinates of the object the run was made from.""" + rotation, translation = calculator.rigid_transform(run.coords, run.aligned_coords) + + def convert(points): + return (np.asarray(points, dtype=float).reshape(-1, 3) - translation) @ rotation + + return convert + + +# --- CGO helpers ----------------------------------------------------------------------------------------- + + +def cgo_cylinder(start, end, radius, color): + return [CYLINDER, *map(float, start), *map(float, end), float(radius), *color, *color] + + +def cgo_circle(points, color, width=3.0): + obj = [LINEWIDTH, float(width), BEGIN, LINE_LOOP, COLOR, *color] + for point in points: + obj += [VERTEX, *map(float, point)] + obj.append(END) + return obj + + +def circle_points(centre, e1, e2, radius, n=90): + angles = np.linspace(0, 2 * math.pi, n, endpoint=False) + return [np.asarray(centre) + radius * (math.cos(a) * np.asarray(e1) + math.sin(a) * np.asarray(e2)) for a in angles] + + +def sterimol_cgo(run): + """CGO list drawing L (blue), Bmin (green) and Bmax (red) of a Sterimol run on the original coordinates.""" + convert = to_original(run) + origin, z_end = convert([[0, 0, 0], [0, 0, run.L]]) + e1, e2 = convert([[1, 0, 0], [0, 1, 0]]) - origin + obj = cgo_cylinder(origin, z_end, 0.08, COLORS["L"]) + vectors = run.sterimol_vectors() or {} + for key in ("bmax", "bmin"): + value = run.Bmax if key == "bmax" else run.Bmin + if value is False: + continue + if key in vectors: + tip = convert([vectors[key]])[0] + obj += cgo_cylinder(origin, tip, 0.08, COLORS[key]) + obj += cgo_circle(circle_points(origin, e1, e2, float(value)), COLORS[key]) + return obj + + +def sphere_cgo(centre, radius, color=COLORS["sphere"], alpha=0.35): + return [ALPHA, float(alpha), COLOR, *color, SPHERE, *map(float, centre), float(radius)] + + +def cone_cgo(apex, axis, half_angle_deg, length, color=COLORS["cone"], alpha=0.35): + """Cone of half angle `half_angle_deg` around `axis` from `apex`, drawn to slant length `length` + (also for angles beyond 90 degrees, where the rim lies behind the apex).""" + axis = np.asarray(axis, dtype=float) + axis = axis / np.linalg.norm(axis) + half = math.radians(half_angle_deg) + rim_centre = np.asarray(apex, dtype=float) + length * math.cos(half) * axis + rim_radius = length * math.sin(half) + return [ALPHA, float(alpha), CONE, *map(float, apex), *map(float, rim_centre), 0.0, float(rim_radius), *color, *color, 1.0, 1.0] + + +def _load(name, obj): + cmd.delete(name) + cmd.load_cgo(obj, name) + + +# --- commands ----------------------------------------------------------------------------------------------- + + +def dbstep_style(): + """Display settings for clean steric figures (white background, soft specular light, orthoscopic view).""" + cmd.bg_color("white") + cmd.set("ray_opaque_background", 0) + cmd.set("specular", 0.25) + cmd.set("spec_power", 300) + cmd.set("spec_reflect", 0.5) + cmd.set("antialias", 1) + cmd.set("orthoscopic", 1) + cmd.set("cgo_transparency", 0.0) + cmd.set("transparency", 0.5) + + +def dbstep_vdw(object_name, radii="bondi", scale=1.0, transparency=0.6): + """Show a translucent van der Waals copy `_vdw` of an object with DBSTEP's radii (bondi, charry-tkatchenko or cpk).""" + structure, ids = _structure(object_name) + values = radii_sets.for_atoms(structure.atoms, structure.coords, str(radii).lower()) * float(scale) + name = "{}_vdw".format(object_name) + cmd.delete(name) + cmd.copy(name, object_name) + cmd.alter(name, "vdw=table.get(index, vdw)", space={"table": {index: float(r) for index, r in zip(ids, values)}}) + cmd.rebuild() + cmd.hide("everything", name) + cmd.show("spheres", name) + cmd.set("sphere_scale", 1.0, name) + cmd.set("sphere_transparency", float(transparency), name) + print(" dbstep: {} van der Waals radii ({}) shown as {}".format(len(values), radii_sets.radii_label(str(radii).lower()), name)) + return dict(zip(ids, values)) + + +def dbstep_sterimol(atom1, atom2, atom3="", radii="bondi", measure="classic", grid=0.05, scale=1.0, noH=0, selection="", name="sterimol"): + """Sterimol L, Bmin and Bmax from atom1 along the atom1-atom2 axis, drawn on the object as CGO `name`. + + `selection` restricts the measured atoms (default: the whole object of atom1); atom3 fixes the roll. + """ + obj = _object_of(atom1) + structure, ids = _structure(obj, selection, (atom1, atom2, atom3)) + i1, i2 = _positions(ids, atom1)[0], _positions(ids, atom2) + i3 = _positions(ids, atom3)[0] if atom3 else False + run = _run(structure, atom1=i1, atom2=i2, atom3=i3, sterimol=True, measure=str(measure), radii=str(radii).lower(), grid=float(grid), scalevdw=float(scale), noH=_flag(noH)) + _load(name, sterimol_cgo(run)) + print(" dbstep: L = {:.2f}, Bmin = {:.2f}, Bmax = {:.2f} Ang ({} radii) -> CGO '{}'".format(run.L, run.Bmin, run.Bmax, radii, name)) + return run + + +def dbstep_vbur(atom1, radius=3.5, radii="bondi", scale=1.0, noH=0, selection="", decompose=1, grid=0.05, name="vbur"): + """Percent buried volume in a sphere of `radius` around atom1, drawn as CGO `name`. + + `selection` restricts the atoms that can fill the sphere (e.g. "polymer and not resi 45" measures the + pocket around residue 45 without waters and without the residue itself). With `decompose` the + B-factor of every residue is set to its contribution and the contributing residues are coloured by it. + """ + obj = _object_of(atom1) + structure, ids = _structure(obj, selection, (atom1,)) + i1 = _positions(ids, atom1)[0] + run = _run(structure, atom1=i1, volume=True, r=float(radius), radii=str(radii).lower(), grid=float(grid), scalevdw=float(scale), noH=_flag(noH), cutoff="auto", decompose=_flag(decompose)) + centre = cmd.get_model(atom1).atom[0].coord + _load(name, sphere_cgo(centre, float(radius))) + if run.contributions: + cmd.alter(obj, "b=0.0") + for label, value in run.contributions.items(): + chain, resi = label.split()[0].split(":", 1) + cmd.alter("{} and chain '{}' and resi {}".format(obj, chain, resi), "b={}".format(float(value))) + cmd.spectrum("b", "white_red", "{} and b > 0".format(obj), minimum=0.0) + top = sorted(run.contributions.items(), key=lambda item: -item[1])[:5] + print(" dbstep: %V_bur = {:.2f} (R = {:.2f} Ang); top contributors: {}".format(run.bur_vol, float(radius), ", ".join("{} {:.1f}".format(k, v) for k, v in top))) + else: + print(" dbstep: %V_bur = {:.2f} (R = {:.2f} Ang) -> CGO '{}'".format(run.bur_vol, float(radius), name)) + return run + + +def dbstep_cone(metal, ligand="", radii="cpk", scale=1.0, selection="", name="cone"): + """Tolman cone angle of the ligand bound to `metal` through the `ligand` atoms (ring atoms or the donor + atom; auto-detected when empty), with the ligand's Sterimol parameters from the metal. Draws the cone, + the metal-centroid distance and the Sterimol axes.""" + obj = _object_of(metal) + structure, ids = _structure(obj, selection, (metal, ligand)) + i1 = _positions(ids, metal)[0] + i2 = _positions(ids, ligand) if ligand else None + kwargs = {"atom2": i2} if i2 else {} + run = _run(structure, atom1=i1, cone=True, radii=str(radii).lower(), scalevdw=float(scale), **kwargs) + apex = np.asarray(cmd.get_model(metal).atom[0].coord) + axis_atoms = [ids[i - 1] for i in run.atom2] + centroid = np.mean([a.coord for a in cmd.get_model("{} and index {}".format(obj, "+".join(str(i) for i in axis_atoms))).atom], axis=0) + length = max(run.metal_centroid, run.L if run.L else 0.0) + _load(name, cone_cgo(apex, centroid - apex, run.cone_angle / 2.0, length) + sterimol_cgo(run)) + cmd.pseudoatom("{}_centroid".format(name), pos=list(map(float, centroid))) + cmd.delete("{}_distance".format(name)) + cmd.distance("{}_distance".format(name), "{}_centroid".format(name), metal) + print(" dbstep: cone angle {:.2f} deg, metal-centroid {:.3f} Ang, ligand of {} atoms (L {:.2f}, Bmin {:.2f}, Bmax {:.2f}) -> CGO '{}'".format( + run.cone_angle, run.metal_centroid, len(run.ligand_atoms), run.L, run.Bmin, run.Bmax, name)) + return run + + +def _expand(files): + names = [] + for pattern in str(files).replace(";", ",").split(","): + pattern = pattern.strip() + if pattern: + matches = sorted(glob.glob(os.path.expanduser(pattern))) + names.extend(matches or [pattern]) + return names + + +def _pdb_string(atoms, coords): + lines = [] + for i, (atom, (x, y, z)) in enumerate(zip(atoms, coords), start=1): + element = "" if atom == "Bq" else atom + lines.append("HETATM{:5d} {:<4} UNK A 1 {:8.3f}{:8.3f}{:8.3f} 1.00 0.00 {:>2}".format(i % 100000, (element or "X")[:4], x, y, z, element[:2])) + return "\n".join(lines + ["END"]) + "\n" + + +def _load_run(name, path, run): + """Load a conformer into PyMOL: the file itself when PyMOL reads the format, else the measured coordinates.""" + cmd.delete(name) + try: + cmd.load(path, name) + if cmd.count_atoms(name) == 0: + raise ValueError("empty") + except Exception: + cmd.delete(name) + cmd.read_pdbstr(_pdb_string(run.atoms, run.coords), name) + + +def dbstep_ensemble(files, atom1, atom2, radii="bondi", temperature=298.15, window="", tag="", vbur=0, radius=3.5, group="ensemble"): + """Boltzmann-weighted Sterimol parameters (and %V_bur with vbur=1) over conformers: one file per + conformer (QM outputs, energies read from them) or one multi-structure file (SDF with energy fields, + multi-frame xyz). `files` is a glob or a comma-separated list; atom1/atom2 are 1-based atom indices. + The conformers are loaded into PyMOL, titled with their population and grouped under `group`.""" + names = _expand(files) + if not names: + raise PluginError("no files match '{}'".format(files)) + kwargs = {"atom1": int(atom1), "atom2": [int(a) for a in str(atom2).split("+") if a] or int(atom2), "sterimol": True, "radii": str(radii).lower(), "quiet": True} + if _flag(vbur): + kwargs.update(volume=True, r=float(radius)) + if len(names) == 1 and trajectory.count_frames(names[0]) > 1: + runs = Dbstep.all_frames(names[0], **kwargs) + else: + runs = [Dbstep.dbstep(f, **kwargs) for f in names] + summary = ensemble.boltzmann_average(runs, tag=tag or None, temperature=float(temperature), window=float(window) if str(window).strip() else None, label=group) + populations = [r.population for r in runs] + print(" {:<28} {:>10} {:>8} {:>8} {:>8} {:>8}".format("structure", "E_rel", "L", "Bmin", "Bmax", "pop")) + for r in runs: + label = r.results[0]["structure"] or r.results[0]["file"] + print(" {:<28} {:>10.2f} {:>8.2f} {:>8.2f} {:>8.2f} {:>8.3f}{}".format(label, r.energy_rel, r.L, r.Bmin, r.Bmax, r.population, "" if r.in_window else " (outside window)")) + row = summary[0] + print(" {:<28} {:>10} {:>8.2f} {:>8.2f} {:>8.2f}".format("Boltzmann average", "", row["L"], row["bmin"], row["bmax"])) + # load into PyMOL: one object per file, or one object with states for a multi-structure file + cmd.delete(group) + if len(names) == 1 and len(runs) > 1: + cmd.load(names[0], group) + for state, r in enumerate(runs, start=1): + cmd.set_title(group, state, "p = {:.3f}".format(r.population)) + else: + best = max(populations) or 1.0 + members = [] + for f, r in zip(names, runs): + member = os.path.basename(f).split(".")[0] + _load_run(member, f, r) + cmd.set_title(member, 1, "p = {:.3f}".format(r.population)) + cmd.set("stick_transparency", 0.85 * (1.0 - r.population / best), member) + members.append(member) + cmd.group(group, " ".join(members)) + return runs, summary + + +def dbstep_conformers(folder, pattern="*.pdb", group="conformers"): + """Load every structure file of a folder into PyMOL under one group (AddConformers of wSterimol).""" + files = sorted(glob.glob(os.path.join(os.path.expanduser(folder), pattern))) + if not files: + raise PluginError("no files matching {} in {}".format(pattern, folder)) + members = [] + for f in files: + member = os.path.basename(f).split(".")[0] + cmd.load(f, member) + members.append(member) + cmd.group(group, " ".join(members)) + print(" dbstep: loaded {} structures into group '{}'".format(len(members), group)) + return members + + +COMMANDS = { + "dbstep_style": dbstep_style, + "dbstep_vdw": dbstep_vdw, + "dbstep_sterimol": dbstep_sterimol, + "dbstep_vbur": dbstep_vbur, + "dbstep_cone": dbstep_cone, + "dbstep_ensemble": dbstep_ensemble, + "dbstep_conformers": dbstep_conformers, +} + + +def register(): + """Register the dbstep_* commands with PyMOL (done on import inside PyMOL).""" + if cmd is None: + raise ImportError("PyMOL is not importable from this Python; install dbstep into PyMOL's Python and import there") + for name, function in COMMANDS.items(): + cmd.extend(name, function) + cmd.auto_arg[0]["dbstep_vdw"] = [cmd.object_sc, "object", ""] + cmd.auto_arg[0]["dbstep_conformers"] = [cmd.object_sc, "folder", ""] + + +def __init_plugin__(app=None): + """Entry point for PyMOL's plugin manager.""" + register() + + +if cmd is not None: + register() + +if __name__ == "__main__": + print("DBSTEP PyMOL plugin: {}".format(os.path.abspath(__file__))) + print("In PyMOL (with dbstep installed in its Python): import dbstep.pymol_plugin") + print("or add to .pymolrc: run {}".format(os.path.abspath(__file__))) + sys.exit(0) diff --git a/dbstep/radii.py b/dbstep/radii.py index 08f5568..1fb3127 100644 --- a/dbstep/radii.py +++ b/dbstep/radii.py @@ -107,19 +107,7 @@ def radii_label(radii_set): return {"charry-tkatchenko": "Charry-Tkatchenko", "cpk": "CPK (Sterimol atom types)"}.get(radii_set, "Bondi") -def atom_radii(mol, options): - """Unscaled VDW radius of every atom of `mol` for options.radii. - - Elements without a CPK type (metals, boron, silicon, ...) fall back to their Bondi radius; - ghost atoms ("Bq") always have radius zero. - """ - atomtypes = np.asarray(mol.ATOMTYPES, dtype=str) - if options.radii != "cpk": - table = radii_table(options.radii) - return np.array([table.get(atom, 2.0) for atom in atomtypes], dtype=float) - types = getattr(mol, "METADATA", {}).get(CPK_TYPE_KEY) - if types is None: # structure parsed without the CPK hook (e.g. built by hand): type it now - types = ["" if t is None else t for t in cpk_types(atomtypes, mol.CARTESIANS)] +def _cpk_radii(atomtypes, types): radii = [] for atom, cpk_t in zip(atomtypes, types): if atom == "Bq": @@ -131,6 +119,29 @@ def atom_radii(mol, options): return np.array(radii, dtype=float) +def for_atoms(atomtypes, coords, radii_set="bondi"): + """Unscaled VDW radius of every atom of a complete structure for the given radii set + (CPK types are derived from the coordinates).""" + atomtypes = np.asarray(atomtypes, dtype=str) + if radii_set != "cpk": + table = radii_table(radii_set) + return np.array([table.get(atom, 2.0) for atom in atomtypes], dtype=float) + return _cpk_radii(atomtypes, ["" if t is None else t for t in cpk_types(atomtypes, coords)]) + + +def atom_radii(mol, options): + """Unscaled VDW radius of every atom of `mol` for options.radii. + + Elements without a CPK type (metals, boron, silicon, ...) fall back to their Bondi radius; + ghost atoms ("Bq") always have radius zero. + """ + atomtypes = np.asarray(mol.ATOMTYPES, dtype=str) + types = getattr(mol, "METADATA", {}).get(CPK_TYPE_KEY) if options.radii == "cpk" else None + if types is None: # not CPK, or a structure parsed without the CPK hook (e.g. built by hand) + return for_atoms(atomtypes, mol.CARTESIANS, options.radii) + return _cpk_radii(atomtypes, types) + + def max_radius(atomtypes, options): """Largest (scaled) radius any of the given atom types can be assigned under options.radii.""" if len(atomtypes) == 0: diff --git a/meta.yaml b/meta.yaml index 5714c7d..1aae167 100644 --- a/meta.yaml +++ b/meta.yaml @@ -1,5 +1,5 @@ {% set name = "dbstep" %} -{% set version = "2.1.0" %} +{% set version = "2.2.0" %} # NOTE: Update version above to match dbstep/__version__ on each release. # The sha256 hash must also be updated from the PyPI sdist. diff --git a/pyproject.toml b/pyproject.toml index d4e9d38..7b594be 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,6 +35,8 @@ dependencies = [ graph2d = ["rdkit", "pandas"] dev = ["pytest", "ruff", "pre-commit"] test = ["pytest"] +# open-source PyMOL from PyPI (pre-release wheels, no 3.14 wheels yet) for the PyMOL plugin and its tests +pymol = ["pymol-open-source>=3.1.0a0; python_version < '3.14'"] [project.urls] Homepage = "https://github.com/patonlab/DBSTEP" diff --git a/tests/test_from_coords.py b/tests/test_from_coords.py new file mode 100644 index 0000000..bae0c94 --- /dev/null +++ b/tests/test_from_coords.py @@ -0,0 +1,56 @@ +"""dbstep on in-memory structures (from_coords / Structure) and the alignment transform.""" + +import numpy as np +import pytest + +from dbstep import Dbstep, calculator, parse_data + + +def read_xyz(path): + lines = open(path).read().splitlines() + atoms = [line.split()[0] for line in lines[2:] if line.strip()] + coords = np.array([[float(v) for v in line.split()[1:4]] for line in lines[2:] if line.strip()]) + return atoms, coords + + +def test_from_coords_matches_file_run(): + atoms, coords = read_xyz("dbstep/data/Et.xyz") + ref = Dbstep.dbstep("dbstep/data/Et.xyz", atom1=2, atom2=5, sterimol=True, volume=True, quiet=True) + mem = Dbstep.from_coords(atoms, coords, name="ethane", atom1=2, atom2=5, sterimol=True, volume=True, quiet=True) + assert (mem.L, mem.Bmin, mem.Bmax, mem.bur_vol) == pytest.approx((ref.L, ref.Bmin, ref.Bmax, ref.bur_vol)) + assert mem.results[0]["file"] == "ethane" and mem.results[0]["path"] == "" + assert mem.n_atoms_total == len(atoms) + + +def test_structure_with_metadata_supports_decompose(): + mol = parse_data.PDBParser("tests/pdb_files/ala5.pdb", "pdb", False, False, 1, [1]) + structure = parse_data.Structure(mol.ATOMTYPES, mol.CARTESIANS, "ala5", mol.METADATA) + ca3 = int(np.flatnonzero((mol.METADATA["resid"] == "A:3") & (mol.METADATA["name"] == "CA"))[0]) + 1 + mem = Dbstep.dbstep(structure, atom1=ca3, volume=True, decompose=True, cutoff="auto", quiet=True) + ref = Dbstep.dbstep("tests/pdb_files/ala5.pdb", residue="A:3", volume=True, decompose=True, quiet=True) + assert mem.bur_vol == pytest.approx(ref.bur_vol, abs=1e-9) + assert mem.contributions == pytest.approx(ref.contributions, abs=1e-9) + + +def test_structure_validation_and_ghosts(): + with pytest.raises(ValueError): + parse_data.Structure(["C", "H"], [[0, 0, 0]]) + # a ghost atom1 occupies nothing but still defines the origin + mol = Dbstep.from_coords(["Bq", "C", "C"], [[0, 0, 0], [0, 0, 1.5], [0, 0, 3.0]], atom1=1, atom2=2, sterimol=True, quiet=True) + assert mol.L == pytest.approx(3.0 + 1.70, abs=1e-6) + + +def test_rigid_transform_recovers_the_alignment(): + atoms, coords = read_xyz("tests/metal_files/RhCpMe5Cl2PMe3.xyz") + run = Dbstep.from_coords(atoms, coords, atom1=17, atom2=[18, 20, 37], atom3=18, sterimol=True, quiet=True) + rotation, translation = calculator.rigid_transform(run.coords, run.aligned_coords) + assert np.allclose(run.coords @ rotation.T + translation, run.aligned_coords, atol=1e-8) + assert np.allclose(rotation @ rotation.T, np.eye(3), atol=1e-10) and np.linalg.det(rotation) == pytest.approx(1.0) + # the aligned origin is atom1 and the aligned z axis points from atom1 towards atom2 + back = (np.array([[0, 0, 0], [0, 0, 1.0]]) - translation) @ rotation + assert np.allclose(back[0], coords[16], atol=1e-8) + direction = coords[[17, 19, 36]].mean(axis=0) - coords[16] + assert np.dot(back[1] - back[0], direction / np.linalg.norm(direction)) == pytest.approx(1.0, abs=1e-8) + vectors = run.sterimol_vectors() + assert set(vectors) == {"bmax", "bmin"} + assert np.hypot(*vectors["bmax"][:2]) == pytest.approx(run.Bmax, abs=2e-3) and np.hypot(*vectors["bmin"][:2]) == pytest.approx(run.Bmin, abs=2e-3) diff --git a/tests/test_pymol_plugin.py b/tests/test_pymol_plugin.py new file mode 100644 index 0000000..157abd6 --- /dev/null +++ b/tests/test_pymol_plugin.py @@ -0,0 +1,110 @@ +"""The PyMOL plugin, run against headless open-source PyMOL (skipped when PyMOL is not installed).""" + +import numpy as np +import pytest + +pymol = pytest.importorskip("pymol") +from pymol import cmd # noqa: E402 + +from dbstep import Dbstep, pymol_plugin as plugin # noqa: E402 + +RH = "tests/metal_files/RhCpMe5Cl2PMe3.xyz" +ALA5 = "tests/pdb_files/ala5.pdb" +ETHER = "tests/sdf_files/ether_conformers.sdf" + + +@pytest.fixture(autouse=True) +def fresh_session(): + cmd.reinitialize() + yield + cmd.reinitialize() + + +def test_commands_are_registered(): + for name in plugin.COMMANDS: + assert name in cmd.keyword + + +def test_sterimol_matches_cli_and_draws_on_the_object(): + cmd.load("dbstep/data/Et.xyz", "et") + run = plugin.dbstep_sterimol("et and id 2", "et and id 5") + ref = Dbstep.dbstep("dbstep/data/Et.xyz", atom1=2, atom2=5, sterimol=True, quiet=True) + assert (run.L, run.Bmin, run.Bmax) == pytest.approx((ref.L, ref.Bmin, ref.Bmax)) + assert "sterimol" in cmd.get_names("all") + # the L cylinder starts at atom1 and points towards atom2 in the object's own frame + obj = plugin.sterimol_cgo(run) + start, end = np.array(obj[1:4]), np.array(obj[4:7]) + xyz = cmd.get_coords("et") + assert np.allclose(start, xyz[1], atol=1e-5) + direction = xyz[4] - xyz[1] + assert np.dot(end - start, direction / np.linalg.norm(direction)) == pytest.approx(run.L, abs=1e-5) + # grid mode and a restricting selection also work + grid = plugin.dbstep_sterimol("et and id 2", "et and id 5", measure="grid", grid=0.1, selection="not id 1", name="grid_sterimol") + assert grid.n_atoms_total == 7 and "grid_sterimol" in cmd.get_names("all") + + +def test_vbur_matches_residue_mode_and_colours_contributions(): + cmd.load(ALA5, "ala") + run = plugin.dbstep_vbur("/ala//A/3/CA", grid=0.1) + ref = Dbstep.dbstep(ALA5, residue="A:3", volume=True, decompose=True, grid=0.1, quiet=True) + assert run.bur_vol == pytest.approx(ref.bur_vol, abs=1e-9) + assert run.contributions == pytest.approx(ref.contributions, abs=1e-9) + assert cmd.get_model("ala and resi 3 and name CA").atom[0].b == pytest.approx(ref.contributions["A:3 ALA"]) + assert cmd.get_model("ala and resi 101").atom[0].b == pytest.approx(ref.contributions["A:101 HOH"]) + assert "vbur" in cmd.get_names("all") + # excluding the residue through the selection reproduces --exclude-self, waters through "polymer" + env = plugin.dbstep_vbur("/ala//A/3/CA", grid=0.1, selection="polymer and not resi 3", decompose=0) + ref_env = Dbstep.dbstep(ALA5, residue="A:3", volume=True, exclude_self=True, nowater=True, nohet=True, grid=0.1, quiet=True) + assert env.bur_vol == pytest.approx(ref_env.bur_vol, abs=1e-9) + assert list(env.atoms).count("Bq") == 1 + + +def test_cone_reproduces_the_legacy_value(): + cmd.load(RH, "rh") + run = plugin.dbstep_cone("rh and elem Rh") + assert run.cone_angle == pytest.approx(173.97, abs=0.01) and run.metal_centroid == pytest.approx(1.833, abs=0.001) + assert {"cone", "cone_centroid", "cone_distance"} <= set(cmd.get_names("all")) + assert np.allclose(cmd.get_coords("cone_centroid")[0], cmd.get_coords("rh and index 3+4+5+24+25").mean(axis=0), atol=1e-4) + phosphine = plugin.dbstep_cone("rh and elem Rh", "rh and elem P", name="pme3") + assert 110 < phosphine.cone_angle < 120 and len(phosphine.ligand_atoms) == 13 + obj = plugin.cone_cgo([0, 0, 0], [0, 0, 1], 60.0, 2.0) + assert obj[6:9] == pytest.approx([0, 0, 1.0]) and obj[10] == pytest.approx(2.0 * np.sin(np.radians(60))) and obj[14:17] == pytest.approx([1.0, 0.7, 0.15]) + + +def test_vdw_copy_carries_cpk_radii(): + cmd.load(RH, "rh") + table = plugin.dbstep_vdw("rh", radii="cpk") + assert "rh_vdw" in cmd.get_object_list() and cmd.count_atoms("rh_vdw") == 41 + radii = {a.index: a.vdw for a in cmd.get_model("rh_vdw").atom} + assert radii[6] == pytest.approx(1.50, abs=1e-3) # methyl carbon (sp3) + assert radii[7] == pytest.approx(1.00, abs=1e-3) and table[7] == pytest.approx(1.0) # H + assert radii[2] == pytest.approx(1.80, abs=1e-3) and radii[17] == pytest.approx(1.40, abs=1e-3) # Cl, P + assert radii[1] == pytest.approx(2.0, abs=1e-3) # Rh has no CPK type: Bondi fallback + plugin.dbstep_vdw("rh", radii="bondi", scale=1.17) + assert {a.index: a.vdw for a in cmd.get_model("rh_vdw").atom}[7] == pytest.approx(1.09 * 1.17, abs=1e-3) + + +def test_ensemble_over_files_and_over_a_multi_structure_file(): + runs, summary = plugin.dbstep_ensemble("tests/qm_files/pentane_*.out.gz", 1, 3, radii="cpk", temperature=298, window=1.0) + assert len(runs) == 9 and summary[0]["L"] + 0.4 == pytest.approx(6.34, abs=0.015) + assert sum(not r.in_window for r in runs) == 2 + assert "ensemble" in cmd.get_names("group_objects") and cmd.count_atoms("pentane_18") == 17 + assert cmd.get_title("pentane_18", 1).startswith("p = 0.27") + runs, summary = plugin.dbstep_ensemble(ETHER, 3, 2, group="ether") + assert len(runs) == 3 and cmd.count_states("ether") == 3 and cmd.get_title("ether", 1).startswith("p = 0.92") + + +def test_conformers_and_style_and_errors(tmp_path): + for i in range(3): + (tmp_path / f"c{i}.xyz").write_text(open("dbstep/data/Et.xyz").read()) + members = plugin.dbstep_conformers(str(tmp_path), pattern="*.xyz") + assert members == ["c0", "c1", "c2"] and "conformers" in cmd.get_names("group_objects") + plugin.dbstep_style() + assert cmd.get_setting_int("orthoscopic") == 1 + with pytest.raises(plugin.PluginError, match="matches no atoms"): + plugin.dbstep_sterimol("c0 and id 99", "c0 and id 2") + with pytest.raises(plugin.PluginError, match="must not include"): + plugin.dbstep_cone("c0 and id 1", "c0 and id 1") + cmd.load("dbstep/data/Me.xyz", "me") + with pytest.raises(plugin.PluginError, match="several objects"): + plugin.dbstep_sterimol("c0 or me", "c0 and id 2") diff --git a/uv.lock b/uv.lock index 902f51c..0d4297a 100644 --- a/uv.lock +++ b/uv.lock @@ -54,7 +54,7 @@ resolution-markers = [ "python_full_version < '3.11'", ] dependencies = [ - { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" } }, + { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, ] sdist = { url = "https://files.pythonhosted.org/packages/66/54/eb9bfc647b19f2009dd5c7f5ec51c4e6ca831725f1aea7a993034f483147/contourpy-1.3.2.tar.gz", hash = "sha256:b6945942715a034c671b7fc54f9588126b0b8bf23db2696e3ca8328f3ff0ab54", size = 13466130, upload-time = "2025-04-15T17:47:53.79Z" } wheels = [ @@ -129,7 +129,7 @@ resolution-markers = [ "python_full_version >= '3.11' and python_full_version < '3.14' and sys_platform != 'emscripten' and sys_platform != 'win32'", ] dependencies = [ - { name = "numpy", version = "2.4.2", source = { registry = "https://pypi.org/simple" } }, + { name = "numpy", version = "2.4.2", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, ] sdist = { url = "https://files.pythonhosted.org/packages/58/01/1253e6698a07380cd31a736d248a3f2a50a7c88779a1813da27503cadc2a/contourpy-1.3.3.tar.gz", hash = "sha256:083e12155b210502d0bca491432bb04d56dc3432f95a979b429f2848c3dbe880", size = 13466174, upload-time = "2025-07-26T12:03:12.549Z" } wheels = [ @@ -237,6 +237,9 @@ graph2d = [ { name = "pandas", version = "3.0.1", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, { name = "rdkit" }, ] +pymol = [ + { name = "pymol-open-source", marker = "python_full_version < '3.14'" }, +] test = [ { name = "pytest" }, ] @@ -256,13 +259,14 @@ requires-dist = [ { name = "numpy" }, { name = "pandas", marker = "extra == 'graph2d'" }, { name = "pre-commit", marker = "extra == 'dev'" }, + { name = "pymol-open-source", marker = "python_full_version < '3.14' and extra == 'pymol'", specifier = ">=3.1.0a0" }, { name = "pytest", marker = "extra == 'dev'" }, { name = "pytest", marker = "extra == 'test'" }, { name = "rdkit", marker = "extra == 'graph2d'" }, { name = "ruff", marker = "extra == 'dev'" }, { name = "scipy" }, ] -provides-extras = ["graph2d", "dev", "test"] +provides-extras = ["graph2d", "dev", "test", "pymol"] [package.metadata.requires-dev] analysis = [ @@ -286,7 +290,7 @@ name = "exceptiongroup" version = "1.3.1" source = { registry = "https://pypi.org/simple" } dependencies = [ - { name = "typing-extensions" }, + { name = "typing-extensions", marker = "python_full_version < '3.11'" }, ] sdist = { url = "https://files.pythonhosted.org/packages/50/79/66800aadf48771f6b62f7eb014e352e5d06856655206165d775e675a02c9/exceptiongroup-1.3.1.tar.gz", hash = "sha256:8b412432c6055b0b7d14c310000ae93352ed6754f70fa8f7c34141f91c4e3219", size = 30371, upload-time = "2025-11-21T23:01:54.787Z" } wheels = [ @@ -738,10 +742,10 @@ resolution-markers = [ "python_full_version < '3.11'", ] dependencies = [ - { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" } }, - { name = "python-dateutil" }, - { name = "pytz" }, - { name = "tzdata" }, + { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "python-dateutil", marker = "python_full_version < '3.11'" }, + { name = "pytz", marker = "python_full_version < '3.11'" }, + { name = "tzdata", marker = "python_full_version < '3.11'" }, ] sdist = { url = "https://files.pythonhosted.org/packages/33/01/d40b85317f86cf08d853a4f495195c73815fdf205eef3993821720274518/pandas-2.3.3.tar.gz", hash = "sha256:e05e1af93b977f7eafa636d043f9f94c7ee3ac81af99c13508215942e64c993b", size = 4495223, upload-time = "2025-09-29T23:34:51.853Z" } wheels = [ @@ -807,9 +811,9 @@ resolution-markers = [ "python_full_version >= '3.11' and python_full_version < '3.14' and sys_platform != 'emscripten' and sys_platform != 'win32'", ] dependencies = [ - { name = "numpy", version = "2.4.2", source = { registry = "https://pypi.org/simple" } }, - { name = "python-dateutil" }, - { name = "tzdata", marker = "sys_platform == 'emscripten' or sys_platform == 'win32'" }, + { name = "numpy", version = "2.4.2", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, + { name = "python-dateutil", marker = "python_full_version >= '3.11'" }, + { name = "tzdata", marker = "(python_full_version >= '3.11' and sys_platform == 'emscripten') or (python_full_version >= '3.11' and sys_platform == 'win32')" }, ] sdist = { url = "https://files.pythonhosted.org/packages/2e/0c/b28ed414f080ee0ad153f848586d61d1878f91689950f037f976ce15f6c8/pandas-3.0.1.tar.gz", hash = "sha256:4186a699674af418f655dbd420ed87f50d56b4cd6603784279d9eef6627823c8", size = 4641901, upload-time = "2026-02-17T22:20:16.434Z" } wheels = [ @@ -1012,6 +1016,33 @@ wheels = [ { url = "https://files.pythonhosted.org/packages/71/46/17f022dd3e953bf20a04a028a21ec746d942f8d2af30fa0f124fa0e6a684/pygments-2.21.0-py3-none-any.whl", hash = "sha256:2363c69b61c4a97c838da3b130dcd6468f4848992b21a82f2a63ec34377137d9", size = 1250147, upload-time = "2026-08-17T08:02:44.912Z" }, ] +[[package]] +name = "pymol-open-source" +version = "3.2.0a0" +source = { registry = "https://pypi.org/simple" } +dependencies = [ + { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, + { name = "numpy", version = "2.4.2", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11' and python_full_version < '3.14'" }, +] +wheels = [ + { url = "https://files.pythonhosted.org/packages/c6/06/ac148be6b62841a1781dc84d1996b137b58ef4817b9db11f57283daaa877/pymol_open_source-3.2.0a0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:615e862f169f6f17b66d20f9edab521abb013de9b1142035fb3d0e63956891bd", size = 7359803, upload-time = "2025-09-09T00:51:04.277Z" }, + { url = "https://files.pythonhosted.org/packages/52/72/804ef7484a87051d437072ca0c930075d32711a0950d1ba9b55888a4dad3/pymol_open_source-3.2.0a0-cp310-cp310-macosx_12_0_x86_64.whl", hash = "sha256:c4d195ff5bfcca8db70979adffb385cc8e97e8117ec13800e93f3822004a6f4f", size = 7626827, upload-time = "2025-08-04T14:39:42.851Z" }, + { url = "https://files.pythonhosted.org/packages/b5/6d/976be1405c1ead4514ddcc2b5306fe64e3630f5dec1dc2cc7e4e85276b2f/pymol_open_source-3.2.0a0-cp310-cp310-manylinux_2_35_x86_64.whl", hash = "sha256:878bcaba7ff18b8bf95707694379d3b85b0b1cd655ef5895c26babb6d2d38954", size = 31358062, upload-time = "2025-08-04T14:39:45.219Z" }, + { url = "https://files.pythonhosted.org/packages/d0/a4/8a369cc7e315b4e48026492254f8251414e255a5e6d8bd9d79481296c2a0/pymol_open_source-3.2.0a0-cp310-cp310-win_amd64.whl", hash = "sha256:302412feef55a4b921220c914c28c6e764e7ab9a02d8b47991a58e658fceb12b", size = 7278967, upload-time = "2025-08-04T14:39:47.888Z" }, + { url = "https://files.pythonhosted.org/packages/52/da/52b5772ecd7a002d9e2072dc5fb31406b2534e58d272e4551e5326d4dc67/pymol_open_source-3.2.0a0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:d7f0f24fdf15a4da37a6b196897ed3fef0ffca36d5404d07a79ae532202286b2", size = 7359805, upload-time = "2025-09-09T00:51:06.108Z" }, + { url = "https://files.pythonhosted.org/packages/7c/4e/617d3e2177e61c68b5e8cf2644ed474fd02f9828448e764003fbba797a59/pymol_open_source-3.2.0a0-cp311-cp311-macosx_12_0_x86_64.whl", hash = "sha256:d694136e0d90a9bbaa702dd412a09e23b84b5819fdfe69db098a177b6d9962b6", size = 7626821, upload-time = "2025-08-04T14:39:49.804Z" }, + { url = "https://files.pythonhosted.org/packages/0c/e8/5d5d400bd04264f2114c76b1843ce0fd17b3a50706d1862c4dec0920c691/pymol_open_source-3.2.0a0-cp311-cp311-manylinux_2_35_x86_64.whl", hash = "sha256:2cb05c1f0584014bf2fa1943388fee297f1e9953b2cdc80945735650822b5164", size = 31358077, upload-time = "2025-08-04T14:39:52.714Z" }, + { url = "https://files.pythonhosted.org/packages/90/81/0bbea8ef94f2b1421f93be624bd69b538cb2bcaa3b7394857cd1284663ae/pymol_open_source-3.2.0a0-cp311-cp311-win_amd64.whl", hash = "sha256:b285e1c7c6165e57f0a137e3dfebd29074d48f14b28d9e386ca8e5dfbe98b8c2", size = 7279013, upload-time = "2025-08-04T14:39:55.053Z" }, + { url = "https://files.pythonhosted.org/packages/98/bc/de878d01a4658d7c7ca986eb2114de2744bcbb0c9452ab554420fcd49877/pymol_open_source-3.2.0a0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:8246536c971c6d6ed426b6b39ce248d5df3635f7954e13044c0ce2743be7d89a", size = 7359871, upload-time = "2025-09-09T00:51:07.943Z" }, + { url = "https://files.pythonhosted.org/packages/c8/8d/718ab51af1b5a8500f2daa3aac353d1da324cb861f58cb8f7707b7517b3f/pymol_open_source-3.2.0a0-cp312-cp312-macosx_12_0_x86_64.whl", hash = "sha256:f7b0e58d5728abfd10c4abf7db88da6dee80a3929b051a3b5eececc717c258b5", size = 7628150, upload-time = "2025-08-04T14:39:56.976Z" }, + { url = "https://files.pythonhosted.org/packages/04/df/82cb271593ce617a6355d330c31500a2f747eca4fadd0e5466ff3a69ae69/pymol_open_source-3.2.0a0-cp312-cp312-manylinux_2_35_x86_64.whl", hash = "sha256:6702a33ca3f23197d1ab026ecbdcc218472c9fcfe92782fdb146efcfc1c82601", size = 31358340, upload-time = "2025-08-04T14:39:59.009Z" }, + { url = "https://files.pythonhosted.org/packages/c5/75/8f07919b9a6e2d374711e33877395f371aafe8d283691969084f54f312fd/pymol_open_source-3.2.0a0-cp312-cp312-win_amd64.whl", hash = "sha256:b3d286c61efe65a8ee6c9974cc2aaeffc7250da3d2616f4cc91bcc824fa05a2a", size = 7279550, upload-time = "2025-08-04T14:40:02.12Z" }, + { url = "https://files.pythonhosted.org/packages/bc/74/536c009dd4a4b054f9485e8f7e3d4b70b143000ed2e3f15d4960f96151c8/pymol_open_source-3.2.0a0-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:f87fd8d0ab9bf15c3c6a63c4e0cdd3b73142d1ac775a1e38cd82f926a07b4465", size = 7359861, upload-time = "2025-09-09T00:51:09.575Z" }, + { url = "https://files.pythonhosted.org/packages/ee/54/491e3b78883647f5c97439d65b52e21c705f68d5fc88613f5f1c90c723d9/pymol_open_source-3.2.0a0-cp313-cp313-macosx_12_0_x86_64.whl", hash = "sha256:6dcd27f0cba0a3557de4672f7d01ecb75f603ebf788180fc34f54c062d8c8616", size = 7628177, upload-time = "2025-08-04T14:40:03.722Z" }, + { url = "https://files.pythonhosted.org/packages/d9/58/d81900193b22758314a74190cf4ea313c055b759e8a66013f4e8bec24622/pymol_open_source-3.2.0a0-cp313-cp313-manylinux_2_35_x86_64.whl", hash = "sha256:7fab3a1d46be170b98e93680a006758a569196eca7d7c920fe9e096b27f4566e", size = 31358095, upload-time = "2025-08-04T14:40:05.609Z" }, + { url = "https://files.pythonhosted.org/packages/59/7c/f6fffe36af45cfef497590ef4d8ab3335e83dc9a8cd289b9c7b39482c6be/pymol_open_source-3.2.0a0-cp313-cp313-win_amd64.whl", hash = "sha256:e080477c2381590f635ce90a29b82ce3be0d4e4cbad83bd66123313c7c7f608d", size = 7279595, upload-time = "2025-08-04T14:40:08.03Z" }, +] + [[package]] name = "pyparsing" version = "3.3.2" @@ -1191,7 +1222,7 @@ resolution-markers = [ "python_full_version < '3.11'", ] dependencies = [ - { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" } }, + { name = "numpy", version = "2.2.6", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version < '3.11'" }, ] sdist = { url = "https://files.pythonhosted.org/packages/0f/37/6964b830433e654ec7485e45a00fc9a27cf868d622838f6b6d9c5ec0d532/scipy-1.15.3.tar.gz", hash = "sha256:eae3cf522bc7df64b42cad3925c876e1b0b6c35c1337c93e12c0f366f55b0eaf", size = 59419214, upload-time = "2025-05-08T16:13:05.955Z" } wheels = [ @@ -1255,7 +1286,7 @@ resolution-markers = [ "python_full_version >= '3.11' and python_full_version < '3.14' and sys_platform != 'emscripten' and sys_platform != 'win32'", ] dependencies = [ - { name = "numpy", version = "2.4.2", source = { registry = "https://pypi.org/simple" } }, + { name = "numpy", version = "2.4.2", source = { registry = "https://pypi.org/simple" }, marker = "python_full_version >= '3.11'" }, ] sdist = { url = "https://files.pythonhosted.org/packages/7a/97/5a3609c4f8d58b039179648e62dd220f89864f56f7357f5d4f45c29eb2cc/scipy-1.17.1.tar.gz", hash = "sha256:95d8e012d8cb8816c226aef832200b1d45109ed4464303e997c5b13122b297c0", size = 30573822, upload-time = "2026-02-23T00:26:24.851Z" } wheels = [