diff --git a/CHANGES/v2.9.md b/CHANGES/v2.9.md
index 3d6306ddd8..dbe80862c3 100644
--- a/CHANGES/v2.9.md
+++ b/CHANGES/v2.9.md
@@ -111,13 +111,14 @@ Changes from version 2.8 which are relevant for users:
     on virtual sites were lost, in MD and in energy minimization. They are now spread on the constructing atoms.
   - GROMACS 2021, 2022, 2023 and 2024 patches: with multiple time stepping, the forces that PLUMED applied
     when it needed the potential energy were lost on the slow steps, and the potential energy passed to PLUMED
-    on the fast steps lacked the slow terms. The forces are now applied to the buffer used by the update,
-    and mdrun stops with an error if PLUMED needs the energy on a fast step.
+    on the fast steps lacked the slow terms. The forces are now applied to the buffer used by the update (and to
+    the one that GROMACS uses for the constraint virial and writes to the trajectory), and mdrun stops with an error
+    if PLUMED needs the energy on a fast step.
   - GROMACS patches: when PLUMED asked to stop the simulation (e.g. with \ref COMMITTOR), multiple simulations
     (e.g. multiple walkers) hung, and with GROMACS 2022-2024 no checkpoint was written at the last step. The request
     now goes through the GROMACS stop handler, so that it reaches all the ranks and the simulations, and mdrun
-    stops at the next step at which GROMACS communicates signals (every nstglobalcomm steps, and at least every 200
-    steps between multiple simulations), writing a checkpoint there.
+    runs one more step after the next step at which GROMACS communicates signals (every nstglobalcomm steps, and
+    between multiple simulations every 200 steps or more), writing a checkpoint at that last step.
   - GROMACS 2020 and 2024 patches: mdrun -rerun with -plumed stopped with an error, as PLUMED was never initialized.
     The rerun is now patched as for the other GROMACS versions.
   - GROMACS patches: in energy minimization with domain decomposition, PLUMED was not told when the atoms were
@@ -126,12 +127,13 @@ Changes from version 2.8 which are relevant for users:
   - GROMACS patches: with -hrex, PLUMED calculated and updated a second time, for the previous step, when GROMACS
     computed the energy of the swapped coordinates: e.g. \ref METAD deposited a second hill, and \ref PRINT printed
     a second line, on the step before an exchange whenever their stride divided it.
-  - GROMACS patches: with shell (polarizable) particles, PLUMED was silently not applied in MD and in reruns.
+  - GROMACS patches: with shell (polarizable) particles or flexible constraints, PLUMED was silently not applied
+    in MD and in reruns.
     mdrun now stops with an error.
-  - GROMACS patches (since the GROMACS 2019 patch): the virial of the forces applied by PLUMED was counted twice,
-    so the pressure was wrong, e.g. in NPT simulations, with any bias that applies forces, unless the PLUMED input
-    needed the potential energy. On the steps when GROMACS computes the virial, PLUMED now runs after do_force(),
-    as it does when it needs the energy.
+  - GROMACS patches (and those for older GROMACS versions, at least since 4.6): the virial of the forces applied
+    by PLUMED was counted twice, so the pressure was wrong, e.g. in NPT simulations, with any bias that applies
+    forces, unless the PLUMED input needed the potential energy. On the steps when GROMACS computes the virial,
+    PLUMED now runs after do_force(), as it does when it needs the energy.
   - GROMACS 2024 patch: with -hrex, the energy of the swapped coordinates was computed with the settings of the
     previous step (often without the energy), so the exchanges were almost never accepted.
   - GROMACS 2024 patch: in expanded ensemble simulations, the `lambda` passed to PLUMED lagged one expanded
diff --git a/patches/gromacs-2020.7.diff/src/gromacs/mdlib/force.cpp b/patches/gromacs-2020.7.diff/src/gromacs/mdlib/force.cpp
index 7349417bfd..627c82980e 100644
--- a/patches/gromacs-2020.7.diff/src/gromacs/mdlib/force.cpp
+++ b/patches/gromacs-2020.7.diff/src/gromacs/mdlib/force.cpp
@@ -76,6 +76,7 @@
 /* PLUMED */
 #include "../../../Plumed.h"
 int    plumedswitch=0;
