Fix UB when generating local indices for constraints
[alexxy/gromacs.git] / src / gromacs / essentialdynamics / edsam.cpp
index d4c9de8e34599f0980ff97637a46a4d1f06979e9..48dafd3d0ab91ccb97649b438728ad917c2b9409 100644 (file)
@@ -904,14 +904,14 @@ static void update_adaption(t_edpar* edi)
 }
 
 
-static void do_single_flood(FILE*            edo,
-                            const rvec       x[],
-                            rvec             force[],
-                            t_edpar*         edi,
-                            int64_t          step,
-                            const matrix     box,
-                            const t_commrec* cr,
-                            gmx_bool         bNS) /* Are we in a neighbor searching step? */
+static void do_single_flood(FILE*                          edo,
+                            gmx::ArrayRef<const gmx::RVec> coords,
+                            gmx::ArrayRef<gmx::RVec>       force,
+                            t_edpar*                       edi,
+                            int64_t                        step,
+                            const matrix                   box,
+                            const t_commrec*               cr,
+                            gmx_bool bNS) /* Are we in a neighbor searching step? */
 {
     int                i;
     matrix             rotmat;   /* rotation matrix */
@@ -932,7 +932,7 @@ static void do_single_flood(FILE*            edo,
                                 buf->shifts_xcoll,
                                 buf->extra_shifts_xcoll,
                                 bNS,
-                                x,
+                                as_rvec_array(coords.data()),
                                 edi->sav.nr,
                                 edi->sav.nr_loc,
                                 edi->sav.anrs_loc,
@@ -948,7 +948,7 @@ static void do_single_flood(FILE*            edo,
                                     buf->shifts_xc_ref,
                                     buf->extra_shifts_xc_ref,
                                     bNS,
-                                    x,
+                                    as_rvec_array(coords.data()),
                                     edi->sref.nr,
                                     edi->sref.nr_loc,
                                     edi->sref.anrs_loc,
@@ -1028,14 +1028,14 @@ static void do_single_flood(FILE*            edo,
 
 
 /* Main flooding routine, called from do_force */
-void do_flood(const t_commrec*  cr,
-              const t_inputrec& ir,
-              const rvec        x[],
-              rvec              force[],
-              gmx_edsam*        ed,
-              const matrix      box,
-              int64_t           step,
-              bool              bNS)
+void do_flood(const t_commrec*               cr,
+              const t_inputrec&              ir,
+              gmx::ArrayRef<const gmx::RVec> coords,
+              gmx::ArrayRef<gmx::RVec>       force,
+              gmx_edsam*                     ed,
+              const matrix                   box,
+              int64_t                        step,
+              bool                           bNS)
 {
     /* Write time to edo, when required. Output the time anyhow since we need
      * it in the output file for ED constraints. */
@@ -1054,7 +1054,7 @@ void do_flood(const t_commrec*  cr,
         /* Call flooding for one matrix */
         if (edi.flood.vecs.neig)
         {
-            do_single_flood(ed->edo, x, force, &edi, step, box, cr, bNS);
+            do_single_flood(ed->edo, coords, force, &edi, step, box, cr, bNS);
         }
     }
 }
@@ -1158,8 +1158,8 @@ static std::unique_ptr<gmx::EssentialDynamics> ed_open(int
                                                        const gmx_output_env_t*     oenv,
                                                        const t_commrec*            cr)
 {
-    auto edHandle = std::make_unique<gmx::EssentialDynamics>();
-    auto ed       = edHandle->getLegacyED();
+    auto  edHandle = std::make_unique<gmx::EssentialDynamics>();
+    auto* ed       = edHandle->getLegacyED();
     /* We want to perform ED (this switch might later be upgraded to EssentialDynamicsType::Flooding) */
     ed->eEDtype = EssentialDynamicsType::EDSampling;
 
@@ -1335,7 +1335,7 @@ static void init_edi(const gmx_mtop_t& mtop, t_edpar* edi)
     {
         if (edi->fitmas)
         {
-            edi->sref.m[i] = mtopGetAtomMass(&mtop, edi->sref.anrs[i], &molb);
+            edi->sref.m[i] = mtopGetAtomMass(mtop, edi->sref.anrs[i], &molb);
         }
         else
         {
@@ -1366,7 +1366,7 @@ static void init_edi(const gmx_mtop_t& mtop, t_edpar* edi)
     snew(edi->sav.m, edi->sav.nr);
     for (i = 0; i < edi->sav.nr; i++)
     {
-        edi->sav.m[i] = mtopGetAtomMass(&mtop, edi->sav.anrs[i], &molb);
+        edi->sav.m[i] = mtopGetAtomMass(mtop, edi->sav.anrs[i], &molb);
         if (edi->pcamas)
         {
             edi->sav.sqrtm[i] = sqrt(edi->sav.m[i]);
@@ -1960,7 +1960,7 @@ static void translate_and_rotate(rvec*  x,        /* The positions to be transla
 
 /* Gets the rms deviation of the positions to the structure s */
 /* fit_to_structure has to be called before calling this routine! */
-static real rmsd_from_structure(rvec* x,           /* The positions under consideration */
+static real rmsd_from_structure(rvec*           x, /* The positions under consideration */
                                 struct gmx_edx* s) /* The structure from which the rmsd shall be computed */
 {
     real rmsd = 0.0;
@@ -2818,8 +2818,8 @@ std::unique_ptr<gmx::EssentialDynamics> init_edsam(const gmx::MDLogger&        m
                     "gmx grompp and the related .mdp options may change also.");
 
     /* Open input and output files, allocate space for ED data structure */
-    auto edHandle = ed_open(mtop.natoms, oh, ediFileName, edoFileName, startingBehavior, oenv, cr);
-    auto ed       = edHandle->getLegacyED();
+    auto  edHandle = ed_open(mtop.natoms, oh, ediFileName, edoFileName, startingBehavior, oenv, cr);
+    auto* ed       = edHandle->getLegacyED();
     GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
     constr->saveEdsamPointer(ed);
 
@@ -3062,14 +3062,15 @@ std::unique_ptr<gmx::EssentialDynamics> init_edsam(const gmx::MDLogger&        m
 
     } /* end of MASTER only section */
 
-    if (PAR(cr))
+    if (haveDDAtomOrdering(*cr))
     {
-        /* Broadcast the essential dynamics / flooding data to all nodes */
+        /* Broadcast the essential dynamics / flooding data to all nodes.
+         * In a single-rank case, only the necessary memory allocation is done. */
         broadcast_ed_data(cr, ed);
     }
     else
     {
-        /* In the single-CPU case, point the local atom numbers pointers to the global
+        /* In the non-DD case, point the local atom numbers pointers to the global
          * one, so that we can use the same notation in serial and parallel case: */
         /* Loop over all ED data sets (usually only one, though) */
         for (auto edi = ed->edpar.begin(); edi != ed->edpar.end(); ++edi)
@@ -3142,7 +3143,13 @@ std::unique_ptr<gmx::EssentialDynamics> init_edsam(const gmx::MDLogger&        m
 }
 
 
-void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[], rvec v[], const matrix box, gmx_edsam* ed)
+void do_edsam(const t_inputrec*        ir,
+              int64_t                  step,
+              const t_commrec*         cr,
+              gmx::ArrayRef<gmx::RVec> coords,
+              gmx::ArrayRef<gmx::RVec> velocities,
+              const matrix             box,
+              gmx_edsam*               ed)
 {
     int    i, edinr, iupdate = 500;
     matrix rotmat;         /* rotation matrix */
@@ -3191,7 +3198,7 @@ void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[]
                                         buf->shifts_xcoll,
                                         buf->extra_shifts_xcoll,
                                         PAR(cr) ? buf->bUpdateShifts : TRUE,
-                                        xs,
+                                        as_rvec_array(coords.data()),
                                         edi.sav.nr,
                                         edi.sav.nr_loc,
                                         edi.sav.anrs_loc,
@@ -3207,7 +3214,7 @@ void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[]
                                             buf->shifts_xc_ref,
                                             buf->extra_shifts_xc_ref,
                                             PAR(cr) ? buf->bUpdateShifts : TRUE,
-                                            xs,
+                                            as_rvec_array(coords.data()),
                                             edi.sref.nr,
                                             edi.sref.nr_loc,
                                             edi.sref.anrs_loc,
@@ -3313,17 +3320,17 @@ void do_edsam(const t_inputrec* ir, int64_t step, const t_commrec* cr, rvec xs[]
                             box, buf->xcoll[edi.sav.c_ind[i]], buf->shifts_xcoll[edi.sav.c_ind[i]], x_unsh);
 
                     /* dx is the ED correction to the positions: */
-                    rvec_sub(x_unsh, xs[edi.sav.anrs_loc[i]], dx);
+                    rvec_sub(x_unsh, coords[edi.sav.anrs_loc[i]], dx);
 
-                    if (v != nullptr)
+                    if (!velocities.empty())
                     {
                         /* dv is the ED correction to the velocity: */
                         svmul(dt_1, dx, dv);
                         /* apply the velocity correction: */
-                        rvec_inc(v[edi.sav.anrs_loc[i]], dv);
+                        rvec_inc(velocities[edi.sav.anrs_loc[i]], dv);
                     }
                     /* Finally apply the position correction due to ED: */
-                    copy_rvec(x_unsh, xs[edi.sav.anrs_loc[i]]);
+                    copy_rvec(x_unsh, coords[edi.sav.anrs_loc[i]]);
                 }
             }
         } /* END of if ( bNeedDoEdsam(edi) ) */