Fix segfault in pull reading
authorJoe Jordan <ejjordan12@gmail.com>
Tue, 13 Oct 2020 12:35:33 +0000 (12:35 +0000)
committerPaul Bauer <paul.bauer.q@gmail.com>
Tue, 13 Oct 2020 12:35:33 +0000 (12:35 +0000)
Triggers asserton now later on :(

Change-Id: I08a390de653e8bbdf23288ed8a28b6dc7ad6305d

21 files changed:
src/gromacs/applied_forces/awh/read_params.cpp
src/gromacs/applied_forces/awh/read_params.h
src/gromacs/fileio/tpxio.cpp
src/gromacs/gmxpreprocess/grompp.cpp
src/gromacs/gmxpreprocess/readir.cpp
src/gromacs/gmxpreprocess/readir.h
src/gromacs/gmxpreprocess/readpull.cpp
src/gromacs/mdlib/constr.cpp
src/gromacs/mdlib/energyoutput.cpp
src/gromacs/mdlib/makeconstraints.h
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdrun/md.cpp
src/gromacs/mdrun/runner.cpp
src/gromacs/mdtypes/inputrec.cpp
src/gromacs/mdtypes/inputrec.h
src/gromacs/mdtypes/pull_params.h
src/gromacs/pulling/pull.cpp
src/gromacs/pulling/pull.h
src/gromacs/pulling/pullutil.cpp
src/gromacs/taskassignment/decidegpuusage.cpp
src/gromacs/tools/check.cpp

index e5945f17d6b06bc7de1bb70dab7bb5d12443a908..5654d57af433467a7c01085b4b1dd099b0fffdd6 100644 (file)
@@ -138,7 +138,7 @@ void checkMtsConsistency(const t_inputrec& inputrec, warninp_t wi)
  */
 void checkPullDimParams(const std::string&   prefix,
                         AwhDimParams*        dimParams,
-                        const pull_params_t* pull_params,
+                        const pull_params_t& pull_params,
                         warninp_t            wi)
 {
     if (dimParams->coordIndex < 0)
@@ -148,18 +148,18 @@ void checkPullDimParams(const std::string&   prefix,
                   "Note that the pull coordinate indexing starts at 1.",
                   prefix.c_str());
     }
-    if (dimParams->coordIndex >= pull_params->ncoord)
+    if (dimParams->coordIndex >= pull_params.ncoord)
     {
         gmx_fatal(FARGS,
                   "The given AWH coordinate index (%d) is larger than the number of pull "
                   "coordinates (%d)",
-                  dimParams->coordIndex + 1, pull_params->ncoord);
+                  dimParams->coordIndex + 1, pull_params.ncoord);
     }
-    if (pull_params->coord[dimParams->coordIndex].rate != 0)
+    if (pull_params.coord[dimParams->coordIndex].rate != 0)
     {
         auto message = formatString(
                 "Setting pull-coord%d-rate (%g) is incompatible with AWH biasing this coordinate",
-                dimParams->coordIndex + 1, pull_params->coord[dimParams->coordIndex].rate);
+                dimParams->coordIndex + 1, pull_params.coord[dimParams->coordIndex].rate);
         warning_error(wi, message);
     }
 
@@ -178,7 +178,7 @@ void checkPullDimParams(const std::string&   prefix,
     }
 
     /* Grid params for each axis */
-    int eGeom = pull_params->coord[dimParams->coordIndex].eGeom;
+    int eGeom = pull_params.coord[dimParams->coordIndex].eGeom;
 
     /* Check that the requested interval is in allowed range */
     if (eGeom == epullgDIST)
@@ -412,7 +412,7 @@ void checkDimParams(const std::string& prefix, AwhDimParams* dimParams, const t_
                       "AWH biasing along a pull dimension is only compatible with COM pulling "
                       "turned on");
         }
