Add trajectory frames: --frames, per-frame runs and CSV time series (proteins/trajectories step 5)

CLAUDE.md

@@ -12,6 +12,7 @@ DBSTEP (DFT-Based Steric Parameters) is a Python package for computing steric pa
   - `calculator.py` — Math/geometry routines (rotations, angles)
   - `sterics.py` — Steric parameter calculations
   - `selection.py` — Atom selection before measurement: radial crop (`--cutoff`) and PDB residue selection (`--residue`, water/het/self filters)
+  - `trajectory.py` — Frames of multi-structure files (`--frames` selection, frame counting)
   - `parse_data.py` — Input file parsing (xyz, sdf/mol, pdb with residue metadata, cube, cclib-supported formats)
   - `constants.py` — Chemical constants (periodic table, Bondi radii, metals)
   - `graph.py` — 2D graph-based steric contribution calculations
@@ -29,8 +30,9 @@ DBSTEP (DFT-Based Steric Parameters) is a Python package for computing steric pa
   - `test_pdb_parser.py` — PDB parsing: columns, element inference, metadata, altlocs, MODEL blocks
   - `test_protein.py` — `--residue` selection: exact equivalence with the XYZ path, crop invariance, water/het/self semantics
   - `test_residue_all.py` — `--residue all`, per-run `results` records and `--csv` output
+  - `test_trajectory.py` — `--frames` parsing, per-frame runs on `ala5_traj.pdb`, multi-frame xyz, CSV time series
   - `cube_files/` — Test cube file fixtures (keep these small; use `benzene_coarse.cube` / `Ne_medium.cube` for new tests)
-  - `pdb_files/` — PDB fixtures: `1a8o.pdb` (real, no H, waters, MSE) and generated `ala5.pdb` (with H, waters, Na); see its README
+  - `pdb_files/` — PDB fixtures: `1a8o.pdb` (real, no H, waters, MSE) and generated `ala5.pdb` (with H, waters, Na) and `ala5_traj.pdb` (10 models, water approaching A:3); see its README
 - `examples/` — Jupyter notebook examples
 - `analysis/` — Standalone analysis scripts and data behind the paper (not part of the package; linted by ruff; extra deps via `uv sync --group analysis`)
 - `reference/` — Reference data

README.md

@@ -103,7 +103,17 @@ With a PDB file you can pick a residue and atoms by name instead of file indices
 * Crystal structures usually lack hydrogens; results differ from a protonated model. `--sambvca` (heavy atoms only, Bondi × 1.17) is a consistent choice for raw PDB files
 
 * `--residue all` runs every polymer residue in turn (waters and ligands are skipped, modified residues such as MSE are included; combine with `--chain`)
-* `--csv results.csv` writes one row per file, structure, residue and radius with the columns `file, structure, residue, atom1, atom2, radius, mol_vol, percent_vbur, percent_sbur, bmin, bmax, L`; this works for any input, not only PDB files
+* `--csv results.csv` writes one row per file, frame, residue and radius with the columns `file, frame, structure, residue, atom1, atom2, radius, mol_vol, percent_vbur, percent_sbur, bmin, bmax, L`; this works for any input, not only PDB files
+
+### Trajectories
+
+Multi-frame `.xyz`, multi-record `.sdf` and multi-MODEL `.pdb` files are trajectories: every frame is measured in turn (the neighbourhood crop is recomputed per frame, so frames cost the same as single structures). `--frames start:stop:stride` selects frames with Python slice rules on the 0-based index, e.g. `--frames 0:1000:10` or `--frames ::5`, and `--csv` collects the time series:
+
+```
+dbstep md_frames.pdb --residue A:45 --vbur --nowater --frames ::10 --csv vbur_A45.csv
+```
+
+From Python, `db.all_frames("md_frames.pdb", frames="::10", residue="A:45", volume=True, nowater=True)` returns one object per frame; the `frame` and `structure` columns of `results` identify it. Binary trajectory formats (DCD, XTC, ...) are not read yet; convert to multi-MODEL PDB or multi-frame XYZ first.
 
 From Python: `db.dbstep("1a8o.pdb", residue="A:186", atom="CA", atom2="CB", volume=True, nowater=True, exclude_self=True)`; the object also exposes `atoms`, `coords` and `metadata` for the atoms that were actually measured, and `results`, the list of row dictionaries that `--csv` writes. `db.all_residues("1a8o.pdb", volume=True, nowater=True)` returns one such object per residue.
 

