Distance restraint force is calculated in an unexpected way in case of negative prefactor

Summary

Forces and for distance restraints are calculated in the file src/gromacs/listed_forces/disre.cpp in the function ta_disres(). The prefactor for the forces is called k0. In an itp file with distance restraints, the ninth entry in any row is a free parameter for the restraint force. Usually, this parameter is supposed to be positive and if this is the case, everything works fine. However, if the parameter is negative, as in this example:

gly10_alpha_small.pdb

grompp.mdp

test.itp

topol.top

topol.tpr

the force calculation is not as in the manual:

https://manual.gromacs.org/documentation/2019/reference-manual/functions/restraints.html#distance-restraints

While negative prefactors to a distance restraint may seem like a nonsensical thing to do, this is something which is of relevance for my use case. At any rate, I think as a matter of principle, the code should do what the manual says even in this case.

Since there seems to be no canonical way of outputting restraint forces, I have checked this bug only by reading the code.

GROMACS version

The releases release-2020, release-2021 and release-2022 are all affected. My version is:

GROMACS version:    2022-beta1-dev-20211028-d0593a5845
GIT SHA1 hash:      d0593a58454b855259dff740cb376a3d43f6c620
Precision:          mixed
Memory model:       64 bit
MPI library:        thread_mpi
OpenMP support:     enabled (GMX_OPENMP_MAX_THREADS = 64)
GPU support:        disabled
SIMD instructions:  AVX2_256
CPU FFT library:    fftw-3.3.8-sse2-avx-avx2-avx2_128
GPU FFT library:    none
RDTSCP usage:       enabled
TNG support:        enabled
Hwloc support:      hwloc-2.1.0
Tracing support:    disabled
C compiler:         /usr/bin/cc GNU 9.3.0
C compiler flags:   -mavx2 -mfma -Wall -Wno-unused -Wunused-value -Wunused-parameter -Wextra -Wno-sign-compare -Wpointer-arith -Wundef -Werror=stringop-truncation -Wno-missing-field-initializers -fexcess-precision=fast -funroll-all-loops -Wno-array-bounds -O3 -DNDEBUG
C++ compiler:       /usr/bin/c++ GNU 9.3.0
C++ compiler flags: -mavx2 -mfma -Wall -Wextra -Wpointer-arith -Wmissing-declarations -Wundef -Wstringop-truncation -Wno-missing-field-initializers -fexcess-precision=fast -funroll-all-loops -Wno-array-bounds -fopenmp -O3 -DNDEBUG

Steps to reproduce

This problem cannot be easily seen from simulation results, however it is obvious from reading the code.

What is the current bug behavior?

In the example above, distance restraint forces will be 1000*(r-r_1) instead of 1000*(r_2−r_1), with r taken from traj_comp.xtc.

What did you expect the correct behavior to be?

Forces should be 1000*(r_2−r_1).

Possible fixes

In the file src/gromacs/listed_forces/disre.cpp replace the following line appears twice (line 564 and 579):

f_scal = std::max(f_scal, fmax_scal);

The bug could be resolved by replacing this line on both occurrences with

if (-f_scal/k0 > up2-up1) {
    f_scal = fmax_scal;
}

(The somewhat awkward term -f_scal/k0 is necessary to account for all the different possibilities for the calculation of f_scal in the function.)