fixed the docs

Author: stefdoerr

doc/source/explanation/atom-selection-language.md

@@ -52,7 +52,6 @@ lookups and element checks:
 | `hydrogen` | All atoms with `element == "H"` |
 | `noh` | All non-hydrogen atoms |
 | `all` | Every atom in the molecule |
-| `none` | No atoms |
 
 ## Per-atom field comparisons
 
@@ -100,10 +99,13 @@ mol.atomselect("beta < 20")
 mol.atomselect("resid 40 to 60")
 mol.atomselect("index 0 to 99")
 
-# Negation with !=
-mol.atomselect("chain != B")
+# Negation: use `not`
+mol.atomselect("not chain B")
 ```
 
+The `!=` operator only applies inside the modulo form (e.g.
+`resid % 2 != 0`); for plain field comparisons use `not`.
+
 ## Boolean composition
 
 Combine selections with `and`, `or`, `not`, and parentheses:

doc/source/explanation/segments-chains-and-bonds.md

@@ -77,9 +77,15 @@ from moleculekit.tools.autosegment import autoSegment
 
 mol_seg = autoSegment(mol)   # returns a modified copy
 import numpy as np
-print(np.unique(mol_seg.segid))   # e.g. ['P0', 'P1', 'W0']
+print(np.unique(mol_seg.segid))   # e.g. ['P0', 'P1', 'P2']
 ```
 
+Each contiguous polypeptide run (and each contiguous run of water residues)
+gets the next `P{i}` segid. Water residues additionally receive `chain = "W"`
+to keep them visually distinct from polymer chains, but their **segid** stays
+on the same `P{i}` numbering sequence as everything else; there is no `W0`
+segid.
+
 The function accepts a `basename` argument to control naming: `basename="P"`
 produces `P0`, `P1`, `P2`, etc. The `fields` argument controls which field(s)
 are written: `("segid",)` (default), `("chain",)`, or `("segid", "chain")`.

doc/source/explanation/system-preparation-pipeline.md

@@ -193,7 +193,7 @@ up the ligand spec automatically.
 ### mutateResidue before systemPrepare
 
 When you want to swap canonical residues (e.g. a point mutation), call
-{py:meth}`~moleculekit.molecule.Molecule.mutateResidue` with `sel` and `new_resname` first. The mutator uses Dunbrack rotamers
+{py:meth}`~moleculekit.molecule.Molecule.mutateResidue` with `sel` and `newres` first. The mutator uses Dunbrack rotamers
 to place the new sidechain, and {py:func}`~moleculekit.tools.preparation.systemPrepare` then protonates and optimizes
 the mutated residue along with the rest of the protein.
 

doc/source/explanation/trajectories-and-frames.md

@@ -51,33 +51,41 @@ The `read` method accepts `frames` and `skip` arguments to load a subset
 without reading the whole file into memory:
 
 ```python
-# Load every 10th frame
+# Load every 10th frame (uniform stride)
 mol.read("run1.xtc", skip=10)
 
-# Load a specific list of frame indices
-mol.read("run1.xtc", frames=range(0, 1000, 10))
-
-# Load only frame 500
+# Load only frame 500 of a single trajectory
 mol.read("run1.xtc", frames=[500])
+
+# Load a specific list of frame indices from a single trajectory
+mol.read("run1.xtc", frames=[[0, 10, 20, 30, 40]])
 ```
 
-`skip=N` applies uniform subsampling; `frames=` takes an arbitrary list of
-indices.
+`skip=N` applies uniform subsampling. `frames=` is **per file**: it must be a
+list with one entry per trajectory in `filename`, and each entry is either an
+int (single frame) or a list of ints (multiple frames). When reading just one
+trajectory this means wrapping your indices in an outer list, as in the third
+example above. For most subsampling use cases `skip=` is simpler.
 
 ### Multi-trajectory loads
 