+int    plumedIsUsed=0; /* unlike plumedswitch, not switched off while PLUMED runs after do_force() */
 plumed plumedmain;
 /* END PLUMED */
 
diff --git a/patches/gromacs-2020.7.diff/src/gromacs/mdlib/sim_util.cpp b/patches/gromacs-2020.7.diff/src/gromacs/mdlib/sim_util.cpp
index ccd84a9a7f..91fadd3e4d 100644
--- a/patches/gromacs-2020.7.diff/src/gromacs/mdlib/sim_util.cpp
+++ b/patches/gromacs-2020.7.diff/src/gromacs/mdlib/sim_util.cpp
@@ -115,7 +115,7 @@
 #include "gromacs/utility/sysinfo.h"
 
 /* PLUMED */
-extern int plumedswitch;
+extern int plumedIsUsed;
 /* END PLUMED */
 
 using gmx::AtomLocality;
@@ -800,7 +800,7 @@ static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec&
             domainWork.haveSpecialForces || domainWork.haveCpuListedForceWork
             || domainWork.haveFreeEnergyWork || simulationWork.useCpuNonbonded || simulationWork.useCpuPme
             || simulationWork.haveEwaldSurfaceContribution || inputrec.nwall > 0
-            || plumedswitch; /* PLUMED: it adds its forces to the CPU force buffer */
+            || plumedIsUsed; /* PLUMED: it adds its forces to the CPU force buffer */
 
     return domainWork;
 }
diff --git a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/legacymdrunoptions.cpp b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
index a6ff208900..1343f73d1b 100644
--- a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
+++ b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
@@ -177,6 +177,7 @@ int LegacyMdrunOptions::updateFromCommandLine(int argc, char** argv, ArrayRef<co
       //
       plumed_cmd(plumedmain,"setPlumedDat",ftp2fn(efDAT,static_cast<int>(filenames.size()), filenames.data()));
       plumedswitch=1;
+      plumedIsUsed=1;
     }
     /* PLUMED HREX*/
     if(getenv("PLUMED_HREX")) plumed_hrex=1;
diff --git a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/legacymdrunoptions.h b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/legacymdrunoptions.h
index 966e5a41a8..eac77318b4 100644
--- a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/legacymdrunoptions.h
+++ b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/legacymdrunoptions.h
@@ -60,6 +60,7 @@
 /* PLUMED */
 #include "../../../Plumed.h"
 extern int    plumedswitch;
+extern int    plumedIsUsed;
 extern plumed plumedmain; 
 /* END PLUMED */
 
