From e6a3a331fb4dcda858d8894555fa74d3ea8aaa04 Mon Sep 17 00:00:00 2001 From: Joe Jordan Date: Tue, 13 Oct 2020 12:35:33 +0000 Subject: [PATCH] Fix segfault in pull reading Triggers asserton now later on :( Change-Id: I08a390de653e8bbdf23288ed8a28b6dc7ad6305d --- .../applied_forces/awh/read_params.cpp | 18 +- src/gromacs/applied_forces/awh/read_params.h | 2 +- src/gromacs/fileio/tpxio.cpp | 69 +++----- src/gromacs/gmxpreprocess/grompp.cpp | 2 +- src/gromacs/gmxpreprocess/readir.cpp | 13 +- src/gromacs/gmxpreprocess/readir.h | 15 +- src/gromacs/gmxpreprocess/readpull.cpp | 166 ++++++++---------- src/gromacs/mdlib/constr.cpp | 2 +- src/gromacs/mdlib/energyoutput.cpp | 2 +- src/gromacs/mdlib/makeconstraints.h | 4 +- src/gromacs/mdlib/sim_util.cpp | 6 +- src/gromacs/mdrun/md.cpp | 2 +- src/gromacs/mdrun/runner.cpp | 2 +- src/gromacs/mdtypes/inputrec.cpp | 77 +++----- src/gromacs/mdtypes/inputrec.h | 4 +- src/gromacs/mdtypes/pull_params.h | 82 ++++----- src/gromacs/pulling/pull.cpp | 72 ++++---- src/gromacs/pulling/pull.h | 6 +- src/gromacs/pulling/pullutil.cpp | 6 +- src/gromacs/taskassignment/decidegpuusage.cpp | 2 +- src/gromacs/tools/check.cpp | 2 +- 21 files changed, 249 insertions(+), 305 deletions(-) diff --git a/src/gromacs/applied_forces/awh/read_params.cpp b/src/gromacs/applied_forces/awh/read_params.cpp index e5945f17d6..5654d57af4 100644 --- a/src/gromacs/applied_forces/awh/read_params.cpp +++ b/src/gromacs/applied_forces/awh/read_params.cpp @@ -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 diff --git a/src/gromacs/applied_forces/awh/read_params.h b/src/gromacs/applied_forces/awh/read_params.h index 7da77f1fdd..bb1cd63fa9 100644 --- a/src/gromacs/applied_forces/awh/read_params.h +++ b/src/gromacs/applied_forces/awh/read_params.h @@ -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, diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index 9ebd98f4b3..9ea2337f86 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -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(); } - do_pull(serializer, ir->pull, file_version, ePullOld); + do_pull(serializer, ir->pull.get(), file_version, ePullOld); } } diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 6a20af125b..46a48b93e4 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -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); } diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index b2d504ca23..2edc32384d 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -45,6 +45,7 @@ #include #include +#include #include #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(); + 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) diff --git a/src/gromacs/gmxpreprocess/readir.h b/src/gromacs/gmxpreprocess/readir.h index c4f0434e43..65050546ff 100644 --- a/src/gromacs/gmxpreprocess/readir.h +++ b/src/gromacs/gmxpreprocess/readir.h @@ -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 read_pullparams(std::vector* 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 pullGroupNames, - const t_blocka* grps, - char** gnames); +void process_pull_groups(gmx::ArrayRef pullGroups, + gmx::ArrayRef 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 pullGroups, + gmx::ArrayRef 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); diff --git a/src/gromacs/gmxpreprocess/readpull.cpp b/src/gromacs/gmxpreprocess/readpull.cpp index 91abefcc45..3ea9de9407 100644 --- a/src/gromacs/gmxpreprocess/readpull.cpp +++ b/src/gromacs/gmxpreprocess/readpull.cpp @@ -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 setupPullGroupWeights(const char* wbuf) { double d; int n; - pg->nweight = 0; + std::vector 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 read_pullparams(std::vector* 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 read_pullparams(std::vector* 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 read_pullparams(std::vector* inp, pull_param std::vector 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 read_pullparams(std::vector* 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 pullGroupNames, - const t_blocka* grps, - char** gnames) +void process_pull_groups(gmx::ArrayRef pullGroups, + gmx::ArrayRef 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 pullGroups, gmx::ArrayRef 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) { diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index 49d6823209..21397f342b 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -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)) { diff --git a/src/gromacs/mdlib/energyoutput.cpp b/src/gromacs/mdlib/energyoutput.cpp index b4bcbf975a..216a49ecf2 100644 --- a/src/gromacs/mdlib/energyoutput.cpp +++ b/src/gromacs/mdlib/energyoutput.cpp @@ -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); diff --git a/src/gromacs/mdlib/makeconstraints.h b/src/gromacs/mdlib/makeconstraints.h index 75c8157372..0b6143a0ce 100644 --- a/src/gromacs/mdlib/makeconstraints.h +++ b/src/gromacs/mdlib/makeconstraints.h @@ -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 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. diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 17a9e44bfa..47b6baac1b 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -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); } diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 96d0666e03..b95eb0b84e 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -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"); diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp index 90275f50a3..9c60986fa3 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -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) { diff --git a/src/gromacs/mdtypes/inputrec.cpp b/src/gromacs/mdtypes/inputrec.cpp index 4f0ec9ca8f..345f6fd5be 100644 --- a/src/gromacs/mdtypes/inputrec.cpp +++ b/src/gromacs/mdtypes/inputrec.cpp @@ -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); } } diff --git a/src/gromacs/mdtypes/inputrec.h b/src/gromacs/mdtypes/inputrec.h index 40ed79ea8d..46cb924592 100644 --- a/src/gromacs/mdtypes/inputrec.h +++ b/src/gromacs/mdtypes/inputrec.h @@ -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; /* 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); diff --git a/src/gromacs/mdtypes/pull_params.h b/src/gromacs/mdtypes/pull_params.h index eae0243089..fbd35e98e7 100644 --- a/src/gromacs/mdtypes/pull_params.h +++ b/src/gromacs/mdtypes/pull_params.h @@ -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. @@ -48,6 +48,10 @@ #ifndef GMX_MDTYPES_PULL_PARAMS_H #define GMX_MDTYPES_PULL_PARAMS_H +#include +#include +#include + #include "gromacs/math/vectypes.h" #include "gromacs/mdtypes/md_enums.h" #include "gromacs/utility/basedefinitions.h" @@ -56,56 +60,54 @@ /*! \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 ind; /**< The global atoms numbers */ + std::vector 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 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 group; /**< groups to pull/restrain/etc/ */ + std::vector coord; /**< the pull coordinates */ +}; /*! \endcond */ diff --git a/src/gromacs/pulling/pull.cpp b/src/gromacs/pulling/pull.cpp index 67cf362e79..3d37732337 100644 --- a/src/gromacs/pulling/pull.cpp +++ b/src/gromacs/pulling/pull.cpp @@ -88,7 +88,7 @@ extern template LocalAtomSet LocalAtomSetManager::add(ArrayRefatomSet.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& 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; } diff --git a/src/gromacs/pulling/pull.h b/src/gromacs/pulling/pull.h index faf4152bbf..2f87228ecb 100644 --- a/src/gromacs/pulling/pull.h +++ b/src/gromacs/pulling/pull.h @@ -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 * diff --git a/src/gromacs/pulling/pullutil.cpp b/src/gromacs/pulling/pullutil.cpp index f3163957de..a1c646c826 100644 --- a/src/gromacs/pulling/pullutil.cpp +++ b/src/gromacs/pulling/pullutil.cpp @@ -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."); diff --git a/src/gromacs/taskassignment/decidegpuusage.cpp b/src/gromacs/taskassignment/decidegpuusage.cpp index 020e6187c1..eb88b85d87 100644 --- a/src/gromacs/taskassignment/decidegpuusage.cpp +++ b/src/gromacs/taskassignment/decidegpuusage.cpp @@ -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"; } diff --git a/src/gromacs/tools/check.cpp b/src/gromacs/tools/check.cpp index 5aa81509a7..dc7d16b573 100644 --- a/src/gromacs/tools/check.cpp +++ b/src/gromacs/tools/check.cpp @@ -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); } -- 2.22.0