Fix residues silently dropped or mismatched by residue identity handling

[NDoering99]

Three bugs in how residues are identified. The first two silently drop or swap residues; the third turns a common user error into an uninterpretable traceback.

1. Residues were located by index arithmetic

DihedralAngles looked residues up positionally:

res = self.traj.residues[res_id - self.first_res_num]

That is only correct when residue numbers form a contiguous, unique run starting at first_res_num. Two common situations break it:

• A numbering gap (an unresolved loop, a spliced construct). Every residue after the gap reads its dihedral from the wrong residue.
• A repeated residue number. Same, plus the tail of the chain is never reached.

The failure is silent. The progress bar fills, every output file is written, and the numbers look plausible.

How it showed up

A CHARMM-GUI system with insertion codes 205A–205O. The PSF format has no insertion-code field, so all 16 residues arrive as resid 205:

311 residues in the topology, 296 distinct residue numbers

MDPath processed 296 of 311 residues, mapped everything past 205 onto a residue 15 positions earlier, and never touched the last 15 residues. Nothing in the output indicated a problem.

Fix

• Residues are looked up through a resid -> position map built from the topology.
• The parallel loop iterates the resids that actually exist rather than range(first_res_num, last_res_num + 1), so a gap no longer queues residues that aren't there.
• Duplicate residue numbers raise a ValueError naming the offending resids and the usual cause, instead of producing corrupt results.

Only amino acids have to be unique. Ligands, ions and waters routinely reuse the protein's numbering and are not graph nodes, so they are excluded from the uniqueness check — mdpath/tests/test_topology.pdb itself has an UNK in chain B numbered 1, alongside ILE 1 in chain A.

2. Force field residue names were rejected as non-protein

All four residue filters used Bio.PDB.Polypeptide.is_aa(), which only knows PDB chemical-component names. It rejects the protonation-state names MD force fields write:

resname	is_aa()	actually
HSD	False	CHARMM δ-protonated His
HSP	False	CHARMM doubly-protonated His
HID, HIE	False	AMBER His
CYX	False	disulfide-bonded Cys
ASH	False	protonated Asp
HSE	True	matched as homoserine, not ε-protonated His

Rejected residues were dropped from the graph skeleton, from residue_CA_coordinates, and from res_num_from_pdb. In the system above this removed 4 histidines from the network — they could not appear in any pathway, with no warning.

AAMAPPING in visualization.py had the same gap, so those residues were also skipped when assigning generic numbers.

Fix

Adds is_amino_acid() in mdpath/src/structure.py, which normalises force field aliases through FORCEFIELD_RESNAME_ALIASES before the is_aa() check, and uses it at all four call sites. AAMAPPING lookups are normalised the same way.

3. A topology with no protein failed with OverflowError

res_num_from_pdb() seeded first_res_num with float("inf") and returned int(first_res_num) unconditionally, so a topology in which nothing was recognised as an amino acid failed with:

OverflowError: cannot convert float infinity to integer

which says nothing about the cause. This is easy to reach by accident: passing a coordinate-only file to -top. MDAnalysis has no residue names to read from a DCD or XTC, so it defaults every residue to UNK, mdpath.py writes that into first_frame.pdb, and nothing downstream is recognised as protein.

Now raises a ValueError naming the file and the likely cause.

Tests

Seven new tests in mdpath/tests/test_structure.py covering force field resname acceptance, non-protein rejection, the duplicate-resid guard, gap-tolerant lookup, non-protein residues reusing protein numbering, and the empty-topology error.

The mock_traj fixture used a dict for traj.residues, which never modelled MDAnalysis' positionally-indexed ResidueGroup — with the old code traj.residues[1 - 1] raised KeyError and was swallowed by the except clause, so the lookup was never actually exercised. Replaced with something indexable.

Full suite: 107 passed, 1 failed. The failure is pre-existing on master — a last-ULP float comparison in test_mutual_information.py::test_hist_range_none_is_exact_no_op (0.6667438140756162 != 0.6667438140756161), unrelated to this change.

Compatibility

No API changes. DihedralAngles.__init__ keeps its signature; num_residues is still accepted and stored. Topologies that were already contiguous and unique behave exactly as before — the only behaviour change for them is that force field residue names are no longer discarded.

🤖 Generated with Claude Code