[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 makes PLUMED check the reference at setup and write a warning to the log when it looks split:

WARNING for action MOLINFO with label @0 : the reference structure in the MOLINFO pdb file may not be 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) ...

It is a warning, not an error, because it also appears for a group of several molecules that are not in contact (a ligand away from the protein, or a membrane, which is continuous across the box), and WHOLEMOLECULES handles those correctly. The message says so.

• 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 warning. 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 give the warning 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).
• basic/rt-molinfo-whole-ligand and basic/rt-molinfo-whole-membrane: trp-cage with a ligand that is closer to a periodic image of the protein than to the protein, and a small bilayer, each with whole references. With all atoms in one ENTITY, the warning is written and the split frames are still made whole as the reference; with one ENTITY per molecule, there is no warning. Without the check, these tests fail (no warning).

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 flagged 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 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.
• 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 reported, 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 reported, 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: 830 performed, 1 error, basic/rt-make-fortran08, which fails the same way without this change in my clang/libc++ build. All 40 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]

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.

[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. @GiovanniBussi, @carlocamilloni, which would you rather have?

[carlocamilloni]

I think a warning may be useful. @GiovanniBussi was also thinking about intermolecular cases if I remember correctly. Like what would happen if the pdb includes a protein and a ligand?

[GiovanniBussi]

@aalhossary I would also like to know what happens if you have a membrane

[aalhossary]

Thanks @carlocamilloni, @GiovanniBussi, good questions. I tested both. The check only looks at the atoms of each group being made whole (each ENTITY), so the answer depends on how the molecules are grouped.

Protein and ligand. I used 3PTB + benzamidine, 3EML + ZM241385 and 2RH1 + carazolol, in cubic and dodecahedral boxes with 10 and 15 Å margins. The ligand was either bound or, in 160 cases per structure, unbound at a random place in the solvent.

• With the ligand in its own ENTITY, or a bound ligand in the protein’s ENTITY, the reference was always accepted.
• With an unbound ligand in the protein’s ENTITY, it was rejected in 63 of 480 cases: exactly the cases where the ligand is closer to a periodic image of the protein than to the protein itself. The error names the ligand atoms. In all 63, WHOLEMOLECULES without the check gave the right result, so these are false alarms.

Membrane. I used a Martini DPPC/cholesterol bilayer, and the all-atom YiiP transporter in POPE/POPG (both from the MDAnalysis test files).

• With the protein, each chain, or each lipid as its own ENTITY, the reference was accepted. Lipids that are not in the group play no role.
• With all lipids, or protein + lipids, as one ENTITY, it was always rejected, naming thousands of atoms. A bilayer is continuous across the box, so as a group it is never whole. Yet WHOLEMOLECULES without the check made every lipid and the protein whole in every frame, also after each lipid had moved by a few Å. So these are false alarms too.

So the check is reliable for a group that is one molecule, but with several molecules in one group it can stop a run that works today. I’ll make it a warning in the log that names the atoms, as Carlo suggested, unless you prefer otherwise.

[GiovanniBussi]

My main doubt about this PR as a whole is that it is very intrusive (lot of code added) but I am not sure I understand the advantage.

• It works well in simple cases. However, with a single molecule, users are not even expected to use the WHOLE flag. The goal is mostly to flag incorrect use of a broken pdb.
• It can give false positives in complexes, which is exactly where the EMST is used

What about a simpler test:

• run over each chain and, for each chain, check that the chain is already whole using the standard "along the chain" convention. If not, give a warning
• expect that the relative positioning of chains provided by the user is correct

I think this should be only a few lines of code. It should give a warning in the example proposed in #1492 What do you think?

[aalhossary]

Thanks @GiovanniBussi, I agree: that is much simpler. I tried it in a quick script on the reference files of the regtests here. For each chain, I walk the atoms in file order and check every step against the minimum image in the box of the CRYST1 record:

• torn.pdb, a case like MOLINFO WHOLE: a reference split by PBC makes WHOLEMOLECULES break whole molecules, silently #1492, is flagged (14 steps cross the boundary), and the whole protein is not.
• The protein and ligand case is not flagged, because the ligand is in its own chain.
• The bilayer, with all its lipids in one chain, is flagged at 5 steps, all between two different lipids. If I check only steps within a residue, or between two protein or nucleic-acid residues (MolDataClass::allowedResidue), those go away and torn.pdb is still flagged.

If that is OK, I will replace the current check with this one. It stays a warning and is skipped when the CRYST1 cell can't be used, and I will keep the regtests.

[GiovanniBussi]

Yes I think this would be fine. I think it can be done with a very short function so that there are not so many modifications to review. Thanks!!!