fix: preserve molecular geometry during manostat scaling

[galjos]

Fixes the MOF-5 stochastic-rescaling manostat distortion reported in review.

What changed:

• Reconstruct each molecule into its current COM image before a manostat box resize, so molecules split by PBC do not carry an old box-length jump through the resize.
• Apply the same reconstruction before Berendsen scaling, because the coordinate invariant is shared even though Berendsen did not show the same large MOF-5 failure.
• Keep stochastic-rescaling velocity updates at the molecular COM level, preserving internal velocities.
• Keep isotropic stochastic mu as length scaling: exp((-compress * deltaP + stochasticFactor) / 3).

Validation:

• Real CUDA/MACE MOF-5 run-02, stochastic/rescale manostat, 50000 steps / 5000 saved frames on a GPU node.
• Clean dev 13a82748: cut mean 1.0943499541 A, min/max 1.0925817400 / 1.0960428055 A; inside mean 1.0918005631 A, min/max 1.0904312459 / 1.0934742641 A; gap 0.0025493910 A.
• Candidate a686892e: cut mean 1.0932610733 A, min/max 1.0922109993 / 1.0972167131 A; inside mean 1.0928381316 A, min/max 1.0907006212 / 1.0970055845 A; gap 0.0004229417 A.
• Both simulations ended normally at final step 660000.

Validator/plots package kept out of this PR tree:

• https://raw.githubusercontent.com/MolarVerse/PQ/pr261-mof5-validator-artifact/pr261-mof5-validator.zip
• SHA256: af944f5e6c411610e704063c03036e1c8f089b193b565c7cc6098d79c4a72e70
• Contains analyze_mof5_ch.py, the comparison PNGs, and CSV/JSON summaries.

Local checks:

• ctest --test-dir build-pr261-publish-static2 -R '^(testMolecule|testManostat)$' --output-on-failure
• git diff --cached --check

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 88.41%. Comparing base (1d02915) to head (28ab3d7).
⚠️ Report is 2 commits behind head on dev.

Additional details and impacted files

@@            Coverage Diff             @@
##              dev     #261      +/-   ##
==========================================
+ Coverage   88.12%   88.41%   +0.28%     
==========================================
  Files         283      283              
  Lines       10874    10910      +36     
  Branches     3376     3381       +5     
==========================================
+ Hits         9583     9646      +63     
+ Misses       1251     1224      -27     
  Partials       40       40              

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[pq-perf-bot]

⚡ Performance (instruction count) — ✅ no regressions

per-benchmark breakdown

benchmark	base Ir	PR Ir	Δ
bondedForces	41.46M	41.46M	-0.00%	✅
boxTransforms	10.67M	10.67M	+0.00%	✅
constraints	11.38M	11.38M	+0.00%	✅
coulombKernel	5.80M	5.80M	+0.00%	✅
forceKernel	15.33M	15.33M	+0.00%	✅
integrator	32.17M	32.17M	+0.00%	✅
kinetics	9.53M	9.53M	+0.00%	✅
linearAlgebra	2.08M	2.08M	+0.00%	✅
nonCoulombPairs	5.48M	5.48M	-0.00%	✅
shiftVector	5.81M	5.81M	+0.00%	✅
virial	11.63M	11.63M	-0.00%	✅

_{Deterministic callgrind instruction counts vs the base branch; gated at ±2%. Not wall-clock.}

[galjos]

Real-PQ reproducer zip:

https://github.com/MolarVerse/PQ/raw/test/pr-261-reproducer/artifacts/pr-261-real-pq-reproducer.zip

Includes native one-step and 1000-step PQ inputs, run_reproducer.sh, long/run_long_reproducer.sh, an ASE wrapping reference (long/ase_reference.py), a plot script (long/plot_long_comparison.py), and captured old/fixed outputs under results/.

Force path: normal PQ MM brute-force path (cell-list = off), shifted Coulomb with zero charges, and a zeroed GUFF self-pair. Forces are intentionally zero; the cases isolate manostat/PBC scaling.

Same one-step input, observed with real PQ:

• pre-fix: FAIL, H2 position (-1.38333333, -0.01666667, 0.0) has x-coordinate outside the [-0.5, 0.5) box range
• this PR: PASS, H2 position (-0.38333333, -0.01666667, 0.0) is inside the box

1000-step input, observed with real PQ and ASE reference:

• pre-fix: FAIL, H2 position (0.71642257, -0.00945294, 0.0) is just outside the final centered box range
• this PR: PASS, H2 position (-0.71640107, -0.00945294, 0.0) matches results/long-ase-reference-final.xyz at written precision
• plot: results/long-h2-x-comparison.png; data: results/long-h2-x-comparison.csv

Run with:
./run_reproducer.sh /path/to/PQ
cd long && ./run_long_reproducer.sh /path/to/PQ
cd long && ./ase_reference.py
cd long && ./plot_long_comparison.py

[ape33]

