Make the coarse-graining deterministic and fix the dummy-bead geometry (#1657) - #1660
Open
amjjbonvin wants to merge 13 commits into
Open
Make the coarse-graining deterministic and fix the dummy-bead geometry (#1657)#1660amjjbonvin wants to merge 13 commits into
amjjbonvin wants to merge 13 commits into
Conversation
The martinize function placed the SCD* dummy beads along a vector drawn from the global `random` module, which is never seeded. Every topocg run therefore produced different CG coordinates for the same input, so a repeated workflow could not be reproduced. martinize() now takes a `seed` and builds its own random.Random, threaded down through map_cg() and add_dummy(). topocg passes the `iniseed` parameter from its defaults.yaml, so the coarse-graining honours the same user-facing seed as the CNS side; iniseed was previously inert for topocg, whose CNS script has nothing stochastic in it. A fresh generator per call keeps a given structure independent of how many were converted before it, and the global random stream is left untouched. caprieval and caprifilter have no iniseed of their own and fall back to libaa2cg.DEFAULT_SEED, kept in sync with the topocg default. While in add_dummy, three related problems: - The bead distances were in nm while everything around them is in angstrom. 0.11 nm is exactly the 1.100 A parent-SCd bond of protein-CG-Martini-2-2.param, and a symmetric pair at 0.14 nm is exactly the 2.800 A SCd-SCd bond, so the values are now 1.1 and 1.4 A. The pair was also placed asymmetrically (+d and -2d), which cannot satisfy both bonds; the beads now straddle the parent at +d and -d. This was not merely cosmetic: generate-topology.cns only minimises if the bonded energy exceeds 10000, so small systems kept the 10x too short geometry. On tests/golden_data/protein.pdb the bonded energy drops from 8000.3 (below the threshold, so never relaxed) to 429.1, with the beads at the force field distances. - The direction was sampled from a cube and normalised, biasing it towards the cube corners. Normalising a Gaussian vector instead is isotropic: over 200k samples the fraction of directions with |z| > 0.9 goes from 0.0615 to 0.1006, against a uniform-on-sphere expectation of 0.1. - The placement comprehension shadowed its loop variable, and only gave the intended result because of that shadowing. Also drop output_cg_restraints() and the never-called extract_groups(). They wrote dna_restraints.def and dna-aa_groups.dat into the current working directory, so each molecule overwrote the previous one's file, and nothing in HADDOCK ever read either. determine_hbonds() is kept: its return value is now unused but it still marks paired bases with bfactor 1, which patch-types-cg-hbond-dna-rna.cns reads to apply the h-bond patches. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
cgtoaa.cns assigns Maxwell velocities at 300 K and runs four rounds of Cartesian MD, but never ran `set seed`, unlike every other CNS recipe that samples. The seed prepared by Python was therefore ignored: CNS fell back to its own default and every replica retraced the same trajectory, so `sampling_factor > 1` burned CPU producing duplicate models. The `+ s_ind` offset in the module was already there to decorrelate the replicas but had no effect. Adding `set seed $seed end` after the log-level block, before the first `maxwell` call, makes the seed effective. `prepare_cns_input` emits `evaluate ($seed=...)` ahead of the recipe, so the symbol is defined by then. Measured with the real CNS binary on e2a-hpr_1GGR_A (1503 atoms) after a full topoaa -> topocg -> cgtoaa chain, running two replicas with seeds 917 and 918: without the fix RMSD 0.000 A, 0 of 1503 atoms moved with the fix RMSD 0.401 A, max 2.814 A, all 1503 atoms moved Only the FILENAME and DATE remarks differed before the fix; the coordinates were bit-identical. Note the effect scales with system size: on a 63-atom stub the two replicas differ by an RMSD of just 0.0009 A, because the short restrained MD cannot diverge much. tests/test_cns_seeding.py guards the invariant behind this bug rather than the single line: any top-level CNS recipe that draws random velocities (`maxwell`) or runs `dynamics` must also run `set seed`. Recipes with no stochastic step are skipped, which is why emscoring and topocg are not flagged for carrying an unused seed. The test fails on the pre-fix cgtoaa.cns. The example config exercises the path end to end with sampling_factor = 5. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Follow-up cleanup on the two preceding commits. No behaviour change: the CG model of a given input is byte-identical before and after, and the dummy beads still land on the force field distances (1.400 / 2.800 / 1.100 A measured on e2a-hpr_1GGR_A). test_cns_seeding checked the wrong scope. It scanned only <module>/cns/<module>.cns, so it skipped flexref and rigidbody reporting "no stochastic step" when both do sample - flexref anneals in five sa_ltad_*.cns includes, rigidbody draws rotations via random() in get_random_rotation.cns - and it never saw topoaa or topocg at all, whose entry point is generate-topology.cns. cgtoaa inlining its MD is the outlier; the convention in this repo is to @@include it, so a later refactor of cgtoaa.cns towards that convention would have flipped the test to SKIP and let the bug back in unnoticed. It now checks per module tree, asks whether any recipe in a module samples and whether any recipe in that module runs `set seed`, and recognises CNS's own random() alongside maxwell and dynamics. Seven modules now assert instead of three; the three remaining skips (emref, emscoring, topocg) are genuinely free of stochastic steps. Still fails on the pre-fix cgtoaa.cns. Discovery comes from `modules_category`/`modules_folder` rather than a filesystem glob and a `__file__` walk, which also drops a pytest generator-parametrize deprecation warning. add_dummy: drop the `while True:` retry guarding a zero-norm Gaussian vector, an unreachable branch no test can cover; hoist the scale factor out of the comprehension; replace `sign = 1 if idx % 2 == 0 else -1` with `enumerate((1, -1)[:n])`, so the 1-or-2 limit is visible rather than implied by a modulo that would silently overwrite SCD1 for n=3. Its `dist=1.1, n=2` defaults are now required arguments: every caller passed both explicitly, and that default pairing matched neither real case (polar is 1.4/2, charged 1.1/1), so it described a 2.2 A pair the force field never asks for. generate_topology's `seed` is likewise required. Its default could only fire from a test, and dropping it removes the topocg -> libaa2cg DEFAULT_SEED import, so the constant no longer travels back into the module that feeds it. DEFAULT_SEED's comment claimed it was kept in sync with the topocg `iniseed` default with nothing enforcing that; a test now asserts it. Tests: hoist the CNS-engine and previous-io fakes into a _prepared_topocg helper, removing about 35 lines duplicated between the two _run tests; use math.dist instead of a hand-rolled euclidean helper that shadowed the `dist=` keyword; state the force field provenance once, in map_cg where the literals live, instead of in three places that must be updated together. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
amjjbonvin
requested review from
AnnaKravchenko,
Comp-era,
VGPReys,
rversin and
rvhonorato
August 8, 2026 08:27
VGPReys
reviewed
Aug 10, 2026
Co-authored-by: Victor Reys <132575181+VGPReys@users.noreply.github.com>
VGPReys
reviewed
Aug 10, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
3 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
What does this PR do and why?
Makes the coarse-graining path reproducible, and fixes the bead geometry it produces along the way (simplification/optimisation of libaa2cg)
1.
libaa2cg.martinize()was non-deterministic.add_dummy()placed theSCD*dummy beads along a direction drawn from the globalrandommodule, which is never seeded. Two runs on the same input therefore produced different CG models.martinize()now takes aseedand uses a dedicatedrandom.Random(seed)instead of the global stream, so it can no longer be perturbed by whatever else drew a random number first in the same interpreter.topocgpasses its owniniseed(fromdefaults.yaml) down throughgenerate_topology(); the default inlibaa2cgis kept in sync with it for the callers that have noiniseedof their own (caprieval,caprifilter). This also givestopocg'siniseedan actual effect — its CNS recipe has no stochastic step, so until now the parameter was inert.2.
cgtoaaignored the seed it was given.prepare_cns_inputwritesevaluate ($seed=...)into every.inp, andcgtoaa/__init__.pyvaries it per replica (iniseed + s_ind), butcgtoaa.cnsnever ranset seed— unlike every other CNS module. It runsmaxwell+dynamics cartesian, so CNS fell back to its own default and every replica retraced the same trajectory:sampling_factor > 1burned CPU for identical models. Added the missingset seed $seed end.3. Dummy-bead distances were in nm, not Å.
add_dummyuseddist=0.11/0.14, i.e. the MARTINI native unit, while HADDOCK/CNS works in Å. The minimisation was correcting this. Now it is correct from the start.4. Direction sampling was biased. The direction came from normalising a vector drawn uniformly from a cube, which concentrates towards the cube's corners. Replaced with a Gaussian vector, which is isotropic under normalisation.
5. A pair of dummy beads was not symmetric about its parent. The
m = 1; ...; m *= -2idiom reused the loop variable, so the second bead landed at-2drather than-d. The pair now straddles the parent bead.6.
extract_groups()reada[0][1]where it meanta[1][1], so both segids came from the same side of the pair. Dead code today (no callers), fixed rather than left as a trap.Also adds
examples/refine-complex/refine-complex-CG-test.cfgexercisingtopocg→cgtoaawithsampling_factor = 5. This example is similar to the protocol used in the Lightdock/HADDOCK membrane docking paper to refine the lightdock modelsHow was this tested?
ruff checkclean;pytest tests/→ 1699 passed, 8 skipped.New regression tests in
tests/test_module_topocg.py:martinizeis reproducible across calls, and changes with a different seed (the atoms that differ are exactly{SCD1, SCD2})martinizeleaves the globalrandomstream untouched, and is independent of the order in which molecules are processedtopocg._run()actually forwardsself.params["iniseed"]togenerate_topologyadd_dummydistances match theSCdforce-field equilibriaadd_dummydirections are isotropictests/test_cns_seeding.pyis new and generic: it walks the module registry and asserts that any module whose CNS recipes containmaxwell,dynamicsorrandom()also runsset seedsomewhere.cgtoaawas the only module failing it.Running twice the new example r
efine-complex-CG-test.cfgresults in exactly the same models at the end of the workflow. Fully reproducible.AI assistance
Claude Code was used to trace the source of the non-determinism, to write the fixes and the regression tests, and to draft this description which was manually edited.
Checklist
CHANGELOG.mdupdated for user-facing changesRelated issues
Issue #1657
Notes for reviewers
topocg/cns/generate-topology.cns:240: theif ($bonded > 10000)guard that triggers an extra minimisation may no longer fire. The dummy beads sitting ~10× too close to their parent were plausibly the dominant bonded-energy term, and with the geometry corrected the initial bonded energy should drop sharply. Worth a look from someone who knows what that threshold was calibrated against.