Problem of position restrains in combination of `refcoord_scaling = com` + walls
Summary
There seems to be a bug related to the position restraints that APPEARED SINCE GROMACS 2025.2.
We observed that, setting a position restrains when refcoord_scaling = com AND walls are set, leads to a very large energy contribution of the position restraint ('Position Rest.') at step 0.
We observed this behaviour for gromacs 2025.3 and 2026.2, while with gromacs 2024.3 and 2023.3 we do not register the same bug.
For this, we suspect it might be related to the changes introduced in #5337 (closed)
Exact steps to reproduce
The attached files represent a minimal example to illustrate the problem. The system is a single water molecule in a big box, using periodic boundary conditions in X and Y and walls in Z. A position restraint is applied to the oxygen atom of the water molecules. Importantly, pressure coupling is turned off but we still have the option 'refcoord_scaling = com' (this should have no effect since there is no pressure coupling, but it seems to be part of the problem).
Run the simulation as follows:
gmx grompp -f params.mdp -c system.gro -r system.gro -p system.top -o run.tpr
gmx mdrun -v -deffnm run
In the output, you will see that the energy contribution of the position restraint ('Position Rest.') at step 0 is quite large:
Step Time
0 0.00000
Energies (kJ/mol)
LJ (SR) Coulomb (SR) Position Rest. Potential Kinetic En.
7.46448e-08 -7.62939e-06 1.22786e+03 1.22786e+03 1.01575e+01
Total Energy Temperature Pressure (bar) Constr. rmsd
1.23802e+03 6.10832e+02 -2.12431e+02 4.04280e-05This seems wrong, since at step 0 the coordinates correspond exactly to the ones from the restraint (note that in the above 'grompp' command we pass as arguments '-c' and '-r' the same file 'system.gro'), and hence the restraint energy should be 0.
If we remove, in the file params.mdp, the line 'refcoord_scaling = com', the position restraint energy drops to 0 (within numerical precision):
Step Time
0 0.00000
Energies (kJ/mol)
LJ (SR) Coulomb (SR) Position Rest. Potential Kinetic En.
7.46448e-08 -7.62939e-06 2.36362e-08 -7.53111e-06 1.01219e+01
Total Energy Temperature Pressure (bar) Constr. rmsd
1.01219e+01 6.08690e+02 -1.20411e-01 3.86510e-05The same happens if we keep 'refcoord_scaling = com' but remove the walls and put 'pbc = xyz':
Step Time
0 0.00000
Energies (kJ/mol)
LJ (SR) Coulomb (SR) Position Rest. Potential Kinetic En.
0.00000e+00 -7.62939e-06 2.36362e-08 -7.60576e-06 1.01218e+01
Total Energy Temperature Pressure (bar) Constr. rmsd
1.01218e+01 8.11586e+02 -1.19836e-01 4.07534e-05Hence, the problem seems to be related to the combination of walls + refcoord_scaling. In the present case, the application of 'refcoord_scaling' does not make sense, but there might be other situations where walls are combined with a barostat, and hence this bug becomes critical.