-        checkPullDimParams(prefix, dimParams, ir->pull, wi);
+        checkPullDimParams(prefix, dimParams, *ir->pull, wi);
     }
     else if (dimParams->eCoordProvider == eawhcoordproviderFREE_ENERGY_LAMBDA)
     {
@@ -1068,7 +1068,7 @@ static void setStateDependentAwhPullDimParams(AwhDimParams*        dimParams,
 }
 
 void setStateDependentAwhParams(AwhParams*           awhParams,
-                                const pull_params_t* pull_params,
+                                const pull_params_t& pull_params,
                                 pull_t*              pull_work,
                                 const matrix         box,
                                 PbcType              pbcType,
@@ -1107,7 +1107,7 @@ void setStateDependentAwhParams(AwhParams*           awhParams,
             AwhDimParams* dimParams = &awhBiasParams->dimParams[d];
             if (dimParams->eCoordProvider == eawhcoordproviderPULL)
             {
-                setStateDependentAwhPullDimParams(dimParams, k, d, pull_params, pull_work, pbc,
+                setStateDependentAwhPullDimParams(dimParams, k, d, &pull_params, pull_work, pbc,
                                                   compressibility, wi);
             }
             else
index 7da77f1fdd49758cd9a96684933666ba0a269b91..bb1cd63fa9cdd5289d33125ff266c1ad16cda8bd 100644 (file)
@@ -96,7 +96,7 @@ void checkAwhParams(const AwhParams* awhParams, const t_inputrec* inputrec, warn
  * \note This function currently relies on the function set_pull_init to have been called.
  */
 void setStateDependentAwhParams(AwhParams*           awhParams,
-                                const pull_params_t* pull_params,
+                                const pull_params_t& pull_params,
                                 pull_t*              pull_work,
                                 const matrix         box,
                                 PbcType              pbcType,
index 9ebd98f4b3e29c87734289ee10b2147f68122351..9ea2337f86727f243dca4cff6997fe0b5e13d248 100644 (file)
@@ -241,20 +241,16 @@ static void do_pullgrp_tpx_pre95(gmx::ISerializer* serializer, t_pull_group* pgr
 {
     rvec tmp;
 
-    serializer->doInt(&pgrp->nat);
-    if (serializer->reading())
-    {
-        snew(pgrp->ind, pgrp->nat);
-    }
-    serializer->doIntArray(pgrp->ind, pgrp->nat);
-    serializer->doInt(&pgrp->nweight);
-    if (serializer->reading())
-    {
-        snew(pgrp->weight, pgrp->nweight);
-    }
-    serializer->doRealArray(pgrp->weight, pgrp->nweight);
+    int numAtoms = pgrp->ind.size();
+    serializer->doInt(&numAtoms);
+    pgrp->ind.resize(numAtoms);
+    serializer->doIntArray(pgrp->ind.data(), numAtoms);
+    int numWeights = pgrp->weight.size();
+    serializer->doInt(&numWeights);
+    pgrp->weight.resize(numWeights);
+    serializer->doRealArray(pgrp->weight.data(), numWeights);
     serializer->doInt(&pgrp->pbcatom);
-    serializer->doRvec(&pcrd->vec);
+    serializer->doRvec(&pcrd->vec.as_vec());
     clear_rvec(pcrd->origin);
     serializer->doRvec(&tmp);
     pcrd->init = tmp[0];
@@ -265,18 +261,14 @@ static void do_pullgrp_tpx_pre95(gmx::ISerializer* serializer, t_pull_group* pgr
 
 static void do_pull_group(gmx::ISerializer* serializer, t_pull_group* pgrp)
 {
-    serializer->doInt(&pgrp->nat);
-    if (serializer->reading())
-    {
-        snew(pgrp->ind, pgrp->nat);
-    }
-    serializer->doIntArray(pgrp->ind, pgrp->nat);
-    serializer->doInt(&pgrp->nweight);
-    if (serializer->reading())
-    {
-        snew(pgrp->weight, pgrp->nweight);
-    }
-    serializer->doRealArray(pgrp->weight, pgrp->nweight);
+    int numAtoms = pgrp->ind.size();
+    serializer->doInt(&numAtoms);
+    pgrp->ind.resize(numAtoms);
+    serializer->doIntArray(pgrp->ind.data(), numAtoms);
+    int numWeights = pgrp->weight.size();
+    serializer->doInt(&numWeights);
+    pgrp->weight.resize(numWeights);
+    serializer->doRealArray(pgrp->weight.data(), numWeights);
     serializer->doInt(&pgrp->pbcatom);
 }
 
@@ -308,14 +300,14 @@ static void do_pull_coord(gmx::ISerializer* serializer,
             }
             else
             {
-                pcrd->externalPotentialProvider = nullptr;
+                pcrd->externalPotentialProvider.clear();
             }
         }
         else
         {
             if (serializer->reading())
             {
-                pcrd->externalPotentialProvider = nullptr;
+                pcrd->externalPotentialProvider.clear();
             }
         }
         /* Note that we try to support adding new geometries without
@@ -326,7 +318,7 @@ static void do_pull_coord(gmx::ISerializer* serializer,
         serializer->doInt(&pcrd->ngroup);
         if (pcrd->ngroup <= c_pullCoordNgroupMax)
         {
-            serializer->doIntArray(pcrd->group, pcrd->ngroup);
+            serializer->doIntArray(pcrd->group.data(), pcrd->ngroup);
         }
         else
         {
@@ -342,7 +334,7 @@ static void do_pull_coord(gmx::ISerializer* serializer,
 
             pcrd->ngroup = 0;
         }
-        serializer->doIvec(&pcrd->dim);
+        serializer->doIvec(&pcrd->dim.as_vec());
     }
     else
     {
@@ -359,7 +351,7 @@ static void do_pull_coord(gmx::ISerializer* serializer,
                 serializer->doInt(&pcrd->group[2]);
                 serializer->doInt(&pcrd->group[3]);
             }
-            serializer->doIvec(&pcrd->dim);
+            serializer->doIvec(&pcrd->dim.as_vec());
         }
         else
         {
@@ -368,8 +360,8 @@ static void do_pull_coord(gmx::ISerializer* serializer,
             copy_ivec(dimOld, pcrd->dim);
         }
     }
-    serializer->doRvec(&pcrd->origin);
-    serializer->doRvec(&pcrd->vec);
+    serializer->doRvec(&pcrd->origin.as_vec());
+    serializer->doRvec(&pcrd->vec.as_vec());
     if (file_version >= tpxv_PullCoordTypeGeom)
     {
         serializer->doBool(&pcrd->bStart);
@@ -775,11 +767,8 @@ static void do_pull(gmx::ISerializer* serializer, pull_params_t* pull, int file_
     {
         pull->bSetPbcRefToPrevStepCOM = FALSE;
     }
-    if (serializer->reading())
-    {
-        snew(pull->group, pull->ngroup);
-        snew(pull->coord, pull->ncoord);
-    }
+    pull->group.resize(pull->ngroup);
+    pull->coord.resize(pull->ncoord);
     if (file_version < 95)
     {
         /* epullgPOS for position pulling, before epullgDIRPBC was removed */
@@ -803,7 +792,7 @@ static void do_pull(gmx::ISerializer* serializer, pull_params_t* pull, int file_
             }
         }
 
-        pull->bPrintCOM = (pull->group[0].nat > 0);
+        pull->bPrintCOM = (!pull->group[0].ind.empty());
     }
     else
     {
@@ -1488,9 +1477,9 @@ static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_v
         {
             if (serializer->reading())
             {
-                snew(ir->pull, 1);
+                ir->pull = std::make_unique<pull_params_t>();
             }
-            do_pull(serializer, ir->pull, file_version, ePullOld);
+            do_pull(serializer, ir->pull.get(), file_version, ePullOld);
         }
     }
 
index 6a20af125b0a1626c976586d1df61cec36014e69..46a48b93e434b724e0038fa7fe0ba875430a501a 100644 (file)
@@ -2304,7 +2304,7 @@ int gmx_grompp(int argc, char* argv[])
             copy_mat(ir->compress, compressibility);
         }
         setStateDependentAwhParams(
-                ir->awhParams, ir->pull, pull, state.box, ir->pbcType, compressibility, &ir->opts,
+                ir->awhParams, *ir->pull, pull, state.box, ir->pbcType, compressibility, &ir->opts,
                 ir->efep != efepNO ? ir->fepvals->all_lambda[efptFEP][ir->fepvals->init_fep_state] : 0,
                 sys, wi);
     }
index b2d504ca235afb9898c598bc6ee2da9fe8e81bf5..2edc32384d8d18298892a67c7b8fd88b918b74aa 100644 (file)
@@ -45,6 +45,7 @@
 #include <cstdlib>
 
 #include <algorithm>
+#include <memory>
 #include <string>
 
 #include "gromacs/applied_forces/awh/read_params.h"
@@ -2209,8 +2210,8 @@ void get_ir(const char*     mdparin,
     ir->bPull = (get_eeenum(&inp, "pull", yesno_names, wi) != 0);
     if (ir->bPull)
     {
-        snew(ir->pull, 1);
-        inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull, wi);
+        ir->pull                        = std::make_unique<pull_params_t>();
+        inputrecStrings->pullGroupNames = read_pullparams(&inp, ir->pull.get(), wi);
 
         if (ir->useMts)
         {
@@ -3072,7 +3073,7 @@ static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
          * belong to different TC or VCM groups it is anyhow difficult
          * to determine the optimal nrdf assignment.
          */
-        pull = ir->pull;
+        pull = ir->pull.get();
 
         for (int i = 0; i < pull->ncoord; i++)
         {
@@ -3089,7 +3090,7 @@ static void calc_nrdf(const gmx_mtop_t* mtop, t_inputrec* ir, char** gnames)
 
                 pgrp = &pull->group[pull->coord[i].group[j]];
 
-                if (pgrp->nat > 0)
+                if (!pgrp->ind.empty())
                 {
                     /* Subtract 1/2 dof from each group */
                     int ai = pgrp->ind[0];
@@ -3756,9 +3757,9 @@ void do_index(const char*                   mdparin,
 
     if (ir->bPull)
     {
-        make_pull_groups(ir->pull, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
+        process_pull_groups(ir->pull->group, inputrecStrings->pullGroupNames, defaultIndexGroups, gnames);
 
-        make_pull_coords(ir->pull);
+        checkPullCoords(ir->pull->group, ir->pull->coord);
     }
 
     if (ir->bRot)
index c4f0434e435fad7100b8af67f8dbc3a8b700c709..65050546ff4003ed05215ab1662986c553a95a47 100644 (file)
@@ -43,6 +43,7 @@
 
 #include "gromacs/fileio/readinp.h"
 #include "gromacs/math/vectypes.h"
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/real.h"
 
 namespace gmx
@@ -59,6 +60,8 @@ struct t_blocka;
 struct t_grpopts;
 struct t_inpfile;
 struct t_inputrec;
+struct t_pull_group;
+struct t_pull_coord;
 struct t_rot;
 struct warninp;
 typedef warninp* warninp_t;
@@ -152,14 +155,14 @@ void do_index(const char*                   mdparin,
 
 std::vector<std::string> read_pullparams(std::vector<t_inpfile>* inp, pull_params_t* pull, warninp_t wi);
 /* Reads the pull parameters, returns a list of the pull group names */
-
-void make_pull_groups(pull_params_t*                   pull,
-                      gmx::ArrayRef<const std::string> pullGroupNames,
-                      const t_blocka*                  grps,
-                      char**                           gnames);
+void process_pull_groups(gmx::ArrayRef<t_pull_group>      pullGroups,
+                         gmx::ArrayRef<const std::string> pullGroupNames,
+                         const t_blocka*                  grps,
+                         char**                           gnames);
 /* Process the pull group parameters after reading the index groups */
 
-void make_pull_coords(pull_params_t* pull);
+void checkPullCoords(gmx::ArrayRef<const t_pull_group> pullGroups,
+                     gmx::ArrayRef<const t_pull_coord> pullCoords);
 /* Process the pull coordinates after reading the pull groups */
 
 pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t* mtop, rvec* x, matrix box, real lambda, warninp_t wi);
index 91abefcc451616d5f6f33b11f3a84eb3ab421eba..3ea9de940767dda272154e3867c9a9e4af57595c 100644 (file)
@@ -55,6 +55,7 @@
 #include "gromacs/pulling/pull.h"
 #include "gromacs/topology/topology.h"
 #include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/futil.h"
@@ -71,21 +72,18 @@ static void string2dvec(const char buf[], dvec nums)
     }
 }
 
-static void init_pull_group(t_pull_group* pg, const char* wbuf)
+static std::vector<real> setupPullGroupWeights(const char* wbuf)
 {
     double d;
     int    n;
 
-    pg->nweight = 0;
+    std::vector<real> weight;
     while (sscanf(wbuf, "%lf %n", &d, &n) == 1)
     {
-        if (pg->nweight % 100 == 0)
-        {
-            srenew(pg->weight, pg->nweight + 100);
-        }
-        pg->weight[pg->nweight++] = d;
+        weight.push_back(d);
         wbuf += n;
     }
+    return weight;
 }
 
 static void process_pull_dim(char* dim_buf, ivec dim, const t_pull_coord* pcrd)
@@ -293,9 +291,6 @@ std::vector<std::string> read_pullparams(std::vector<t_inpfile>* inp, pull_param
     char provider[STRLEN], groups[STRLEN], dim_buf[STRLEN];
     char wbuf[STRLEN], origin_buf[STRLEN], vec_buf[STRLEN];
 
-    t_pull_group* pgrp;
-    t_pull_coord* pcrd;
-
     /* read pull parameters */
     printStringNoNewline(inp, "Cylinder radius for dynamic reaction force groups (nm)");
     pull->cylinder_r     = get_ereal(inp, "pull-cylinder-r", 1.5, wi);
@@ -325,10 +320,6 @@ std::vector<std::string> read_pullparams(std::vector<t_inpfile>* inp, pull_param
         gmx_fatal(FARGS, "pull-ncoords should be >= 1");
     }
 
-    snew(pull->group, pull->ngroup);
-
-    snew(pull->coord, pull->ncoord);
-
     /* pull group options */
     printStringNoNewline(inp, "Group and coordinate parameters");
 
@@ -336,61 +327,64 @@ std::vector<std::string> read_pullparams(std::vector<t_inpfile>* inp, pull_param
     std::vector<std::string> pullGroups(pull->ngroup);
     char                     readBuffer[STRLEN];
     /* Group 0 is the absolute reference, we don't read anything for 0 */
+    pull->group.emplace_back(t_pull_group());
     for (int groupNum = 1; groupNum < pull->ngroup; groupNum++)
     {
-        pgrp = &pull->group[groupNum];
+        t_pull_group pullGroup; //= &pull->group[groupNum];
         sprintf(buf, "pull-group%d-name", groupNum);
         setStringEntry(inp, buf, readBuffer, "");
         pullGroups[groupNum] = readBuffer;
         sprintf(buf, "pull-group%d-weights", groupNum);
         setStringEntry(inp, buf, wbuf, "");
         sprintf(buf, "pull-group%d-pbcatom", groupNum);
-        pgrp->pbcatom = get_eint(inp, buf, 0, wi);
+        pullGroup.pbcatom = get_eint(inp, buf, 0, wi);
 
         /* Initialize the pull group */
-        init_pull_group(pgrp, wbuf);
+        pullGroup.weight = setupPullGroupWeights(wbuf);
+        pull->group.emplace_back(pullGroup);
     }
 
     /* Read the pull coordinates */
     for (int coordNum = 1; coordNum < pull->ncoord + 1; coordNum++)
     {
-        pcrd = &pull->coord[coordNum - 1];
+        t_pull_coord pullCoord; // = &pull->coord[coordNum - 1];
         sprintf(buf, "pull-coord%d-type", coordNum);
-        pcrd->eType = get_eeenum(inp, buf, epull_names, wi);
+        pullCoord.eType = get_eeenum(inp, buf, epull_names, wi);
         sprintf(buf, "pull-coord%d-potential-provider", coordNum);
         setStringEntry(inp, buf, provider, "");
-        pcrd->externalPotentialProvider = gmx_strdup(provider);
+        pullCoord.externalPotentialProvider = gmx_strdup(provider);
         sprintf(buf, "pull-coord%d-geometry", coordNum);
-        pcrd->eGeom = get_eeenum(inp, buf, epullg_names, wi);
+        pullCoord.eGeom = get_eeenum(inp, buf, epullg_names, wi);
         sprintf(buf, "pull-coord%d-groups", coordNum);
         setStringEntry(inp, buf, groups, "");
 
-        switch (pcrd->eGeom)
+        switch (pullCoord.eGeom)
         {
-            case epullgDIHEDRAL: pcrd->ngroup = 6; break;
+            case epullgDIHEDRAL: pullCoord.ngroup = 6; break;
             case epullgDIRRELATIVE:
-            case epullgANGLE: pcrd->ngroup = 4; break;
-            default: pcrd->ngroup = 2; break;
+            case epullgANGLE: pullCoord.ngroup = 4; break;
+            default: pullCoord.ngroup = 2; break;
         }
 
-        nscan = sscanf(groups, "%d %d %d %d %d %d %d", &pcrd->group[0], &pcrd->group[1],
-                       &pcrd->group[2], &pcrd->group[3], &pcrd->group[4], &pcrd->group[5], &idum);
-        if (nscan != pcrd->ngroup)
+        nscan = sscanf(groups, "%d %d %d %d %d %d %d", &pullCoord.group[0], &pullCoord.group[1],
+                       &pullCoord.group[2], &pullCoord.group[3], &pullCoord.group[4],
+                       &pullCoord.group[5], &idum);
+        if (nscan != pullCoord.ngroup)
         {
             auto message =
                     gmx::formatString("%s should contain %d pull group indices with geometry %s",
-                                      buf, pcrd->ngroup, epullg_names[pcrd->eGeom]);
+                                      buf, pullCoord.ngroup, epullg_names[pullCoord.eGeom]);
             set_warning_line(wi, nullptr, -1);
             warning_error(wi, message);
         }
-        for (int g = 0; g < pcrd->ngroup; g++)
+        for (int g = 0; g < pullCoord.ngroup; g++)
         {
-            if (pcrd->group[g] < 0 || pcrd->group[g] >= pull->ngroup)
+            if (pullCoord.group[g] < 0 || pullCoord.group[g] >= pull->ngroup)
             {
                 /* Quit with a fatal error to avoid invalid memory access */
                 gmx_fatal(FARGS,
                           "%s contains an invalid pull group %d, you should have %d <= group <= %d",
-                          buf, pcrd->group[g], 0, pull->ngroup - 1);
+                          buf, pullCoord.group[g], 0, pull->ngroup - 1);
             }
         }
 
@@ -401,121 +395,114 @@ std::vector<std::string> read_pullparams(std::vector<t_inpfile>* inp, pull_param
         sprintf(buf, "pull-coord%d-vec", coordNum);
         setStringEntry(inp, buf, vec_buf, "0.0 0.0 0.0");
         sprintf(buf, "pull-coord%d-start", coordNum);
-        pcrd->bStart = (get_eeenum(inp, buf, yesno_names, wi) != 0);
+        pullCoord.bStart = (get_eeenum(inp, buf, yesno_names, wi) != 0);
         sprintf(buf, "pull-coord%d-init", coordNum);
-        pcrd->init = get_ereal(inp, buf, 0.0, wi);
+        pullCoord.init = get_ereal(inp, buf, 0.0, wi);
         sprintf(buf, "pull-coord%d-rate", coordNum);
-        pcrd->rate = get_ereal(inp, buf, 0.0, wi);
+        pullCoord.rate = get_ereal(inp, buf, 0.0, wi);
         sprintf(buf, "pull-coord%d-k", coordNum);
-        pcrd->k = get_ereal(inp, buf, 0.0, wi);
+        pullCoord.k = get_ereal(inp, buf, 0.0, wi);
         sprintf(buf, "pull-coord%d-kB", coordNum);
-        pcrd->kB = get_ereal(inp, buf, pcrd->k, wi);
+        pullCoord.kB = get_ereal(inp, buf, pullCoord.k, wi);
 
         /* Initialize the pull coordinate */
-        init_pull_coord(pcrd, coordNum, dim_buf, origin_buf, vec_buf, wi);
+        init_pull_coord(&pullCoord, coordNum, dim_buf, origin_buf, vec_buf, wi);
+        pull->coord.emplace_back(pullCoord);
     }
 
     return pullGroups;
 }
 
-void make_pull_groups(pull_params_t*                   pull,
-                      gmx::ArrayRef<const std::string> pullGroupNames,
-                      const t_blocka*                  grps,
-                      char**                           gnames)
+void process_pull_groups(gmx::ArrayRef<t_pull_group>      pullGroups,
+                         gmx::ArrayRef<const std::string> pullGroupNames,
+                         const t_blocka*                  grps,
+                         char**                           gnames)
 {
-    int           g, ig = -1, i;
-    t_pull_group* pgrp;
-
     /* Absolute reference group (might not be used) is special */
-    pgrp                = &pull->group[0];
-    pgrp->nat           = 0;
-    pgrp->pbcatom       = -1;
-    pgrp->pbcatom_input = -1;
+    pullGroups.front().pbcatom       = -1;
+    pullGroups.front().pbcatom_input = -1;
 
-    for (g = 1; g < pull->ngroup; g++)
+    // Skip pull group 0 here, as that is the absolute reference
+    for (int g = 1; g < int(pullGroups.size()); g++)
     {
-        pgrp                = &pull->group[g];
-        pgrp->pbcatom_input = pgrp->pbcatom;
+        auto& pullGroup = pullGroups[g];
 
         if (pullGroupNames[g].empty())
         {
             gmx_fatal(FARGS, "Pull option pull_group%d required by grompp has not been set.", g);
         }
 
-        ig        = search_string(pullGroupNames[g].c_str(), grps->nr, gnames);
-        pgrp->nat = grps->index[ig + 1] - grps->index[ig];
+        int ig                = search_string(pullGroupNames[g].c_str(), grps->nr, gnames);
+        int numPullGroupAtoms = grps->index[ig + 1] - grps->index[ig];
 
-        fprintf(stderr, "Pull group %d '%s' has %d atoms\n", g, pullGroupNames[g].c_str(), pgrp->nat);
+        fprintf(stderr, "Pull group %d '%s' has %d atoms\n", g, pullGroupNames[g].c_str(), numPullGroupAtoms);
 
-        if (pgrp->nat == 0)
+        if (numPullGroupAtoms == 0)
         {
             gmx_fatal(FARGS, "Pull group %d '%s' is empty", g, pullGroupNames[g].c_str());
         }
 
-        snew(pgrp->ind, pgrp->nat);
-        for (i = 0; i < pgrp->nat; i++)
+        for (int i = 0; i < numPullGroupAtoms; i++)
         {
-            pgrp->ind[i] = grps->a[grps->index[ig] + i];
+            pullGroup.ind.push_back(grps->a[grps->index[ig] + i]);
         }
 
-        if (pgrp->nweight > 0 && pgrp->nweight != pgrp->nat)
+        if (!pullGroup.weight.empty() && pullGroup.weight.size() != pullGroup.ind.size())
         {
             gmx_fatal(FARGS,
-                      "Number of weights (%d) for pull group %d '%s' does not match the number of "
-                      "atoms (%d)",
-                      pgrp->nweight, g, pullGroupNames[g].c_str(), pgrp->nat);
+                      "Number of weights (%ld) for pull group %d '%s' does not match the number of "
+                      "atoms (%ld)",
+                      gmx::ssize(pullGroup.weight), g, pullGroupNames[g].c_str(),
+                      gmx::ssize(pullGroup.ind));
         }
 
-        if (pgrp->nat == 1)
+        if (pullGroup.ind.size() == 1)
         {
             /* No pbc is required for this group */
-            pgrp->pbcatom = -1;
+            pullGroup.pbcatom = -1;
         }
         else
         {
-            if (pgrp->pbcatom > 0)
+            if (pullGroup.pbcatom > 0)
             {
-                pgrp->pbcatom -= 1;
+                pullGroup.pbcatom -= 1;
             }
-            else if (pgrp->pbcatom == 0)
+            else if (pullGroup.pbcatom == 0)
             {
-                pgrp->pbcatom = pgrp->ind[(pgrp->nat - 1) / 2];
+                pullGroup.pbcatom = pullGroup.ind[(pullGroup.ind.size() - 1) / 2];
             }
             else
             {
                 /* Use cosine weighting */
-                pgrp->pbcatom = -1;
+                pullGroup.pbcatom = -1;
             }
         }
     }
 }
 
-void make_pull_coords(pull_params_t* pull)
+void checkPullCoords(gmx::ArrayRef<const t_pull_group> pullGroups, gmx::ArrayRef<const t_pull_coord> pullCoords)
 {
-    int           c;
-    t_pull_coord* pcrd;
-
-    for (c = 0; c < pull->ncoord; c++)
+    for (int c = 0; c < int(pullCoords.size()); c++)
     {
-        pcrd = &pull->coord[c];
+        t_pull_coord pcrd = pullCoords[c];
 
-        if (pcrd->group[0] < 0 || pcrd->group[0] >= pull->ngroup || pcrd->group[1] < 0
-            || pcrd->group[1] >= pull->ngroup)
+        if (pcrd.group[0] < 0 || pcrd.group[0] >= int(pullGroups.size()) || pcrd.group[1] < 0
+            || pcrd.group[1] >= int(pullGroups.size()))
         {
             gmx_fatal(FARGS,
                       "Pull group index in pull-coord%d-groups out of range, should be between %d "
                       "and %d",
-                      c + 1, 0, pull->ngroup + 1);
+                      c + 1, 0, int(pullGroups.size()) + 1);
         }
 
-        if (pcrd->group[0] == pcrd->group[1])
+        if (pcrd.group[0] == pcrd.group[1])
         {
             gmx_fatal(FARGS, "Identical pull group indices in pull-coord%d-groups", c + 1);
         }
 
-        if (pcrd->eGeom == epullgCYL)
+        if (pcrd.eGeom == epullgCYL)
         {
-            if (pull->group[pcrd->group[0]].nweight > 0)
+            if (!pullGroups[pcrd.group[0]].weight.empty())
             {
                 gmx_fatal(
                         FARGS,
@@ -527,13 +514,12 @@ void make_pull_coords(pull_params_t* pull)
 
 pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t* mtop, rvec* x, matrix box, real lambda, warninp_t wi)
 {
-    pull_params_t* pull;
-    pull_t*        pull_work;
-    t_pbc          pbc;
-    int            c;
-    double         t_start;
+    pull_t* pull_work;
+    t_pbc   pbc;
+    int     c;
+    double  t_start;
 
-    pull = ir->pull;
+    pull_params_t*           pull = ir->pull.get();
     gmx::LocalAtomSetManager atomSets;
     pull_work    = init_pull(nullptr, pull, ir, mtop, nullptr, &atomSets, lambda);
     auto mdAtoms = gmx::makeMDAtoms(nullptr, *mtop, *ir, false);
@@ -617,8 +603,8 @@ pull_t* set_pull_init(t_inputrec* ir, const gmx_mtop_t* mtop, rvec* x, matrix bo
 
         pgrp0 = &pull->group[pcrd->group[0]];
         pgrp1 = &pull->group[pcrd->group[1]];
-        fprintf(stderr, "%8d  %8d  %8d\n", pcrd->group[0], pgrp0->nat, pgrp0->pbcatom + 1);
-        fprintf(stderr, "%8d  %8d  %8d ", pcrd->group[1], pgrp1->nat, pgrp1->pbcatom + 1);
+        fprintf(stderr, "%8d  %8zu  %8d\n", pcrd->group[0], pgrp0->ind.size(), pgrp0->pbcatom + 1);
+        fprintf(stderr, "%8d  %8zu  %8d ", pcrd->group[1], pgrp1->ind.size(), pgrp1->pbcatom + 1);
 
         if (pcrd->bStart)
         {
index 49d682320906e12184951599239dd2528767fb26..21397f342b3dedc28d3fd458da2e1235255bca42 100644 (file)
@@ -672,7 +672,7 @@ bool Constraints::Impl::apply(bool                      bLog,
 
     if (econq == ConstraintVariable::Positions)
     {
-        if (ir.bPull && pull_have_constraint(pull_work))
+        if (ir.bPull && pull_have_constraint(*pull_work))
         {
             if (EI_DYNAMICS(ir.eI))
             {
index b4bcbf975adbc67e57408b16e38798c38da82e7e..216a49ecf2c35340482dd1bcc20f9f666299db83 100644 (file)
@@ -239,7 +239,7 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
     bEner_[F_DISPCORR]   = (ir->eDispCorr != edispcNO);
     bEner_[F_DISRESVIOL] = (gmx_mtop_ftype_count(mtop, F_DISRES) > 0);
     bEner_[F_ORIRESDEV]  = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
-    bEner_[F_COM_PULL]   = ((ir->bPull && pull_have_potential(pull_work)) || ir->bRot);
+    bEner_[F_COM_PULL]   = ((ir->bPull && pull_have_potential(*pull_work)) || ir->bRot);
 
     MdModulesEnergyOutputToDensityFittingRequestChecker mdModulesAddOutputToDensityFittingFieldRequest;
     mdModulesNotifier.simulationSetupNotifications_.notify(&mdModulesAddOutputToDensityFittingFieldRequest);
index 75c81573729ba619852b82ed27a533c112d7ea82..0b6143a0ce2ed984ecdf7cf5456c4196e43f5459 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
@@ -106,7 +106,7 @@ std::unique_ptr<Constraints> makeConstraints(const gmx_mtop_t& mtop,
     GMX_RELEASE_ASSERT(!ir.bPull || pull_work != nullptr,
                        "When COM pulling is active, it must be initialized before constraints are "
                        "initialized");
-    bool doPullingWithConstraints = ir.bPull && pull_have_constraint(pull_work);
+    bool doPullingWithConstraints = ir.bPull && pull_have_constraint(*pull_work);
     if (numConstraints + numSettles == 0 && !doPullingWithConstraints && !doEssentialDynamics)
     {
         // No work, so don't make a Constraints object.
index 17a9e44bfa2b41640631b40e7930fdef6823a57e..47b6baac1bf0ef2fddddae85a59d91230af4e361 100644 (file)
@@ -548,7 +548,7 @@ static bool haveSpecialForces(const t_inputrec&          inputrec,
 {
 
     return ((computeForces && forceProviders.hasForceProvider()) || // forceProviders
-            (inputrec.bPull && pull_have_potential(pull_work)) ||   // pull
+            (inputrec.bPull && pull_have_potential(*pull_work)) ||  // pull
             inputrec.bRot ||                                        // enforced rotation
             (ed != nullptr) ||                                      // flooding
             (inputrec.bIMD && computeForces));                      // IMD
@@ -623,7 +623,7 @@ static void computeSpecialForces(FILE*                          fplog,
         forceProviders->calculateForces(forceProviderInput, &forceProviderOutput);
     }
 
-    if (inputrec->bPull && pull_have_potential(pull_work))
+    if (inputrec->bPull && pull_have_potential(*pull_work))
     {
         const int mtsLevel = forceGroupMtsLevel(inputrec->mtsLevels, gmx::MtsForceGroups::Pull);
         if (mtsLevel == 0 || stepWork.computeSlowForces)
@@ -1642,7 +1642,7 @@ void do_force(FILE*                               fplog,
 
     ForceOutputs* forceOutNonbonded = nonbondedAtMtsLevel1 ? forceOutMtsLevel1 : &forceOutMtsLevel0;
 
-    if (inputrec->bPull && pull_have_constraint(pull_work))
+    if (inputrec->bPull && pull_have_constraint(*pull_work))
     {
         clear_pull_forces(pull_work);
     }
index 96d0666e03ed29fbfaad97d38f6235e40691f753..b95eb0b84e3667eca582cd070297a1625a8388b0 100644 (file)
@@ -386,7 +386,7 @@ void gmx::LegacySimulator::do_md()
                            "Virtual sites are not supported with the GPU update.\n");
         GMX_RELEASE_ASSERT(ed == nullptr,
                            "Essential dynamics is not supported with the GPU update.\n");
-        GMX_RELEASE_ASSERT(!ir->bPull || !pull_have_constraint(ir->pull),
+        GMX_RELEASE_ASSERT(!ir->bPull || !pull_have_constraint(*ir->pull),
                            "Constraints pulling is not supported with the GPU update.\n");
         GMX_RELEASE_ASSERT(fcdata.orires->nr == 0,
                            "Orientation restraints are not supported with the GPU update.\n");
index 90275f50a3605c562a59e205c8c6ad41f81eae91..9c60986fa382906011474ac80df20468189795db 100644 (file)
@@ -1630,7 +1630,7 @@ int Mdrunner::mdrunner()
         if (inputrec->bPull)
         {
             /* Initialize pull code */
-            pull_work = init_pull(fplog, inputrec->pull, inputrec.get(), &mtop, cr, &atomSets,
+            pull_work = init_pull(fplog, inputrec->pull.get(), inputrec.get(), &mtop, cr, &atomSets,
                                   inputrec->fepvals->init_lambda);
             if (inputrec->pull->bXOutAverage || inputrec->pull->bFOutAverage)
             {
index 4f0ec9ca8feb62887dc2e51971ce07ae6f48e846..345f6fd5bef8b2833a43ce8a05c8c41eb7ff013d 100644 (file)
@@ -275,28 +275,6 @@ gmx_bool ir_vdw_might_be_zero_at_cutoff(const t_inputrec* ir)
     return (ir_vdw_is_zero_at_cutoff(ir) || ir->vdwtype == evdwUSER);
 }
 
-static void done_pull_group(t_pull_group* pgrp)
-{
-    if (pgrp->nat > 0)
-    {
-        sfree(pgrp->ind);
-        sfree(pgrp->weight);
-    }
-}
-
-static void done_pull_params(pull_params_t* pull)
-{
-    int i;
-
-    for (i = 0; i < pull->ngroup + 1; i++)
-    {
-        done_pull_group(pull->group);
-    }
-
-    sfree(pull->group);
-    sfree(pull->coord);
-}
-
 static void done_lambdas(t_lambda* fep)
 {
     if (fep->n_lambda > 0)
@@ -349,11 +327,6 @@ void done_inputrec(t_inputrec* ir)
     sfree(ir->expandedvals);
     sfree(ir->simtempvals);
 
-    if (ir->pull)
-    {
-        done_pull_params(ir->pull);
-        sfree(ir->pull);
-    }
     done_t_rot(ir->rot);
     delete ir->params;
 }
@@ -486,8 +459,8 @@ static void pr_pull_group(FILE* fp, int indent, int g, const t_pull_group* pgrp)
     pr_indent(fp, indent);
     fprintf(fp, "pull-group %d:\n", g);
     indent += 2;
-    pr_ivec_block(fp, indent, "atom", pgrp->ind, pgrp->nat, TRUE);
-    pr_rvec(fp, indent, "weight", pgrp->weight, pgrp->nweight, TRUE);
+    pr_ivec_block(fp, indent, "atom", pgrp->ind.data(), pgrp->ind.size(), TRUE);
+    pr_rvec(fp, indent, "weight", pgrp->weight.data(), pgrp->weight.size(), TRUE);
     PI("pbcatom", pgrp->pbcatom);
 }
 
@@ -500,7 +473,7 @@ static void pr_pull_coord(FILE* fp, int indent, int c, const t_pull_coord* pcrd)
     PS("type", EPULLTYPE(pcrd->eType));
     if (pcrd->eType == epullEXTERNAL)
     {
-        PS("potential-provider", pcrd->externalPotentialProvider);
+        PS("potential-provider", pcrd->externalPotentialProvider.c_str());
     }
     PS("geometry", EPULLGEOM(pcrd->eGeom));
     for (g = 0; g < pcrd->ngroup; g++)
@@ -629,29 +602,29 @@ static void pr_fepvals(FILE* fp, int indent, const t_lambda* fep, gmx_bool bMDPf
     PS("dhdl-derivatives", DHDLDERIVATIVESTYPE(fep->dhdl_derivatives));
 };
 
-static void pr_pull(FILE* fp, int indent, const pull_params_t* pull)
+static void pr_pull(FILE* fp, int indent, const pull_params_t& pull)
 {
     int g;
 
-    PR("pull-cylinder-r", pull->cylinder_r);
-    PR("pull-constr-tol", pull->constr_tol);
-    PS("pull-print-COM", EBOOL(pull->bPrintCOM));
-    PS("pull-print-ref-value", EBOOL(pull->bPrintRefValue));
-    PS("pull-print-components", EBOOL(pull->bPrintComp));
-    PI("pull-nstxout", pull->nstxout);
-    PI("pull-nstfout", pull->nstfout);
-    PS("pull-pbc-ref-prev-step-com", EBOOL(pull->bSetPbcRefToPrevStepCOM));
-    PS("pull-xout-average", EBOOL(pull->bXOutAverage));
-    PS("pull-fout-average", EBOOL(pull->bFOutAverage));
-    PI("pull-ngroups", pull->ngroup);
-    for (g = 0; g < pull->ngroup; g++)
+    PR("pull-cylinder-r", pull.cylinder_r);
+    PR("pull-constr-tol", pull.constr_tol);
+    PS("pull-print-COM", EBOOL(pull.bPrintCOM));
+    PS("pull-print-ref-value", EBOOL(pull.bPrintRefValue));
+    PS("pull-print-components", EBOOL(pull.bPrintComp));
+    PI("pull-nstxout", pull.nstxout);
+    PI("pull-nstfout", pull.nstfout);
+    PS("pull-pbc-ref-prev-step-com", EBOOL(pull.bSetPbcRefToPrevStepCOM));
+    PS("pull-xout-average", EBOOL(pull.bXOutAverage));
+    PS("pull-fout-average", EBOOL(pull.bFOutAverage));
+    PI("pull-ngroups", pull.ngroup);
+    for (g = 0; g < pull.ngroup; g++)
     {
-        pr_pull_group(fp, indent, g, &pull->group[g]);
+        pr_pull_group(fp, indent, g, &pull.group[g]);
     }
-    PI("pull-ncoords", pull->ncoord);
-    for (g = 0; g < pull->ncoord; g++)
+    PI("pull-ncoords", pull.ncoord);
+    for (g = 0; g < pull.ncoord; g++)
     {
-        pr_pull_coord(fp, indent, g, &pull->coord[g]);
+        pr_pull_coord(fp, indent, g, &pull.coord[g]);
     }
 }
 
@@ -995,7 +968,7 @@ void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir,
         PS("pull", EBOOL(ir->bPull));
         if (ir->bPull)
         {
-            pr_pull(fp, indent, ir->pull);
+            pr_pull(fp, indent, *ir->pull);
         }
 
         /* AWH BIASING */
@@ -1476,14 +1449,12 @@ void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real f
     gmx::compareKeyValueTrees(&writer, *ir1->params, *ir2->params, ftol, abstol);
 }
 
-void comp_pull_AB(FILE* fp, pull_params_t* pull, real ftol, real abstol)
+void comp_pull_AB(FILE* fp, const pull_params_t& pull, real ftol, real abstol)
 {
-    int i;
-
-    for (i = 0; i < pull->ncoord; i++)
+    for (int i = 0; i < pull.ncoord; i++)
     {
         fprintf(fp, "comparing pull coord %d\n", i);
-        cmp_real(fp, "pull-coord->k", -1, pull->coord[i].k, pull->coord[i].kB, ftol, abstol);
+        cmp_real(fp, "pull-coord->k", -1, pull.coord[i].k, pull.coord[i].kB, ftol, abstol);
     }
 }
 
index 40ed79ea8de4909d66637d16e4b412c3c827e3fa..46cb924592cacf8d4fc62a9f2773747cb065ce35 100644 (file)
@@ -515,7 +515,7 @@ struct t_inputrec // NOLINT (clang-analyzer-optin.performance.Padding)
     //! Do we do COM pulling?
     gmx_bool bPull;
     //! The data for center of mass pulling
-    pull_params_t* pull;
+    std::unique_ptr<pull_params_t> pull;
 
     /* AWH bias data */
     //! Whether to use AWH biasing for PMF calculations
@@ -620,7 +620,7 @@ void pr_inputrec(FILE* fp, int indent, const char* title, const t_inputrec* ir,
 
 void cmp_inputrec(FILE* fp, const t_inputrec* ir1, const t_inputrec* ir2, real ftol, real abstol);
 
-void comp_pull_AB(FILE* fp, pull_params_t* pull, real ftol, real abstol);
+void comp_pull_AB(FILE* fp, const pull_params_t& pull, real ftol, real abstol);
 
 
 gmx_bool inputrecDeform(const t_inputrec* ir);
index eae024308941025c74229387ea123e86735bc07f..fbd35e98e799d1a0b0ec789d072274e6b7356f80 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2015,2016,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2015,2016,2018,2019,2020, by the GROMACS development team, led by
  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
  * and including many others, as listed in the AUTHORS file in the
  * top-level source directory and at http://www.gromacs.org.
 #ifndef GMX_MDTYPES_PULL_PARAMS_H
 #define GMX_MDTYPES_PULL_PARAMS_H
 
+#include <array>
+#include <string>
+#include <vector>
+
 #include "gromacs/math/vectypes.h"
 #include "gromacs/mdtypes/md_enums.h"
 #include "gromacs/utility/basedefinitions.h"
 /*! \cond INTERNAL */
 
 /*! \brief Struct that defines a pull group */
-typedef struct
+struct t_pull_group
 {
-    int   nat;     /**< Number of atoms in the pull group */
-    int*  ind;     /**< The global atoms numbers */
-    int   nweight; /**< The number of weights (0 or nat) */
-    real* weight;  /**< Weights (use all 1 when weight==NULL) */
-    int   pbcatom; /**< The reference atom for pbc (global number) */
-    int   pbcatom_input; /**< The reference atom for pbc (global number) as specified in the input parameters */
-} t_pull_group;
+    std::vector<int>  ind;     /**< The global atoms numbers */
+    std::vector<real> weight;  /**< Weights (use all 1 when weight==NULL) */
+    int               pbcatom; /**< The reference atom for pbc (global number) */
+    int               pbcatom_input; /**< The reference atom for pbc (global number) as specified in the input parameters */
+};
 
 /*! Maximum number of pull groups that can be used in a pull coordinate */
 static const int c_pullCoordNgroupMax = 6;
 
 /*! \brief Struct that defines a pull coordinate */
-typedef struct
+struct t_pull_coord
 {
-    int      eType; /**< The pull type: umbrella, constraint, ... */
-    char*    externalPotentialProvider; /**< Name of the module providing the external potential, only used with eType==epullEXTERNAL */
-    int      eGeom;                       /**< The pull geometry */
-    int      ngroup;                      /**< The number of groups, depends on eGeom */
-    int      group[c_pullCoordNgroupMax]; /**< The pull groups: indices into the group arrays in pull_t and pull_params_t, ngroup indices are used */
-    ivec     dim;                         /**< Used to select components for constraint */
-    rvec     origin;                      /**< The origin for the absolute reference */
-    rvec     vec;                         /**< The pull vector, direction or position */
-    gmx_bool bStart;                      /**< Set init based on the initial structure */
-    real     init;                        /**< Initial reference displacement (nm) or (deg) */
-    real     rate;                        /**< Rate of motion (nm/ps) or (deg/ps) */
-    real     k; /**< Force constant (kJ/(mol nm^2) or kJ/(mol rad^2) for umbrella pull type, or kJ/(mol nm) or kJ/(mol rad) for constant force pull type */
-    real     kB; /**< Force constant for state B */
-} t_pull_coord;
+    int                                   eType; /**< The pull type: umbrella, constraint, ... */
+    std::string                           externalPotentialProvider; /**< Name of the module providing the external potential, only used with eType==epullEXTERNAL */
+    int                                   eGeom;  /**< The pull geometry */
+    int                                   ngroup; /**< The number of groups, depends on eGeom */
+    std::array<int, c_pullCoordNgroupMax> group; /**< The pull groups: indices into the group arrays in pull_t and pull_params_t, ngroup indices are used */
+    gmx::IVec                             dim;   /**< Used to select components for constraint */
+    gmx::RVec                             origin; /**< The origin for the absolute reference */
+    gmx::RVec                             vec;    /**< The pull vector, direction or position */
+    bool                                  bStart; /**< Set init based on the initial structure */
+    real                                  init; /**< Initial reference displacement (nm) or (deg) */
+    real                                  rate; /**< Rate of motion (nm/ps) or (deg/ps) */
+    real                                  k; /**< Force constant (kJ/(mol nm^2) or kJ/(mol rad^2) for umbrella pull type, or kJ/(mol nm) or kJ/(mol rad) for constant force pull type */
+    real                                  kB; /**< Force constant for state B */
+};
 
 /*! \brief Struct containing all pull parameters */
-typedef struct pull_params_t
+struct pull_params_t
 {
-    int      ngroup;         /**< Number of pull groups */
-    int      ncoord;         /**< Number of pull coordinates */
-    real     cylinder_r;     /**< Radius of cylinder for dynamic COM (nm) */
-    real     constr_tol;     /**< Absolute tolerance for constraints in (nm) */
-    gmx_bool bPrintCOM;      /**< Print coordinates of COM for each coord */
-    gmx_bool bPrintRefValue; /**< Print the reference value for each coord */
-    gmx_bool bPrintComp; /**< Print cartesian components for each coord with geometry=distance */
-    gmx_bool bSetPbcRefToPrevStepCOM; /**< Use the COM of each group from the previous step as reference */
-    int      nstxout;                 /**< Output interval for pull x */
-    int      nstfout;                 /**< Output interval for pull f */
-    bool     bXOutAverage; /**< Write the average coordinate during the output interval */
-    bool     bFOutAverage; /**< Write the average force during the output interval */
+    int  ngroup;         /**< Number of pull groups */
+    int  ncoord;         /**< Number of pull coordinates */
+    real cylinder_r;     /**< Radius of cylinder for dynamic COM (nm) */
+    real constr_tol;     /**< Absolute tolerance for constraints in (nm) */
+    bool bPrintCOM;      /**< Print coordinates of COM for each coord */
+    bool bPrintRefValue; /**< Print the reference value for each coord */
+    bool bPrintComp;     /**< Print cartesian components for each coord with geometry=distance */
+    bool bSetPbcRefToPrevStepCOM; /**< Use the COM of each group from the previous step as reference */
+    int  nstxout;                 /**< Output interval for pull x */
+    int  nstfout;                 /**< Output interval for pull f */
+    bool bXOutAverage;            /**< Write the average coordinate during the output interval */
+    bool bFOutAverage;            /**< Write the average force during the output interval */
 
-    t_pull_group* group; /**< groups to pull/restrain/etc/ */
-    t_pull_coord* coord; /**< the pull coordinates */
-} pull_params_t;
+    std::vector<t_pull_group> group; /**< groups to pull/restrain/etc/ */
+    std::vector<t_pull_coord> coord; /**< the pull coordinates */
+};
 
 /*! \endcond */
 
index 67cf362e79100c5a1b87802201b240503a891a6c..3d3773233716aab3c617ba7c7006646de3efbdcf 100644 (file)
@@ -88,7 +88,7 @@ extern template LocalAtomSet LocalAtomSetManager::add<void, void>(ArrayRef<const
 
 static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
 {
-    if (params.nat <= 1)
+    if (params.ind.size() <= 1)
     {
         /* no PBC required */
         return epgrppbcNONE;
@@ -208,7 +208,7 @@ static void apply_forces_grp(const pull_group_work_t* pgrp,
 {
     auto localAtomIndices = pgrp->atomSet.localIndex();
 
-    if (pgrp->params.nat == 1 && pgrp->atomSet.numAtomsLocal() == 1)
+    if (pgrp->params.ind.size() == 1 && pgrp->atomSet.numAtomsLocal() == 1)
     {
         /* Only one atom and our rank has this atom: we can skip
          * the mass weighting, which means that this code also works
@@ -356,7 +356,7 @@ static void apply_forces_coord(struct pull_t*               pull,
             apply_forces_vec_torque(pull, &pcrd, masses, f);
         }
 
-        if (pull->group[pcrd.params.group[0]].params.nat > 0)
+        if (!pull->group[pcrd.params.group[0]].params.ind.empty())
         {
             apply_forces_grp(&pull->group[pcrd.params.group[0]], masses, forces.force01, -1, f,
                              pull->nthreads);
@@ -457,7 +457,7 @@ static void low_get_pull_coord_dr(const struct pull_t*     pull,
     const pull_group_work_t* pgrp0 = &pull->group[pcrd->params.group[0]];
 
     /* Only the first group can be an absolute reference, in that case nat=0 */
-    if (pgrp0->params.nat == 0)
+    if (pgrp0->params.ind.empty())
     {
         for (int m = 0; m < DIM; m++)
         {
@@ -1380,15 +1380,15 @@ void register_external_pull_potential(struct pull_t* pull, int coord_index, cons
                 provider, coord_index + 1, epull_names[pcrd->params.eType], epull_names[epullEXTERNAL]);
     }
 
-    GMX_RELEASE_ASSERT(pcrd->params.externalPotentialProvider != nullptr,
+    GMX_RELEASE_ASSERT(!pcrd->params.externalPotentialProvider.empty(),
                        "The external potential provider string for a pull coordinate is NULL");
 
-    if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider) != 0)
+    if (gmx_strcasecmp(provider, pcrd->params.externalPotentialProvider.c_str()) != 0)
     {
         gmx_fatal(FARGS,
                   "Module '%s' attempted to register an external potential for pull coordinate %d "
                   "which expects the external potential to be provided by a module named '%s'",
-                  provider, coord_index + 1, pcrd->params.externalPotentialProvider);
+                  provider, coord_index + 1, pcrd->params.externalPotentialProvider.c_str());
     }
 
     /* Lock to avoid (extremely unlikely) simultaneous reading and writing of
@@ -1435,7 +1435,7 @@ static void check_external_potential_registration(const struct pull_t* pull)
                   "pull coordinates. The first coordinate without provider is number %zu, which "
                   "expects a module named '%s' to provide the external potential.",
                   pull->numUnregisteredExternalPotentials, c + 1,
-                  pull->coord[c].params.externalPotentialProvider);
+                  pull->coord[c].params.externalPotentialProvider.c_str());
     }
 }
 
@@ -1746,7 +1746,8 @@ static void init_pull_group_index(FILE*              fplog,
      * But we still want to have the correct mass-weighted COMs.
      * So we store the real masses in the weights.
      */
-    const bool setWeights = (pg->params.nweight > 0 || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
+    const bool setWeights =
+            (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD);
 
     /* In parallel, store we need to extract localWeights from weights at DD time */
     std::vector<real>& weights = ((cr && PAR(cr)) ? pg->globalWeights : pg->localWeights);
@@ -1759,7 +1760,7 @@ static void init_pull_group_index(FILE*              fplog,
     double wmass   = 0;
     double wwmass  = 0;
     int    molb    = 0;
-    for (int i = 0; i < pg->params.nat; i++)
+    for (int i = 0; i < int(pg->params.ind.size()); i++)
     {
         int ii = pg->params.ind[i];
         if (bConstraint && ir->opts.nFreeze)
@@ -1784,7 +1785,7 @@ static void init_pull_group_index(FILE*              fplog,
             m = (1 - lambda) * atom.m + lambda * atom.mB;
         }
         real w;
-        if (pg->params.nweight > 0)
+        if (!pg->params.weight.empty())
         {
             w = pg->params.weight[i];
         }
@@ -1834,7 +1835,7 @@ static void init_pull_group_index(FILE*              fplog,
         /* We can have single atom groups with zero mass with potential pulling
          * without cosine weighting.
          */
-        if (pg->params.nat == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
+        if (pg->params.ind.size() == 1 && !bConstraint && pg->epgrppbc != epgrppbcCOS)
         {
             /* With one atom the mass doesn't matter */
             wwmass = 1;
@@ -1842,13 +1843,13 @@ static void init_pull_group_index(FILE*              fplog,
         else
         {
             gmx_fatal(FARGS, "The total%s mass of pull group %d is zero",
-                      pg->params.weight ? " weighted" : "", g);
+                      !pg->params.weight.empty() ? " weighted" : "", g);
         }
     }
     if (fplog)
     {
-        fprintf(fplog, "Pull group %d: %5d atoms, mass %9.3f", g, pg->params.nat, tmass);
-        if (pg->params.weight || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
+        fprintf(fplog, "Pull group %d: %5zu atoms, mass %9.3f", g, pg->params.ind.size(), tmass);
+        if (!pg->params.weight.empty() || EI_ENERGY_MINIMIZATION(ir->eI) || ir->eI == eiBD)
         {
             fprintf(fplog, ", weighted mass %9.3f", wmass * wmass / wwmass);
         }
@@ -1869,7 +1870,7 @@ static void init_pull_group_index(FILE*              fplog,
         int ndim = 0;
         for (int d = 0; d < DIM; d++)
         {
-            ndim += pulldim_con[d] * pg->params.nat;
+            ndim += pulldim_con[d] * pg->params.ind.size();
         }
         if (fplog && nfrozen > 0 && nfrozen < ndim)
         {
@@ -1899,15 +1900,10 @@ struct pull_t* init_pull(FILE*                     fplog,
 
     /* Copy the pull parameters */
     pull->params = *pull_params;
-    /* Avoid pointer copies */
-    pull->params.group = nullptr;
-    pull->params.coord = nullptr;
 
     for (int i = 0; i < pull_params->ngroup; ++i)
     {
-        pull->group.emplace_back(pull_params->group[i],
-                                 atomSets->add({ pull_params->group[i].ind,
-                                                 pull_params->group[i].ind + pull_params->group[i].nat }),
+        pull->group.emplace_back(pull_params->group[i], atomSets->add(pull_params->group[i].ind),
                                  pull_params->bSetPbcRefToPrevStepCOM);
     }
 
@@ -1932,7 +1928,7 @@ struct pull_t* init_pull(FILE*                     fplog,
     pull->bXOutAverage = pull_params->bXOutAverage;
     pull->bFOutAverage = pull_params->bFOutAverage;
 
-    GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0,
+    GMX_RELEASE_ASSERT(pull->group[0].params.ind.empty(),
                        "pull group 0 is an absolute reference group and should not contain atoms");
 
     pull->numCoordinatesWithExternalPotential = 0;
@@ -2076,8 +2072,8 @@ struct pull_t* init_pull(FILE*                     fplog,
         bAbs = FALSE;
         for (const pull_coord_work_t& coord : pull->coord)
         {
-            if (pull->group[coord.params.group[0]].params.nat == 0
-                || pull->group[coord.params.group[1]].params.nat == 0)
+            if (pull->group[coord.params.group[0]].params.ind.empty()
+                || pull->group[coord.params.group[1]].params.ind.empty())
             {
                 bAbs = TRUE;
             }
@@ -2106,7 +2102,7 @@ struct pull_t* init_pull(FILE*                     fplog,
         // Don't include the reference group 0 in loop
         for (size_t g = 1; g < pull->group.size(); g++)
         {
-            if (pull->group[g].params.nat > 1 && pull->group[g].params.pbcatom < 0)
+            if (pull->group[g].params.ind.size() > 1 && pull->group[g].params.pbcatom < 0)
             {
                 /* We are using cosine weighting */
                 fprintf(fplog, "Cosine weighting is used for group %zu\n", g);
@@ -2126,7 +2122,7 @@ struct pull_t* init_pull(FILE*                     fplog,
         pull_group_work_t* pgrp;
 
         pgrp = &pull->group[g];
-        if (pgrp->params.nat > 0)
+        if (!pgrp->params.ind.empty())
         {
             /* There is an issue when a group is used in multiple coordinates
              * and constraints are applied in different dimensions with atoms
@@ -2182,7 +2178,7 @@ struct pull_t* init_pull(FILE*                     fplog,
             {
                 case epgrppbcREFAT: pull->bRefAt = TRUE; break;
                 case epgrppbcCOS:
-                    if (pgrp->params.weight != nullptr)
+                    if (!pgrp->params.weight.empty())
                     {
                         gmx_fatal(FARGS,
                                   "Pull groups can not have relative weights and cosine weighting "
@@ -2225,7 +2221,7 @@ struct pull_t* init_pull(FILE*                     fplog,
         {
             if (coord.params.eGeom == epullgCYL)
             {
-                if (pull->group[coord.params.group[0]].params.nat == 0)
+                if (pull->group[coord.params.group[0]].params.ind.empty())
                 {
                     gmx_fatal(FARGS,
                               "A cylinder pull group is not supported when using absolute "
@@ -2352,25 +2348,21 @@ void finish_pull(struct pull_t* pull)
     destroy_pull(pull);
 }
 
-gmx_bool pull_have_potential(const struct pull_t* pull)
+bool pull_have_potential(const pull_t& pull)
 {
-    return pull->bPotential;
+    return pull.bPotential;
 }
 
-gmx_bool pull_have_constraint(const struct pull_t* pull)
+bool pull_have_constraint(const pull_t& pull)
 {
-    return pull->bConstraint;
+    return pull.bConstraint;
 }
 
-bool pull_have_constraint(const pull_params_t* pullParameters)
+bool pull_have_constraint(const pull_params_t& pullParameters)
 {
-    if (pullParameters == nullptr)
+    for (int c = 0; c < pullParameters.ncoord; c++)
     {
-        return false;
-    }
-    for (int c = 0; c < pullParameters->ncoord; c++)
-    {
-        if (pullParameters->coord[c].eType == epullCONSTRAINT)
+        if (pullParameters.coord[c].eType == epullCONSTRAINT)
         {
             return true;
         }
index faf4152bbf6cdf3ff2a2e9b2934bbc9c2f1a891f..2f87228ecb929a2edfe1c7d9827cabb7c9f6d4bc 100644 (file)
@@ -329,20 +329,20 @@ bool pullCheckPbcWithinGroup(const pull_t&                  pull,
  *
  * \param[in] pull     The pull data structure.
  */
-gmx_bool pull_have_potential(const struct pull_t* pull);
+bool pull_have_potential(const pull_t& pull);
 
 
 /*! \brief Returns if we have pull coordinates with constraint pulling.
  *
  * \param[in] pull     The pull data structure.
  */
-gmx_bool pull_have_constraint(const struct pull_t* pull);
+bool pull_have_constraint(const pull_t& pull);
 
 /*! \brief Returns if inputrec has pull coordinates with constraint pulling.
  *
  * \param[in] pullParameters  Pulling input parameters from input record.
  */
-bool pull_have_constraint(const pull_params_t* pullParameters);
+bool pull_have_constraint(const pull_params_t& pullParameters);
 
 /*! \brief Returns the maxing distance for pulling
  *
index f3163957deb92e91f292cd82b02cdfef3796cfc2..a1c646c826ead2e847862b71783b96de6ff48b2a 100644 (file)
@@ -610,7 +610,7 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc
                  * Note that with constraint pulling the mass does matter, but
                  * in that case a check group mass != 0 has been done before.
                  */
-                if (pgrp->params.nat == 1 && pgrp->atomSet.numAtomsLocal() == 1
+                if (pgrp->params.ind.size() == 1 && pgrp->atomSet.numAtomsLocal() == 1
                     && masses[pgrp->atomSet.localIndex()[0]] == 0)
                 {
                     GMX_ASSERT(xp == nullptr,
@@ -724,7 +724,7 @@ void pull_calc_coms(const t_commrec* cr, pull_t* pull, const real* masses, t_pbc
         pgrp = &pull->group[g];
         if (pgrp->needToCalcCom)
         {
-            GMX_ASSERT(pgrp->params.nat > 0,
+            GMX_ASSERT(!pgrp->params.ind.empty(),
                        "Normal pull groups should have atoms, only group 0, which should have "
                        "bCalcCom=FALSE has nat=0");
 
@@ -1049,7 +1049,7 @@ void initPullComFromPrevStep(const t_commrec* cr, pull_t* pull, const real* mass
 
         if (pgrp->needToCalcCom && pgrp->epgrppbc == epgrppbcPREVSTEPCOM)
         {
-            GMX_ASSERT(pgrp->params.nat > 1,
+            GMX_ASSERT(pgrp->params.ind.size() > 1,
                        "Groups with no atoms, or only one atom, should not "
                        "use the COM from the previous step as reference.");
 
index 020e6187c1803268c0991549b8f02a2b0a369df2..eb88b85d87d8be21e11694a5df8abff2995b6c48 100644 (file)
@@ -640,7 +640,7 @@ bool decideWhetherToUseGpuForUpdate(const bool                     isDomainDecom
     {
         errorMessage += "Essential dynamics is not supported.\n";
     }
-    if (inputrec.bPull && pull_have_constraint(inputrec.pull))
+    if (inputrec.bPull && pull_have_constraint(*inputrec.pull))
     {
         errorMessage += "Constraints pulling is not supported.\n";
     }
index 5aa81509a78d80ed329d127eb340d7fb0b3366d7..dc7d16b57372ea83e8f53aa3a009a8bd48fc40a7 100644 (file)
@@ -125,7 +125,7 @@ static void comp_tpx(const char* fn1, const char* fn2, gmx_bool bRMSD, real ftol
         {
             if (ir[0]->bPull)
             {
-                comp_pull_AB(stdout, ir[0]->pull, ftol, abstol);
+                comp_pull_AB(stdout, *ir[0]->pull, ftol, abstol);
             }
             compareMtopAB(stdout, mtop[0], ftol, abstol);
         }