Calculation of distance restraint energies is faulty
Summary
Forces and energies for distance restraints are calculated in the file src/gromacs/listed_forces/disre.cpp in the function ta_disres(). While forces are calculated as described in the manual https://manual.gromacs.org/documentation/2019/reference-manual/functions/restraints.html#distance-restraints the calculation of the energy ignores the parameter up2. This means that if a distance restraint is strongly violated, the energy can be vastly overestimated. Since this does not affect the dynamics (forces are calculated correctly) this bug is not terrible. Nevertheless it should be fixed.
I have checked this bug both by numerical experiments as well as by reading the code.
GROMACS version
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 simple example system should take ~30 seconds to run.
What is the current bug behavior?
In the example above, distance restraint energies from ener.edr will be 500*(r-r_1)^2 instead of 500*(r_2−r_1)(2r−r_2−r_1), with r taken from traj_comp.xtc.
What did you expect the correct behavior to be?
Energies should be 500*(r_2−r_1)(2r−r_2−r_1).
Possible fixes
In the file src/gromacs/listed_forces/disre.cpp replace the line (515)
vtot += 0.5 * k0 * gmx::square(tav_viol) * pairFac;
with
if (tav_viol > up2-up1) {
vtot += 0.5 * k0 * (up2-up1) * (2*tav_viol+up1-up2) * pairFac;
} else {
vtot += 0.5 * k0 * gmx::square(tav_viol) * pairFac;
}