Fix MTS constraint dH/dlambda contribution

With MTS the slow force contributions to dH/dlambda were
multiplied by the MTS factor.

Fixes #5060

Author: berkhess

docs/release-notes/2023/2023.5.rst

@@ -63,6 +63,17 @@ for this in grompp was lost. Now both grompp and mdrun will exit with an error.
 :issue:`3775`
 :issue:`5053`
 
+Correct dH/dlambda with MTS and perturbed constraints
+"""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+When doing free-energy calculations with multiple time stepping and perturbed
+constraint lengths, the slow force contributions to dH/dlambda were incorrectly
+multiplied by the MTS factor. The effect of this is usually rather small
+when only integrating the PME grid part less frequently, but can be substantial
+when the non-bonded pair interactions are not integrated at every step.
+
+:issue:`5060`
+
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 

src/gromacs/mdrun/md.cpp

@@ -1588,12 +1588,14 @@ void gmx::LegacySimulator::do_md()
                 else
                 {
                     /* With multiple time stepping we need to do an additional normal
-                     * update step to obtain the virial, as the actual MTS integration
+                     * update step to obtain the virial and dH/dl, as the actual MTS integration
                      * using an acceleration where the slow forces are multiplied by mtsFactor.
                      * Using that acceleration would result in a virial with the slow
                      * force contribution would be a factor mtsFactor too large.
                      */
-                    if (simulationWork.useMts && bCalcVir && constr != nullptr)
+                    const bool separateVirialConstraining =
+                            (simulationWork.useMts && (bCalcVir || computeDHDL) && constr != nullptr);
+                    if (separateVirialConstraining)
                     {
                         upd.update_for_constraint_virial(*ir,
                                                          md->homenr,
@@ -1643,8 +1645,8 @@ void gmx::LegacySimulator::do_md()
                                           step,
                                           state,
                                           upd.xp()->arrayRefWithPadding(),
-                                          &dvdl_constr,
-                                          bCalcVir && !simulationWork.useMts,
+                                          separateVirialConstraining ? nullptr : &dvdl_constr,
+                                          bCalcVir && !separateVirialConstraining,
                                           shake_vir);
 
                     upd.update_sd_second_half(*ir,