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-05

This 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-05

The 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-05

Hence, 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.

position_restraints_water.itp

params.mdp

system.gro

system.top

Edited by Emanuele Bosoni