diff --git a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/md.cpp b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/md.cpp
index 2ab4e4b4c0..5ab4b4ca57 100644
--- a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/md.cpp
+++ b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/md.cpp
@@ -703,8 +703,8 @@ void gmx::LegacySimulator::do_md()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
@@ -750,7 +750,7 @@ void gmx::LegacySimulator::do_md()
                 });
       }
       if(vsite){
-        /* The forces that PLUMED applies on virtual sites in the energy path are spread below */
+        /* When PLUMED runs after do_force(), the forces it applies on virtual sites are spread below */
         plumedIsVirtualSite.assign(top_global->natoms,false);
         for(const AtomProxy atomP : AtomRange(*top_global)){
           plumedIsVirtualSite[atomP.globalAtomNumber()]=(atomP.atom().ptype==eptVSite);
diff --git a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/rerun.cpp b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/rerun.cpp
index 0137808e67..3a5b7282a3 100644
--- a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/rerun.cpp
+++ b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/rerun.cpp
@@ -418,8 +418,8 @@ void gmx::LegacySimulator::do_rerun()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
diff --git a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/runner.cpp b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/runner.cpp
index 4e3353c4ca..909616ec97 100644
--- a/patches/gromacs-2020.7.diff/src/gromacs/mdrun/runner.cpp
+++ b/patches/gromacs-2020.7.diff/src/gromacs/mdrun/runner.cpp
@@ -1599,7 +1599,7 @@ int Mdrunner::mdrunner()
           if(useModularSimulator) gmx_fatal(FARGS, "PLUMED is not yet compatible with GROMACS new modular simulator");
 #if GMX_THREAD_MPI
           /* all the thread-MPI ranks would share the same PLUMED object */
-          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: use GROMACS built with MPI");
+          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: run with -ntmpi 1, or use GROMACS built with MPI");
 #endif
           /* detect plumed API version */
           int pversion=0;
diff --git a/patches/gromacs-2021.7.diff/src/gromacs/mdlib/sim_util.cpp b/patches/gromacs-2021.7.diff/src/gromacs/mdlib/sim_util.cpp
index c4e6803e06..a5d5d3612a 100644
--- a/patches/gromacs-2021.7.diff/src/gromacs/mdlib/sim_util.cpp
+++ b/patches/gromacs-2021.7.diff/src/gromacs/mdlib/sim_util.cpp
@@ -138,6 +138,7 @@ using gmx::StepWorkload;
 /* PLUMED */
 #include "../../../Plumed.h"
 extern int    plumedswitch;
+extern int    plumedIsUsed;
 extern plumed plumedmain;
 /* END PLUMED */
 
@@ -882,7 +883,7 @@ static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec&
             domainWork.haveSpecialForces || domainWork.haveCpuListedForceWork
             || domainWork.haveFreeEnergyWork || simulationWork.useCpuNonbonded || simulationWork.useCpuPme
             || simulationWork.haveEwaldSurfaceContribution || inputrec.nwall > 0
-            || plumedswitch; /* PLUMED: it adds its forces to the CPU force buffer */
+            || plumedIsUsed; /* PLUMED: it adds its forces to the CPU force buffer */
 
     return domainWork;
 }
diff --git a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/legacymdrunoptions.cpp b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
index 8706772915..9376c0333e 100644
--- a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
+++ b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
@@ -178,6 +178,7 @@ int LegacyMdrunOptions::updateFromCommandLine(int argc, char** argv, ArrayRef<co
       //
       plumed_cmd(plumedmain,"setPlumedDat",ftp2fn(efDAT,static_cast<int>(filenames.size()), filenames.data()));
       plumedswitch=1;
+      plumedIsUsed=1;
     }
     /* PLUMED HREX*/
     if(getenv("PLUMED_HREX")) plumed_hrex=1;
diff --git a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/legacymdrunoptions.h b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/legacymdrunoptions.h
index 13ee9b89f9..f27c702446 100644
--- a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/legacymdrunoptions.h
+++ b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/legacymdrunoptions.h
@@ -60,6 +60,7 @@
 /* PLUMED */
 #include "../../../Plumed.h"
 extern int    plumedswitch;
+extern int    plumedIsUsed;
 extern plumed plumedmain; 
 /* END PLUMED */
 
diff --git a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/md.cpp b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/md.cpp
index 9b6b813eab..4933e54bb2 100644
--- a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/md.cpp
+++ b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/md.cpp
@@ -720,8 +720,8 @@ void gmx::LegacySimulator::do_md()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
@@ -767,7 +767,7 @@ void gmx::LegacySimulator::do_md()
                 });
       }
       if(vsite){
-        /* The forces that PLUMED applies on virtual sites in the energy path are spread below */
+        /* When PLUMED runs after do_force(), the forces it applies on virtual sites are spread below */
         plumedIsVirtualSite.assign(top_global->natoms,false);
         for(const AtomProxy atomP : AtomRange(*top_global)){
           plumedIsVirtualSite[atomP.globalAtomNumber()]=(atomP.atom().ptype==eptVSite);
@@ -1206,10 +1206,13 @@ void gmx::LegacySimulator::do_md()
             if(plumedswitch){
               if(plumedAfterForces){
                 ArrayRef<RVec> forceForPlumed=f.view().force();
+                std::vector<RVec> forceBeforePlumed;
                 if(fr->useMts && step % ir->mtsLevels[1].stepFactor == 0){
-                  /* On the slow steps the update uses the combined MTS force buffer */
+                  /* On the slow steps the update uses the combined MTS force buffer, and the constraint virial
+                     the normal one: PLUMED's forces are added to the latter below */
                   forceForPlumed=f.view().forceMtsCombined();
                   plumed_cmd(plumedmain,"setForces",&forceForPlumed[0][0]);
+                  forceBeforePlumed.assign(forceForPlumed.begin(),forceForPlumed.begin()+mdatoms->homenr);
                 }
                 if(plumedNeedsEnergy){
                   msmul(force_vir,2.0,plumed_vir);
@@ -1230,6 +1233,10 @@ void gmx::LegacySimulator::do_md()
                                       VirtualSitesHandler::VirialHandling::NonLinear,{},vsiteVir,nrnb,state->box,wcycle);
                   m_add(force_vir,vsiteVir,force_vir);
                 }
+                if(!forceBeforePlumed.empty()){
+                  ArrayRef<RVec> forceNormal=f.view().force();
+                  for(int i=0;i<mdatoms->homenr;i++) forceNormal[i]+=forceForPlumed[i]-forceBeforePlumed[i];
+                }
               } else {
                 msmul(plumed_vir,0.5,plumed_vir);
                 m_add(force_vir,plumed_vir,force_vir);
diff --git a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/rerun.cpp b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/rerun.cpp
index 28ff670f9d..4f594f36fb 100644
--- a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/rerun.cpp
+++ b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/rerun.cpp
@@ -396,8 +396,8 @@ void gmx::LegacySimulator::do_rerun()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
diff --git a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/runner.cpp b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/runner.cpp
index 7e9461d029..a6ef5cee92 100644
--- a/patches/gromacs-2021.7.diff/src/gromacs/mdrun/runner.cpp
+++ b/patches/gromacs-2021.7.diff/src/gromacs/mdrun/runner.cpp
@@ -175,6 +175,7 @@
 /* PLUMED */
 #include "../../../Plumed.h"
 int    plumedswitch;
+int    plumedIsUsed=0; /* unlike plumedswitch, not switched off while PLUMED runs after do_force() */
 plumed plumedmain; 
 /* END PLUMED */
 
@@ -1769,7 +1770,7 @@ int Mdrunner::mdrunner()
           if(useModularSimulator) gmx_fatal(FARGS, "PLUMED is not yet compatible with GROMACS new modular simulator");
 #if GMX_THREAD_MPI
           /* all the thread-MPI ranks would share the same PLUMED object */
-          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: use GROMACS built with MPI");
+          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: run with -ntmpi 1, or use GROMACS built with MPI");
 #endif
           /* detect plumed API version */
           int pversion=0;
diff --git a/patches/gromacs-2022.5.diff/src/gromacs/mdlib/sim_util.cpp b/patches/gromacs-2022.5.diff/src/gromacs/mdlib/sim_util.cpp
index b1b4607852..b9a15038c8 100644
--- a/patches/gromacs-2022.5.diff/src/gromacs/mdlib/sim_util.cpp
+++ b/patches/gromacs-2022.5.diff/src/gromacs/mdlib/sim_util.cpp
@@ -137,6 +137,7 @@ using gmx::StepWorkload;
 /* PLUMED */
 #include "../../../Plumed.h"
 extern int    plumedswitch;
+extern int    plumedIsUsed;
 extern plumed plumedmain;
 /* END PLUMED */
 
@@ -939,7 +940,7 @@ static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec&
             domainWork.haveSpecialForces || domainWork.haveCpuListedForceWork
             || domainWork.haveFreeEnergyWork || simulationWork.useCpuNonbonded || simulationWork.useCpuPme
             || simulationWork.haveEwaldSurfaceContribution || inputrec.nwall > 0
-            || plumedswitch; /* PLUMED: it adds its forces to the CPU force buffer */
+            || plumedIsUsed; /* PLUMED: it adds its forces to the CPU force buffer */
     domainWork.haveLocalForceContribInCpuBuffer =
             domainWork.haveCpuLocalForceWork || simulationWork.havePpDomainDecomposition;
     domainWork.haveNonLocalForceContribInCpuBuffer =
diff --git a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/legacymdrunoptions.cpp b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
index 024919022a..014dd57b63 100644
--- a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
+++ b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
@@ -185,6 +185,7 @@ int LegacyMdrunOptions::updateFromCommandLine(int argc, char** argv, ArrayRef<co
       //
       plumed_cmd(plumedmain,"setPlumedDat",ftp2fn(efDAT,static_cast<int>(filenames.size()), filenames.data()));
       plumedswitch=1;
+      plumedIsUsed=1;
     }
     /* PLUMED HREX*/
     if(getenv("PLUMED_HREX")) plumed_hrex=1;
diff --git a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/legacymdrunoptions.h b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/legacymdrunoptions.h
index 26a9b7d9c0..4668ccba91 100644
--- a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/legacymdrunoptions.h
+++ b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/legacymdrunoptions.h
@@ -57,6 +57,7 @@
 /* PLUMED */
 #include "../../../Plumed.h"
 extern int    plumedswitch;
+extern int    plumedIsUsed;
 extern plumed plumedmain; 
 /* END PLUMED */
 
diff --git a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/md.cpp b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/md.cpp
index 73d102ea32..0b742a8d69 100644
--- a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/md.cpp
+++ b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/md.cpp
@@ -823,8 +823,8 @@ void gmx::LegacySimulator::do_md()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
@@ -870,7 +870,7 @@ void gmx::LegacySimulator::do_md()
                 });
       }
       if(vsite){
-        /* The forces that PLUMED applies on virtual sites in the energy path are spread below */
+        /* When PLUMED runs after do_force(), the forces it applies on virtual sites are spread below */
         plumedIsVirtualSite.assign(top_global.natoms,false);
         for(const AtomProxy atomP : AtomRange(top_global)){
           plumedIsVirtualSite[atomP.globalAtomNumber()]=(atomP.atom().ptype==ParticleType::VSite);
@@ -1530,10 +1530,13 @@ void gmx::LegacySimulator::do_md()
             if(plumedswitch){
               if(plumedAfterForces){
                 ArrayRef<RVec> forceForPlumed=f.view().force();
+                std::vector<RVec> forceBeforePlumed;
                 if(simulationWork.useMts && step % ir->mtsLevels[1].stepFactor == 0){
-                  /* On the slow steps the update uses the combined MTS force buffer */
+                  /* On the slow steps the update uses the combined MTS force buffer, and the constraint virial
+                     the normal one: PLUMED's forces are added to the latter below */
                   forceForPlumed=f.view().forceMtsCombined();
                   plumed_cmd(plumedmain,"setForces",&forceForPlumed[0][0]);
+                  forceBeforePlumed.assign(forceForPlumed.begin(),forceForPlumed.begin()+md->homenr);
                 }
                 if(plumedNeedsEnergy){
                   msmul(force_vir,2.0,plumed_vir);
@@ -1554,6 +1557,10 @@ void gmx::LegacySimulator::do_md()
                                       VirtualSitesHandler::VirialHandling::NonLinear,{},vsiteVir,nrnb,state->box,wcycle);
                   m_add(force_vir,vsiteVir,force_vir);
                 }
+                if(!forceBeforePlumed.empty()){
+                  ArrayRef<RVec> forceNormal=f.view().force();
+                  for(int i=0;i<md->homenr;i++) forceNormal[i]+=forceForPlumed[i]-forceBeforePlumed[i];
+                }
               } else {
                 msmul(plumed_vir,0.5,plumed_vir);
                 m_add(force_vir,plumed_vir,force_vir);
diff --git a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/rerun.cpp b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/rerun.cpp
index 6f4e2e1fd2..6e8ad48b2b 100644
--- a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/rerun.cpp
+++ b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/rerun.cpp
@@ -444,8 +444,8 @@ void gmx::LegacySimulator::do_rerun()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
diff --git a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/runner.cpp b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/runner.cpp
index f0e8228936..8cebdeec35 100644
--- a/patches/gromacs-2022.5.diff/src/gromacs/mdrun/runner.cpp
+++ b/patches/gromacs-2022.5.diff/src/gromacs/mdrun/runner.cpp
@@ -178,6 +178,7 @@
 /* PLUMED */
 #include "../../../Plumed.h"
 int    plumedswitch;
+int    plumedIsUsed=0; /* unlike plumedswitch, not switched off while PLUMED runs after do_force() */
 plumed plumedmain; 
 /* END PLUMED */
 
@@ -2133,7 +2134,7 @@ int Mdrunner::mdrunner()
           if(useModularSimulator) gmx_fatal(FARGS, "PLUMED is not yet compatible with GROMACS new modular simulator");
 #if GMX_THREAD_MPI
           /* all the thread-MPI ranks would share the same PLUMED object */
-          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: use GROMACS built with MPI");
+          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: run with -ntmpi 1, or use GROMACS built with MPI");
 #endif
           /* detect plumed API version */
           int pversion=0;
diff --git a/patches/gromacs-2023.5.diff/src/gromacs/mdlib/sim_util.cpp b/patches/gromacs-2023.5.diff/src/gromacs/mdlib/sim_util.cpp
index 9aec5daf46..35f252fdd6 100644
--- a/patches/gromacs-2023.5.diff/src/gromacs/mdlib/sim_util.cpp
+++ b/patches/gromacs-2023.5.diff/src/gromacs/mdlib/sim_util.cpp
@@ -137,6 +137,7 @@ using gmx::StepWorkload;
 /* PLUMED */
 #include "../../../Plumed.h"
 extern int    plumedswitch;
+extern int    plumedIsUsed;
 extern plumed plumedmain;
 /* END PLUMED */
 
@@ -968,7 +969,7 @@ static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec&
             domainWork.haveSpecialForces || domainWork.haveCpuListedForceWork
             || domainWork.haveFreeEnergyWork || simulationWork.useCpuNonbonded || simulationWork.useCpuPme
             || simulationWork.haveEwaldSurfaceContribution || inputrec.nwall > 0
-            || plumedswitch; /* PLUMED: it adds its forces to the CPU force buffer */
+            || plumedIsUsed; /* PLUMED: it adds its forces to the CPU force buffer */
     domainWork.haveCpuNonLocalForceWork = domainWork.haveCpuBondedWork || domainWork.haveFreeEnergyWork;
     domainWork.haveLocalForceContribInCpuBuffer =
             domainWork.haveCpuLocalForceWork || simulationWork.havePpDomainDecomposition;
diff --git a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/legacymdrunoptions.cpp b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
index 024919022a..014dd57b63 100644
--- a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
+++ b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/legacymdrunoptions.cpp
@@ -185,6 +185,7 @@ int LegacyMdrunOptions::updateFromCommandLine(int argc, char** argv, ArrayRef<co
       //
       plumed_cmd(plumedmain,"setPlumedDat",ftp2fn(efDAT,static_cast<int>(filenames.size()), filenames.data()));
       plumedswitch=1;
+      plumedIsUsed=1;
     }
     /* PLUMED HREX*/
     if(getenv("PLUMED_HREX")) plumed_hrex=1;
diff --git a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/legacymdrunoptions.h b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/legacymdrunoptions.h
index cf894001b5..6cd6e3a993 100644
--- a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/legacymdrunoptions.h
+++ b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/legacymdrunoptions.h
@@ -57,6 +57,7 @@
 /* PLUMED */
 #include "../../../Plumed.h"
 extern int    plumedswitch;
+extern int    plumedIsUsed;
 extern plumed plumedmain; 
 /* END PLUMED */
 
diff --git a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/md.cpp b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/md.cpp
index 5ba9ee34e0..80ac6fcb09 100644
--- a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/md.cpp
+++ b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/md.cpp
@@ -813,8 +813,8 @@ void gmx::LegacySimulator::do_md()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
@@ -860,7 +860,7 @@ void gmx::LegacySimulator::do_md()
                 });
       }
       if(vsite){
-        /* The forces that PLUMED applies on virtual sites in the energy path are spread below */
+        /* When PLUMED runs after do_force(), the forces it applies on virtual sites are spread below */
         plumedIsVirtualSite.assign(top_global.natoms,false);
         for(const AtomProxy atomP : AtomRange(top_global)){
           plumedIsVirtualSite[atomP.globalAtomNumber()]=(atomP.atom().ptype==ParticleType::VSite);
@@ -1566,10 +1566,13 @@ void gmx::LegacySimulator::do_md()
             if(plumedswitch){
               if(plumedAfterForces){
                 ArrayRef<RVec> forceForPlumed=f.view().force();
+                std::vector<RVec> forceBeforePlumed;
                 if(simulationWork.useMts && step % ir->mtsLevels[1].stepFactor == 0){
-                  /* On the slow steps the update uses the combined MTS force buffer */
+                  /* On the slow steps the update uses the combined MTS force buffer, and the constraint virial
+                     the normal one: PLUMED's forces are added to the latter below */
                   forceForPlumed=f.view().forceMtsCombined();
                   plumed_cmd(plumedmain,"setForces",&forceForPlumed[0][0]);
+                  forceBeforePlumed.assign(forceForPlumed.begin(),forceForPlumed.begin()+md->homenr);
                 }
                 if(plumedNeedsEnergy){
                   msmul(force_vir,2.0,plumed_vir);
@@ -1590,6 +1593,10 @@ void gmx::LegacySimulator::do_md()
                                       VirtualSitesHandler::VirialHandling::NonLinear,{},vsiteVir,nrnb,state->box,wcycle);
                   m_add(force_vir,vsiteVir,force_vir);
                 }
+                if(!forceBeforePlumed.empty()){
+                  ArrayRef<RVec> forceNormal=f.view().force();
+                  for(int i=0;i<md->homenr;i++) forceNormal[i]+=forceForPlumed[i]-forceBeforePlumed[i];
+                }
               } else {
                 msmul(plumed_vir,0.5,plumed_vir);
                 m_add(force_vir,plumed_vir,force_vir);
diff --git a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/rerun.cpp b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/rerun.cpp
index db03341b4e..af848243fc 100644
--- a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/rerun.cpp
+++ b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/rerun.cpp
@@ -439,8 +439,8 @@ void gmx::LegacySimulator::do_rerun()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
diff --git a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/runner.cpp b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/runner.cpp
index 5342fc6868..757989a849 100644
--- a/patches/gromacs-2023.5.diff/src/gromacs/mdrun/runner.cpp
+++ b/patches/gromacs-2023.5.diff/src/gromacs/mdrun/runner.cpp
@@ -183,6 +183,7 @@
 /* PLUMED */
 #include "../../../Plumed.h"
 int    plumedswitch;
+int    plumedIsUsed=0; /* unlike plumedswitch, not switched off while PLUMED runs after do_force() */
 plumed plumedmain; 
 /* END PLUMED */
 
@@ -2247,7 +2248,7 @@ int Mdrunner::mdrunner()
           if(useModularSimulator) gmx_fatal(FARGS, "PLUMED is not yet compatible with GROMACS new modular simulator");
 #if GMX_THREAD_MPI
           /* all the thread-MPI ranks would share the same PLUMED object */
-          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: use GROMACS built with MPI");
+          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: run with -ntmpi 1, or use GROMACS built with MPI");
 #endif
           /* detect plumed API version */
           int pversion=0;
diff --git a/patches/gromacs-2024.3.diff/src/gromacs/mdrun/md.cpp b/patches/gromacs-2024.3.diff/src/gromacs/mdrun/md.cpp
index ee049c81f2..b803972755 100644
--- a/patches/gromacs-2024.3.diff/src/gromacs/mdrun/md.cpp
+++ b/patches/gromacs-2024.3.diff/src/gromacs/mdrun/md.cpp
@@ -818,8 +818,8 @@ void gmx::LegacySimulator::do_md()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
@@ -865,7 +865,7 @@ void gmx::LegacySimulator::do_md()
                 });
       }
       if(virtualSites_){
-        /* The forces that PLUMED applies on virtual sites in the energy path are spread below */
+        /* When PLUMED runs after do_force(), the forces it applies on virtual sites are spread below */
         plumedIsVirtualSite.assign(topGlobal_.natoms,false);
         for(const AtomProxy atomP : AtomRange(topGlobal_)){
           plumedIsVirtualSite[atomP.globalAtomNumber()]=(atomP.atom().ptype==ParticleType::VSite);
@@ -1617,10 +1617,13 @@ void gmx::LegacySimulator::do_md()
             if(plumedswitch){
               if(plumedAfterForces){
                 ArrayRef<RVec> forceForPlumed=f.view().force();
+                std::vector<RVec> forceBeforePlumed;
                 if(simulationWork.useMts && step % ir->mtsLevels[1].stepFactor == 0){
-                  /* On the slow steps the update uses the combined MTS force buffer */
+                  /* On the slow steps the update uses the combined MTS force buffer, and the constraint virial
+                     the normal one: PLUMED's forces are added to the latter below */
                   forceForPlumed=f.view().forceMtsCombined();
                   plumed_cmd(plumedmain,"setForces",&forceForPlumed[0][0]);
+                  forceBeforePlumed.assign(forceForPlumed.begin(),forceForPlumed.begin()+md->homenr);
                 }
                 if(plumedNeedsEnergy){
                   msmul(force_vir,2.0,plumed_vir);
@@ -1641,6 +1644,10 @@ void gmx::LegacySimulator::do_md()
                                               VirtualSitesHandler::VirialHandling::NonLinear,{},vsiteVir,nrnb_,state_->box,wallCycleCounters_);
                   m_add(force_vir,vsiteVir,force_vir);
                 }
+                if(!forceBeforePlumed.empty()){
+                  ArrayRef<RVec> forceNormal=f.view().force();
+                  for(int i=0;i<md->homenr;i++) forceNormal[i]+=forceForPlumed[i]-forceBeforePlumed[i];
+                }
               } else {
                 msmul(plumed_vir,0.5,plumed_vir);
                 m_add(force_vir,plumed_vir,force_vir);
diff --git a/patches/gromacs-2024.3.diff/src/gromacs/mdrun/rerun.cpp b/patches/gromacs-2024.3.diff/src/gromacs/mdrun/rerun.cpp
index ca09f7d512..0af63e39a3 100644
--- a/patches/gromacs-2024.3.diff/src/gromacs/mdrun/rerun.cpp
+++ b/patches/gromacs-2024.3.diff/src/gromacs/mdrun/rerun.cpp
@@ -443,8 +443,8 @@ void gmx::LegacySimulator::do_rerun()
     if(plumedswitch){
       if(shellfc){
         /* relax_shell_flexcon() calls do_force() without setting PLUMED up for the step */
-        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles with this version of GROMACS: "
-                        "PLUMED would not be applied.");
+        gmx_fatal(FARGS,"PLUMED cannot be used with shell (polarizable) particles or flexible constraints with this version "
+                        "of GROMACS: PLUMED would not be applied.");
       }
       /* detect plumed API version */
       int pversion=0;
diff --git a/patches/gromacs-2024.3.diff/src/gromacs/mdrun/runner.cpp b/patches/gromacs-2024.3.diff/src/gromacs/mdrun/runner.cpp
index 59e72da05f..70a1906f05 100644
--- a/patches/gromacs-2024.3.diff/src/gromacs/mdrun/runner.cpp
+++ b/patches/gromacs-2024.3.diff/src/gromacs/mdrun/runner.cpp
@@ -2330,7 +2330,7 @@ int Mdrunner::mdrunner()
           if(useModularSimulator) gmx_fatal(FARGS, "PLUMED is not yet compatible with GROMACS new modular simulator");
 #if GMX_THREAD_MPI
           /* all the thread-MPI ranks would share the same PLUMED object */
-          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: use GROMACS built with MPI");
+          if(PAR(cr)) gmx_fatal(FARGS, "PLUMED cannot be used with more than one thread-MPI rank: run with -ntmpi 1, or use GROMACS built with MPI");
 #endif
           /* detect plumed API version */
           int pversion=0;