MOLINFO WHOLE: check that the reference is not split by the periodic boundary (#1492)

[aalhossary]

Description

Closes #1492.

MOLINFO ... WHOLE is taken on trust. If the reference pdb is in fact split by the periodic boundary, WHOLEMOLECULES splits whole molecules in every frame, silently (details and a reproducer in #1492). This PR turns that silent wrong result into an error at setup:

the reference structure in the MOLINFO pdb file is not whole: 12 of its 862 atoms here are closer to a periodic
image of the others, in the box of its CRYST1 record, than to the others themselves. Atoms: 163 (LYS10) 164 (LYS10) ...

• Tree::buildTree() (used by WHOLEMOLECULES with MOLINFO ... WHOLE) now calls checkReferenceIsWhole() for its group of atoms. When the reference has a usable box, the check builds the tree a second time on minimum-image distances in that box and walks it over the reference: a whole reference comes back unchanged, wrapped atoms come back moved by a box vector, and those are the atoms named in the error. The tree itself is unchanged.

• A usable box (GenericMolInfo::getReferenceBox()) is a CRYST1 record with space group P 1, which is what simulation programs write. No CRYST1, a degenerate cell, the 1 Å placeholder cube, or a crystallographic cell (any other space group) skip the check, so those references behave exactly as before. PDB now keeps the space group for this.

• Speed: a linear-time grid search (hasCloseImage()) first checks whether any atom comes closer to a periodic image of another atom than the longest link of the tree. If none does, the second tree is provably the same as the first, and it is skipped. This is the usual case for a whole reference. Extra setup time (optimized build, best of 3, rhombic dodecahedron / triclinic / cubic cells):

  atoms	whole driver run without the check	extra time for the check	extra if the second tree were always built
  1,000	0.02 s	< 0.001 s	0.01–0.02 s
  25,000	1.4 s	< 0.01 s	3–7 s
  50,000	5.5 s	< 0.02 s	13–29 s

• The check runs per group, not on the whole MOLINFO file: on a solvated reference, the solvent always spans the box, so a file-wide check would flag it.

Target release

I would like my code to appear in release master (v2.11)

Type of contribution

• changes to code or doc authored by PLUMED developers, or additions of code in the core or within the default modules
• changes to a module not authored by you
• new module contribution or edit of a module authored by you

Copyright

• I agree to transfer the copyright of the code I have written to the PLUMED developers or to the author of the code I am modifying.

Tests

• I added a new regtest or modified an existing regtest to validate my changes.
• I verified that all regtests are passed successfully on GitHub Actions.

New regtests, on trp-cage (1L2Y, model 1):

• basic/rt-molinfo-whole-check reads a whole and a split copy with MOLINFO ... WHOLE and reports which are rejected and which atoms are named. Copies without a usable CRYST1 (none, placeholder, zero, crystallographic space group) are accepted. Without the check, this test fails.
• basic/rt-molinfo-whole-repair makes five split frames whole with a whole reference (rhombic dodecahedron, truncated octahedron, triclinic and cubic cells; atom- and residue-wise wrapping).

Beyond the regtests, the check and the repair were tested against ground truth outside PLUMED:

• Synthetic references (~117,000 cases): 19 structures (proteins with and without hydrogens, complexes, ligands, DNA, tRNA, collagen, membrane proteins) × selections (all atoms, heavy atoms, backbone, Cα, every 5th/10th Cα) × 5 cell shapes × margins from 20 Å down to 1 Å × placements (random, centroid at corners/faces/edges, rim of the protein 0.05–4 Å across a face, edge or corner of the brick or of the compact cell) × wrapping (atom-wise brick or compact, residue-wise). The truth is known by construction. The check rejected every split reference (including those with 1–3 atoms across) and named exactly the right atoms, accepted every whole one, and the repair with the whole reference was right in every case, whenever the group's own periodic images are farther than its longest tree link. That held for every all-atom, heavy-atom, backbone and Cα selection. The only exceptions were very sparse selections (every 5th/10th Cα) in tight boxes, where the geometry is ambiguous; there, WHOLEMOLECULES is wrong in more than half of the cases even with a whole reference, and the error message names this cause too.
• gmx trjconv -pbc whole (bond graph from the tpr) on 6 systems × 60 split frames: identical to PLUMED's repair in all 360.
• --debug-dd on 2 and 4 ranks and --multi 2: output byte-identical to serial; a split reference raises the error on every rank, without hanging.
• Real files: references written by GROMACS pass, and the split GB3 reference of MOLINFO WHOLE: a reference split by PBC makes WHOLEMOLECULES break whole molecules, silently #1492 is rejected, naming its 12 atoms. Crystal structures from the PDB are not checked (crystallographic space group); relabelled as P 1, 3 of 13 whole asymmetric units were rejected, each time because the group held more than one unbonded part (two chains, or protein and lipids) and one part was closer to a crystal neighbour; this is why other space groups are skipped.
• Locally (debug build), the full suite: 817 performed, 1 error, basic/rt-make-fortran08, which fails the same way without this change in my clang/libc++ build. All 37 regtests that use WHOLE or WHOLEMOLECULES pass.

Prepared with the help of Claude (Anthropic).

[GiovanniBussi]

Thanks @aalhossary !

I am not sure this is a good idea. I see that the test is only done when using WHOLE - which indeed is the only case where the coordinates in the PDB are used to reconstruct images. However, having a unique way of making the coordinates whole, or of checking if they are whole, looks difficult to me. In a way, this is the reason to ask user to provide coordinates.

[aalhossary]

Thanks @GiovanniBussi. I agree that making the reference whole is the user’s job, and the PR does not try to do it, nor to define one way of doing it. It only checks the reference against the box written in the same file, and only when there is one (CRYST1, space group P 1).
The question it asks is narrow: Are some atoms closer to a periodic image of the rest of the group than to the rest itself? If so, WHOLEMOLECULES with this reference will split the group in every frame, as in #1492. Then it stops and names those atoms.
It never changes coordinates, and with no CRYST1, a placeholder or crystallographic cell, or without WHOLE, nothing changes.

BTW, An ambiguous case does exist: a very sparse group (e.g. every 10th Cα) in a tight box can have its own periodic images closer than its longest tree link. There, the check may reject a whole reference, but WHOLEMOLECULES is unreliable there anyway (wrong in 882 of the 1568 such cases I tested), and the message says the atoms may be too sparse for the box. For all-atom, heavy-atom, backbone and Cα groups I found no false positives.

If you prefer, I can make it a warning in the log instead of an error, or put it behind a keyword. Which would you rather have?

[aalhossary]

The full test report for this PR, for reference: how the check and the repair were tested against ground truth.

Summary. 89,780 synthetic references with known truth (117,710 runs, one set repeated with the fast path). In the 88,016 cases where the geometry is unambiguous (see D ≤ L below), every split reference was rejected, every whole one accepted, the displaced atoms named exactly, and the frame repaired correctly. The repair agrees with gmx trjconv -pbc whole on 360/360 frames. Extra setup time is < 0.02 s at 50,000 atoms.

Synthetic references. Each case takes a whole structure, places it in a cell, wraps it, and runs PLUMED twice: (A) the result as the MOLINFO WHOLE reference, where the check must reject it exactly when atoms were wrapped and name those atoms; (B) the wrapped coordinates as a frame, repaired with the true whole structure as reference. 19 structures (1AKE 1B27 1BNA 1C3W 1CAG 1EHZ 1L2Y 1UBQ 2K39 2LZM 2RH1 3EML 3PTB 4AKE 4HHB 6LU7 7K3G, GB3 and ubiquitin from MD); cubic, triclinic, rhombic dodecahedron (square and hexagonal) and truncated octahedron cells; centroid at corners, faces and edges, rim of the molecule 0.05–4 Å across a face, edge or corner, or random; atom-wise (brick or compact cell) or residue-wise wrapping.

test set	selections	margins	cases	split	whole	result
all atoms, 19 structures	all	20, 6, 2.5 Å	27,930	26,619	1,311	27,930 correct
same, fast path, new seed	all	20, 6, 2.5 Å	27,930	26,541	1,389	27,930 correct
selections, 6 structures	heavy, backbone, Cα, every 5th / 10th Cα, guide atoms	20, 6, 2.5 Å	45,570	36,007	9,563	43,904 / 43,904 with D > L; 98 / 1,666 with D ≤ L
stress, 10 structures	all, Cα	1 Å	9,800	8,831	969	9,800 correct
residue-wise wrapping, 12 structures	all, Cα	20, 6, 2.5 Å	6,480	6,284	196	6,480 correct

2,186 of the all-atom cases had only 1–3 atoms across the boundary; all were caught. Without the check, WHOLEMOLECULES with the split reference broke a whole frame in 24,834 of the 26,619 split all-atom cases.

The ambiguous case, D ≤ L. L is the longest link of the group's nearest-neighbour tree; D is the shortest distance between an atom and a periodic image of another atom of the same whole group. When D > L, the tree in the box equals the plain tree, so the check is exact. When D ≤ L, an image is as close as a real neighbour, and "whole" is not defined by the coordinates alone. This happened only for every 5th / 10th Cα in tight boxes: 1,666 of the 45,570 selection cases, 1,568 of them answered wrongly. WHOLEMOLECULES without the check was itself wrong in 882 of those 1,568. The error message names this cause ("or when these atoms are too sparse for that box").

Outside tools and real files.

test	result
gmx trjconv -pbc whole (bond graph from the tpr), 6 systems × 60 split frames	360 / 360 identical
--debug-dd on 2 and 4 ranks, --multi 2 × 2 ranks (4HHB, GB3, ubiquitin)	byte-identical to serial
split reference in the same parallel runs	error on every rank, no hang
real files with their own CRYST1 (GROMACS references, PDB entries as deposited, crystal waters, cells relabelled P 1), 90 files	77 as expected, 13 explained below
nm and Å units, references missing hydrogens or residues, several ENTITY	correct

The 13 real-file cases where the verdict differs from "is the molecule whole":

reference	check	reading
split GB3 with no CRYST1, the 1 Å placeholder, zero lengths, all zero (4)	skipped	no usable box; behaves as before
whole GB3 in a cubic cell of 20, 10 or 2 Å (3)	rejected	the cell is smaller than the protein, so it cannot be the box it was written in
ubiquitin MD, protein and solvent as one ENTITY	rejected	solvent fills the box, so the group is never whole; this is why the check is per group
1AKE, 1B27, 1C3W crystal cell relabelled P 1 (3 of 13 relabelled)	rejected	the group holds unbonded parts (two chains, or protein and lipids), and one is closer to a crystal neighbour; with one ENTITY per chain, all pass
1EHZ relabelled P 1, tRNA only	rejected	64 atoms closer to a crystal neighbour
4AKE, deposited in P 1, with crystal waters in the group	rejected (8 atoms)	the waters are closer to the neighbouring cell; without them it passes

The relabelled cells are why only P 1 (or a blank space group) is checked. The split GB3 reference of #1492 is rejected with its 12 atoms named, also with its cell 2% larger or smaller (an NPT frame) and with a blank space group.

Cost (optimized serial build, best of 3):

atoms	whole driver run without the check	extra for the check	extra if the second tree were always built
1,000	0.02 s	< 0.001 s	0.01–0.02 s
25,000	1.4 s	< 0.01 s	3–7 s
50,000	5.5 s	< 0.02 s	13–29 s

Regtests. Full local suite (debug build, clang/libc++): 817 run; the one failure, basic/rt-make-fortran08, fails the same way without this change in that toolchain. All 37 regtests that use WHOLE or WHOLEMOLECULES pass.