Reimplement constant acceleration groups

[alexxy/gromacs.git] / src / gromacs / mdrun / shellfc.cpp

diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp
index 73c14ac3b5453c08d824dd5b0b2a1f83a69e2505..7657596888178a4f4b6e4cb6ebaab1a9b8cd84af 100644 (file)
--- a/src/gromacs/mdrun/shellfc.cpp
+++ b/src/gromacs/mdrun/shellfc.cpp
@@ -62,7 +62,6 @@
 #include "gromacs/mdlib/force.h"
 #include "gromacs/mdlib/force_flags.h"
 #include "gromacs/mdlib/gmx_omp_nthreads.h"
-#include "gromacs/mdlib/mdatoms.h"
 #include "gromacs/mdlib/vsite.h"
 #include "gromacs/mdtypes/commrec.h"
 #include "gromacs/mdtypes/enerdata.h"
@@ -288,7 +287,7 @@ gmx_shellfc_t* init_shell_flexcon(FILE*             fplog,
                   "not supported in combination with shell particles.\nPlease make a new tpr file.",
                   nstcalcenergy);
     }
-    if (usingDomainDecomposition)
+    if (nshell > 0 && usingDomainDecomposition)
     {
         gmx_fatal(
                 FARGS,
@@ -536,12 +535,12 @@ gmx_shellfc_t* init_shell_flexcon(FILE*             fplog,
     return shfc;
 }
 
-void gmx::make_local_shells(const t_commrec* cr, const t_mdatoms* md, gmx_shellfc_t* shfc)
+void gmx::make_local_shells(const t_commrec* cr, const t_mdatoms& md, gmx_shellfc_t* shfc)
 {
     int           a0, a1;
     gmx_domdec_t* dd = nullptr;
 
-    if (DOMAINDECOMP(cr))
+    if (haveDDAtomOrdering(*cr))
     {
         dd = cr->dd;
         a0 = 0;
@@ -559,9 +558,10 @@ void gmx::make_local_shells(const t_commrec* cr, const t_mdatoms* md, gmx_shellf
 
     std::vector<t_shell>& shells = shfc->shells;
     shells.clear();
+    auto* ptype = md.ptype;
     for (int i = a0; i < a1; i++)
     {
-        if (md->ptype[i] == ParticleType::Shell)
+        if (ptype[i] == ParticleType::Shell)
         {
             if (dd)
             {
@@ -819,7 +819,7 @@ static void init_adir(gmx_shellfc_t*            shfc,
                       const t_commrec*          cr,
                       int                       dd_ac1,
                       int64_t                   step,
-                      const t_mdatoms*          md,
+                      const t_mdatoms&          md,
                       int                       end,
                       ArrayRefWithPadding<RVec> xOld,
                       ArrayRef<RVec>            x_init,
@@ -833,7 +833,7 @@ static void init_adir(gmx_shellfc_t*            shfc,
     double dt, w_dt;
     int    n, d;
 
-    if (DOMAINDECOMP(cr))
+    if (haveDDAtomOrdering(*cr))
     {
         n = dd_ac1;
     }
@@ -848,14 +848,14 @@ static void init_adir(gmx_shellfc_t*            shfc,
     rvec* x_old = as_rvec_array(xOld.paddedArrayRef().data());
     rvec* x     = as_rvec_array(xCurrent.paddedArrayRef().data());
 
-    const ParticleType* ptype = md->ptype;
+    auto* ptype   = md.ptype;
+    auto  invmass = gmx::arrayRefFromArray(md.invmass, md.nr);
+    dt            = ir->delta_t;
 
-    dt = ir->delta_t;
-
-    /* Does NOT work with freeze groups (yet) */
+    /* Does NOT work with freeze or acceleration groups (yet) */
     for (n = 0; n < end; n++)
     {
-        w_dt = md->invmass[n] * dt;
+        w_dt = invmass[n] * dt;
 
         for (d = 0; d < DIM; d++)
         {
@@ -910,7 +910,7 @@ static void init_adir(gmx_shellfc_t*            shfc,
         for (d = 0; d < DIM; d++)
         {
             xnew[n][d] = -(2 * x[n][d] - xnold[n][d] - xnew[n][d]) / gmx::square(dt)
-                         - f[n][d] * md->invmass[n];
+                         - f[n][d] * invmass[n];
         }
         clear_rvec(acc_dir[n]);
     }
@@ -955,9 +955,10 @@ void relax_shell_flexcon(FILE*                         fplog,
                          const history_t*              hist,
                          gmx::ForceBuffersView*        f,
                          tensor                        force_vir,
-                         const t_mdatoms*              md,
+                         const t_mdatoms&              md,
+                         CpuPpLongRangeNonbondeds*     longRangeNonbondeds,
                          t_nrnb*                       nrnb,
-                         gmx_wallcycle_t               wcycle,
+                         gmx_wallcycle*                wcycle,
                          gmx_shellfc_t*                shfc,
                          t_forcerec*                   fr,
                          gmx::MdrunScheduleWorkload*   runScheduleWork,
@@ -971,7 +972,7 @@ void relax_shell_flexcon(FILE*                         fplog,
     real dum = 0;
     char sbuf[22];
     int  nat, dd_ac0, dd_ac1 = 0, i;
-    int  homenr = md->homenr, end = homenr;
+    int  homenr = md.homenr, end = homenr;
     int  d, Min = 0, count = 0;
 #define Try (1 - Min) /* At start Try = 1 */
 
@@ -982,7 +983,7 @@ void relax_shell_flexcon(FILE*                         fplog,
     ArrayRef<t_shell> shells       = shfc->shells;
     const int         nflexcon     = shfc->nflexcon;
 
-    if (DOMAINDECOMP(cr))
+    if (haveDDAtomOrdering(*cr))
     {
         nat = dd_natoms_vsite(*cr->dd);
         if (nflexcon > 0)
@@ -1018,14 +1019,14 @@ void relax_shell_flexcon(FILE*                         fplog,
     ArrayRef<RVec> x = xPadded.unpaddedArrayRef();
     ArrayRef<RVec> v = vPadded.unpaddedArrayRef();
 
-    if (bDoNS && inputrec->pbcType != PbcType::No && !DOMAINDECOMP(cr))
+    if (bDoNS && inputrec->pbcType != PbcType::No && !haveDDAtomOrdering(*cr))
     {
         /* This is the only time where the coordinates are used
          * before do_force is called, which normally puts all
          * charge groups in the box.
          */
         put_atoms_in_box_omp(
-                fr->pbcType, box, x.subArray(0, md->homenr), gmx_omp_nthreads_get(emntDefault));
+                fr->pbcType, box, x.subArray(0, md.homenr), gmx_omp_nthreads_get(ModuleMultiThread::Default));
     }
 
     if (nflexcon)
@@ -1042,13 +1043,13 @@ void relax_shell_flexcon(FILE*                         fplog,
         }
     }
 
+    auto massT = gmx::arrayRefFromArray(md.massT, md.nr);
     /* Do a prediction of the shell positions, when appropriate.
      * Without velocities (EM, NM, BD) we only do initial prediction.
      */
     if (shfc->predictShells && !bCont && (EI_STATE_VELOCITY(inputrec->eI) || bInit))
     {
-        predict_shells(
-                fplog, x, v, inputrec->delta_t, shells, gmx::arrayRefFromArray(md->massT, md->homenr), bInit);
+        predict_shells(fplog, x, v, inputrec->delta_t, shells, massT, bInit);
     }
 
     /* Calculate the forces first time around */
@@ -1075,7 +1076,7 @@ void relax_shell_flexcon(FILE*                         fplog,
              hist,
              &forceViewInit,
              force_vir,
-             md,
+             &md,
              enerd,
              lambda,
              fr,
@@ -1084,6 +1085,7 @@ void relax_shell_flexcon(FILE*                         fplog,
              mu_tot,
              t,
              nullptr,
+             longRangeNonbondeds,
              (bDoNS ? GMX_FORCE_NS : 0) | shellfc_flags,
              ddBalanceRegionHandler);
 
@@ -1109,7 +1111,7 @@ void relax_shell_flexcon(FILE*                         fplog,
 
         for (i = 0; i < end; i++)
         {
-            sf_dir += md->massT[i] * norm2(shfc->acc_dir[i]);
+            sf_dir += massT[i] * norm2(shfc->acc_dir[i]);
         }
     }
     accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals.get());
@@ -1124,7 +1126,7 @@ void relax_shell_flexcon(FILE*                         fplog,
 
     if (gmx_debug_at)
     {
-        pr_rvecs(debug, 0, "force0", as_rvec_array(force[Min].data()), md->nr);
+        pr_rvecs(debug, 0, "force0", as_rvec_array(force[Min].data()), md.nr);
     }
 
     if (!shells.empty() || nflexcon > 0)
@@ -1215,7 +1217,7 @@ void relax_shell_flexcon(FILE*                         fplog,
                  hist,
                  &forceViewTry,
                  force_vir,
-                 md,
+                 &md,
                  enerd,
                  lambda,
                  fr,
@@ -1224,6 +1226,7 @@ void relax_shell_flexcon(FILE*                         fplog,
                  mu_tot,
                  t,
                  nullptr,
+                 longRangeNonbondeds,
                  shellfc_flags,
                  ddBalanceRegionHandler);
         accumulatePotentialEnergies(enerd, lambda, inputrec->fepvals.get());
@@ -1255,7 +1258,7 @@ void relax_shell_flexcon(FILE*                         fplog,
             ArrayRef<const RVec> acc_dir = shfc->acc_dir;
             for (i = 0; i < end; i++)
             {
-                sf_dir += md->massT[i] * norm2(acc_dir[i]);
+                sf_dir += massT[i] * norm2(acc_dir[i]);
             }
         }