GROMACS 2020-2024 patches: follow-up of #1505 (GPU buffer operations with 2020-2023, constraint virial with MTS, messages and CHANGES)

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

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 */
 

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;
 }

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;

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 */
 

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);

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;

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;

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;
 }

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;

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 */
 

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);

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;

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;

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 =

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;

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 */
 

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);

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;

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;

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;

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;

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 */