-Pass a list of trajectory files to the `Molecule` constructor, or call
-`read(..., append=True)` repeatedly. Frames are concatenated in order:
+Pass topology and trajectories as a single list to the `Molecule` constructor,
+or call `read(..., append=True)` repeatedly. Frames are concatenated in order:
 
 ```python
-mol = Molecule("system.prmtop", ["run1.xtc", "run2.xtc", "run3.xtc"])
+mol = Molecule(["system.prmtop", "run1.xtc", "run2.xtc", "run3.xtc"])
 # Equivalent to:
 mol = Molecule("system.prmtop")
 mol.read("run1.xtc")
 mol.read("run2.xtc", append=True)
 mol.read("run3.xtc", append=True)
 ```
 
+The second positional argument of {py:class}`~moleculekit.molecule.Molecule`
+is `name=`, **not** a trajectory file — passing trajectories there assigns
+them to `mol.name` and loads nothing. Always use a single list (as above) or
+the explicit `read(..., append=True)` form.
+
 `mol.numFrames` is the total number of frames from all trajectories. Frame
 ordering mirrors the order the files were passed.
 
@@ -123,9 +131,9 @@ precision, what they store, and typical file size:
 
 | Format | Extension | Precision | Box | Velocities | Notes |
 |---|---|---|---|---|---|
-| XTC | `.xtc` | Reduced (lossy ~0.001 Å) | Yes | No | GROMACS native; very compact |
+| XTC | `.xtc` | Reduced (lossy XDR fixed-point, ~0.01 Å) | Yes | No | GROMACS native; very compact |
 | DCD | `.dcd` | Full `float32` | Varies | No | NAMD/CHARMM native |
-| TRR | `.trr` | Full `float64` | Yes | Yes | GROMACS native; larger files |
+| TRR | `.trr` | Full (file stores up to `float64`; moleculekit holds `coords` as `float32`) | Yes | Yes | GROMACS native; larger files |
 | NetCDF | `.nc`, `.ncdf` | Full `float32` | Yes | Optional | AMBER native |
 
 XTC is the most common format encountered in practice. Its compressed storage

doc/source/howto/compute-dihedrals.md

@@ -34,7 +34,7 @@ print(f"psi(17) after : {np.degrees(mol.getDihedral(psi17)):.1f}°")
 | Method / argument | What it does |
 |---|---|
 | `mol.getDihedral(atom_quad)` | Returns the dihedral angle in **radians** for the four atom indices in `atom_quad`. Operates on the current `mol.frame`. |
-| `mol.setDihedral(atom_quad, radians, bonds=None, guessBonds=False)` | Rotates the downstream half of the molecule around the central bond of the dihedral so the angle becomes `radians`. For a chain of modifications pass `bonds=mol._getBonds()` once to keep the bond table from being re-guessed on every call. |
+| `mol.setDihedral(atom_quad, radians, bonds=None, guessBonds=False)` | Rotates the downstream half of the molecule around the central bond of the dihedral so the angle becomes `radians`. For a chain of modifications pass `bonds=mol.bonds` (or a pre-built bond array) once to keep the bond table from being re-built on every call. |
 
 ## Common variations
 