This does not fix the problem at hand. I have attached example input of a MOF-5 simulation, where calculating the shake distances for the second run shows the issue. All C-H shake distances are associated with linker molecules. The shake distances of linkers cut by the box are systematically longer than those of linkers that are fully inside the box (see attached graph: shake.pdf). This is only the case when using the stochastic rescaling manostat. Using the exact same input, but switching to the Berendsen manostat fixes this issue (comparison graph will follow).
Here the example input, including a script for calculating the average distances:
stoch_res.zip

[galjos]

@ape33 you were right: the wrapping fix was not the stochastic-only bug.

The MOF-5 case points to the extra velocity scaling in StochasticRescalingManostat. That path scaled every atom velocity by inverse(mu), while positions are scaled by molecule COM. In c5b8c822 this now scales only the molecule COM velocity and keeps internal relative velocities unchanged.

So: PBC wrapping stays as a separate shared fix, but this velocity change is the actual fix for your stochastic-rescaling case. CI is green, and the new tests cover that invariant.

[ape33]

@ape33 you were right: the wrapping fix was not the stochastic-only bug.

The MOF-5 case points to the extra velocity scaling in StochasticRescalingManostat. That path scaled every atom velocity by inverse(mu), while positions are scaled by molecule COM. In c5b8c822 this now scales only the molecule COM velocity and keeps internal relative velocities unchanged.

So: PBC wrapping stays as a separate shared fix, but this velocity change is the actual fix for your stochastic-rescaling case. CI is green, and the new tests cover that invariant.

@galjos
Does the MOF-5 example I have provided produce consistent shake distances?

[ape33]

@galjos
the issue persists, even after the newest commits

[ape33]

@galjos
using the exact same input, but switching to the Berendsen manostat results in the following shake distances with expected behavior
shake.pdf

[galjos]

I replaced the earlier velocity-only approach with the molecule-reconstruction fix and reran the real MOF-5 CUDA/MACE validation.

50000 steps / 5000 saved frames, stochastic/rescale manostat:

• clean dev 13a82748: cut mean 1.0943499541 A, inside mean 1.0918005631 A, gap 0.0025493910 A
• candidate a686892e: cut mean 1.0932610733 A, inside mean 1.0928381316 A, gap 0.0004229417 A

The validator package and comparison PNG/CSV/JSON outputs are intentionally kept out of this PR tree and are available as a separate zip:
https://raw.githubusercontent.com/MolarVerse/PQ/pr261-mof5-validator-artifact/pr261-mof5-validator.zip

SHA256: af944f5e6c411610e704063c03036e1c8f089b193b565c7cc6098d79c4a72e70

[galjos]

@ape33 I moved the validator/plots out of the PR tree so they will not land in dev or main.

Download:
https://raw.githubusercontent.com/MolarVerse/PQ/pr261-mof5-validator-artifact/pr261-mof5-validator.zip

SHA256: af944f5e6c411610e704063c03036e1c8f089b193b565c7cc6098d79c4a72e70

The zip contains analyze_mof5_ch.py, the clean-dev vs candidate comparison PNGs, and the CSV/JSON summaries from the 50000-step CUDA/MACE MOF-5 validation.

[ape33]

please do not merge yet, I would like to run tests first

[ape33]

Nice, these updates seem to fix the issue!
However, is the scaling of the velocities based on the COM necessary and even meaningful? For the velocity rescaling thermostat for example also an atom-based scaling is employed.
On this I would like to hear the opinion of @97gamjak

[galjos]

Replying to ape33's review: #261 (review)

Good question. I think the meaningfulness comes from matching the velocity update to the coordinate DOF used by the barostat. After this PR, the manostat scales positions via molecule COMs to preserve constrained geometry, so the reciprocal stochastic velocity scaling should act on the molecule COM velocity as well. Atom-wise scaling would change intramolecular relative velocities while the coordinate update intentionally leaves them unchanged.

This is also how I read the comparison with other engines: GROMACS C-rescale scales atom velocities, but its coordinate scaling path is atom-based too; OpenMM/Amber/LAMMPS use molecule/rigid-body centers when rigid molecular geometry is what is being preserved. The velocity-rescaling thermostat is different because it intentionally acts on atomic thermal velocities.

[ape33]

Replying to ape33's review: #261 (review)

Good question. I think the meaningfulness comes from matching the velocity update to the coordinate DOF used by the barostat. After this PR, the manostat scales positions via molecule COMs to preserve constrained geometry, so the reciprocal stochastic velocity scaling should act on the molecule COM velocity as well. Atom-wise scaling would change intramolecular relative velocities while the coordinate update intentionally leaves them unchanged.

This is also how I read the comparison with other engines: GROMACS C-rescale scales atom velocities, but its coordinate scaling path is atom-based too; OpenMM/Amber/LAMMPS use molecule/rigid-body centers when rigid molecular geometry is what is being preserved. The velocity-rescaling thermostat is different because it intentionally acts on atomic thermal velocities.

Sounds reasonable to me, thank you for the detailed explanation!
Still I would like to hear @97gamjak's opinion before merging the PR.

[galjos]

@ape33 me too