From 999ebbb1b20c2b57c21d4a4b78c5780403848afd Mon Sep 17 00:00:00 2001 From: Paul Bauer Date: Wed, 13 Feb 2019 18:46:15 +0100 Subject: [PATCH] Refactor t_molinfo to MoleculeInformation class Refs #2833 Change-Id: I66a6afc8ade636fb3002492e6c4e1d30ae51019a --- src/gromacs/gmxpreprocess/convparm.cpp | 10 +- src/gromacs/gmxpreprocess/convparm.h | 7 +- src/gromacs/gmxpreprocess/grompp.cpp | 221 ++++++++++++------------ src/gromacs/gmxpreprocess/grompp_impl.h | 59 +++++-- src/gromacs/gmxpreprocess/tomorse.cpp | 66 +++---- src/gromacs/gmxpreprocess/tomorse.h | 6 +- src/gromacs/gmxpreprocess/topio.cpp | 116 ++++++------- src/gromacs/gmxpreprocess/topio.h | 39 +++-- src/gromacs/gmxpreprocess/toppush.cpp | 49 +++--- src/gromacs/gmxpreprocess/toppush.h | 11 +- src/gromacs/gmxpreprocess/toputil.cpp | 23 --- src/gromacs/gmxpreprocess/toputil.h | 6 +- 12 files changed, 320 insertions(+), 293 deletions(-) diff --git a/src/gromacs/gmxpreprocess/convparm.cpp b/src/gromacs/gmxpreprocess/convparm.cpp index df6889d52a..1ea48876c9 100644 --- a/src/gromacs/gmxpreprocess/convparm.cpp +++ b/src/gromacs/gmxpreprocess/convparm.cpp @@ -497,7 +497,7 @@ static void append_interaction(InteractionList *ilist, } } -static void enter_function(t_params *p, t_functype ftype, int comb, real reppow, +static void enter_function(const t_params *p, t_functype ftype, int comb, real reppow, gmx_ffparams_t *ffparams, InteractionList *il, bool bNB, bool bAppend) { @@ -520,7 +520,8 @@ static void enter_function(t_params *p, t_functype ftype, int comb, real reppow, } void convert_params(int atnr, t_params nbtypes[], - t_molinfo *mi, t_molinfo *intermolecular_interactions, + gmx::ArrayRef mi, + const MoleculeInformation *intermolecular_interactions, int comb, double reppow, real fudgeQQ, gmx_mtop_t *mtop) { @@ -528,7 +529,6 @@ void convert_params(int atnr, t_params nbtypes[], unsigned long flags; gmx_ffparams_t *ffp; gmx_moltype_t *molt; - t_params *plist; ffp = &mtop->ffparams; ffp->atnr = atnr; @@ -548,7 +548,7 @@ void convert_params(int atnr, t_params nbtypes[], { molt->ilist[i].iatoms.clear(); - plist = mi[mt].plist; + const t_params *plist = mi[mt].plist; flags = interaction_function[i].flags; if ((i != F_LJ) && (i != F_BHAM) && ((flags & IF_BOND) || @@ -572,7 +572,7 @@ void convert_params(int atnr, t_params nbtypes[], { (*mtop->intermolecular_ilist)[i].iatoms.clear(); - plist = intermolecular_interactions->plist; + const t_params *plist = intermolecular_interactions->plist; if (plist[i].nr > 0) { diff --git a/src/gromacs/gmxpreprocess/convparm.h b/src/gromacs/gmxpreprocess/convparm.h index 496d866d3d..2b61d6fcfd 100644 --- a/src/gromacs/gmxpreprocess/convparm.h +++ b/src/gromacs/gmxpreprocess/convparm.h @@ -38,15 +38,16 @@ #ifndef GMX_GMXPREPROCESS_CONVPARM_H #define GMX_GMXPREPROCESS_CONVPARM_H +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/real.h" struct gmx_mtop_t; -struct t_molinfo; +struct MoleculeInformation; struct t_params; void convert_params(int atnr, t_params nbtypes[], - t_molinfo *mi, - t_molinfo *intermolecular_interactions, + gmx::ArrayRef mi, + const MoleculeInformation *intermolecular_interactions, int comb, double reppow, real fudgeQQ, gmx_mtop_t *mtop); diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 08ce54c750..c7d7ffa6eb 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -104,16 +104,37 @@ #include "gromacs/utility/smalloc.h" #include "gromacs/utility/snprintf.h" -static int rm_interactions(int ifunc, int nrmols, t_molinfo mols[]) +void MoleculeInformation::initMolInfo() { - int i, n; + init_plist(plist); + init_block(&cgs); + init_block(&mols); + init_blocka(&excls); + init_t_atoms(&atoms, 0, FALSE); +} + +void MoleculeInformation::partialCleanUp() +{ + done_block(&mols); + done_plist(plist); +} - n = 0; +void MoleculeInformation::fullCleanUp() +{ + done_atom (&atoms); + done_block(&cgs); + done_block(&mols); + done_plist(plist); +} + +static int rm_interactions(int ifunc, gmx::ArrayRef mols) +{ + int n = 0; /* For all the molecule types */ - for (i = 0; i < nrmols; i++) + for (auto &mol : mols) { - n += mols[i].plist[ifunc].nr; - mols[i].plist[ifunc].nr = 0; + n += mol.plist[ifunc].nr; + mol.plist[ifunc].nr = 0; } return n; } @@ -406,7 +427,7 @@ static void check_shells_inputrec(gmx_mtop_t *mtop, /* TODO Decide whether this function can be consolidated with * gmx_mtop_ftype_count */ -static int nint_ftype(gmx_mtop_t *mtop, t_molinfo *mi, int ftype) +static int nint_ftype(gmx_mtop_t *mtop, gmx::ArrayRef mi, int ftype) { int nint = 0; for (const gmx_molblock_t &molb : mtop->molblock) @@ -421,73 +442,63 @@ static int nint_ftype(gmx_mtop_t *mtop, t_molinfo *mi, int ftype) * in the order of use in the molblocks, * unused molecule types are deleted. */ -static void renumber_moltypes(gmx_mtop_t *sys, - int *nmolinfo, t_molinfo **molinfo) +static void renumber_moltypes(gmx_mtop_t *sys, + std::vector *molinfo) { - int *order, norder; - t_molinfo *minew; - snew(order, *nmolinfo); - norder = 0; + std::vector order; for (gmx_molblock_t &molblock : sys->molblock) { - int i; - for (i = 0; i < norder; i++) + const auto found = std::find_if(order.begin(), order.end(), + [&molblock](const auto &entry) + { return molblock.type == entry; }); + if (found == order.end()) { - if (order[i] == molblock.type) - { - break; - } + /* This type did not occur yet, add it */ + order.push_back(molblock.type); + molblock.type = order.size() - 1; } - if (i == norder) + else { - /* This type did not occur yet, add it */ - order[norder] = molblock.type; - /* Renumber the moltype in the topology */ - norder++; + molblock.type = std::distance(order.begin(), found); } - molblock.type = i; } /* We still need to reorder the molinfo structs */ - snew(minew, norder); - for (int mi = 0; mi < *nmolinfo; mi++) + std::vector minew(order.size()); + int index = 0; + for (auto &mi : *molinfo) { - int i; - for (i = 0; i < norder; i++) + const auto found = std::find(order.begin(), order.end(), index); + if (found != order.end()) { - if (order[i] == mi) - { - break; - } - } - if (i == norder) - { - done_mi(&(*molinfo)[mi]); + int pos = std::distance(order.begin(), found); + minew[pos] = mi; } else { - minew[i] = (*molinfo)[mi]; + // Need to manually clean up memory .... + mi.fullCleanUp(); } + index++; } - sfree(order); - sfree(*molinfo); - *nmolinfo = norder; *molinfo = minew; } -static void molinfo2mtop(int nmi, t_molinfo *mi, gmx_mtop_t *mtop) +static void molinfo2mtop(gmx::ArrayRef mi, gmx_mtop_t *mtop) { - mtop->moltype.resize(nmi); - for (int m = 0; m < nmi; m++) + mtop->moltype.resize(mi.size()); + int pos = 0; + for (const auto &mol : mi) { - gmx_moltype_t &molt = mtop->moltype[m]; - molt.name = mi[m].name; - molt.atoms = mi[m].atoms; + gmx_moltype_t &molt = mtop->moltype[pos]; + molt.name = mol.name; + molt.atoms = mol.atoms; /* ilists are copied later */ - molt.cgs = mi[m].cgs; - molt.excls = mi[m].excls; + molt.cgs = mol.cgs; + molt.excls = mol.excls; + pos++; } } @@ -496,23 +507,23 @@ new_status(const char *topfile, const char *topppfile, const char *confin, t_gromppopts *opts, t_inputrec *ir, gmx_bool bZero, bool bGenVel, bool bVerbose, t_state *state, gpp_atomtype *atype, gmx_mtop_t *sys, - int *nmi, t_molinfo **mi, t_molinfo **intermolecular_interactions, + std::vector *mi, + std::unique_ptr *intermolecular_interactions, t_params plist[], int *comb, double *reppow, real *fudgeQQ, gmx_bool bMorse, warninp *wi) { - t_molinfo *molinfo = nullptr; - std::vector molblock; - int i, nrmols, nmismatch; - bool ffParametrizedWithHBondConstraints; - char buf[STRLEN]; - char warn_buf[STRLEN]; + std::vector molblock; + int i, nrmols, nmismatch; + bool ffParametrizedWithHBondConstraints; + char buf[STRLEN]; + char warn_buf[STRLEN]; /* TOPOLOGY processing */ sys->name = do_top(bVerbose, topfile, topppfile, opts, bZero, &(sys->symtab), plist, comb, reppow, fudgeQQ, - atype, &nrmols, &molinfo, intermolecular_interactions, + atype, &nrmols, mi, intermolecular_interactions, ir, &molblock, &ffParametrizedWithHBondConstraints, @@ -534,23 +545,23 @@ new_status(const char *topfile, const char *topppfile, const char *confin, /* Add a new molblock to the topology */ sys->molblock.push_back(molb); } - sys->natoms += molb.nmol*molinfo[sys->molblock.back().type].atoms.nr; + sys->natoms += molb.nmol*(*mi)[sys->molblock.back().type].atoms.nr; } if (sys->molblock.empty()) { gmx_fatal(FARGS, "No molecules were defined in the system"); } - renumber_moltypes(sys, &nrmols, &molinfo); + renumber_moltypes(sys, mi); if (bMorse) { - convert_harmonics(nrmols, molinfo, atype); + convert_harmonics(*mi, atype); } if (ir->eDisre == edrNone) { - i = rm_interactions(F_DISRES, nrmols, molinfo); + i = rm_interactions(F_DISRES, *mi); if (i > 0) { set_warning_line(wi, "unknown", -1); @@ -560,7 +571,7 @@ new_status(const char *topfile, const char *topppfile, const char *confin, } if (!opts->bOrire) { - i = rm_interactions(F_ORIRES, nrmols, molinfo); + i = rm_interactions(F_ORIRES, *mi); if (i > 0) { set_warning_line(wi, "unknown", -1); @@ -570,7 +581,7 @@ new_status(const char *topfile, const char *topppfile, const char *confin, } /* Copy structures from msys to sys */ - molinfo2mtop(nrmols, molinfo, sys); + molinfo2mtop(*mi, sys); gmx_mtop_finalize(sys); @@ -654,9 +665,9 @@ new_status(const char *topfile, const char *topppfile, const char *confin, fprintf(stderr, "double-checking input for internal consistency...\n"); } { - bool bHasNormalConstraints = 0 < (nint_ftype(sys, molinfo, F_CONSTR) + - nint_ftype(sys, molinfo, F_CONSTRNC)); - bool bHasAnyConstraints = bHasNormalConstraints || 0 < nint_ftype(sys, molinfo, F_SETTLE); + bool bHasNormalConstraints = 0 < (nint_ftype(sys, *mi, F_CONSTR) + + nint_ftype(sys, *mi, F_CONSTRNC)); + bool bHasAnyConstraints = bHasNormalConstraints || 0 < nint_ftype(sys, *mi, F_SETTLE); double_check(ir, state->box, bHasNormalConstraints, bHasAnyConstraints, @@ -686,9 +697,6 @@ new_status(const char *topfile, const char *topppfile, const char *confin, stop_cm(stdout, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array()); sfree(mass); } - - *nmi = nrmols; - *mi = molinfo; } static void copy_state(const char *slog, t_trxframe *fr, @@ -811,7 +819,9 @@ static void cont_status(const char *slog, const char *ener, } } -static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB, +static void read_posres(gmx_mtop_t *mtop, + gmx::ArrayRef molinfo, + gmx_bool bTopB, const char *fn, int rc_scaling, int ePBC, rvec com, @@ -826,7 +836,6 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB, int natoms, npbcdim = 0; char warn_buf[STRLEN]; int a, i, ai, j, k, nat_molb; - t_params *pr, *prfb; t_atom *atom; snew(top, 1); @@ -861,8 +870,8 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB, for (gmx_molblock_t &molb : mtop->molblock) { nat_molb = molb.nmol*mtop->moltype[molb.type].atoms.nr; - pr = &(molinfo[molb.type].plist[F_POSRES]); - prfb = &(molinfo[molb.type].plist[F_FBPOSRES]); + const t_params *pr = &(molinfo[molb.type].plist[F_POSRES]); + const t_params *prfb = &(molinfo[molb.type].plist[F_FBPOSRES]); if (pr->nr > 0 || prfb->nr > 0) { atom = mtop->moltype[molb.type].atoms.atom; @@ -988,13 +997,14 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB, sfree(hadAtom); } -static void gen_posres(gmx_mtop_t *mtop, t_molinfo *mi, +static void gen_posres(gmx_mtop_t *mtop, + gmx::ArrayRef mi, const char *fnA, const char *fnB, int rc_scaling, int ePBC, rvec com, rvec comB, warninp *wi) { - read_posres (mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi); + read_posres(mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi); /* It is safer to simply read the b-state posres rather than trying * to be smart and copy the positions. */ @@ -1022,7 +1032,7 @@ static void set_wall_atomtype(gpp_atomtype *at, t_gromppopts *opts, } } -static int nrdf_internal(t_atoms *atoms) +static int nrdf_internal(const t_atoms *atoms) { int i, nmass, nrdf; @@ -1192,17 +1202,18 @@ static void init_cmap_grid(gmx_cmap_t *cmap_grid, int ngrid, int grid_spacing) } -static int count_constraints(const gmx_mtop_t *mtop, t_molinfo *mi, warninp *wi) +static int count_constraints(const gmx_mtop_t *mtop, + gmx::ArrayRef mi, + warninp *wi) { int count, count_mol, i; - t_params *plist; char buf[STRLEN]; count = 0; for (const gmx_molblock_t &molb : mtop->molblock) { count_mol = 0; - plist = mi[molb.type].plist; + const t_params *plist = mi[molb.type].plist; for (i = 0; i < F_NRE; i++) { @@ -1580,7 +1591,7 @@ static void set_verlet_buffer(const gmx_mtop_t *mtop, int gmx_grompp(int argc, char *argv[]) { - const char *desc[] = { + const char *desc[] = { "[THISMODULE] (the gromacs preprocessor)", "reads a molecular topology file, checks the validity of the", "file, expands the topology from a molecular description to an atomic", @@ -1678,24 +1689,24 @@ int gmx_grompp(int argc, char *argv[]) "interpret the output messages before attempting to bypass them with", "this option." }; - t_gromppopts *opts; - int nmi; - t_molinfo *mi, *intermolecular_interactions; - gpp_atomtype *atype; - int nvsite, comb; - t_params *plist; - real fudgeQQ; - double reppow; - const char *mdparin; - int ntype; - bool bNeedVel, bGenVel; - gmx_bool have_atomnumber; - gmx_output_env_t *oenv; - gmx_bool bVerbose = FALSE; - warninp *wi; - char warn_buf[STRLEN]; - - t_filenm fnm[] = { + t_gromppopts *opts; + std::vector mi; + std::unique_ptr intermolecular_interactions; + gpp_atomtype *atype; + int nvsite, comb; + t_params *plist; + real fudgeQQ; + double reppow; + const char *mdparin; + int ntype; + bool bNeedVel, bGenVel; + gmx_bool have_atomnumber; + gmx_output_env_t *oenv; + gmx_bool bVerbose = FALSE; + warninp *wi; + char warn_buf[STRLEN]; + + t_filenm fnm[] = { { efMDP, nullptr, nullptr, ffREAD }, { efMDP, "-po", "mdout", ffWRITE }, { efSTX, "-c", nullptr, ffREAD }, @@ -1806,7 +1817,7 @@ int gmx_grompp(int argc, char *argv[]) t_state state; new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm), opts, ir, bZero, bGenVel, bVerbose, &state, - atype, &sys, &nmi, &mi, &intermolecular_interactions, + atype, &sys, &mi, &intermolecular_interactions, plist, &comb, &reppow, &fudgeQQ, opts->bMorse, wi); @@ -1842,7 +1853,7 @@ int gmx_grompp(int argc, char *argv[]) gmx_mtop_remove_chargegroups(&sys); } - if (count_constraints(&sys, mi, wi) && (ir->eConstrAlg == econtSHAKE)) + if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE)) { if (ir->eI == eiCG || ir->eI == eiLBFGS) { @@ -1976,7 +1987,7 @@ int gmx_grompp(int argc, char *argv[]) } ntype = get_atomtype_ntypes(atype); - convert_params(ntype, plist, mi, intermolecular_interactions, + convert_params(ntype, plist, mi, intermolecular_interactions.get(), comb, reppow, fudgeQQ, &sys); if (debug) @@ -2314,14 +2325,12 @@ int gmx_grompp(int argc, char *argv[]) sfree(opts); done_plist(plist); sfree(plist); - for (int i = 0; i < nmi; i++) + for (auto &mol : mi) { // Some of the contents of molinfo have been stolen, so - // done_mi can't be called. - done_block(&mi[i].mols); - done_plist(mi[i].plist); + // fullCleanUp can't be called. + mol.partialCleanUp(); } - sfree(mi); done_atomtype(atype); done_inputrec_strings(); output_env_done(oenv); diff --git a/src/gromacs/gmxpreprocess/grompp_impl.h b/src/gromacs/gmxpreprocess/grompp_impl.h index 62470454d9..8c6489a15f 100644 --- a/src/gromacs/gmxpreprocess/grompp_impl.h +++ b/src/gromacs/gmxpreprocess/grompp_impl.h @@ -94,17 +94,56 @@ struct t_excls int *e; /* The excluded atoms */ }; -struct t_molinfo + +/*! \libinternal \brief + * Holds the molecule information during preprocessing. + */ +struct MoleculeInformation { - char **name; - int nrexcl; /* Number of exclusions per atom */ - bool excl_set; /* Have exclusions been generated? */ - bool bProcessed; /* Has the mol been processed */ - t_atoms atoms; /* Atoms */ - t_block cgs; /* Charge groups */ - t_block mols; /* Molecules */ - t_blocka excls; /* Exclusions */ - t_params plist[F_NRE]; /* Parameters in old style */ + //! Name of the molecule. + char **name = nullptr; + //!Number of exclusions per atom. + int nrexcl = 0; + //! Have exclusions been generated?. + bool excl_set = false; + //! Has the mol been processed. + bool bProcessed = false; + //! Atoms in the moelcule. + t_atoms atoms; + //! Charge groups in the molecule + t_block cgs; + //! Molecules separated in datastructure. + t_block mols; + //! Exclusions in the molecule. + t_blocka excls; + //! Parameters in old style. + t_params plist[F_NRE]; + + /*! \brief + * Initializer. + * + * This should be removed as soon as the underlying datastructures + * have been cleaned up to use proper initialization and can be copy + * constructed. + */ + void initMolInfo(); + + /*! \brief + * Partial clean up function. + * + * Should be removed once this datastructure actually owns all its own memory and + * elements of it are not stolen by other structures and properly copy constructed + * or moved. + * Cleans up the mols and plist datastructures but not cgs and excls. + */ + void partialCleanUp(); + + /*! \brief + * Full clean up function. + * + * Should be removed once the destructor can always do this. + */ + void fullCleanUp(); }; struct t_mols diff --git a/src/gromacs/gmxpreprocess/tomorse.cpp b/src/gromacs/gmxpreprocess/tomorse.cpp index 90b2962746..9bd2c09262 100644 --- a/src/gromacs/gmxpreprocess/tomorse.cpp +++ b/src/gromacs/gmxpreprocess/tomorse.cpp @@ -187,14 +187,11 @@ static real search_e_diss(int n2m, t_2morse t2m[], char *ai, char *aj) } } -void convert_harmonics(int nrmols, t_molinfo mols[], gpp_atomtype *atype) +void convert_harmonics(gmx::ArrayRef mols, gpp_atomtype *atype) { int n2m; t_2morse *t2m; - int i, j, k, last, ni, nj; - int nrharm, nrmorse, bb; - real edis, kb, b0, beta; bool *bRemoveHarm; /* First get the data */ @@ -206,61 +203,63 @@ void convert_harmonics(int nrmols, t_molinfo mols[], gpp_atomtype *atype) } /* For all the molecule types */ - for (i = 0; (i < nrmols); i++) + int i = 0; + for (auto &mol : mols) { /* Check how many morse and harmonic BONDSs there are, increase size of * morse with the number of harmonics */ - nrmorse = mols[i].plist[F_MORSE].nr; + int nrmorse = mol.plist[F_MORSE].nr; - for (bb = 0; (bb < F_NRE); bb++) + for (int bb = 0; (bb < F_NRE); bb++) { if ((interaction_function[bb].flags & IF_BTYPE) && (bb != F_MORSE)) { - nrharm = mols[i].plist[bb].nr; - pr_alloc(nrharm, &(mols[i].plist[F_MORSE])); + int nrharm = mol.plist[bb].nr; + pr_alloc(nrharm, &(mol.plist[F_MORSE])); snew(bRemoveHarm, nrharm); /* Now loop over the harmonics, trying to convert them */ - for (j = 0; (j < nrharm); j++) + for (int j = 0; (j < nrharm); j++) { - ni = mols[i].plist[bb].param[j].ai(); - nj = mols[i].plist[bb].param[j].aj(); - edis = + int ni = mol.plist[bb].param[j].ai(); + int nj = mol.plist[bb].param[j].aj(); + real edis = search_e_diss(n2m, t2m, - get_atomtype_name(mols[i].atoms.atom[ni].type, atype), - get_atomtype_name(mols[i].atoms.atom[nj].type, atype)); + get_atomtype_name(mol.atoms.atom[ni].type, atype), + get_atomtype_name(mol.atoms.atom[nj].type, atype)); if (edis != 0) { - bRemoveHarm[j] = TRUE; - b0 = mols[i].plist[bb].param[j].c[0]; - kb = mols[i].plist[bb].param[j].c[1]; - beta = std::sqrt(kb/(2*edis)); - mols[i].plist[F_MORSE].param[nrmorse].a[0] = ni; - mols[i].plist[F_MORSE].param[nrmorse].a[1] = nj; - mols[i].plist[F_MORSE].param[nrmorse].c[0] = b0; - mols[i].plist[F_MORSE].param[nrmorse].c[1] = edis; - mols[i].plist[F_MORSE].param[nrmorse].c[2] = beta; + bRemoveHarm[j] = true; + real b0 = mol.plist[bb].param[j].c[0]; + real kb = mol.plist[bb].param[j].c[1]; + real beta = std::sqrt(kb/(2*edis)); + mol.plist[F_MORSE].param[nrmorse].a[0] = ni; + mol.plist[F_MORSE].param[nrmorse].a[1] = nj; + mol.plist[F_MORSE].param[nrmorse].c[0] = b0; + mol.plist[F_MORSE].param[nrmorse].c[1] = edis; + mol.plist[F_MORSE].param[nrmorse].c[2] = beta; nrmorse++; } } - mols[i].plist[F_MORSE].nr = nrmorse; + mol.plist[F_MORSE].nr = nrmorse; /* Now remove the harmonics */ - for (j = last = 0; (j < nrharm); j++) + int last = 0; + for (int j = 0; (j < nrharm); j++) { if (!bRemoveHarm[j]) { /* Copy it to the last position */ - for (k = 0; (k < MAXATOMLIST); k++) + for (int k = 0; (k < MAXATOMLIST); k++) { - mols[i].plist[bb].param[last].a[k] = - mols[i].plist[bb].param[j].a[k]; + mol.plist[bb].param[last].a[k] = + mol.plist[bb].param[j].a[k]; } - for (k = 0; (k < MAXFORCEPARAM); k++) + for (int k = 0; (k < MAXFORCEPARAM); k++) { - mols[i].plist[bb].param[last].c[k] = - mols[i].plist[bb].param[j].c[k]; + mol.plist[bb].param[last].c[k] = + mol.plist[bb].param[j].c[k]; } last++; } @@ -268,9 +267,10 @@ void convert_harmonics(int nrmols, t_molinfo mols[], gpp_atomtype *atype) sfree(bRemoveHarm); fprintf(stderr, "Converted %d out of %d %s to morse bonds for mol %d\n", nrharm-last, nrharm, interaction_function[bb].name, i); - mols[i].plist[bb].nr = last; + mol.plist[bb].nr = last; } } + i++; } sfree(t2m); } diff --git a/src/gromacs/gmxpreprocess/tomorse.h b/src/gromacs/gmxpreprocess/tomorse.h index 778e477cbc..9d2f510493 100644 --- a/src/gromacs/gmxpreprocess/tomorse.h +++ b/src/gromacs/gmxpreprocess/tomorse.h @@ -38,9 +38,11 @@ #ifndef GMX_GMXPREPROCESS_TOMORSE_H #define GMX_GMXPREPROCESS_TOMORSE_H +#include "gromacs/utility/arrayref.h" + struct gpp_atomtype; -struct t_molinfo; +struct MoleculeInformation; -void convert_harmonics(int nrmols, t_molinfo mols[], gpp_atomtype *atype); +void convert_harmonics(gmx::ArrayRef mols, gpp_atomtype *atype); #endif diff --git a/src/gromacs/gmxpreprocess/topio.cpp b/src/gromacs/gmxpreprocess/topio.cpp index 97a5ff7f99..5878613420 100644 --- a/src/gromacs/gmxpreprocess/topio.cpp +++ b/src/gromacs/gmxpreprocess/topio.cpp @@ -47,6 +47,7 @@ #include #include +#include #include #include @@ -349,9 +350,9 @@ static char ** cpp_opts(const char *define, const char *include, } -static void make_atoms_sys(const std::vector &molblock, - const t_molinfo *molinfo, - t_atoms *atoms) +static void make_atoms_sys(gmx::ArrayRef molblock, + gmx::ArrayRef molinfo, + t_atoms *atoms) { atoms->nr = 0; atoms->atom = nullptr; @@ -378,8 +379,8 @@ static char **read_topol(const char *infile, const char *outfile, t_symtab *symtab, gpp_atomtype *atype, int *nrmols, - t_molinfo **molinfo, - t_molinfo **intermolecular_interactions, + std::vector *molinfo, + std::unique_ptr *intermolecular_interactions, t_params plist[], int *combination_rule, double *reppow, @@ -392,32 +393,32 @@ static char **read_topol(const char *infile, const char *outfile, bool usingFullRangeElectrostatics, warninp *wi) { - FILE *out; - int i, sl, nb_funct; - char *pline = nullptr, **title = nullptr; - char line[STRLEN], errbuf[256], comb_str[256], nb_str[256]; - char genpairs[32]; - char *dirstr, *dummy2; - int nrcopies, nmol, nscan, ncombs, ncopy; - double fLJ, fQQ, fPOW; - t_molinfo *mi0 = nullptr; - DirStack *DS; - Directive d, newd; - t_nbparam **nbparam, **pair; - gmx::ExclusionBlocks *exclusionBlocks; - real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */ - bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B; - double qt = 0, qBt = 0; /* total charge */ - gpp_bond_atomtype *batype; - int lastcg = -1; - int dcatt = -1, nmol_couple; + FILE *out; + int sl, nb_funct; + char *pline = nullptr, **title = nullptr; + char line[STRLEN], errbuf[256], comb_str[256], nb_str[256]; + char genpairs[32]; + char *dirstr, *dummy2; + int nrcopies, nmol, nscan, ncombs, ncopy; + double fLJ, fQQ, fPOW; + MoleculeInformation *mi0 = nullptr; + DirStack *DS; + Directive d, newd; + t_nbparam **nbparam, **pair; + gmx::ExclusionBlocks *exclusionBlocks; + real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */ + bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B; + double qt = 0, qBt = 0; /* total charge */ + gpp_bond_atomtype *batype; + int lastcg = -1; + int dcatt = -1, nmol_couple; /* File handling variables */ - int status; - bool done; - gmx_cpp_t handle; - char *tmp_line = nullptr; - char warn_buf[STRLEN]; - const char *floating_point_arithmetic_tip = + int status; + bool done; + gmx_cpp_t handle; + char *tmp_line = nullptr; + char warn_buf[STRLEN]; + const char *floating_point_arithmetic_tip = "Total charge should normally be an integer. See\n" "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n" "for discussion on how close it should be to an integer.\n"; @@ -452,8 +453,6 @@ static char **read_topol(const char *infile, const char *outfile, *reppow = 12.0; /* Default value for repulsion power */ - *intermolecular_interactions = nullptr; - /* Init the number of CMAP torsion angles and grid spacing */ plist[F_CMAP].grid_spacing = 0; plist[F_CMAP].nc = 0; @@ -578,9 +577,9 @@ static char **read_topol(const char *infile, const char *outfile, * to process the intermolecular interactions * by making a "molecule" of the size of the system. */ - snew(*intermolecular_interactions, 1); - init_molinfo(*intermolecular_interactions); - mi0 = *intermolecular_interactions; + *intermolecular_interactions = std::make_unique( ); + mi0 = intermolecular_interactions->get(); + mi0->initMolInfo(); make_atoms_sys(*molblock, *molinfo, &mi0->atoms); } @@ -723,10 +722,11 @@ static char **read_topol(const char *infile, const char *outfile, bReadMolType = TRUE; } - push_molt(symtab, &nmol, molinfo, pline, wi); + push_molt(symtab, molinfo, pline, wi); + nmol = molinfo->size(); srenew(exclusionBlocks, nmol); exclusionBlocks[nmol-1].nr = 0; - mi0 = &((*molinfo)[nmol-1]); + mi0 = &molinfo->back(); mi0->atoms.haveMass = TRUE; mi0->atoms.haveCharge = TRUE; mi0->atoms.haveType = TRUE; @@ -791,7 +791,7 @@ static char **read_topol(const char *infile, const char *outfile, int whichmol; bool bCouple; - push_mol(nmol, *molinfo, pline, &whichmol, &nrcopies, wi); + push_mol(*molinfo, pline, &whichmol, &nrcopies, wi); mi0 = &((*molinfo)[whichmol]); molblock->resize(molblock->size() + 1); molblock->back().type = whichmol; @@ -924,7 +924,7 @@ static char **read_topol(const char *infile, const char *outfile, } DS_Done (&DS); - for (i = 0; i < nmol; i++) + for (int i = 0; i < gmx::index(molinfo->size()); i++) { gmx::doneExclusionBlocks(&(exclusionBlocks[i])); } @@ -934,7 +934,7 @@ static char **read_topol(const char *infile, const char *outfile, if (*intermolecular_interactions != nullptr) { - sfree(mi0->atoms.atom); + sfree(intermolecular_interactions->get()->atoms.atom); } *nrmols = nmol; @@ -942,24 +942,24 @@ static char **read_topol(const char *infile, const char *outfile, return title; } -char **do_top(bool bVerbose, - const char *topfile, - const char *topppfile, - t_gromppopts *opts, - bool bZero, - t_symtab *symtab, - t_params plist[], - int *combination_rule, - double *repulsion_power, - real *fudgeQQ, - gpp_atomtype *atype, - int *nrmols, - t_molinfo **molinfo, - t_molinfo **intermolecular_interactions, - const t_inputrec *ir, - std::vector *molblock, - bool *ffParametrizedWithHBondConstraints, - warninp *wi) +char **do_top(bool bVerbose, + const char *topfile, + const char *topppfile, + t_gromppopts *opts, + bool bZero, + t_symtab *symtab, + t_params plist[], + int *combination_rule, + double *repulsion_power, + real *fudgeQQ, + gpp_atomtype *atype, + int *nrmols, + std::vector *molinfo, + std::unique_ptr *intermolecular_interactions, + const t_inputrec *ir, + std::vector *molblock, + bool *ffParametrizedWithHBondConstraints, + warninp *wi) { /* Tmpfile might contain a long path */ const char *tmpfile; diff --git a/src/gromacs/gmxpreprocess/topio.h b/src/gromacs/gmxpreprocess/topio.h index c3ec5d5ff3..c71d5b7a38 100644 --- a/src/gromacs/gmxpreprocess/topio.h +++ b/src/gromacs/gmxpreprocess/topio.h @@ -38,6 +38,7 @@ #ifndef GMX_GMXPREPROCESS_TOPIO_H #define GMX_GMXPREPROCESS_TOPIO_H +#include #include #include "gromacs/utility/real.h" @@ -47,7 +48,7 @@ struct gmx_mtop_t; struct gpp_atomtype; struct t_gromppopts; struct t_inputrec; -struct t_molinfo; +struct MoleculeInformation; struct t_params; struct t_symtab; struct warninp; @@ -57,24 +58,24 @@ typedef warninp *warninp_t; double check_mol(const gmx_mtop_t *mtop, warninp_t wi); /* Check mass and charge */ -char **do_top(bool bVerbose, - const char *topfile, - const char *topppfile, - t_gromppopts *opts, - bool bZero, - t_symtab *symtab, - t_params plist[], - int *combination_rule, - double *repulsion_power, - real *fudgeQQ, - gpp_atomtype *atype, - int *nrmols, - t_molinfo **molinfo, - t_molinfo **intermolecular_interactions, - const t_inputrec *ir, - std::vector *molblock, - bool *ffParametrizedWithHBondConstraints, - warninp_t wi); +char **do_top(bool bVerbose, + const char *topfile, + const char *topppfile, + t_gromppopts *opts, + bool bZero, + t_symtab *symtab, + t_params plist[], + int *combination_rule, + double *repulsion_power, + real *fudgeQQ, + gpp_atomtype *atype, + int *nrmols, + std::vector *molinfo, + std::unique_ptr *intermolecular_interactions, + const t_inputrec *ir, + std::vector *molblock, + bool *ffParametrizedWithHBondConstraints, + warninp_t wi); /* This routine expects sys->molt[m].ilist to be of size F_NRE and ordered. */ void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi, GmxQmmmMode qmmmMode); diff --git a/src/gromacs/gmxpreprocess/toppush.cpp b/src/gromacs/gmxpreprocess/toppush.cpp index abf5722ddf..738a0afb9a 100644 --- a/src/gromacs/gmxpreprocess/toppush.cpp +++ b/src/gromacs/gmxpreprocess/toppush.cpp @@ -1437,12 +1437,13 @@ void push_atom(t_symtab *symtab, t_block *cgs, typeB == type ? ctype : ctypeB, mB, qB, wi); } -void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line, - warninp *wi) +void push_molt(t_symtab *symtab, + std::vector *mol, + char *line, + warninp *wi) { char type[STRLEN]; - int nrexcl, i; - t_molinfo *newmol; + int nrexcl; if ((sscanf(line, "%s%d", type, &nrexcl)) != 2) { @@ -1450,26 +1451,22 @@ void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line, } /* Test if this moleculetype overwrites another */ - i = 0; - while (i < *nmol) + const auto found = std::find_if(mol->begin(), mol->end(), + [&type](const auto &m) + { return strcmp(*(m.name), type) == 0; }); + if (found != mol->end()) { - if (strcmp(*((*mol)[i].name), type) == 0) - { - auto message = gmx::formatString("moleculetype %s is redefined", type); - warning_error_and_exit(wi, message, FARGS); - } - i++; + auto message = gmx::formatString("moleculetype %s is redefined", type); + warning_error_and_exit(wi, message, FARGS); } - (*nmol)++; - srenew(*mol, *nmol); - newmol = &((*mol)[*nmol-1]); - init_molinfo(newmol); + mol->emplace_back(); + mol->back().initMolInfo(); /* Fill in the values */ - newmol->name = put_symtab(symtab, type); - newmol->nrexcl = nrexcl; - newmol->excl_set = FALSE; + mol->back().name = put_symtab(symtab, type); + mol->back().nrexcl = nrexcl; + mol->back().excl_set = false; } static bool default_nb_params(int ftype, t_params bt[], t_atoms *at, @@ -2371,7 +2368,7 @@ void push_vsitesn(Directive d, t_params bond[], sfree(weight); } -void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol, +void push_mol(gmx::ArrayRef mols, char *pline, int *whichmol, int *nrcopies, warninp *wi) { @@ -2393,18 +2390,20 @@ void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol, int nrci = 0; int matchci = -1; int matchcs = -1; - for (int i = 0; i < nrmols; i++) + int i = 0; + for (const auto &mol : mols) { - if (strcmp(type, *(mols[i].name)) == 0) + if (strcmp(type, *(mol.name)) == 0) { nrcs++; matchcs = i; } - if (gmx_strcasecmp(type, *(mols[i].name)) == 0) + if (gmx_strcasecmp(type, *(mol.name)) == 0) { nrci++; matchci = i; } + i++; } if (nrcs == 1) @@ -2576,7 +2575,7 @@ static void convert_pairs_to_pairsQ(t_params *plist, plist[F_LJ14].param = nullptr; } -static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp, warninp *wi) +static void generate_LJCpairsNB(MoleculeInformation *mol, int nb_funct, t_params *nbp, warninp *wi) { int n, ntype, i, j, k; t_atom *atom; @@ -2696,7 +2695,7 @@ static void decouple_atoms(t_atoms *atoms, int atomtype_decouple, } } -void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ, +void convert_moltype_couple(MoleculeInformation *mol, int atomtype_decouple, real fudgeQQ, int couple_lam0, int couple_lam1, bool bCoupleIntra, int nb_funct, t_params *nbp, warninp *wi) diff --git a/src/gromacs/gmxpreprocess/toppush.h b/src/gromacs/gmxpreprocess/toppush.h index ee76fe3f00..aa8d8b9bdb 100644 --- a/src/gromacs/gmxpreprocess/toppush.h +++ b/src/gromacs/gmxpreprocess/toppush.h @@ -38,6 +38,9 @@ #ifndef GMX_GMXPREPROCESS_TOPPUSH_H #define GMX_GMXPREPROCESS_TOPPUSH_H +#include + +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/real.h" enum class Directive : int; @@ -45,7 +48,7 @@ struct gpp_atomtype; struct gpp_bond_atomtype; struct t_atoms; struct t_block; -struct t_molinfo; +struct MoleculeInformation; struct t_nbparam; struct t_param; struct t_params; @@ -104,11 +107,11 @@ void push_vsitesn(Directive d, t_params bond[], t_atoms *at, char *line, warninp *wi); -void push_mol(int nrmols, t_molinfo mols[], char *pline, +void push_mol(gmx::ArrayRef mols, char *pline, int *whichmol, int *nrcopies, warninp *wi); -void push_molt(struct t_symtab *symtab, int *nmol, t_molinfo **mol, char *line, +void push_molt(struct t_symtab *symtab, std::vector *mol, char *line, warninp *wi); void push_excl(char *line, gmx::ExclusionBlocks *b2, warninp *wi); @@ -123,7 +126,7 @@ int add_atomtype_decoupled(struct t_symtab *symtab, gpp_atomtype *at, * Returns the atom type number. */ -void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, +void convert_moltype_couple(MoleculeInformation *mol, int atomtype_decouple, real fudgeQQ, int couple_lam0, int couple_lam1, bool bCoupleIntra, diff --git a/src/gromacs/gmxpreprocess/toputil.cpp b/src/gromacs/gmxpreprocess/toputil.cpp index c5a69f6014..ae1a78e340 100644 --- a/src/gromacs/gmxpreprocess/toputil.cpp +++ b/src/gromacs/gmxpreprocess/toputil.cpp @@ -178,29 +178,6 @@ void add_param_to_list(t_params *list, t_param *b) list->nr++; } - -void init_molinfo(t_molinfo *mol) -{ - mol->nrexcl = 0; - mol->excl_set = FALSE; - mol->bProcessed = FALSE; - init_plist(mol->plist); - init_block(&mol->cgs); - init_block(&mol->mols); - init_blocka(&mol->excls); - init_t_atoms(&mol->atoms, 0, FALSE); -} - -/* FREEING MEMORY */ - -void done_mi(t_molinfo *mi) -{ - done_atom (&(mi->atoms)); - done_block(&(mi->cgs)); - done_block(&(mi->mols)); - done_plist(mi->plist); -} - /* PRINTING STRUCTURES */ static void print_bt(FILE *out, Directive d, gpp_atomtype *at, diff --git a/src/gromacs/gmxpreprocess/toputil.h b/src/gromacs/gmxpreprocess/toputil.h index e581c81272..5664d606b1 100644 --- a/src/gromacs/gmxpreprocess/toputil.h +++ b/src/gromacs/gmxpreprocess/toputil.h @@ -45,7 +45,7 @@ struct gpp_atomtype; struct t_atoms; struct t_blocka; struct t_excls; -struct t_molinfo; +struct MoleculeInformation; struct t_param; struct t_params; @@ -66,10 +66,6 @@ void add_param_to_list(t_params *list, t_param *b); void init_plist(t_params plist[]); void done_plist(t_params *plist); -void init_molinfo(t_molinfo *mol); - -/* FREE */ -void done_mi(t_molinfo *mi); /* PRINTING */ -- 2.22.0