Merge branch release-2021 into merge-2021-into-master
[alexxy/gromacs.git] / src / gromacs / mdlib / stat.cpp
index 7c41c480c83ac85b39bd97b41e9e0d6406eee185..7bfdc6eebed13b327e9af3436c84b1f1320c24b1 100644 (file)
@@ -49,7 +49,6 @@
 #include "gromacs/gmxlib/network.h"
 #include "gromacs/math/utilities.h"
 #include "gromacs/math/vec.h"
-#include "gromacs/mdlib/constr.h"
 #include "gromacs/mdlib/md_support.h"
 #include "gromacs/mdlib/rbin.h"
 #include "gromacs/mdlib/tgroup.h"
@@ -141,46 +140,38 @@ static int filter_enerdterm(const real* afrom, gmx_bool bToBuffer, real* ato, gm
     return to;
 }
 
-void global_stat(const gmx_global_stat*  gs,
-                 const t_commrec*        cr,
-                 gmx_enerdata_t*         enerd,
-                 tensor                  fvir,
-                 tensor                  svir,
-                 const t_inputrec*       inputrec,
-                 gmx_ekindata_t*         ekind,
-                 const gmx::Constraints* constr,
-                 t_vcm*                  vcm,
-                 int                     nsig,
-                 real*                   sig,
-                 int*                    totalNumberOfBondedInteractions,
-                 bool                    bSumEkinhOld,
-                 int                     flags)
+void global_stat(const gmx_global_stat& gs,
+                 const t_commrec*       cr,
+                 gmx_enerdata_t*        enerd,
+                 tensor                 fvir,
+                 tensor                 svir,
+                 const t_inputrec&      inputrec,
+                 gmx_ekindata_t*        ekind,
+                 gmx::ArrayRef<real>    constraintsRmsdData,
+                 t_vcm*                 vcm,
+                 gmx::ArrayRef<real>    sig,
+                 bool                   bSumEkinhOld,
+                 int                    flags)
 /* instead of current system, gmx_booleans for summing virial, kinetic energy, and other terms */
 {
-    t_bin* rb;
-    int *  itc0, *itc1;
-    int    ie = 0, ifv = 0, isv = 0, irmsd = 0;
+    int ie = 0, ifv = 0, isv = 0, irmsd = 0;
     int idedl = 0, idedlo = 0, idvdll = 0, idvdlnl = 0, iepl = 0, icm = 0, imass = 0, ica = 0, inb = 0;
-    int      isig = -1;
-    int      icj = -1, ici = -1, icx = -1;
-    int      inn[egNR];
-    real     copyenerd[F_NRE];
-    int      nener, j;
-    double   nb;
-    gmx_bool bVV, bTemp, bEner, bPres, bConstrVir, bEkinAveVel, bReadEkin;
-    bool checkNumberOfBondedInteractions = (flags & CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS) != 0;
+    int isig = -1;
+    int icj = -1, ici = -1, icx = -1;
 
-    bVV         = EI_VV(inputrec->eI);
-    bTemp       = ((flags & CGLO_TEMPERATURE) != 0);
-    bEner       = ((flags & CGLO_ENERGY) != 0);
-    bPres       = ((flags & CGLO_PRESSURE) != 0);
-    bConstrVir  = ((flags & CGLO_CONSTRAINT) != 0);
-    bEkinAveVel = (inputrec->eI == eiVV || (inputrec->eI == eiVVAK && bPres));
-    bReadEkin   = ((flags & CGLO_READEKIN) != 0);
+    bool checkNumberOfBondedInteractions = (flags & CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS) != 0;
+    bool bVV                             = EI_VV(inputrec.eI);
+    bool bTemp                           = ((flags & CGLO_TEMPERATURE) != 0);
+    bool bEner                           = ((flags & CGLO_ENERGY) != 0);
+    bool bPres                           = ((flags & CGLO_PRESSURE) != 0);
+    bool bConstrVir                      = ((flags & CGLO_CONSTRAINT) != 0);
+    bool bEkinAveVel                     = (inputrec.eI == IntegrationAlgorithm::VV
+                        || (inputrec.eI == IntegrationAlgorithm::VVAK && bPres));
+    bool bReadEkin                       = ((flags & CGLO_READEKIN) != 0);
 
-    rb   = gs->rb;
-    itc0 = gs->itc0;
-    itc1 = gs->itc1;
+    t_bin* rb   = gs.rb;
+    int*   itc0 = gs.itc0;
+    int*   itc1 = gs.itc1;
 
 
     reset_bin(rb);
@@ -193,7 +184,8 @@ void global_stat(const gmx_global_stat*  gs,
        communicated and summed when they need to be, to avoid repeating
        the sums and overcounting. */
 
-    nener = filter_enerdterm(enerd->term, TRUE, copyenerd, bTemp, bPres, bEner);
+    std::array<real, F_NRE> copyenerd;
+    int nener = filter_enerdterm(enerd->term.data(), TRUE, copyenerd.data(), bTemp, bPres, bEner);
 
     /* First, the data that needs to be communicated with velocity verlet every time
        This is just the constraint virial.*/
@@ -207,7 +199,7 @@ void global_stat(const gmx_global_stat*  gs,
     {
         if (ekind)
         {
-            for (j = 0; (j < inputrec->opts.ngtc); j++)
+            for (int j = 0; (j < inputrec.opts.ngtc); j++)
             {
                 if (bSumEkinhOld)
                 {
@@ -240,30 +232,26 @@ void global_stat(const gmx_global_stat*  gs,
         ifv = add_binr(rb, DIM * DIM, fvir[0]);
     }
 
-    gmx::ArrayRef<real> rmsdData;
+    gmx::EnumerationArray<NonBondedEnergyTerms, int> inn;
     if (bEner)
     {
-        ie = add_binr(rb, nener, copyenerd);
-        if (constr)
+        ie = add_binr(rb, nener, copyenerd.data());
+        if (!constraintsRmsdData.empty())
         {
-            rmsdData = constr->rmsdData();
-            if (!rmsdData.empty())
-            {
-                irmsd = add_binr(rb, 2, rmsdData.data());
-            }
+            irmsd = add_binr(rb, 2, constraintsRmsdData.data());
         }
-
-        for (j = 0; (j < egNR); j++)
+        for (auto key : gmx::keysOf(inn))
         {
-            inn[j] = add_binr(rb, enerd->grpp.nener, enerd->grpp.ener[j].data());
+            inn[key] = add_binr(rb, enerd->grpp.nener, enerd->grpp.energyGroupPairTerms[key].data());
         }
-        if (inputrec->efep != efepNO)
+        if (inputrec.efep != FreeEnergyPerturbationType::No)
         {
-            idvdll  = add_bind(rb, efptNR, enerd->dvdl_lin);
-            idvdlnl = add_bind(rb, efptNR, enerd->dvdl_nonlin);
+            idvdll  = add_bind(rb, enerd->dvdl_lin);
+            idvdlnl = add_bind(rb, enerd->dvdl_nonlin);
             if (enerd->foreignLambdaTerms.numLambdas() > 0)
             {
-                iepl = add_bind(rb, enerd->foreignLambdaTerms.energies().size(),
+                iepl = add_bind(rb,
+                                enerd->foreignLambdaTerms.energies().size(),
                                 enerd->foreignLambdaTerms.energies().data());
             }
         }
@@ -273,7 +261,7 @@ void global_stat(const gmx_global_stat*  gs,
     {
         icm   = add_binr(rb, DIM * vcm->nr, vcm->group_p[0]);
         imass = add_binr(rb, vcm->nr, vcm->group_mass.data());
-        if (vcm->mode == ecmANGULAR)
+        if (vcm->mode == ComRemovalAlgorithm::Angular)
         {
             icj = add_binr(rb, DIM * vcm->nr, vcm->group_j[0]);
             icx = add_binr(rb, DIM * vcm->nr, vcm->group_x[0]);
@@ -281,21 +269,20 @@ void global_stat(const gmx_global_stat*  gs,
         }
     }
 
+    double nb;
     if (checkNumberOfBondedInteractions)
     {
-        nb  = cr->dd->nbonded_local;
+        GMX_RELEASE_ASSERT(DOMAINDECOMP(cr),
+                           "No need to check number of bonded interactions when not using domain "
+                           "decomposition");
+        nb  = numBondedInteractions(*cr->dd);
         inb = add_bind(rb, 1, &nb);
     }
-    if (nsig > 0)
+    if (!sig.empty())
     {
-        isig = add_binr(rb, nsig, sig);
+        isig = add_binr(rb, sig);
     }
 
-    /* Global sum it all */
-    if (debug)
-    {
-        fprintf(debug, "Summing %d energies\n", rb->maxreal);
-    }
     sum_bin(rb, cr);
 
     /* Extract all the data locally */
@@ -310,7 +297,7 @@ void global_stat(const gmx_global_stat*  gs,
     {
         if (ekind)
         {
-            for (j = 0; (j < inputrec->opts.ngtc); j++)
+            for (int j = 0; (j < inputrec.opts.ngtc); j++)
             {
                 if (bSumEkinhOld)
                 {
@@ -343,35 +330,36 @@ void global_stat(const gmx_global_stat*  gs,
 
     if (bEner)
     {
-        extract_binr(rb, ie, nener, copyenerd);
-        if (!rmsdData.empty())
+        extract_binr(rb, ie, nener, copyenerd.data());
+        if (!constraintsRmsdData.empty())
         {
-            extract_binr(rb, irmsd, rmsdData);
+            extract_binr(rb, irmsd, constraintsRmsdData);
         }
-
-        for (j = 0; (j < egNR); j++)
+        for (auto key : gmx::keysOf(inn))
         {
-            extract_binr(rb, inn[j], enerd->grpp.nener, enerd->grpp.ener[j].data());
+            extract_binr(rb, inn[key], enerd->grpp.nener, enerd->grpp.energyGroupPairTerms[key].data());
         }
-        if (inputrec->efep != efepNO)
+        if (inputrec.efep != FreeEnergyPerturbationType::No)
         {
-            extract_bind(rb, idvdll, efptNR, enerd->dvdl_lin);
-            extract_bind(rb, idvdlnl, efptNR, enerd->dvdl_nonlin);
+            extract_bind(rb, idvdll, enerd->dvdl_lin);
+            extract_bind(rb, idvdlnl, enerd->dvdl_nonlin);
             if (enerd->foreignLambdaTerms.numLambdas() > 0)
             {
-                extract_bind(rb, iepl, enerd->foreignLambdaTerms.energies().size(),
+                extract_bind(rb,
+                             iepl,
+                             enerd->foreignLambdaTerms.energies().size(),
                              enerd->foreignLambdaTerms.energies().data());
             }
         }
 
-        filter_enerdterm(copyenerd, FALSE, enerd->term, bTemp, bPres, bEner);
+        filter_enerdterm(copyenerd.data(), FALSE, enerd->term.data(), bTemp, bPres, bEner);
     }
 
     if (vcm)
     {
         extract_binr(rb, icm, DIM * vcm->nr, vcm->group_p[0]);
         extract_binr(rb, imass, vcm->nr, vcm->group_mass.data());
-        if (vcm->mode == ecmANGULAR)
+        if (vcm->mode == ComRemovalAlgorithm::Angular)
         {
             extract_binr(rb, icj, DIM * vcm->nr, vcm->group_j[0]);
             extract_binr(rb, icx, DIM * vcm->nr, vcm->group_x[0]);
@@ -382,11 +370,14 @@ void global_stat(const gmx_global_stat*  gs,
     if (checkNumberOfBondedInteractions)
     {
         extract_bind(rb, inb, 1, &nb);
-        *totalNumberOfBondedInteractions = gmx::roundToInt(nb);
+        GMX_RELEASE_ASSERT(DOMAINDECOMP(cr),
+                           "No need to check number of bonded interactions when not using domain "
+                           "decomposition");
+        setNumberOfBondedInteractionsOverAllDomains(*cr->dd, gmx::roundToInt(nb));
     }
 
-    if (nsig > 0)
+    if (!sig.empty())
     {
-        extract_binr(rb, isig, nsig, sig);
+        extract_binr(rb, isig, sig);
     }
 }