dbstep/Dbstep.py

@@ -7,7 +7,7 @@
 import numpy as np
 from optparse import OptionParser
 
-from dbstep import sterics, parse_data, calculator, writer, selection
+from dbstep import sterics, parse_data, calculator, writer, selection, trajectory
 from dbstep.constants import periodic_table, bondi, charry_tkatchenko, metals
 
 class dbstep:
@@ -20,8 +20,8 @@ class dbstep:
 			atom1, atom2 (the reference atom indices as given in the input file),
 			cutoff, n_atoms_total, n_atoms_kept (radial crop, see --cutoff),
 			residue_label (PDB residue mode, e.g. "A:45 LEU"),
-			results (list of dicts, one per radius: file, structure, residue, atom1, atom2, radius,
-				mol_vol, percent_vbur, percent_sbur, bmin, bmax, L; also what --csv writes),
+			results (list of dicts, one per radius: file, frame, structure, residue, atom1, atom2,
+				radius, mol_vol, percent_vbur, percent_sbur, bmin, bmax, L; also what --csv writes),
 			atoms, coords, spec_atoms (the structure actually measured, after selection and crop,
 				in input orientation; spec_atoms are 1-indexed into it)
 
@@ -407,6 +407,7 @@ def _compute(self, mol, file, options, origin, occ_grid, occ_vol, point_tree, gr
 			# Record the result row (also written by --csv)
 			self.results.append({
 				"file": plain_name,
+				"frame": structure_idx if structure_idx is not None else "",
 				"structure": structure_label,
 				"residue": self.residue_label or "",
 				"atom1": self.atom1,
@@ -558,6 +559,7 @@ def set_options(kwargs):
 		"self_only": ["self_only", False],
 		"chain": ["chain", False],
 		"csv": ["csv", False],
+		"frames": ["frames", False],
 		"measure": ["measure", "classic"],
 		"pos": ["pos", False],
 		"dp": ["dp", 2],
@@ -613,6 +615,41 @@ def all_residues(file, **kwargs):
 	return runs
 
 
+def run_file(file, options):
+	"""All runs for one input file: every selected frame (see --frames) times every residue for --residue all.
+
+	Returns:
+		list of dbstep objects in frame, then residue, order
+	"""
+	runs = []
+	previous = options.structure
+	try:
+		for frame in trajectory.frame_indices(file, options):
+			options.structure = frame
+			if options.residue and str(options.residue).lower() == "all":
+				runs.extend(all_residues(file, options=options))
+			else:
+				runs.append(dbstep(file, options=options))
+	finally:
+		options.structure = previous
+	return runs
+
+
+def all_frames(file, frames=None, **kwargs):
+	"""Run every selected frame of a multi-structure file (multi-frame xyz, multi-record sdf, multi-MODEL pdb).
+
+	`frames` follows --frames ('start:stop:stride', 0-based Python slice rules, or a list of indices;
+	default all). Accepts the same keyword arguments as dbstep(), including residue="all".
+
+	Returns:
+		list of dbstep objects, one per frame (times residues for residue="all")
+	"""
+	options = kwargs["options"] if "options" in kwargs else set_options(kwargs)
+	if frames is not None:
+		options.frames = frames
+	return run_file(file, options)
+
+
 def main():
 	files = []
 	# get command line inputs. Use -h to list all possible arguments and default values
@@ -635,6 +672,7 @@ def main():
 	parser.add_option("--chain", dest="chain", action="store", help="[PDB] Keep only this chain (plus the selected residue)", default=False, metavar="chain")
 	parser.add_option("--exclude-self", dest="exclude_self", action="store_true", help="[PDB] The selected residue occupies no volume: measure only its environment", default=False)
 	parser.add_option("--self-only", dest="self_only", action="store_true", help="[PDB] Keep only the selected residue (same as measuring it extracted to its own file)", default=False)
+	parser.add_option("--frames", dest="frames", action="store", help="Frames of a multi-structure file to run, 0-based with Python slice rules: start:stop:stride, e.g. 0:1000:10, ::5, or a single index (default: all)", default=False, metavar="frames")
 	parser.add_option("--csv", dest="csv", action="store", help="Write all result rows (one per file/structure/residue/radius) to this CSV file", default=False, metavar="file")
 	parser.add_option("--cutoff", dest="cutoff", action="store", help="Ignore atoms farther than this distance (Angstrom) from atom1; 'auto' keeps exactly the atoms that can occupy the buried-volume sphere. Keeps the grid small for large systems (default: off)", default=False, metavar="cutoff")
 	parser.add_option("--isoval", dest="isoval", action="store", help="Density isovalue cutoff (default: 0.0016)", type="float", default=0.0016, metavar="isoval")
@@ -728,6 +766,8 @@ def main():
 					print("   PDB residue mode: measuring residue {}{}".format(options.residue, ", excluding " + " and ".join(filters) if filters else ""))
 					if options.self_only:
 						print("   Only the selected residue is kept (--self-only)")
+				if options.frames:
+					print("   Frames {} of each multi-structure file will be run".format(options.frames))
 				if options.cutoff or options.residue:
 					cutoff = options.cutoff if options.cutoff else "auto"
 					print("   Atoms farther than {} from atom1 are ignored (--cutoff)".format("the auto cutoff" if str(cutoff).lower() == "auto" else "{} Angstrom".format(cutoff)))
@@ -750,18 +790,8 @@ def main():
 			vec_df[numeric_cols] = vec_df[numeric_cols].round(2)
 			vec_df.to_csv(os.path.splitext(file)[0] + "_2d_output.csv", index=False)
 		else:
-			# Multi-structure files (multi-xyz, multi-sdf, multi-MODEL pdb): one run per structure
-			_, ext = os.path.splitext(file)
-			counters = {".xyz": parse_data.get_xyz_structures, ".sdf": parse_data.get_sdf_structures, ".mol": parse_data.get_sdf_structures,
-						".pdb": parse_data.get_pdb_models, ".ent": parse_data.get_pdb_models}
-			n_structures = len(counters[ext](file)) if ext in counters else 1
-			for idx in (range(n_structures) if n_structures > 1 else [None]):
-				options.structure = idx
-				if options.residue and str(options.residue).lower() == "all":
-					runs.extend(all_residues(file, options=options))
-				else:
-					runs.append(dbstep(file, options=options))
-			options.structure = None
+			# Multi-structure files (multi-xyz, multi-sdf, multi-MODEL pdb): one run per selected frame
+			runs.extend(run_file(file, options))
 
 	if dbstep._column_width and not options.quiet:
 		print("   " + "-" * (dbstep._column_width - 3))

dbstep/trajectory.py

@@ -0,0 +1,80 @@
+# -*- coding: UTF-8 -*-
+import os
+import sys
+
+from dbstep import parse_data
+
+"""
+trajectory
+
+Frames of multi-structure inputs: multi-frame XYZ, multi-record SDF and multi-MODEL PDB files.
+Each frame is measured by its own dbstep run (with the crop recomputed, since neighbours move),
+and --frames selects which frames to run.
+"""
+
+# functions returning the structure boundaries of a multi-structure file, by extension
+FRAME_COUNTERS = {
+	".xyz": parse_data.get_xyz_structures,
+	".sdf": parse_data.get_sdf_structures,
+	".mol": parse_data.get_sdf_structures,
+	".pdb": parse_data.get_pdb_models,
+	".ent": parse_data.get_pdb_models,
+}
+
+
+def count_frames(file):
+	"""Number of structures (frames) in a file; 1 for single-structure formats."""
+	if not isinstance(file, str):
+		return 1
+	_, ext = os.path.splitext(file)
+	if ext not in FRAME_COUNTERS:
+		return 1
+	return max(1, len(FRAME_COUNTERS[ext](file)))
+
+
+def parse_frames(spec, n_frames):
+	"""Frame indices selected by a --frames value.
+
+	The value uses Python slice semantics on the 0-based frame index: 'start:stop:stride' with any
+	part omitted ('::10' every tenth frame, '100:' from frame 100 on, ':50' the first fifty), a
+	single index, negative indices from the end, or a comma-separated list of these. A list or
+	tuple of integers is accepted from Python. Unset means every frame.
+
+	Returns:
+		list of 0-based frame indices in the order given
+	"""
+	if spec in (False, None, ""):
+		return list(range(n_frames))
+	tokens = [str(t).strip() for t in spec] if isinstance(spec, (list, tuple)) else [t.strip() for t in str(spec).split(",")]
+	indices = []
+	for token in tokens:
+		if not token:
+			continue
+		try:
+			if ":" in token:
+				parts = token.split(":")
+				if len(parts) > 3:
+					raise ValueError
+				indices.extend(range(n_frames)[slice(*[int(p) if p.strip() else None for p in parts])])
+			else:
+				indices.append(int(token))
+		except ValueError:
+			sys.exit(f"   Can't read frames '{spec}'. Use start:stop:stride (0-based, Python slice rules), e.g. --frames 0:1000:10")
+	selected = []
+	for index in indices:
+		frame = index + n_frames if index < 0 else index
+		if not 0 <= frame < n_frames:
+			sys.exit(f"   Frame {index} is out of range (file has {n_frames} frame{'s' if n_frames != 1 else ''})")
+		selected.append(frame)
+	if not selected:
+		sys.exit(f"   --frames {spec} selects no frames")
+	return selected
+
+
+def frame_indices(file, options):
+	"""Frames to run for a file: a list of 0-based indices, or [None] for a single-structure file without --frames."""
+	n_frames = count_frames(file)
+	frames = getattr(options, "frames", False)
+	if n_frames == 1 and not frames:
+		return [None]
+	return parse_frames(frames, n_frames)

dbstep/writer.py

@@ -277,7 +277,7 @@ def xyz_export(file, mol):
 	log.close()
 
 
-CSV_COLUMNS = ["file", "structure", "residue", "atom1", "atom2", "radius", "mol_vol", "percent_vbur", "percent_sbur", "bmin", "bmax", "L"]
+CSV_COLUMNS = ["file", "frame", "structure", "residue", "atom1", "atom2", "radius", "mol_vol", "percent_vbur", "percent_sbur", "bmin", "bmax", "L"]
 
 
 def csv_export(path, rows):

docs/plans/2.0-proteins-and-trajectories.md

@@ -201,8 +201,8 @@ Each step is one PR, merged in order. Estimates are working days for one person.
 | 1 | Radial crop (merged, #44) | `selection.crop`, `--cutoff`, `Mol_Vol(cut)` header, tests on existing XYZ data (crop vs no crop identical for small molecules; crop bounds the grid for a synthetic 5,000-atom cluster) | — | 1 |
 | 2 | PDB parser (merged, #45) | `PDBParser`, `get_pdb_models`, metadata arrays, `DataParser.keep`, parser unit tests, fix the README claim about PDB | — | 1.5 |
 | 3 | Residue selection (merged, #46) | `selection.py` name resolution and masks, `--residue/--atom/--atom2/--nowater/--nohet/--exclude-self/--self-only/--chain`, auto cutoff for PDB, residue label in output, exact-equivalence and semantics tests, `1crn.pdb` + `ala5.pdb` fixtures | 1, 2 | 3 |
-| 4 | `--residue all` and CSV output (PR opened) | residue loop, `writer.csv_export`, `--csv` | 3 | 1 |
-| 5 | Trajectories, native | multi-MODEL PDB and multi-frame XYZ frame loop, `--frames`, per-frame crop, `ala5_traj.pdb` tests | 3, 4 | 1.5 |
+| 4 | `--residue all` and CSV output (merged, #47) | residue loop, `writer.csv_export`, `--csv` | 3 | 1 |
+| 5 | Trajectories, native (PR opened) | multi-MODEL PDB and multi-frame XYZ frame loop, `--frames`, per-frame crop, `ala5_traj.pdb` tests | 3, 4 | 1.5 |
 | 6 | Trajectories, MDAnalysis | `traj` extra, `MDAnalysisFrames`, `--topology`, extra CI job | 5 | 1.5 |
 | 7 | 2.0 breaking changes | Python 3.10+, argparse migration, `from_rdkit()`, list-valued `exclude`, deprecation notes | can start any time, merge before release | 2 |
 | 8 | Docs and release | README section "Proteins and trajectories" with worked examples on `1crn.pdb`, CHANGELOG, version 2.0.0, release | all | 1 |

tests/pdb_files/README.md

@@ -6,3 +6,6 @@
 - `ala5.pdb` — penta-alanine with all hydrogens (53 atoms), two waters and a sodium ion, generated
   by `make_ala5.py` with RDKit (embedding seed 7, MMFF-optimised). Small enough that a grid over the
   whole molecule is feasible, which the crop-invariance tests rely on.
+- `ala5_traj.pdb` — 10 MODELs derived from `ala5.pdb` by `make_ala5_traj.py`: frame k is translated by
+  0.7·k Å along x, and the first water moves from 4.8 to 3.0 Å from CA of residue 3, so %V_bur around
+  A:3 rises monotonically with the frame index unless `--nowater` is used.