@@ -59,7 +59,7 @@ angles = np.array([
 
 - Both `getDihedral` and `dihedralAngle` return angles in **radians**; convert with `np.degrees(angle)` if you want degrees.
 - `mol.getDihedral` operates on the current `mol.frame` — set `mol.frame = i` first if you want a specific frame.
-- `mol.setDihedral` rotates the downstream side of the dihedral in place; the upstream side is held fixed. If the topology is ambiguous (the rotation would split the molecule wrong), the call may fail — guard against this by passing an explicit `bonds=` array.
+- `mol.setDihedral` rotates the downstream side of the dihedral in place; the upstream side is held fixed. If the topology is ambiguous (the rotation would split the molecule wrong), the call may fail — guard against this by passing an explicit `bonds=mol.bonds` array.
 - For computing many dihedrals at once across a trajectory, prefer {py:class}`~moleculekit.projections.metricdihedral.MetricDihedral` from `moleculekit.projections.metricdihedral` — it batches the work efficiently.
 
 ## See also

doc/source/howto/compute-protein-ligand-interactions.md

@@ -8,57 +8,48 @@ Detect hydrogen bonds, π–π stacking, cation–π, and σ-hole interactions b
 
 ```python
 import numpy as np
-from rdkit.Chem import ChemicalFeatures
-from rdkit import RDConfig
-import os
 from moleculekit.molecule import Molecule
-from moleculekit.interactions.interactions import hbonds_calculate
+from moleculekit.smallmol.smallmol import SmallMol
+from moleculekit.interactions.interactions import (
+    hbonds_calculate,
+    get_ligand_donors_acceptors,
+)
 
 mol = Molecule("3PTB")
 mol.templateResidueFromSmiles(
     mol.resname == "BEN",
-    smiles="NC(=N)c1ccccc1",
+    smiles="[NH2+]=C(N)c1ccccc1",
     addHs=True,
 )
 
-# Convert the ligand to RDKit, then ask RDKit which atoms are H-bond
-# donors and acceptors.
-lig = mol.copy(sel="resname BEN")
-rdlig = lig.toRDKitMol(sanitize=True)
-
-fdef = os.path.join(RDConfig.RDDataDir, "BaseFeatures.fdef")
-factory = ChemicalFeatures.BuildFeatureFactory(fdef)
-features = factory.GetFeaturesForMol(rdlig)
-
+# Build a SmallMol view of the ligand, then ask the library helper for
+# donor/acceptor indices in the parent Molecule's atom numbering.
+lig = SmallMol(mol.copy(sel="resname BEN"))
 start_idx = int(mol.atomselect("resname BEN", indexes=True)[0])
-acceptors = [idx + start_idx for f in features if f.GetFamily() == "Acceptor"
-             for idx in f.GetAtomIds()]
-donor_pairs = [(idx + start_idx, idx + start_idx)
-               for f in features if f.GetFamily() == "Donor"
-               for idx in f.GetAtomIds()]
+donors, acceptors = get_ligand_donors_acceptors(lig, start_idx=start_idx)
 
 # For non-periodic structures hbonds_calculate needs a zero box.
 mol.box = np.zeros((3, mol.numFrames), dtype=np.float32)
-hbonds = hbonds_calculate(mol, np.array(donor_pairs), np.array(acceptors), sel2="protein")
+hbonds = hbonds_calculate(mol, donors, acceptors, sel2="protein")
 print(f"H-bonds in frame 0: {len(hbonds[0])}")
 ```
 
 ## Parameters that matter
 
 | Function/parameter | What it does |
 |---|---|
-| `toRDKitMol(sanitize=True)` | Convert the ligand selection to an RDKit Mol for feature detection. |
-| `RDConfig.RDDataDir / "BaseFeatures.fdef"` | RDKit's standard donor/acceptor feature definition file. |
-| `start_idx` | Offset of the first ligand atom in the parent `Molecule`; required when the ligand does not start at atom 0. |
-| `hbonds_calculate(mol, donors, acceptors, sel2=...)` | Returns a list (one entry per frame) of H-bond arrays. |
+| {py:class}`~moleculekit.smallmol.smallmol.SmallMol` | Lightweight RDKit-backed view of the ligand used by the helper to walk the bond graph. |
+| `start_idx` | Offset of the first ligand atom in the parent `Molecule`; the helper adds it so the returned indices reference the parent. |
+| `get_ligand_donors_acceptors(smol, start_idx=...)` | Returns `(donors, acceptors)`. `donors` is a `(N, 2)` array of `[heavy_atom_idx, hydrogen_idx]` pairs; `acceptors` is a 1D array of heavy-atom indices. |
+| `hbonds_calculate(mol, donors, acceptors, sel2=...)` | Returns a list (one entry per frame) of H-bonds; each entry is a list of `(donor_heavy, donor_h, acceptor)` triples. |
 
 ## Common variations
 
 ```python
-# Visualise H-bonds in VMD (requires VMD on PATH)
-from moleculekit.interactions.interactions import view_hbonds
-
-view_hbonds(mol, hbonds)
+# Convert the per-frame list of triples into a single numpy array of donors,
+# Hs, and acceptors for frame 0.
+frame0 = np.array(hbonds[0])   # shape (n_hbonds, 3)
+print("donor heavy atoms, Hs, acceptors:\n", frame0)
 ```
 
 ```python
@@ -75,9 +66,9 @@ rings = get_protein_rings(mol)
 
 ## Gotchas
 
-- Reasonable hydrogens must be present — run {py:func}`~moleculekit.tools.preparation.systemPrepare` before this step if the input lacks explicit H atoms.
-- Feature detection uses RDKit through {py:meth}`~moleculekit.molecule.Molecule.toRDKitMol`; make sure RDKit is installed.
-- `start_idx` is critical when the ligand atoms do not start at index 0 in the parent molecule.
+- Reasonable hydrogens must be present — run {py:func}`~moleculekit.tools.preparation.systemPrepare` (or pass `addHs=True` to {py:meth}`~moleculekit.molecule.Molecule.templateResidueFromSmiles`) before this step.
+- {py:func}`~moleculekit.interactions.interactions.get_ligand_donors_acceptors` returns donor *pairs* `[heavy, H]`, **not** lone heavy atoms. Building your own pairs as `(heavy, heavy)` silently produces wrong H-bond results because `hbonds_calculate` reads column 1 as the hydrogen index.
+- `start_idx` is critical when the ligand atoms do not start at index 0 in the parent molecule — the helper offsets every returned index by `start_idx`.
 - For non-periodic structures set `mol.box` to a zero array before calling `hbonds_calculate`.
 
 ## See also

doc/source/howto/compute-rmsd-rmsf.md

@@ -33,11 +33,14 @@ print(rmsd)
 ## Common variations
 
 ```python
+import numpy as np
 # RMSF per Cα atom over a trajectory
 from moleculekit.projections.metricfluctuation import MetricFluctuation
 
-rmsf = MetricFluctuation("protein and name CA").project(mol)
-# rmsf shape: (n_frames, n_CA); for a single per-atom value over the trajectory use rmsf.mean(axis=0)
+# .project returns squared per-atom displacement (shape (n_frames, n_CA),
+# units Ų). RMSF = sqrt(mean over frames of the squared displacement).
+sq_disp = MetricFluctuation("protein and name CA").project(mol)
+rmsf = np.sqrt(sq_disp.mean(axis=0))   # shape (n_CA,), units Å
 ```
 
 ```python
@@ -51,7 +54,7 @@ rmsd_aligned = molRMSD(mol, ref, ca_idx, ca_idx)
 - {py:func}`~moleculekit.util.molRMSD` computes raw (non-mass-weighted) RMSD; align the structures first if you want to remove rigid-body motion.
 - RMSF is per-atom and depends on the choice of reference (mean structure); the result is sensitive to conformational heterogeneity.
 - Both `mol` and `refmol` must have the same number of selected atoms.
-- {py:class}`~moleculekit.projections.metricfluctuation.MetricFluctuation` returns an `(n_frames, n_atoms)` array — use `rmsf.mean(axis=0)` to get a per-atom fluctuation averaged over all frames.
+- {py:class}`~moleculekit.projections.metricfluctuation.MetricFluctuation` returns **squared** per-atom displacement (`(n_frames, n_atoms)`, units Ų). For RMSF you must take `np.sqrt(sq_disp.mean(axis=0))` — `sq_disp.mean(axis=0)` alone is mean-squared-fluctuation, not RMSF.
 
 ## See also
 

doc/source/howto/fetch-from-rcsb-and-opm.md

@@ -47,14 +47,24 @@ mol, thickness = get_opm_pdb("1BL8")
 from moleculekit.opm import align_to_opm
 
 mol = Molecule("my_structure.pdb")
-mol_opm = align_to_opm(mol, maxalignments=3)
+results = align_to_opm(mol, maxalignments=3)
+
+# results is a list, one entry per OPM hit. Each entry has the OPM PDB ID,
+# the membrane thickness, and a list of high-scoring sequence pairs (HSPs)
+# whose aligned_mol is `mol` re-imaged into the OPM frame.
+for hit in results:
+    print(hit["pdbid"], "thickness:", hit["thickness"])
+    for hsp in hit["hsps"]:
+        aligned_mol = hsp["aligned_mol"]   # Molecule, oriented in the OPM membrane frame
+        print(f"  TM={hsp['TM-Score']:.2f}  RMSD={hsp['Common RMSD']:.2f} Å")
 ```
 
 ## Gotchas
 
 - RCSB downloads respect the server rate limits; avoid hammering the API in tight loops.
 - Set the `LOCAL_PDB_REPO` environment variable to a local PDB mirror directory to avoid repeated network downloads.
-- OPM membership requires a known PDB ID or a successful BLAST sequence alignment; structures absent from OPM will raise or return `None`.
+- OPM membership requires a known PDB ID or a successful BLAST sequence alignment; when nothing matches, {py:func}`~moleculekit.opm.align_to_opm` returns an **empty list** (not `None`).
+- {py:func}`~moleculekit.opm.align_to_opm` returns a `list[dict]` — the aligned `Molecule` objects live under `hit["hsps"][j]["aligned_mol"]`, not at the top level.
 - {py:func}`~moleculekit.opm.get_opm_pdb` with `keep=False` (default) strips the dummy membrane atoms that OPM adds; pass `keep=True` if you need them for visualization.
 
 ## See also

doc/source/howto/read-a-trajectory.md

@@ -26,9 +26,15 @@ print(mol.numFrames)
 ## Common variations
 
 ```python
-# Load topology first, then read a specific frame range
+# Read every 10th frame (uniform stride)
 mol = Molecule("topology.psf")
-mol.read("trajectory.xtc", frames=list(range(0, 500)))
+mol.read("trajectory.xtc", skip=10)
+```
+
+```python
+# Read a single specific frame
+mol = Molecule("topology.psf")
+mol.read("trajectory.xtc", frames=[500])
 ```
 
 ```python
@@ -41,9 +47,10 @@ print(mol.numFrames)
 
 ## Gotchas
 
-- All frames are loaded into memory at once; for very long trajectories use `frames=` to read a slice.
+- All frames are loaded into memory at once; for very long trajectories use `skip` or `frames=` to read a subset.
 - The atom count of the trajectory must exactly match the topology — a mismatch raises an error.
-- `skip` and `frames` can be combined: `frames` selects the subset to consider, then `skip` further subsamples within that subset.
+- `frames=` is **per file**: when reading a single trajectory it must be a length-1 list whose single element is either an int or a list of ints (e.g. `frames=[[0, 10, 20]]`). For uniform subsampling prefer `skip=`.
+- `skip` strides the final merged coordinate array; it applies independently of `frames` and is the simplest way to subsample a long trajectory.
 - When `append=False` (default) a second `mol.read(traj)` call replaces the existing coordinates rather than extending them.
 
 ## See also

doc/source/howto/view-a-molecule.md

@@ -46,7 +46,7 @@ mol.view(viewer="pymol")
 When `viewer=` is not given, moleculekit picks a backend in this order:
 
 1. `MOLECULEKIT_VIEWER` environment variable (e.g. `export MOLECULEKIT_VIEWER=molstar`).
-2. `moleculekit.config["viewer"]` runtime setting.
+2. The runtime setting from {py:func}`moleculekit.config.config` (e.g. `from moleculekit.config import config; config(viewer="molstar")`).
 3. VMD found on `$PATH`.
 4. PyMOL found on `$PATH`.
 5. Molstar fallback (always available).