*/
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)
"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);
}
}
/* 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)
"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)
{
}
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,
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
* \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,
{
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];
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);
}
}
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
serializer->doInt(&pcrd->ngroup);
if (pcrd->ngroup <= c_pullCoordNgroupMax)
{
- serializer->doIntArray(pcrd->group, pcrd->ngroup);
+ serializer->doIntArray(pcrd->group.data(), pcrd->ngroup);
}
else
{
pcrd->ngroup = 0;
}
- serializer->doIvec(&pcrd->dim);
+ serializer->doIvec(&pcrd->dim.as_vec());
}
else
{
serializer->doInt(&pcrd->group[2]);
serializer->doInt(&pcrd->group[3]);
}
- serializer->doIvec(&pcrd->dim);
+ serializer->doIvec(&pcrd->dim.as_vec());
}
else
{
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);
{
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 */
}
}
- pull->bPrintCOM = (pull->group[0].nat > 0);
+ pull->bPrintCOM = (!pull->group[0].ind.empty());
}
else
{
{
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);
}
}
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);
}
#include <cstdlib>
#include <algorithm>
+#include <memory>
#include <string>
#include "gromacs/applied_forces/awh/read_params.h"
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)
{
* 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++)
{
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];
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)
#include "gromacs/fileio/readinp.h"
#include "gromacs/math/vectypes.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/real.h"
namespace gmx
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;
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);
#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"
}
}
-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)
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);
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");
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);
}
}
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,
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);
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)
{
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))
{
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);
/*
* 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.
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.
{
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
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)
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);
}
"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");
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)
{
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)
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;
}
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);
}
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++)
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]);
}
}
PS("pull", EBOOL(ir->bPull));
if (ir->bPull)
{
- pr_pull(fp, indent, ir->pull);
+ pr_pull(fp, indent, *ir->pull);
}
/* AWH BIASING */
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);
}
}
//! 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
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);
/*
* 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 */
static int groupPbcFromParams(const t_pull_group& params, bool setPbcRefToPrevStepCOM)
{
- if (params.nat <= 1)
+ if (params.ind.size() <= 1)
{
/* no PBC required */
return epgrppbcNONE;
{
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
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);
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++)
{
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
"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());
}
}
* 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);
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)
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];
}
/* 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;
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);
}
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)
{
/* 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);
}
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;
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;
}
// 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);
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
{
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 "
{
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 "
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;
}
*
* \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
*
* 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,
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");
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.");
{
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";
}
{
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);
}