[jamesmkrieger]

ProDy cannot currently round-trip a CHARMM PSF: parsePSF followed by writePSF produces a file with a different topology from the one it read. Two independent causes, plus a third that hides one of them.

1. A CMAP cross-term is eight atom indices, not four

A !NCRTERM record holds the two coupled dihedrals of a CMAP term — eight indices, one term per line. From msbA's PSF:

        20        22        24        31        22        24        31        33

that is φ (20-22-24-31) and ψ (22-24-31-33) of one residue. parsePSF read count=n_crossterms*4, which consumed only the first 571 of the file's 1142 records and split each one into two unrelated 4-tuples — filling the array to the expected count with terms that do not exist. writePSF then wrote those 4-tuples two to a line, so the output disagreed with its own !NCRTERM count.

Fixed in the parser, the writer, setCrossterms (which now requires the eight columns and documents what they are), and evalCrossterms, which now takes the term width from the array instead of assuming four, so its neighbour map covers all eight atoms.

2. Sorting the indices inside a term changes which term it is

The atom order within an angle, dihedral, improper or cross-term is what identifies it: an angle's middle index is its vertex, a torsion's four are a sequence, a cross-term's eight are two dihedrals. All four setters began with .sort(1). msbA's PSF gives its first angles as

         2         1         3         2         1         4         2         1         5

— all three about vertex 1 — and they were stored as 1 2 3, 1 2 4, 1 2 5, so Angle.getVectors, documented as returning the bond vectors that originate from the central atom, took atom 2 as the vertex of all three.

Donors and acceptors are ordered pairs too (!NDON is a donor and its hydrogen, !NACC an acceptor and its antecedent), so they are no longer sorted or deduplicated either. The term-level argsort chains go with the sorts: nothing downstream needs a canonical term order, and keeping the order read is what makes the round-trip exact.

setBonds and setNBExclusions keep their within-pair sort deliberately — those pairs really are symmetric, so sorting each is what collapses (i, j) and (j, i) in the np.unique that follows, and for bonds it is also what guarantees i < j in the '%d %d' keys of _bondIndex that Bond looks up. Their row-level argsort chains are dropped as dead code, since np.unique(..., axis=0) already sorts rows.

3. !NCRTERM was dropped when it directly follows the exclusions

The loop reading !NNB stops on the first line containing ! — which has already consumed the next section's header. CHARMM puts !NGRP and !NUMLP there, so real files survived, but a file that goes straight to !NCRTERM — including one written by writePSF — had the scan start past the header, and every cross-term was silently dropped with no error. The line the loop stopped on is now checked before scanning on.

Testing

• New prody/tests/trajectory/test_psffile.py and a hand-made ten-atom fixture topology.psf with every section populated, laid out the way CHARMM writes one, with no term's indices already sorted. 5 of the 7 tests fail without these fixes.
• prody/tests/atomic, trajectory and proteins: 701 passed. The 7 failures are pre-existing and environmental (no mmtf, no Hpb, network fetches) and are unaffected.
• Real-world check on a 720,343-atom CHARMM-GUI system: an OpenMM System built from the original PSF and from the ProDy round-tripped PSF give identical potential energy, −9,820,146.3 kJ/mol, with the same 704,223 constraints.

Not addressed here: writePSF still omits !NGRP, !NUMLP and the per-atom !NNB pointer array. OpenMM tolerates that; CHARMM and NAMD may not.

🤖 Generated with Claude Code

[jamesmkrieger]

the sporadic insty parallel bug is fixed in PR #2277

[karolamik13]

Looks good. I tested on a few PSF files.

Code:
from prody import *
import numpy as np

name = 'test.psf'
p = parsePSF(name)
writePSF(name[:-4]+'_out.psf', p)
p2 = parsePSF(name[:-4]+'_out.psf')

for x in ['bonds', 'angles', 'dihedrals', 'impropers', 'donors', 'acceptors', 'crossterms']:
print(x, np.array_equal(getattr(p, '' + x), getattr(p2, '' + x)))

print(p._crossterms.shape)
print(p2._crossterms.shape)

Output:
bonds True
angles True
dihedrals True
impropers True
donors True
acceptors True
crossterms True
(1012, 8)
(1012, 8)