PyMOL plugin, in-memory structure input, and the 2.2.0 release

.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:

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

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

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`.
 

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().
 

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

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

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