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** [gly10_alpha_small.pdb](/uploads/ccdf5fabef0eb0cd98ccdf68f83ac9b0/gly10_alpha_small.pdb) [test.itp](/uploads/ff35b193bf40cb0545a7cad3102a14e7/test.itp) [topol.top](/uploads/e903878391df1a0479870d1be9d7213e/topol.top) [topol.tpr](/uploads/d12e11d61fe5a7917bb4c9906a0b76b4/topol.tpr) [grompp.mdp](/uploads/f9e700615bbfcb7d571942290fc342db/grompp.mdp) 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)*(2*r−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)*(2*r−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; }
issue