From f51890ea2850ab56e5b6c8f73959b6b756ad6b29 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Thu, 13 Sep 2018 01:29:27 +0200 Subject: [PATCH] Add C++ version of t_ilist gmx_moltype_t in gmx_mtop_t now uses a C++ version of t_ilist called InteractionList. The interface of IteractionList and t_ilist match such that templated code can work with either. Thus we do not need to change all code that uses t_idef now as well. Change-Id: I90b829dae874f38d52ec0ceb538342313b2f8dd9 --- src/gromacs/domdec/domdec.cpp | 8 +- src/gromacs/domdec/domdec_constraints.cpp | 53 ++++----- src/gromacs/domdec/domdec_topology.cpp | 62 ++++++---- src/gromacs/fileio/tngio.cpp | 88 ++++++-------- src/gromacs/fileio/tpxio.cpp | 51 ++++---- src/gromacs/gmxana/gmx_nmeig.cpp | 2 +- src/gromacs/gmxpreprocess/convparm.cpp | 32 ++--- src/gromacs/gmxpreprocess/grompp.cpp | 34 +++--- src/gromacs/gmxpreprocess/readir.cpp | 39 +++--- src/gromacs/gmxpreprocess/topio.cpp | 36 +++--- src/gromacs/gmxpreprocess/vsite_parm.cpp | 15 +-- src/gromacs/listed-forces/disre.cpp | 18 +-- src/gromacs/listed-forces/orires.cpp | 14 +-- src/gromacs/mdlib/broadcaststructs.cpp | 23 ++-- src/gromacs/mdlib/calc_verletbuf.cpp | 14 +-- src/gromacs/mdlib/constr.cpp | 139 +++++++++++++--------- src/gromacs/mdlib/constr.h | 35 +++++- src/gromacs/mdlib/constraintrange.cpp | 20 ++-- src/gromacs/mdlib/forcerec.cpp | 13 +- src/gromacs/mdlib/lincs.cpp | 56 ++++----- src/gromacs/mdlib/membed.cpp | 4 +- src/gromacs/mdlib/perf_est.cpp | 2 +- src/gromacs/mdlib/settle.cpp | 12 +- src/gromacs/mdlib/shellfc.cpp | 5 +- src/gromacs/mdlib/sim_util.cpp | 4 +- src/gromacs/mdlib/tests/settle.cpp | 19 ++- src/gromacs/mdlib/vsite.cpp | 45 ++++--- src/gromacs/pbcutil/mshift.cpp | 80 +++++++------ src/gromacs/pbcutil/mshift.h | 11 +- src/gromacs/tools/convert_tpr.cpp | 7 +- src/gromacs/topology/idef.cpp | 39 +++--- src/gromacs/topology/idef.h | 50 +++++++- src/gromacs/topology/mtop_util.cpp | 50 ++++---- src/gromacs/topology/mtop_util.h | 13 +- src/gromacs/topology/topology.cpp | 16 +-- src/gromacs/topology/topology.h | 34 +++--- src/gromacs/topology/topsort.cpp | 9 +- 37 files changed, 617 insertions(+), 535 deletions(-) diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index d7ab17c0fa..3e1bd81c8c 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -3226,9 +3226,9 @@ static real *get_slb_frac(const gmx::MDLogger &mdlog, static int multi_body_bondeds_count(const gmx_mtop_t *mtop) { - int n, nmol, ftype; - gmx_mtop_ilistloop_t iloop; - const t_ilist *il; + int n, nmol, ftype; + gmx_mtop_ilistloop_t iloop; + const InteractionList *il; n = 0; iloop = gmx_mtop_ilistloop_init(mtop); @@ -3239,7 +3239,7 @@ static int multi_body_bondeds_count(const gmx_mtop_t *mtop) if ((interaction_function[ftype].flags & IF_BOND) && NRAL(ftype) > 2) { - n += nmol*il[ftype].nr/(1 + NRAL(ftype)); + n += nmol*il[ftype].size()/(1 + NRAL(ftype)); } } } diff --git a/src/gromacs/domdec/domdec_constraints.cpp b/src/gromacs/domdec/domdec_constraints.cpp index ef4ef0765f..4e2d634dbf 100644 --- a/src/gromacs/domdec/domdec_constraints.cpp +++ b/src/gromacs/domdec/domdec_constraints.cpp @@ -135,7 +135,8 @@ void dd_clear_local_constraint_indices(gmx_domdec_t *dd) /*! \brief Walks over the constraints out from the local atoms into the non-local atoms and adds them to a list */ static void walk_out(int con, int con_offset, int a, int offset, int nrec, - int ncon1, const t_iatom *ia1, const t_iatom *ia2, + gmx::ArrayRef ia1, + gmx::ArrayRef ia2, const t_blocka *at2con, const gmx_ga2la_t &ga2la, gmx_bool bHomeConnect, gmx_domdec_constraints_t *dc, @@ -157,7 +158,7 @@ static void walk_out(int con, int con_offset, int a, int offset, int nrec, il_local->nalloc = over_alloc_dd(il_local->nr+3); srenew(il_local->iatoms, il_local->nalloc); } - iap = constr_iatomptr(ncon1, ia1, ia2, con); + iap = constr_iatomptr(ia1, ia2, con); il_local->iatoms[il_local->nr++] = iap[0]; a1_gl = offset + iap[1]; a2_gl = offset + iap[2]; @@ -200,7 +201,7 @@ static void walk_out(int con, int con_offset, int a, int offset, int nrec, if (coni != con) { /* Walk further */ - iap = constr_iatomptr(ncon1, ia1, ia2, coni); + iap = constr_iatomptr(ia1, ia2, coni); if (a == iap[1]) { b = iap[2]; @@ -212,7 +213,7 @@ static void walk_out(int con, int con_offset, int a, int offset, int nrec, if (!ga2la.findHome(offset + b)) { walk_out(coni, con_offset, b, offset, nrec-1, - ncon1, ia1, ia2, at2con, + ia1, ia2, at2con, ga2la, FALSE, dc, dcc, il_local, ireq); } } @@ -224,7 +225,7 @@ static void walk_out(int con, int con_offset, int a, int offset, int nrec, static void atoms_to_settles(gmx_domdec_t *dd, const gmx_mtop_t *mtop, const int *cginfo, - const int *const*at2settle_mt, + gmx::ArrayRef < const std::vector < int>> at2settle_mt, int cg_start, int cg_end, t_ilist *ils_local, std::vector *ireq) @@ -248,13 +249,13 @@ static void atoms_to_settles(gmx_domdec_t *dd, if (settle >= 0) { - int offset = a_gl - a_mol; + int offset = a_gl - a_mol; - t_iatom *ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms; + const int *ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms.data(); - int a_gls[3]; - gmx_bool bAssign = FALSE; - int nlocal = 0; + int a_gls[3]; + gmx_bool bAssign = FALSE; + int nlocal = 0; for (int sa = 0; sa < nral; sa++) { int a_glsa = offset + ia1[settle*(1+nral)+1+sa]; @@ -312,8 +313,6 @@ static void atoms_to_constraints(gmx_domdec_t *dd, std::vector *ireq) { const t_blocka *at2con; - int ncon1; - t_iatom *ia1, *ia2, *iap; int b_lo, offset, b_mol, i, con, con_offset; gmx_domdec_constraints_t *dc = dd->constraints; @@ -336,12 +335,10 @@ static void atoms_to_constraints(gmx_domdec_t *dd, int molnr, a_mol; mtopGetMolblockIndex(mtop, a_gl, &mb, &molnr, &a_mol); - const gmx_molblock_t *molb = &mtop->molblock[mb]; + const gmx_molblock_t &molb = mtop->molblock[mb]; - ncon1 = mtop->moltype[molb->type].ilist[F_CONSTR].nr/NRAL(F_SETTLE); - - ia1 = mtop->moltype[molb->type].ilist[F_CONSTR].iatoms; - ia2 = mtop->moltype[molb->type].ilist[F_CONSTRNC].iatoms; + gmx::ArrayRef ia1 = mtop->moltype[molb.type].ilist[F_CONSTR].iatoms; + gmx::ArrayRef ia2 = mtop->moltype[molb.type].ilist[F_CONSTRNC].iatoms; /* Calculate the global constraint number offset for the molecule. * This is only required for the global index to make sure @@ -352,11 +349,11 @@ static void atoms_to_constraints(gmx_domdec_t *dd, /* The global atom number offset for this molecule */ offset = a_gl - a_mol; - at2con = &at2con_mt[molb->type]; + at2con = &at2con_mt[molb.type]; for (i = at2con->index[a_mol]; i < at2con->index[a_mol+1]; i++) { - con = at2con->a[i]; - iap = constr_iatomptr(ncon1, ia1, ia2, con); + con = at2con->a[i]; + const int *iap = constr_iatomptr(ia1, ia2, con); if (a_mol == iap[1]) { b_mol = iap[2]; @@ -394,7 +391,7 @@ static void atoms_to_constraints(gmx_domdec_t *dd, * after this first call. */ walk_out(con, con_offset, b_mol, offset, nrec, - ncon1, ia1, ia2, at2con, + ia1, ia2, at2con, ga2la, TRUE, dc, dcc, ilc_local, ireq); } } @@ -424,7 +421,6 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start, t_ilist *ilc_local, *ils_local; std::vector *ireq; gmx::ArrayRef at2con_mt; - const int *const*at2settle_mt; gmx::HashedMap *ga2la_specat; int at_end, i, j; t_iatom *iap; @@ -463,6 +459,8 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start, ireq = nullptr; } + gmx::ArrayRef < const std::vector < int>> at2settle_mt; + /* When settle works inside charge groups, we assigned them already */ if (dd->bInterCGsettles) { // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr? @@ -470,13 +468,8 @@ int dd_make_local_constraints(gmx_domdec_t *dd, int at_start, at2settle_mt = constr->atom2settle_moltype(); ils_local->nr = 0; } - else - { - /* Settle works inside charge groups, we assigned them already */ - at2settle_mt = nullptr; - } - if (at2settle_mt == nullptr) + if (at2settle_mt.empty()) { atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec, ilc_local, ireq); @@ -643,8 +636,8 @@ void init_domdec_constraints(gmx_domdec_t *dd, molb = &mtop->molblock[mb]; dc->molb_con_offset[mb] = ncon; dc->molb_ncon_mol[mb] = - mtop->moltype[molb->type].ilist[F_CONSTR].nr/3 + - mtop->moltype[molb->type].ilist[F_CONSTRNC].nr/3; + mtop->moltype[molb->type].ilist[F_CONSTR].size()/3 + + mtop->moltype[molb->type].ilist[F_CONSTRNC].size()/3; ncon += molb->nmol*dc->molb_ncon_mol[mb]; } diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp index d0cde4d366..55498e12c4 100644 --- a/src/gromacs/domdec/domdec_topology.cpp +++ b/src/gromacs/domdec/domdec_topology.cpp @@ -493,7 +493,8 @@ static void count_excls(const t_block *cgs, const t_blocka *excls, } /*! \brief Run the reverse ilist generation and store it in r_il when \p bAssign = TRUE */ -static int low_make_reverse_ilist(const t_ilist *il_mt, const t_atom *atom, +static int low_make_reverse_ilist(const InteractionList *il_mt, + const t_atom *atom, gmx::ArrayRef < const std::vector < int>> vsitePbc, int *count, gmx_bool bConstr, gmx_bool bSettle, @@ -503,12 +504,9 @@ static int low_make_reverse_ilist(const t_ilist *il_mt, const t_atom *atom, gmx_bool bLinkToAllAtoms, gmx_bool bAssign) { - int ftype, nral, i, j, nlink, link; - const t_ilist *il; - const t_iatom *ia; + int ftype, j, nlink, link; int a; int nint; - gmx_bool bVSite; nint = 0; for (ftype = 0; ftype < F_NRE; ftype++) @@ -517,12 +515,12 @@ static int low_make_reverse_ilist(const t_ilist *il_mt, const t_atom *atom, (bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) || (bSettle && ftype == F_SETTLE)) { - bVSite = ((interaction_function[ftype].flags & IF_VSITE) != 0u); - nral = NRAL(ftype); - il = &il_mt[ftype]; - for (i = 0; i < il->nr; i += 1+nral) + const bool bVSite = ((interaction_function[ftype].flags & IF_VSITE) != 0u); + const int nral = NRAL(ftype); + const auto &il = il_mt[ftype]; + for (int i = 0; i < il.size(); i += 1+nral) { - ia = il->iatoms + i; + const int* ia = il.iatoms.data() + i; if (bLinkToAllAtoms) { if (bVSite) @@ -599,7 +597,7 @@ static int low_make_reverse_ilist(const t_ilist *il_mt, const t_atom *atom, } /*! \brief Make the reverse ilist: a list of bonded interactions linked to atoms */ -static int make_reverse_ilist(const t_ilist *ilist, +static int make_reverse_ilist(const InteractionList *ilist, const t_atoms *atoms, gmx::ArrayRef < const std::vector < int>> vsitePbc, gmx_bool bConstr, gmx_bool bSettle, @@ -700,8 +698,11 @@ static gmx_reverse_top_t make_reverse_top(const gmx_mtop_t *mtop, gmx_bool bFE, atoms_global.nr = mtop->natoms; atoms_global.atom = nullptr; /* Only used with virtual sites */ + GMX_RELEASE_ASSERT(mtop->intermolecular_ilist.get(), "We should have an ilist when intermolecular interactions are on"); + *nint += - make_reverse_ilist(mtop->intermolecular_ilist, &atoms_global, + make_reverse_ilist(mtop->intermolecular_ilist->data(), + &atoms_global, gmx::EmptyArrayRef(), rt.bConstr, rt.bSettle, rt.bBCheck, FALSE, &rt.ril_intermol); @@ -2277,7 +2278,10 @@ t_blocka *make_charge_group_links(const gmx_mtop_t *mtop, gmx_domdec_t *dd, atoms.nr = mtop->natoms; atoms.atom = nullptr; - make_reverse_ilist(mtop->intermolecular_ilist, &atoms, + GMX_RELEASE_ASSERT(mtop->intermolecular_ilist.get(), "We should have an ilist when intermolecular interactions are on"); + + make_reverse_ilist(mtop->intermolecular_ilist->data(), + &atoms, gmx::EmptyArrayRef(), FALSE, FALSE, FALSE, TRUE, &ril_intermol); } @@ -2454,11 +2458,11 @@ static void bonded_cg_distance_mol(const gmx_moltype_t *molt, { if (dd_check_ftype(ftype, bBCheck, FALSE, FALSE)) { - const t_ilist *il = &molt->ilist[ftype]; + const auto *il = &molt->ilist[ftype]; int nral = NRAL(ftype); if (nral > 1) { - for (int i = 0; i < il->nr; i += 1+nral) + for (int i = 0; i < il->size(); i += 1+nral) { for (int ai = 0; ai < nral; ai++) { @@ -2503,7 +2507,7 @@ static void bonded_cg_distance_mol(const gmx_moltype_t *molt, } /*! \brief Set the distance, function type and atom indices for the longest atom distance involved in intermolecular interactions for two-body and multi-body bonded interactions */ -static void bonded_distance_intermol(const t_ilist *ilists_intermol, +static void bonded_distance_intermol(const InteractionList *ilists_intermol, gmx_bool bBCheck, const rvec *x, int ePBC, const matrix box, bonded_distance_t *bd_2b, @@ -2517,13 +2521,13 @@ static void bonded_distance_intermol(const t_ilist *ilists_intermol, { if (dd_check_ftype(ftype, bBCheck, FALSE, FALSE)) { - const t_ilist *il = &ilists_intermol[ftype]; + const auto *il = &ilists_intermol[ftype]; int nral = NRAL(ftype); /* No nral>1 check here, since intermol interactions always * have nral>=2 (and the code is also correct for nral=1). */ - for (int i = 0; i < il->nr; i += 1+nral) + for (int i = 0; i < il->size(); i += 1+nral) { for (int ai = 0; ai < nral; ai++) { @@ -2557,7 +2561,7 @@ static bool moltypeHasVsite(const gmx_moltype_t &molt) for (int i = 0; i < F_NRE; i++) { if ((interaction_function[i].flags & IF_VSITE) && - molt.ilist[i].nr > 0) + molt.ilist[i].size() > 0) { hasVsite = true; } @@ -2601,8 +2605,19 @@ static void get_cgcm_mol(const gmx_moltype_t *molt, if (moltypeHasVsite(*molt)) { + /* Convert to old, deprecated format */ + t_ilist ilist[F_NRE]; + for (int ftype = 0; ftype < F_NRE; ftype++) + { + if (interaction_function[ftype].flags & IF_VSITE) + { + ilist[ftype].nr = molt->ilist[ftype].size(); + ilist[ftype].iatoms = const_cast(molt->ilist[ftype].iatoms.data()); + } + } + construct_vsites(nullptr, xs, 0.0, nullptr, - ffparams->iparams, molt->ilist, + ffparams->iparams, ilist, epbcNONE, TRUE, nullptr, nullptr); } @@ -2639,8 +2654,7 @@ void dd_bonded_cg_distance(const gmx::MDLogger &mdlog, { if (ir->ePBC != epbcNONE) { - mk_graph_ilist(nullptr, molt.ilist, 0, molt.atoms.nr, FALSE, FALSE, - &graph); + mk_graph_moltype(molt, &graph); } std::vector at2cg = make_at2cg(molt.cgs); @@ -2685,7 +2699,9 @@ void dd_bonded_cg_distance(const gmx::MDLogger &mdlog, gmx_fatal(FARGS, "The combination of intermolecular interactions, charge groups and domain decomposition is not supported. Use cutoff-scheme=Verlet (which removes the charge groups) or run without domain decomposition."); } - bonded_distance_intermol(mtop->intermolecular_ilist, + GMX_RELEASE_ASSERT(mtop->intermolecular_ilist.get(), "We should have an ilist when intermolecular interactions are on"); + + bonded_distance_intermol(mtop->intermolecular_ilist->data(), bBCheck, x, ir->ePBC, box, &bd_2b, &bd_mb); diff --git a/src/gromacs/fileio/tngio.cpp b/src/gromacs/fileio/tngio.cpp index f6fd09c6ef..2fa0e2d898 100644 --- a/src/gromacs/fileio/tngio.cpp +++ b/src/gromacs/fileio/tngio.cpp @@ -293,7 +293,6 @@ void gmx_tng_add_mtop(gmx_tng_trajectory_t gmx_tng, int j; std::vector atomCharges; std::vector atomMasses; - const t_ilist *ilist; tng_bond_t tngBond; char datatype; @@ -334,29 +333,23 @@ void gmx_tng_add_mtop(gmx_tng_trajectory_t gmx_tng, { if (IS_CHEMBOND(i)) { - ilist = &molType->ilist[i]; - if (ilist) + const InteractionList &ilist = molType->ilist[i]; + j = 1; + while (j < ilist.size()) { - j = 1; - while (j < ilist->nr) - { - tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+1], &tngBond); - j += 3; - } + tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 1], &tngBond); + j += 3; } } } /* Settle is described using three atoms */ - ilist = &molType->ilist[F_SETTLE]; - if (ilist) + const InteractionList &ilist = molType->ilist[F_SETTLE]; + j = 1; + while (j < ilist.size()) { - j = 1; - while (j < ilist->nr) - { - tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+1], &tngBond); - tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+2], &tngBond); - j += 4; - } + tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 1], &tngBond); + tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 2], &tngBond); + j += 4; } /* First copy atom charges and masses, first atom by atom and then copy the memory for the molecule instances. * FIXME: Atom B state data should also be written to TNG (v 2.0?) */ @@ -618,7 +611,6 @@ static void add_selection_groups(gmx_tng_trajectory_t gmx_tng, const t_atoms *atoms; const t_atom *at; const t_resinfo *resInfo; - const t_ilist *ilist; int nameIndex; int atom_offset = 0; tng_molecule_t mol, iterMol; @@ -709,51 +701,41 @@ static void add_selection_groups(gmx_tng_trajectory_t gmx_tng, { if (IS_CHEMBOND(k)) { - ilist = &molType.ilist[k]; - if (ilist) + const InteractionList &ilist = molType.ilist[k]; + for (int l = 1; l < ilist.size(); l += 3) { - int l = 1; - while (l < ilist->nr) + int atom1, atom2; + atom1 = ilist.iatoms[l] + atom_offset; + atom2 = ilist.iatoms[l + 1] + atom_offset; + if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0 && + getGroupType(&mtop->groups, egcCompressedX, atom2) == 0) { - int atom1, atom2; - atom1 = ilist->iatoms[l] + atom_offset; - atom2 = ilist->iatoms[l+1] + atom_offset; - if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0 && - getGroupType(&mtop->groups, egcCompressedX, atom2) == 0) - { - tng_molecule_bond_add(tng, mol, ilist->iatoms[l], - ilist->iatoms[l+1], &tngBond); - } - l += 3; + tng_molecule_bond_add(tng, mol, ilist.iatoms[l], + ilist.iatoms[l + 1], &tngBond); } } } } /* Settle is described using three atoms */ - ilist = &molType.ilist[F_SETTLE]; - if (ilist) + const InteractionList &ilist = molType.ilist[F_SETTLE]; + for (int l = 1; l < ilist.size(); l += 4) { - int l = 1; - while (l < ilist->nr) + int atom1, atom2, atom3; + atom1 = ilist.iatoms[l] + atom_offset; + atom2 = ilist.iatoms[l + 1] + atom_offset; + atom3 = ilist.iatoms[l + 2] + atom_offset; + if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0) { - int atom1, atom2, atom3; - atom1 = ilist->iatoms[l] + atom_offset; - atom2 = ilist->iatoms[l+1] + atom_offset; - atom3 = ilist->iatoms[l+2] + atom_offset; - if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0) + if (getGroupType(&mtop->groups, egcCompressedX, atom2) == 0) { - if (getGroupType(&mtop->groups, egcCompressedX, atom2) == 0) - { - tng_molecule_bond_add(tng, mol, atom1, - atom2, &tngBond); - } - if (getGroupType(&mtop->groups, egcCompressedX, atom3) == 0) - { - tng_molecule_bond_add(tng, mol, atom1, - atom3, &tngBond); - } + tng_molecule_bond_add(tng, mol, atom1, + atom2, &tngBond); + } + if (getGroupType(&mtop->groups, egcCompressedX, atom3) == 0) + { + tng_molecule_bond_add(tng, mol, atom1, + atom3, &tngBond); } - l += 4; } } } diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index 3c79f26038..fa44d4ed05 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -47,6 +47,7 @@ #include #include +#include "gromacs/compat/make_unique.h" #include "gromacs/fileio/filetypes.h" #include "gromacs/fileio/gmxfio.h" #include "gromacs/fileio/gmxfio-xdr.h" @@ -1963,14 +1964,15 @@ static void do_iparams(t_fileio *fio, t_functype ftype, t_iparams *iparams, } } -static void do_ilist(t_fileio *fio, t_ilist *ilist, gmx_bool bRead) +static void do_ilist(t_fileio *fio, InteractionList *ilist, gmx_bool bRead) { - gmx_fio_do_int(fio, ilist->nr); + int nr = ilist->size(); + gmx_fio_do_int(fio, nr); if (bRead) { - snew(ilist->iatoms, ilist->nr); + ilist->iatoms.resize(nr); } - gmx_fio_ndo_int(fio, ilist->iatoms, ilist->nr); + gmx_fio_ndo_int(fio, ilist->iatoms.data(), ilist->size()); } static void do_ffparams(t_fileio *fio, gmx_ffparams_t *ffparams, @@ -2022,23 +2024,22 @@ static void do_ffparams(t_fileio *fio, gmx_ffparams_t *ffparams, } } -static void add_settle_atoms(t_ilist *ilist) +static void add_settle_atoms(InteractionList *ilist) { int i; /* Settle used to only store the first atom: add the other two */ - srenew(ilist->iatoms, 2*ilist->nr); - for (i = ilist->nr/2-1; i >= 0; i--) + ilist->iatoms.resize(2*ilist->size()); + for (i = ilist->size()/4 - 1; i >= 0; i--) { ilist->iatoms[4*i+0] = ilist->iatoms[2*i+0]; ilist->iatoms[4*i+1] = ilist->iatoms[2*i+1]; ilist->iatoms[4*i+2] = ilist->iatoms[2*i+1] + 1; ilist->iatoms[4*i+3] = ilist->iatoms[2*i+1] + 2; } - ilist->nr = 2*ilist->nr; } -static void do_ilists(t_fileio *fio, t_ilist *ilist, gmx_bool bRead, +static void do_ilists(t_fileio *fio, InteractionList *ilist, gmx_bool bRead, int file_version) { int j; @@ -2059,13 +2060,12 @@ static void do_ilists(t_fileio *fio, t_ilist *ilist, gmx_bool bRead, } if (bClear) { - ilist[j].nr = 0; - ilist[j].iatoms = nullptr; + ilist[j].iatoms.clear(); } else { do_ilist(fio, &ilist[j], bRead); - if (file_version < 78 && j == F_SETTLE && ilist[j].nr > 0) + if (file_version < 78 && j == F_SETTLE && ilist[j].size() > 0) { add_settle_atoms(&ilist[j]); } @@ -2454,27 +2454,26 @@ static void do_molblock(t_fileio *fio, gmx_molblock_t *molb, static void set_disres_npair(gmx_mtop_t *mtop) { - t_iparams *ip; - gmx_mtop_ilistloop_t iloop; - const t_ilist *ilist, *il; - int nmol, i, npair; - t_iatom *a; + t_iparams *ip; + gmx_mtop_ilistloop_t iloop; + const InteractionList *ilist; + int nmol; ip = mtop->ffparams.iparams; iloop = gmx_mtop_ilistloop_init(mtop); while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol)) { - il = &ilist[F_DISRES]; + const InteractionList &il = ilist[F_DISRES]; - if (il->nr > 0) + if (il.size() > 0) { - a = il->iatoms; - npair = 0; - for (i = 0; i < il->nr; i += 3) + gmx::ArrayRef a = il.iatoms; + int npair = 0; + for (int i = 0; i < il.size(); i += 3) { npair++; - if (i+3 == il->nr || ip[a[i]].disres.label != ip[a[i+3]].disres.label) + if (i+3 == il.size() || ip[a[i]].disres.label != ip[a[i+3]].disres.label) { ip[a[i]].disres.npair = npair; npair = 0; @@ -2528,9 +2527,9 @@ static void do_mtop(t_fileio *fio, gmx_mtop_t *mtop, gmx_bool bRead, { if (bRead) { - snew(mtop->intermolecular_ilist, F_NRE); + mtop->intermolecular_ilist = gmx::compat::make_unique(); } - do_ilists(fio, mtop->intermolecular_ilist, bRead, file_version); + do_ilists(fio, mtop->intermolecular_ilist->data(), bRead, file_version); } } else @@ -2907,7 +2906,7 @@ static int do_tpx(t_fileio *fio, gmx_bool bRead, { if (fileVersion < 57) { - if (mtop->moltype[0].ilist[F_DISRES].nr > 0) + if (mtop->moltype[0].ilist[F_DISRES].size() > 0) { ir->eDisre = edrSimple; } diff --git a/src/gromacs/gmxana/gmx_nmeig.cpp b/src/gromacs/gmxana/gmx_nmeig.cpp index fb2f9bf262..d5533e4665 100644 --- a/src/gromacs/gmxana/gmx_nmeig.cpp +++ b/src/gromacs/gmxana/gmx_nmeig.cpp @@ -107,7 +107,7 @@ static size_t get_nharm_mt(const gmx_moltype_t *mt) for (i = 0; (i < asize(harm_func)); i++) { ft = harm_func[i]; - nh += mt->ilist[ft].nr/(interaction_function[ft].nratoms+1); + nh += mt->ilist[ft].size()/(interaction_function[ft].nratoms+1); } return nh; } diff --git a/src/gromacs/gmxpreprocess/convparm.cpp b/src/gromacs/gmxpreprocess/convparm.cpp index 525aa064e6..a5f3f05ccc 100644 --- a/src/gromacs/gmxpreprocess/convparm.cpp +++ b/src/gromacs/gmxpreprocess/convparm.cpp @@ -43,6 +43,7 @@ #include #include +#include "gromacs/compat/make_unique.h" #include "gromacs/gmxpreprocess/gpp_atomtype.h" #include "gromacs/gmxpreprocess/topio.h" #include "gromacs/gmxpreprocess/toputil.h" @@ -482,27 +483,22 @@ static int enter_params(gmx_ffparams_t *ffparams, t_functype ftype, return type; } -static void append_interaction(t_ilist *ilist, +static void append_interaction(InteractionList *ilist, int type, int nral, const int a[MAXATOMLIST]) { - int i, where1; - - where1 = ilist->nr; - ilist->nr += nral+1; - - ilist->iatoms[where1++] = type; - for (i = 0; (i < nral); i++) + ilist->iatoms.push_back(type); + for (int i = 0; (i < nral); i++) { - ilist->iatoms[where1++] = a[i]; + ilist->iatoms.push_back(a[i]); } } static void enter_function(t_params *p, t_functype ftype, int comb, real reppow, - gmx_ffparams_t *ffparams, t_ilist *il, + gmx_ffparams_t *ffparams, InteractionList *il, int *maxtypes, bool bNB, bool bAppend) { - int k, type, nr, nral, delta, start; + int k, type, nr, nral, start; start = ffparams->ntypes; nr = p->nr; @@ -521,8 +517,6 @@ static void enter_function(t_params *p, t_functype ftype, int comb, real reppow, { assert(il); nral = NRAL(ftype); - delta = nr*(nral+1); - srenew(il->iatoms, il->nr+delta); append_interaction(il, type, nral, p->param[k].a); } } @@ -558,8 +552,7 @@ void convert_params(int atnr, t_params nbtypes[], molt = &mtop->moltype[mt]; for (i = 0; (i < F_NRE); i++) { - molt->ilist[i].nr = 0; - molt->ilist[i].iatoms = nullptr; + molt->ilist[i].iatoms.clear(); plist = mi[mt].plist; @@ -579,12 +572,11 @@ void convert_params(int atnr, t_params nbtypes[], if (intermolecular_interactions != nullptr) { /* Process the intermolecular interaction list */ - snew(mtop->intermolecular_ilist, F_NRE); + mtop->intermolecular_ilist = gmx::compat::make_unique(); for (i = 0; (i < F_NRE); i++) { - mtop->intermolecular_ilist[i].nr = 0; - mtop->intermolecular_ilist[i].iatoms = nullptr; + (*mtop->intermolecular_ilist)[i].iatoms.clear(); plist = intermolecular_interactions->plist; @@ -611,7 +603,7 @@ void convert_params(int atnr, t_params nbtypes[], else { enter_function(&(plist[i]), static_cast(i), comb, reppow, - ffp, &mtop->intermolecular_ilist[i], + ffp, &(*mtop->intermolecular_ilist)[i], &maxtypes, FALSE, FALSE); mtop->bIntermolecularInteractions = TRUE; @@ -621,7 +613,7 @@ void convert_params(int atnr, t_params nbtypes[], if (!mtop->bIntermolecularInteractions) { - sfree(mtop->intermolecular_ilist); + mtop->intermolecular_ilist.reset(nullptr); } } diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 338de9a82b..bd036ef254 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -264,10 +264,10 @@ static void check_bonds_timestep(const gmx_mtop_t *mtop, double dt, warninp_t wi const gmx_moltype_t *w_moltype = nullptr; for (const gmx_moltype_t &moltype : mtop->moltype) { - const t_atom *atom = moltype.atoms.atom; - const t_ilist *ilist = moltype.ilist; - const t_ilist *ilc = &ilist[F_CONSTR]; - const t_ilist *ils = &ilist[F_SETTLE]; + const t_atom *atom = moltype.atoms.atom; + const InteractionList *ilist = moltype.ilist; + const InteractionList *ilc = &ilist[F_CONSTR]; + const InteractionList *ils = &ilist[F_SETTLE]; for (ftype = 0; ftype < F_NRE; ftype++) { if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC)) @@ -275,8 +275,8 @@ static void check_bonds_timestep(const gmx_mtop_t *mtop, double dt, warninp_t wi continue; } - const t_ilist *ilb = &ilist[ftype]; - for (i = 0; i < ilb->nr; i += 3) + const InteractionList *ilb = &ilist[ftype]; + for (i = 0; i < ilb->size(); i += 3) { fc = ip[ilb->iatoms[i]].harmonic.krA; re = ip[ilb->iatoms[i]].harmonic.rA; @@ -305,7 +305,7 @@ static void check_bonds_timestep(const gmx_mtop_t *mtop, double dt, warninp_t wi if (period2 < limit2) { bFound = false; - for (j = 0; j < ilc->nr; j += 3) + for (j = 0; j < ilc->size(); j += 3) { if ((ilc->iatoms[j+1] == a1 && ilc->iatoms[j+2] == a2) || (ilc->iatoms[j+1] == a2 && ilc->iatoms[j+2] == a1)) @@ -313,7 +313,7 @@ static void check_bonds_timestep(const gmx_mtop_t *mtop, double dt, warninp_t wi bFound = true; } } - for (j = 0; j < ils->nr; j += 4) + for (j = 0; j < ils->size(); j += 4) { if ((a1 == ils->iatoms[j+1] || a1 == ils->iatoms[j+2] || a1 == ils->iatoms[j+3]) && (a2 == ils->iatoms[j+1] || a2 == ils->iatoms[j+2] || a2 == ils->iatoms[j+3])) @@ -1318,10 +1318,10 @@ static void checkForUnboundAtoms(const gmx_moltype_t *molt, (interaction_function[ftype].flags & IF_CONSTRAINT) || ftype == F_SETTLE) { - const t_ilist *il = &molt->ilist[ftype]; - int nral = NRAL(ftype); + const InteractionList *il = &molt->ilist[ftype]; + const int nral = NRAL(ftype); - for (int i = 0; i < il->nr; i += 1 + nral) + for (int i = 0; i < il->size(); i += 1 + nral) { for (int j = 0; j < nral; j++) { @@ -1377,7 +1377,8 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t *molt, const t_iparams *iparams, real massFactorThreshold) { - if (molt->ilist[F_CONSTR].nr == 0 && molt->ilist[F_CONSTRNC].nr == 0) + if (molt->ilist[F_CONSTR].size() == 0 && + molt->ilist[F_CONSTRNC].size() == 0) { return false; } @@ -1385,8 +1386,7 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t *molt, const t_atom * atom = molt->atoms.atom; t_blocka atomToConstraints = - gmx::make_at2con(molt->atoms.nr, - molt->ilist, iparams, + gmx::make_at2con(*molt, iparams, gmx::FlexibleConstraintTreatment::Exclude); bool haveDecoupledMode = false; @@ -1394,9 +1394,9 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t *molt, { if (interaction_function[ftype].flags & IF_ATYPE) { - const int nral = NRAL(ftype); - const t_ilist *il = &molt->ilist[ftype]; - for (int i = 0; i < il->nr; i += 1 + nral) + const int nral = NRAL(ftype); + const InteractionList *il = &molt->ilist[ftype]; + for (int i = 0; i < il->size(); i += 1 + nral) { /* Here we check for the mass difference between the atoms * at both ends of the angle, that the atoms at the ends have diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp index 62885e2ebd..b27b707f1b 100644 --- a/src/gromacs/gmxpreprocess/readir.cpp +++ b/src/gromacs/gmxpreprocess/readir.cpp @@ -2760,7 +2760,6 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) const gmx_groups_t *groups; pull_params_t *pull; int natoms, ai, aj, i, j, d, g, imin, jmin; - t_iatom *ia; int *nrdf2, *na_vcm, na_tot; double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub; ivec *dof_vcm; @@ -2834,8 +2833,8 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) { for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++) { - ia = molt.ilist[ftype].iatoms; - for (i = 0; i < molt.ilist[ftype].nr; ) + gmx::ArrayRef ia = molt.ilist[ftype].iatoms; + for (i = 0; i < molt.ilist[ftype].size(); ) { /* Subtract degrees of freedom for the constraints, * if the particles still have degrees of freedom left. @@ -2844,12 +2843,12 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) * contribute to the constraints the degrees of freedom do not * change. */ - ai = as + ia[1]; - aj = as + ia[2]; - if (((atom[ia[1]].ptype == eptNucleus) || - (atom[ia[1]].ptype == eptAtom)) && - ((atom[ia[2]].ptype == eptNucleus) || - (atom[ia[2]].ptype == eptAtom))) + ai = as + ia[i + 1]; + aj = as + ia[i + 2]; + if (((atom[ia[i + 1]].ptype == eptNucleus) || + (atom[ia[i + 1]].ptype == eptAtom)) && + ((atom[ia[i + 2]].ptype == eptNucleus) || + (atom[ia[i + 2]].ptype == eptAtom))) { if (nrdf2[ai] > 0) { @@ -2876,23 +2875,21 @@ static void calc_nrdf(const gmx_mtop_t *mtop, t_inputrec *ir, char **gnames) nrdf_vcm[getGroupType(groups, egcVCM, ai)] -= 0.5*imin; nrdf_vcm[getGroupType(groups, egcVCM, aj)] -= 0.5*jmin; } - ia += interaction_function[ftype].nratoms+1; i += interaction_function[ftype].nratoms+1; } } - ia = molt.ilist[F_SETTLE].iatoms; - for (i = 0; i < molt.ilist[F_SETTLE].nr; ) + gmx::ArrayRef ia = molt.ilist[F_SETTLE].iatoms; + for (i = 0; i < molt.ilist[F_SETTLE].size(); ) { /* Subtract 1 dof from every atom in the SETTLE */ for (j = 0; j < 3; j++) { - ai = as + ia[1+j]; + ai = as + ia[i + 1 + j]; imin = std::min(2, nrdf2[ai]); nrdf2[ai] -= imin; nrdf_tc [getGroupType(groups, egcTC, ai)] -= 0.5*imin; nrdf_vcm[getGroupType(groups, egcVCM, ai)] -= 0.5*imin; } - ia += 4; i += 4; } as += molt.atoms.nr; @@ -3721,11 +3718,11 @@ static bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys, bool posres_only, ivec AbsRef) { - int d, g, i; - gmx_mtop_ilistloop_t iloop; - const t_ilist *ilist; - int nmol; - t_iparams *pr; + int d, g, i; + gmx_mtop_ilistloop_t iloop; + const InteractionList *ilist; + int nmol; + t_iparams *pr; clear_ivec(AbsRef); @@ -3756,7 +3753,7 @@ static bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys, if (nmol > 0 && (AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0)) { - for (i = 0; i < ilist[F_POSRES].nr; i += 2) + for (i = 0; i < ilist[F_POSRES].size(); i += 2) { pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]]; for (d = 0; d < DIM; d++) @@ -3767,7 +3764,7 @@ static bool absolute_reference(t_inputrec *ir, gmx_mtop_t *sys, } } } - for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2) + for (i = 0; i < ilist[F_FBPOSRES].size(); i += 2) { /* Check for flat-bottom posres */ pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]]; diff --git a/src/gromacs/gmxpreprocess/topio.cpp b/src/gromacs/gmxpreprocess/topio.cpp index 755e4f9bee..a43bb17147 100644 --- a/src/gromacs/gmxpreprocess/topio.cpp +++ b/src/gromacs/gmxpreprocess/topio.cpp @@ -1087,12 +1087,12 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr */ int ftype_connbond = 0; int ind_connbond = 0; - if (molt->ilist[F_CONNBONDS].nr != 0) + if (molt->ilist[F_CONNBONDS].size() != 0) { fprintf(stderr, "nr. of CONNBONDS present already: %d\n", - molt->ilist[F_CONNBONDS].nr/3); + molt->ilist[F_CONNBONDS].size()/3); ftype_connbond = molt->ilist[F_CONNBONDS].iatoms[0]; - ind_connbond = molt->ilist[F_CONNBONDS].nr; + ind_connbond = molt->ilist[F_CONNBONDS].size(); } /* now we delete all bonded interactions, except the ones describing * a chemical bond. These are converted to CONNBONDS @@ -1106,7 +1106,7 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr } int nratoms = interaction_function[ftype].nratoms; int j = 0; - while (j < molt->ilist[ftype].nr) + while (j < molt->ilist[ftype].size()) { bool bexcl; @@ -1124,11 +1124,11 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr */ if (bexcl && IS_CHEMBOND(ftype)) { - srenew(molt->ilist[F_CONNBONDS].iatoms, ind_connbond + 3); - molt->ilist[F_CONNBONDS].nr += 3; - molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = ftype_connbond; - molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = a1; - molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = a2; + InteractionList &ilist = molt->ilist[F_CONNBONDS]; + ilist.iatoms.resize(ind_connbond + 3); + ilist.iatoms[ind_connbond++] = ftype_connbond; + ilist.iatoms[ind_connbond++] = a1; + ilist.iatoms[ind_connbond++] = a2; } } else @@ -1160,11 +1160,12 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr /* since the interaction involves QM atoms, these should be * removed from the MM ilist */ - molt->ilist[ftype].nr -= (nratoms+1); - for (int l = j; l < molt->ilist[ftype].nr; l++) + InteractionList &ilist = molt->ilist[ftype]; + for (int k = j; k < ilist.size() - (nratoms + 1); k++) { - molt->ilist[ftype].iatoms[l] = molt->ilist[ftype].iatoms[l+(nratoms+1)]; + ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)]; } + ilist.iatoms.resize(ilist.size() - (nratoms + 1)); } else { @@ -1183,7 +1184,7 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr if (IS_CHEMBOND(i)) { int j = 0; - while (j < molt->ilist[i].nr) + while (j < molt->ilist[i].size()) { int a1 = molt->ilist[i].iatoms[j+1]; int a2 = molt->ilist[i].iatoms[j+2]; @@ -1269,7 +1270,7 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr { int nratoms = interaction_function[i].nratoms; int j = 0; - while (j < molt->ilist[i].nr) + while (j < molt->ilist[i].size()) { int a1 = molt->ilist[i].iatoms[j+1]; int a2 = molt->ilist[i].iatoms[j+2]; @@ -1281,11 +1282,12 @@ static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *gr /* since the interaction involves QM atoms, these should be * removed from the MM ilist */ - molt->ilist[i].nr -= (nratoms+1); - for (int k = j; k < molt->ilist[i].nr; k++) + InteractionList &ilist = molt->ilist[i]; + for (int k = j; k < ilist.size() - (nratoms + 1); k++) { - molt->ilist[i].iatoms[k] = molt->ilist[i].iatoms[k+(nratoms+1)]; + ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)]; } + ilist.iatoms.resize(ilist.size() - (nratoms + 1)); } else { diff --git a/src/gromacs/gmxpreprocess/vsite_parm.cpp b/src/gromacs/gmxpreprocess/vsite_parm.cpp index b441ed1a9a..5caf15877d 100644 --- a/src/gromacs/gmxpreprocess/vsite_parm.cpp +++ b/src/gromacs/gmxpreprocess/vsite_parm.cpp @@ -54,6 +54,7 @@ #include "gromacs/mdtypes/md_enums.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/topology.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/cstringutil.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/gmxassert.h" @@ -893,9 +894,6 @@ int set_vsites(bool bVerbose, t_atoms *atoms, gpp_atomtype_t atype, void set_vsites_ptype(bool bVerbose, gmx_moltype_t *molt) { int ftype, i; - int nra, nrd; - t_ilist *il; - t_iatom *ia, avsite; if (bVerbose) { @@ -903,12 +901,12 @@ void set_vsites_ptype(bool bVerbose, gmx_moltype_t *molt) } for (ftype = 0; ftype < F_NRE; ftype++) { - il = &molt->ilist[ftype]; + InteractionList *il = &molt->ilist[ftype]; if (interaction_function[ftype].flags & IF_VSITE) { - nra = interaction_function[ftype].nratoms; - nrd = il->nr; - ia = il->iatoms; + const int nra = interaction_function[ftype].nratoms; + const int nrd = il->size(); + gmx::ArrayRef ia = il->iatoms; if (debug && nrd) { @@ -919,11 +917,10 @@ void set_vsites_ptype(bool bVerbose, gmx_moltype_t *molt) for (i = 0; (i < nrd); ) { /* The virtual site */ - avsite = ia[1]; + int avsite = ia[i + 1]; molt->atoms.atom[avsite].ptype = eptVSite; i += nra+1; - ia += nra+1; } } } diff --git a/src/gromacs/listed-forces/disre.cpp b/src/gromacs/listed-forces/disre.cpp index a1b677fe54..d43cdb268b 100644 --- a/src/gromacs/listed-forces/disre.cpp +++ b/src/gromacs/listed-forces/disre.cpp @@ -71,13 +71,13 @@ void init_disres(FILE *fplog, const gmx_mtop_t *mtop, const gmx_multisim_t *ms, t_fcdata *fcd, t_state *state, gmx_bool bIsREMD) { - int fa, nmol, npair, np; - t_disresdata *dd; - history_t *hist; - gmx_mtop_ilistloop_t iloop; - const t_ilist *il; - char *ptr; - int type_min, type_max; + int fa, nmol, npair, np; + t_disresdata *dd; + history_t *hist; + gmx_mtop_ilistloop_t iloop; + const InteractionList *il; + char *ptr; + int type_min, type_max; dd = &(fcd->disres); @@ -132,13 +132,13 @@ void init_disres(FILE *fplog, const gmx_mtop_t *mtop, iloop = gmx_mtop_ilistloop_init(mtop); while (gmx_mtop_ilistloop_next(iloop, &il, &nmol)) { - if (nmol > 1 && il[F_DISRES].nr > 0 && ir->eDisre != edrEnsemble) + if (nmol > 1 && il[F_DISRES].size() > 0 && ir->eDisre != edrEnsemble) { gmx_fatal(FARGS, "NMR distance restraints with multiple copies of the same molecule are currently only supported with ensemble averaging. If you just want to restrain distances between atom pairs using a flat-bottomed potential, use a restraint potential (bonds type 10) instead."); } np = 0; - for (fa = 0; fa < il[F_DISRES].nr; fa += 3) + for (fa = 0; fa < il[F_DISRES].size(); fa += 3) { int type; diff --git a/src/gromacs/listed-forces/orires.cpp b/src/gromacs/listed-forces/orires.cpp index 90e4c9e0af..f76c52f1ae 100644 --- a/src/gromacs/listed-forces/orires.cpp +++ b/src/gromacs/listed-forces/orires.cpp @@ -108,12 +108,12 @@ void init_orires(FILE *fplog, od->eig = nullptr; od->v = nullptr; - int *nr_ex = nullptr; - int typeMin = INT_MAX; - int typeMax = 0; - gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop); - const t_ilist *il; - int nmol; + int *nr_ex = nullptr; + int typeMin = INT_MAX; + int typeMax = 0; + gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop); + const InteractionList *il; + int nmol; while (gmx_mtop_ilistloop_next(iloop, &il, &nmol)) { if (nmol > 1) @@ -121,7 +121,7 @@ void init_orires(FILE *fplog, gmx_fatal(FARGS, "Found %d copies of a molecule with orientation restrains while the current code only supports a single copy. If you want to ensemble average, run multiple copies of the system using the multi-sim feature of mdrun.", nmol); } - for (int i = 0; i < il[F_ORIRES].nr; i += 3) + for (int i = 0; i < il[F_ORIRES].size(); i += 3) { int type = il[F_ORIRES].iatoms[i]; int ex = mtop->ffparams.iparams[type].orires.ex; diff --git a/src/gromacs/mdlib/broadcaststructs.cpp b/src/gromacs/mdlib/broadcaststructs.cpp index 29157e7ff8..fb015cc9e3 100644 --- a/src/gromacs/mdlib/broadcaststructs.cpp +++ b/src/gromacs/mdlib/broadcaststructs.cpp @@ -41,6 +41,7 @@ #include +#include "gromacs/compat/make_unique.h" #include "gromacs/gmxlib/network.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/mdrun.h" @@ -308,7 +309,7 @@ void broadcastStateWithoutDynamics(const t_commrec *cr, t_state *state) } } -static void bc_ilists(const t_commrec *cr, t_ilist *ilist) +static void bc_ilists(const t_commrec *cr, InteractionList *ilist) { int ftype; @@ -317,11 +318,12 @@ static void bc_ilists(const t_commrec *cr, t_ilist *ilist) { for (ftype = 0; ftype < F_NRE; ftype++) { - if (ilist[ftype].nr > 0) + if (ilist[ftype].size() > 0) { block_bc(cr, ftype); - block_bc(cr, ilist[ftype].nr); - nblock_bc(cr, ilist[ftype].nr, ilist[ftype].iatoms); + int nr = ilist[ftype].size(); + block_bc(cr, nr); + nblock_bc(cr, nr, ilist[ftype].iatoms.data()); } } ftype = -1; @@ -331,16 +333,17 @@ static void bc_ilists(const t_commrec *cr, t_ilist *ilist) { for (ftype = 0; ftype < F_NRE; ftype++) { - ilist[ftype].nr = 0; + ilist[ftype].iatoms.clear(); } do { block_bc(cr, ftype); if (ftype >= 0) { - block_bc(cr, ilist[ftype].nr); - snew_bc(cr, ilist[ftype].iatoms, ilist[ftype].nr); - nblock_bc(cr, ilist[ftype].nr, ilist[ftype].iatoms); + int nr; + block_bc(cr, nr); + ilist[ftype].iatoms.resize(nr); + nblock_bc(cr, nr, ilist[ftype].iatoms.data()); } } while (ftype >= 0); @@ -816,8 +819,8 @@ void bcast_ir_mtop(const t_commrec *cr, t_inputrec *inputrec, gmx_mtop_t *mtop) block_bc(cr, mtop->bIntermolecularInteractions); if (mtop->bIntermolecularInteractions) { - snew_bc(cr, mtop->intermolecular_ilist, F_NRE); - bc_ilists(cr, mtop->intermolecular_ilist); + mtop->intermolecular_ilist = gmx::compat::make_unique(); + bc_ilists(cr, mtop->intermolecular_ilist->data()); } int nmolblock = mtop->molblock.size(); diff --git a/src/gromacs/mdlib/calc_verletbuf.cpp b/src/gromacs/mdlib/calc_verletbuf.cpp index f0a48401fc..80a7f19f7b 100644 --- a/src/gromacs/mdlib/calc_verletbuf.cpp +++ b/src/gromacs/mdlib/calc_verletbuf.cpp @@ -226,7 +226,6 @@ static void get_vsite_masses(const gmx_moltype_t *moltype, int *n_nonlin_vsite) { int ft, i; - const t_ilist *il; *n_nonlin_vsite = 0; @@ -235,9 +234,9 @@ static void get_vsite_masses(const gmx_moltype_t *moltype, { if (IS_VSITE(ft)) { - il = &moltype->ilist[ft]; + const InteractionList *il = &moltype->ilist[ft]; - for (i = 0; i < il->nr; i += 1+NRAL(ft)) + for (i = 0; i < il->size(); i += 1+NRAL(ft)) { const t_iparams *ip; real inv_mass, coeff, m_aj; @@ -356,7 +355,6 @@ static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop, verletbuf_atomtype_t *att; int natt; int ft, i, a1, a2, a3, a; - const t_ilist *il; const t_iparams *ip; atom_nonbonded_kinetic_prop_t *prop; real *vsite_m; @@ -385,9 +383,9 @@ static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop, for (ft = F_CONSTR; ft <= F_CONSTRNC; ft++) { - il = &moltype.ilist[ft]; + const InteractionList *il = &moltype.ilist[ft]; - for (i = 0; i < il->nr; i += 1+NRAL(ft)) + for (i = 0; i < il->size(); i += 1+NRAL(ft)) { ip = &mtop->ffparams.iparams[il->iatoms[i]]; a1 = il->iatoms[i+1]; @@ -407,9 +405,9 @@ static void get_verlet_buffer_atomtypes(const gmx_mtop_t *mtop, } } - il = &moltype.ilist[F_SETTLE]; + const InteractionList *il = &moltype.ilist[F_SETTLE]; - for (i = 0; i < il->nr; i += 1+NRAL(F_SETTLE)) + for (i = 0; i < il->size(); i += 1+NRAL(F_SETTLE)) { ip = &mtop->ffparams.iparams[il->iatoms[i]]; a1 = il->iatoms[i+1]; diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index 82b722ef02..128e68fbcd 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -133,10 +133,8 @@ class Constraints::Impl int nflexcon = 0; //! A list of atoms to constraints for each moleculetype. std::vector at2con_mt; - //! The size of at2settle = number of moltypes - int n_at2settle_mt = 0; - //! A list of atoms to settles. - int **at2settle_mt = nullptr; + //! A list of atoms to settles for each moleculetype + std::vector < std::vector < int>> at2settle_mt; //! Whether any SETTLES cross charge-group boundaries. bool bInterCGsettles = false; //! LINCS data. @@ -741,10 +739,24 @@ FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegra } } -t_blocka make_at2con(int numAtoms, - const t_ilist *ilists, - const t_iparams *iparams, - FlexibleConstraintTreatment flexibleConstraintTreatment) +/*! \brief Returns a block struct to go from atoms to constraints + * + * The block struct will contain constraint indices with lower indices + * directly matching the order in F_CONSTR and higher indices matching + * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR. + * + * \param[in] numAtoms The number of atoms to construct the list for + * \param[in] ilists The interaction lists, size F_NRE + * \param[in] iparams Interaction parameters, can be null when flexibleConstraintTreatment=Include + * \param[in] flexibleConstraintTreatment The flexible constraint treatment, see enum above + * \returns a block struct with all constraints for each atom + */ +template +static t_blocka +makeAtomsToConstraintsList(int numAtoms, + const T *ilists, + const t_iparams *iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment) { GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr, "With flexible constraint detection we need valid iparams"); @@ -752,9 +764,9 @@ t_blocka make_at2con(int numAtoms, for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++) { - const t_ilist &ilist = ilists[ftype]; - const int stride = 1 + NRAL(ftype); - for (int i = 0; i < ilist.nr; i += stride) + const T &ilist = ilists[ftype]; + const int stride = 1 + NRAL(ftype); + for (int i = 0; i < ilist.size(); i += stride) { if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || !isConstraintFlexible(iparams, ilist.iatoms[i])) @@ -788,9 +800,9 @@ t_blocka make_at2con(int numAtoms, int numConstraints = 0; for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++) { - const t_ilist &ilist = ilists[ftype]; - const int stride = 1 + NRAL(ftype); - for (int i = 0; i < ilist.nr; i += stride) + const T &ilist = ilists[ftype]; + const int stride = 1 + NRAL(ftype); + for (int i = 0; i < ilist.size(); i += stride) { if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || !isConstraintFlexible(iparams, ilist.iatoms[i])) @@ -808,49 +820,65 @@ t_blocka make_at2con(int numAtoms, return at2con; } -int countFlexibleConstraints(const t_ilist *ilist, - const t_iparams *iparams) +t_blocka make_at2con(int numAtoms, + const t_ilist *ilist, + const t_iparams *iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment) +{ + return makeAtomsToConstraintsList(numAtoms, ilist, iparams, flexibleConstraintTreatment); +} + +t_blocka make_at2con(const gmx_moltype_t &moltype, + const t_iparams *iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment) +{ + return makeAtomsToConstraintsList(moltype.atoms.nr, moltype.ilist, iparams, flexibleConstraintTreatment); +} + +//! Return the number of flexible constraints in the \c ilist and \c iparams. +template +static int +countFlexibleConstraintsTemplate(const T *ilist, + const t_iparams *iparams) { int nflexcon = 0; for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++) { const int numIatomsPerConstraint = 3; - int ncon = ilist[ftype].nr / numIatomsPerConstraint; - t_iatom *ia = ilist[ftype].iatoms; - for (int con = 0; con < ncon; con++) + for (int i = 0; i < ilist[ftype].size(); i += numIatomsPerConstraint) { - if (iparams[ia[0]].constr.dA == 0 && - iparams[ia[0]].constr.dB == 0) + const int type = ilist[ftype].iatoms[i]; + if (iparams[type].constr.dA == 0 && + iparams[type].constr.dB == 0) { nflexcon++; } - ia += numIatomsPerConstraint; } } return nflexcon; } -//! Returns the index of the settle to which each atom belongs. -static int *make_at2settle(int natoms, const t_ilist *ilist) +int countFlexibleConstraints(const t_ilist *ilist, + const t_iparams *iparams) { - int *at2s; - int a, stride, s; + return countFlexibleConstraintsTemplate(ilist, iparams); +} - snew(at2s, natoms); +//! Returns the index of the settle to which each atom belongs. +static std::vector make_at2settle(int natoms, + const InteractionList &ilist) +{ /* Set all to no settle */ - for (a = 0; a < natoms; a++) - { - at2s[a] = -1; - } + std::vector at2s(natoms, -1); - stride = 1 + NRAL(F_SETTLE); + const int stride = 1 + NRAL(F_SETTLE); - for (s = 0; s < ilist->nr; s += stride) + for (int s = 0; s < ilist.size(); s += stride) { - at2s[ilist->iatoms[s+1]] = s/stride; - at2s[ilist->iatoms[s+2]] = s/stride; - at2s[ilist->iatoms[s+3]] = s/stride; + at2s[ilist.iatoms[s + 1]] = s/stride; + at2s[ilist.iatoms[s + 2]] = s/stride; + at2s[ilist.iatoms[s + 3]] = s/stride; } return at2s; @@ -918,8 +946,7 @@ makeAtomToConstraintMappings(const gmx_mtop_t &mtop, mapping.reserve(mtop.moltype.size()); for (const gmx_moltype_t &moltype : mtop.moltype) { - mapping.push_back(make_at2con(moltype.atoms.nr, - moltype.ilist, + mapping.push_back(make_at2con(moltype, mtop.ffparams.iparams, flexibleConstraintTreatment)); } @@ -986,7 +1013,8 @@ Constraints::Impl::Impl(const gmx_mtop_t &mtop_p, for (const gmx_molblock_t &molblock : mtop.molblock) { - int count = countFlexibleConstraints(mtop.moltype[molblock.type].ilist, + int count = + countFlexibleConstraintsTemplate(mtop.moltype[molblock.type].ilist, mtop.ffparams.iparams); nflexcon += molblock.nmol*count; } @@ -1050,13 +1078,10 @@ Constraints::Impl::Impl(const gmx_mtop_t &mtop_p, settled = settle_init(mtop); /* Make an atom to settle index for use in domain decomposition */ - n_at2settle_mt = mtop.moltype.size(); - snew(at2settle_mt, n_at2settle_mt); for (size_t mt = 0; mt < mtop.moltype.size(); mt++) { - at2settle_mt[mt] = - make_at2settle(mtop.moltype[mt].atoms.nr, - &mtop.moltype[mt].ilist[F_SETTLE]); + at2settle_mt.emplace_back(make_at2settle(mtop.moltype[mt].atoms.nr, + mtop.moltype[mt].ilist[F_SETTLE])); } /* Allocate thread-local work arrays */ @@ -1111,7 +1136,7 @@ Constraints::atom2constraints_moltype() const return impl_->at2con_mt; } -int *const* Constraints::atom2settle_moltype() const +ArrayRef < const std::vector < int>> Constraints::atom2settle_moltype() const { return impl_->at2settle_mt; } @@ -1121,7 +1146,6 @@ bool inter_charge_group_constraints(const gmx_mtop_t &mtop) { const gmx_moltype_t *molt; const t_block *cgs; - const t_ilist *il; int *at2cg, cg, a, ftype, i; bool bInterCG; @@ -1130,9 +1154,9 @@ bool inter_charge_group_constraints(const gmx_mtop_t &mtop) { molt = &mtop.moltype[mtop.molblock[mb].type]; - if (molt->ilist[F_CONSTR].nr > 0 || - molt->ilist[F_CONSTRNC].nr > 0 || - molt->ilist[F_SETTLE].nr > 0) + if (molt->ilist[F_CONSTR].size() > 0 || + molt->ilist[F_CONSTRNC].size() > 0 || + molt->ilist[F_SETTLE].size() > 0) { cgs = &molt->cgs; snew(at2cg, molt->atoms.nr); @@ -1146,10 +1170,10 @@ bool inter_charge_group_constraints(const gmx_mtop_t &mtop) for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++) { - il = &molt->ilist[ftype]; - for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(ftype)) + const InteractionList &il = molt->ilist[ftype]; + for (i = 0; i < il.size() && !bInterCG; i += 1+NRAL(ftype)) { - if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]]) + if (at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 2]]) { bInterCG = TRUE; } @@ -1167,7 +1191,6 @@ bool inter_charge_group_settles(const gmx_mtop_t &mtop) { const gmx_moltype_t *molt; const t_block *cgs; - const t_ilist *il; int *at2cg, cg, a, ftype, i; bool bInterCG; @@ -1176,7 +1199,7 @@ bool inter_charge_group_settles(const gmx_mtop_t &mtop) { molt = &mtop.moltype[mtop.molblock[mb].type]; - if (molt->ilist[F_SETTLE].nr > 0) + if (molt->ilist[F_SETTLE].size() > 0) { cgs = &molt->cgs; snew(at2cg, molt->atoms.nr); @@ -1190,11 +1213,11 @@ bool inter_charge_group_settles(const gmx_mtop_t &mtop) for (ftype = F_SETTLE; ftype <= F_SETTLE; ftype++) { - il = &molt->ilist[ftype]; - for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(F_SETTLE)) + const InteractionList &il = molt->ilist[ftype]; + for (i = 0; i < il.size() && !bInterCG; i += 1+NRAL(F_SETTLE)) { - if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]] || - at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+3]]) + if (at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 2]] || + at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 3]]) { bInterCG = TRUE; } diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h index f1c5b1879b..f2f768dc7b 100644 --- a/src/gromacs/mdlib/constr.h +++ b/src/gromacs/mdlib/constr.h @@ -58,6 +58,7 @@ struct gmx_edsam; struct gmx_localtop_t; +struct gmx_moltype_t; struct gmx_mtop_t; struct gmx_multisim_t; struct gmx_wallcycle; @@ -177,7 +178,7 @@ class Constraints //! Getter for use by domain decomposition. const ArrayRef atom2constraints_moltype() const; //! Getter for use by domain decomposition. - int *const* atom2settle_moltype() const; + ArrayRef < const std::vector < int>> atom2settle_moltype() const; /*! \brief Return the data for reduction for determining * constraint RMS relative deviations, or an empty ArrayRef @@ -223,6 +224,21 @@ enum class FlexibleConstraintTreatment FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator); +/*! \brief Returns a block struct to go from atoms to constraints + * + * The block struct will contain constraint indices with lower indices + * directly matching the order in F_CONSTR and higher indices matching + * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR. + * + * \param[in] moltype The molecule data + * \param[in] iparams Interaction parameters, can be null when flexibleConstraintTreatment=Include + * \param[in] flexibleConstraintTreatment The flexible constraint treatment, see enum above + * \returns a block struct with all constraints for each atom + */ +t_blocka make_at2con(const gmx_moltype_t &moltype, + const t_iparams *iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment); + /*! \brief Returns a block struct to go from atoms to constraints * * The block struct will contain constraint indices with lower indices @@ -247,10 +263,23 @@ const t_blocka *atom2constraints_moltype(const Constraints *constr); int countFlexibleConstraints(const t_ilist *ilist, const t_iparams *iparams); -/*! \brief Macro for getting the constraint iatoms for a constraint number con +/*! \brief Returns the constraint iatoms for a constraint number con * which comes from a list where F_CONSTR and F_CONSTRNC constraints * are concatenated. */ -#define constr_iatomptr(nconstr, iatom_constr, iatom_constrnc, con) ((con) < (nconstr) ? (iatom_constr)+(con)*3 : (iatom_constrnc)+((con)-(nconstr))*3) +inline const int * +constr_iatomptr(gmx::ArrayRef iatom_constr, + gmx::ArrayRef iatom_constrnc, + int con) +{ + if (con*3 < iatom_constr.size()) + { + return iatom_constr.data() + con*3; + } + else + { + return iatom_constrnc.data() + con*3 - iatom_constr.size(); + } +}; /*! \brief Returns whether there are inter charge group constraints */ bool inter_charge_group_constraints(const gmx_mtop_t &mtop); diff --git a/src/gromacs/mdlib/constraintrange.cpp b/src/gromacs/mdlib/constraintrange.cpp index 9e8da0e0d6..9e1931a8ba 100644 --- a/src/gromacs/mdlib/constraintrange.cpp +++ b/src/gromacs/mdlib/constraintrange.cpp @@ -58,24 +58,20 @@ namespace gmx //! Recursing function to help find all adjacent constraints. static void constr_recur(const t_blocka *at2con, - const t_ilist *ilist, const t_iparams *iparams, + const InteractionList *ilist, const t_iparams *iparams, gmx_bool bTopB, int at, int depth, int nc, int *path, real r0, real r1, real *r2max, int *count) { - int ncon1; - t_iatom *ia1, *ia2; int c, con, a1; gmx_bool bUse; - t_iatom *ia; real len, rn0, rn1; (*count)++; - ncon1 = ilist[F_CONSTR].nr/3; - ia1 = ilist[F_CONSTR].iatoms; - ia2 = ilist[F_CONSTRNC].iatoms; + gmx::ArrayRef ia1 = ilist[F_CONSTR].iatoms; + gmx::ArrayRef ia2 = ilist[F_CONSTRNC].iatoms; /* Loop over all constraints connected to this atom */ for (c = at2con->index[at]; c < at2con->index[at+1]; c++) @@ -92,7 +88,7 @@ static void constr_recur(const t_blocka *at2con, } if (bUse) { - ia = constr_iatomptr(ncon1, ia1, ia2, con); + const int *ia = constr_iatomptr(ia1, ia2, con); /* Flexible constraints currently have length 0, which is incorrect */ if (!bTopB) { @@ -124,7 +120,7 @@ static void constr_recur(const t_blocka *at2con, { fprintf(debug, " %d %5.3f", path[a1], - iparams[constr_iatomptr(ncon1, ia1, ia2, con)[0]].constr.dA); + iparams[constr_iatomptr(ia1, ia2, con)[0]].constr.dA); } fprintf(debug, " %d %5.3f\n", con, len); } @@ -163,15 +159,15 @@ static real constr_r_max_moltype(const gmx_moltype_t *molt, t_blocka at2con; real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1; - if (molt->ilist[F_CONSTR].nr == 0 && - molt->ilist[F_CONSTRNC].nr == 0) + if (molt->ilist[F_CONSTR].size() == 0 && + molt->ilist[F_CONSTRNC].size() == 0) { return 0; } natoms = molt->atoms.nr; - at2con = make_at2con(natoms, molt->ilist, iparams, + at2con = make_at2con(*molt, iparams, flexibleConstraintTreatment(EI_DYNAMICS(ir->eI))); snew(path, 1+ir->nProjOrder); for (at = 0; at < 1+ir->nProjOrder; at++) diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index 7ff7ac628e..68ccf35436 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -707,7 +707,7 @@ static cginfo_mb_t *init_cginfo_mb(FILE *fplog, const gmx_mtop_t *mtop, int nral; nral = NRAL(ftype); - for (ia = 0; ia < molt->ilist[ftype].nr; ia += 1+nral) + for (ia = 0; ia < molt->ilist[ftype].size(); ia += 1+nral) { int a; @@ -1371,8 +1371,7 @@ static void make_nbf_tables(FILE *fp, static void count_tables(int ftype1, int ftype2, const gmx_mtop_t *mtop, int *ncount, int **count) { - const t_ilist *il; - int ftype, stride, i, j, tabnr; + int ftype, i, j, tabnr; // Loop over all moleculetypes for (const gmx_moltype_t &molt : mtop->moltype) @@ -1383,13 +1382,13 @@ static void count_tables(int ftype1, int ftype2, const gmx_mtop_t *mtop, // If the current interaction type is one of the types whose tables we're trying to count... if (ftype == ftype1 || ftype == ftype2) { - il = &molt.ilist[ftype]; - stride = 1 + NRAL(ftype); + const InteractionList &il = molt.ilist[ftype]; + const int stride = 1 + NRAL(ftype); // ... and there are actually some interactions for this type - for (i = 0; i < il->nr; i += stride) + for (i = 0; i < il.size(); i += stride) { // Find out which table index the user wanted - tabnr = mtop->ffparams.iparams[il->iatoms[i]].tab.table; + tabnr = mtop->ffparams.iparams[il.iatoms[i]].tab.table; if (tabnr < 0) { gmx_fatal(FARGS, "A bonded table number is smaller than 0: %d\n", tabnr); diff --git a/src/gromacs/mdlib/lincs.cpp b/src/gromacs/mdlib/lincs.cpp index d8b1d96b62..e37b56ae78 100644 --- a/src/gromacs/mdlib/lincs.cpp +++ b/src/gromacs/mdlib/lincs.cpp @@ -1385,36 +1385,34 @@ static void set_lincs_matrix(Lincs *li, real *invmass, real lambda) } //! Finds all triangles of atoms that share constraints to a central atom. -static int count_triangle_constraints(const t_ilist *ilist, - const t_blocka &at2con) +static int count_triangle_constraints(const InteractionList *ilist, + const t_blocka &at2con) { int ncon1, ncon_tot; - int c0, a00, a01, n1, c1, a10, a11, ac1, n2, c2, a20, a21; + int c0, n1, c1, ac1, n2, c2; int ncon_triangle; - bool bTriangle; - t_iatom *ia1, *ia2, *iap; - ncon1 = ilist[F_CONSTR].nr/3; - ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3; + ncon1 = ilist[F_CONSTR].size()/3; + ncon_tot = ncon1 + ilist[F_CONSTRNC].size()/3; - ia1 = ilist[F_CONSTR].iatoms; - ia2 = ilist[F_CONSTRNC].iatoms; + gmx::ArrayRef ia1 = ilist[F_CONSTR].iatoms; + gmx::ArrayRef ia2 = ilist[F_CONSTRNC].iatoms; ncon_triangle = 0; for (c0 = 0; c0 < ncon_tot; c0++) { - bTriangle = FALSE; - iap = constr_iatomptr(ncon1, ia1, ia2, c0); - a00 = iap[1]; - a01 = iap[2]; + bool bTriangle = FALSE; + const int *iap = constr_iatomptr(ia1, ia2, c0); + const int a00 = iap[1]; + const int a01 = iap[2]; for (n1 = at2con.index[a01]; n1 < at2con.index[a01+1]; n1++) { c1 = at2con.a[n1]; if (c1 != c0) { - iap = constr_iatomptr(ncon1, ia1, ia2, c1); - a10 = iap[1]; - a11 = iap[2]; + const int *iap = constr_iatomptr(ia1, ia2, c1); + const int a10 = iap[1]; + const int a11 = iap[2]; if (a10 == a01) { ac1 = a11; @@ -1428,9 +1426,9 @@ static int count_triangle_constraints(const t_ilist *ilist, c2 = at2con.a[n2]; if (c2 != c0 && c2 != c1) { - iap = constr_iatomptr(ncon1, ia1, ia2, c2); - a20 = iap[1]; - a21 = iap[2]; + const int *iap = constr_iatomptr(ia1, ia2, c2); + const int a20 = iap[1]; + const int a21 = iap[2]; if (a20 == a00 || a21 == a00) { bTriangle = TRUE; @@ -1449,26 +1447,24 @@ static int count_triangle_constraints(const t_ilist *ilist, } //! Finds sequences of sequential constraints. -static bool more_than_two_sequential_constraints(const t_ilist *ilist, - const t_blocka &at2con) +static bool more_than_two_sequential_constraints(const InteractionList *ilist, + const t_blocka &at2con) { - t_iatom *ia1, *ia2, *iap; int ncon1, ncon_tot, c; - int a1, a2; bool bMoreThanTwoSequentialConstraints; - ncon1 = ilist[F_CONSTR].nr/3; - ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3; + ncon1 = ilist[F_CONSTR].size()/3; + ncon_tot = ncon1 + ilist[F_CONSTRNC].size()/3; - ia1 = ilist[F_CONSTR].iatoms; - ia2 = ilist[F_CONSTRNC].iatoms; + gmx::ArrayRef ia1 = ilist[F_CONSTR].iatoms; + gmx::ArrayRef ia2 = ilist[F_CONSTRNC].iatoms; bMoreThanTwoSequentialConstraints = FALSE; for (c = 0; c < ncon_tot && !bMoreThanTwoSequentialConstraints; c++) { - iap = constr_iatomptr(ncon1, ia1, ia2, c); - a1 = iap[1]; - a2 = iap[2]; + const int *iap = constr_iatomptr(ia1, ia2, c); + const int a1 = iap[1]; + const int a2 = iap[2]; /* Check if this constraint has constraints connected at both atoms */ if (at2con.index[a1+1] - at2con.index[a1] > 1 && at2con.index[a2+1] - at2con.index[a2] > 1) diff --git a/src/gromacs/mdlib/membed.cpp b/src/gromacs/mdlib/membed.cpp index 70a2132969..dcb6206b71 100644 --- a/src/gromacs/mdlib/membed.cpp +++ b/src/gromacs/mdlib/membed.cpp @@ -867,12 +867,12 @@ static int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop) { for (j = 0; j < F_LJ; j++) { - mtop->moltype[i].ilist[j].nr = 0; + mtop->moltype[i].ilist[j].iatoms.clear(); } for (j = F_POSRES; j <= F_VSITEN; j++) { - mtop->moltype[i].ilist[j].nr = 0; + mtop->moltype[i].ilist[j].iatoms.clear(); } } } diff --git a/src/gromacs/mdlib/perf_est.cpp b/src/gromacs/mdlib/perf_est.cpp index b2a871f4fc..4705fca790 100644 --- a/src/gromacs/mdlib/perf_est.cpp +++ b/src/gromacs/mdlib/perf_est.cpp @@ -246,7 +246,7 @@ void count_bonded_distances(const gmx_mtop_t *mtop, const t_inputrec *ir, nd_c = NRAL(ftype) - 1; break; } - nbonds = molb.nmol*molt->ilist[ftype].nr/(1 + NRAL(ftype)); + nbonds = molb.nmol*molt->ilist[ftype].size()/(1 + NRAL(ftype)); ndtot_c += nbonds*nd_c; ndtot_simd += nbonds*nd_simd; } diff --git a/src/gromacs/mdlib/settle.cpp b/src/gromacs/mdlib/settle.cpp index ae85ab9e6d..1d96943314 100644 --- a/src/gromacs/mdlib/settle.cpp +++ b/src/gromacs/mdlib/settle.cpp @@ -181,14 +181,14 @@ static void settleparam_init(settleparam_t *p, settledata *settle_init(const gmx_mtop_t &mtop) { /* Check that we have only one settle type */ - int settle_type = -1; - gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop); - const t_ilist *ilist; - int nmol; - const int nral1 = 1 + NRAL(F_SETTLE); + int settle_type = -1; + gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop); + const InteractionList *ilist; + int nmol; + const int nral1 = 1 + NRAL(F_SETTLE); while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol)) { - for (int i = 0; i < ilist[F_SETTLE].nr; i += nral1) + for (int i = 0; i < ilist[F_SETTLE].size(); i += nral1) { if (settle_type == -1) { diff --git a/src/gromacs/mdlib/shellfc.cpp b/src/gromacs/mdlib/shellfc.cpp index 667c958150..e0b3fefb67 100644 --- a/src/gromacs/mdlib/shellfc.cpp +++ b/src/gromacs/mdlib/shellfc.cpp @@ -312,7 +312,6 @@ gmx_shellfc_t *init_shell_flexcon(FILE *fplog, int aS, aN = 0; /* Shell and nucleus */ int bondtypes[] = { F_BONDS, F_HARMONIC, F_CUBICBONDS, F_POLARIZATION, F_ANHARM_POL, F_WATER_POL }; #define NBT asize(bondtypes) - t_iatom *ia; gmx_mtop_atomloop_all_t aloop; const gmx_ffparams_t *ffparams; @@ -406,8 +405,8 @@ gmx_shellfc_t *init_shell_flexcon(FILE *fplog, { for (j = 0; (j < NBT); j++) { - ia = molt->ilist[bondtypes[j]].iatoms; - for (i = 0; (i < molt->ilist[bondtypes[j]].nr); ) + const int *ia = molt->ilist[bondtypes[j]].iatoms.data(); + for (i = 0; (i < molt->ilist[bondtypes[j]].size()); ) { type = ia[0]; ftype = ffparams->functype[type]; diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index 042bb6ceab..df4366d86c 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -2641,9 +2641,7 @@ static void low_do_pbc_mtop(FILE *fplog, int ePBC, const matrix box, } else { - /* Pass NULL iso fplog to avoid graph prints for each molecule type */ - mk_graph_ilist(nullptr, moltype.ilist, - 0, moltype.atoms.nr, FALSE, FALSE, graph); + mk_graph_moltype(moltype, graph); for (mol = 0; mol < molb.nmol; mol++) { diff --git a/src/gromacs/mdlib/tests/settle.cpp b/src/gromacs/mdlib/tests/settle.cpp index 6fd711ec93..dafed2cc8d 100644 --- a/src/gromacs/mdlib/tests/settle.cpp +++ b/src/gromacs/mdlib/tests/settle.cpp @@ -200,18 +200,14 @@ TEST_P(SettleTest, SatisfiesConstraints) mtop.moltype.resize(1); mtop.molblock.resize(1); mtop.molblock[0].type = 0; - const int settleStride = 1 + atomsPerSettle; - int *iatoms; - snew(iatoms, numSettles*settleStride); + std::vector &iatoms = mtop.moltype[0].ilist[F_SETTLE].iatoms; for (int i = 0; i < numSettles; ++i) { - iatoms[i*settleStride + 0] = settleType; - iatoms[i*settleStride + 1] = i*atomsPerSettle + 0; - iatoms[i*settleStride + 2] = i*atomsPerSettle + 1; - iatoms[i*settleStride + 3] = i*atomsPerSettle + 2; + iatoms.push_back(settleType); + iatoms.push_back(i*atomsPerSettle + 0); + iatoms.push_back(i*atomsPerSettle + 1); + iatoms.push_back(i*atomsPerSettle + 2); } - mtop.moltype[0].ilist[F_SETTLE].iatoms = iatoms; - mtop.moltype[0].ilist[F_SETTLE].nr = numSettles*settleStride; // Set up the SETTLE parameters. mtop.ffparams.ntypes = 1; @@ -239,8 +235,9 @@ TEST_P(SettleTest, SatisfiesConstraints) mdatoms.homenr = numSettles * atomsPerSettle; // Finally make the settle data structures - settledata *settled = settle_init(mtop); - settle_set_constraints(settled, &mtop.moltype[0].ilist[F_SETTLE], mdatoms); + settledata *settled = settle_init(mtop); + const t_ilist ilist = { mtop.moltype[0].ilist[F_SETTLE].size(), 0, mtop.moltype[0].ilist[F_SETTLE].iatoms.data(), 0 }; + settle_set_constraints(settled, &ilist, mdatoms); // Copy the original positions from the array of doubles to a vector of reals std::vector startingPositions(std::begin(g_positions), std::end(g_positions)); diff --git a/src/gromacs/mdlib/vsite.cpp b/src/gromacs/mdlib/vsite.cpp index db5cb53fff..476b137dfa 100644 --- a/src/gromacs/mdlib/vsite.cpp +++ b/src/gromacs/mdlib/vsite.cpp @@ -164,12 +164,13 @@ struct VsiteThread * * \param[in] ilist The interaction list */ -static int vsiteIlistNrCount(const t_ilist *ilist) +template +static int vsiteIlistNrCount(const T *ilist) { int nr = 0; for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++) { - nr += ilist[ftype].nr; + nr += ilist[ftype].size(); } return nr; @@ -769,9 +770,16 @@ void constructVsitesGlobal(const gmx_mtop_t &mtop, int atomOffset = mtop.moleculeBlockIndices[mb].globalAtomStart; for (int mol = 0; mol < molb.nmol; mol++) { + t_ilist ilist[F_NRE]; + for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++) + { + ilist[ftype].nr = molt.ilist[ftype].size(); + ilist[ftype].iatoms = const_cast(molt.ilist[ftype].iatoms.data()); + } + construct_vsites(nullptr, as_rvec_array(x.data()) + atomOffset, 0.0, nullptr, - mtop.ffparams.iparams, molt.ilist, + mtop.ffparams.iparams, ilist, epbcNONE, TRUE, nullptr, nullptr); atomOffset += molt.atoms.nr; } @@ -1846,15 +1854,14 @@ int count_intercg_vsites(const gmx_mtop_t *mtop) std::vector a2cg = atom2cg(molt.cgs); for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++) { - int nral = NRAL(ftype); - const t_ilist &il = molt.ilist[ftype]; - const t_iatom *ia = il.iatoms; - for (int i = 0; i < il.nr; i += 1 + nral) + const int nral = NRAL(ftype); + const InteractionList &il = molt.ilist[ftype]; + for (int i = 0; i < il.size(); i += 1 + nral) { - int cg = a2cg[ia[1+i]]; + int cg = a2cg[il.iatoms[1 + i]]; for (int a = 1; a < nral; a++) { - if (a2cg[ia[1+a]] != cg) + if (a2cg[il.iatoms[1 + a]] != cg) { n_intercg_vsite += molb.nmol; break; @@ -1867,8 +1874,9 @@ int count_intercg_vsites(const gmx_mtop_t *mtop) return n_intercg_vsite; } +template static VsitePbc -get_vsite_pbc(const t_iparams *iparams, const t_ilist *ilist, +get_vsite_pbc(const t_iparams *iparams, const T *ilist, const t_atom *atom, const t_mdatoms *md, const t_block &cgs) { @@ -1893,17 +1901,16 @@ get_vsite_pbc(const t_iparams *iparams, const t_ilist *ilist, { { // TODO remove me int nral = NRAL(ftype); - const t_ilist *il = &ilist[ftype]; - const t_iatom *ia = il->iatoms; + const T &il = ilist[ftype]; std::vector &vsite_pbc_f = vsite_pbc[ftype - F_VSITE2]; - vsite_pbc_f.resize(il->nr/(1 + nral)); + vsite_pbc_f.resize(il.size()/(1 + nral)); int i = 0; - while (i < il->nr) + while (i < il.size()) { int vsi = i/(1 + nral); - t_iatom vsite = ia[i+1]; + t_iatom vsite = il.iatoms[i + 1]; int cg_v = a2cg[vsite]; /* A value of -2 signals that this vsite and its contructing * atoms are all within the same cg, so no pbc is required. @@ -1913,10 +1920,10 @@ get_vsite_pbc(const t_iparams *iparams, const t_ilist *ilist, int nc3 = 0; if (ftype == F_VSITEN) { - nc3 = 3*iparams[ia[i]].vsiten.n; + nc3 = 3*iparams[il.iatoms[i]].vsiten.n; for (int j = 0; j < nc3; j += 3) { - if (a2cg[ia[i+j+2]] != cg_v) + if (a2cg[il.iatoms[i + j + 2]] != cg_v) { vsite_pbc_f[vsi] = -1; } @@ -1926,7 +1933,7 @@ get_vsite_pbc(const t_iparams *iparams, const t_ilist *ilist, { for (int a = 1; a < nral; a++) { - if (a2cg[ia[i+1+a]] != cg_v) + if (a2cg[il.iatoms[i + 1 + a]] != cg_v) { vsite_pbc_f[vsi] = -1; } @@ -1953,7 +1960,7 @@ get_vsite_pbc(const t_iparams *iparams, const t_ilist *ilist, */ vsite_pbc_f[vsi] = vsite; } - else if (cg_v != a2cg[ia[1+i+1]]) + else if (cg_v != a2cg[il.iatoms[1 + i + 1]]) { /* This vsite has a different charge group index * than it's first constructing atom diff --git a/src/gromacs/pbcutil/mshift.cpp b/src/gromacs/pbcutil/mshift.cpp index ef86515707..cc410c7452 100644 --- a/src/gromacs/pbcutil/mshift.cpp +++ b/src/gromacs/pbcutil/mshift.cpp @@ -46,6 +46,7 @@ #include "gromacs/pbcutil/pbc.h" #include "gromacs/topology/idef.h" #include "gromacs/topology/ifunc.h" +#include "gromacs/topology/topology.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" #include "gromacs/utility/strconvert.h" @@ -93,26 +94,25 @@ static void add_gbond(t_graph *g, int a0, int a1) * When a non-null part array is supplied with part indices for each atom, * edges are only added when atoms have a different part index. */ -static bool mk_igraph(t_graph *g, int ftype, const t_ilist *il, +template +static bool mk_igraph(t_graph *g, int ftype, const T &il, int at_start, int at_end, const int *part) { - t_iatom *ia; int i, j, np; int end; bool addedEdge = false; - end = il->nr; - ia = il->iatoms; + end = il.size(); i = 0; while (i < end) { np = interaction_function[ftype].nratoms; - if (np > 1 && ia[1] >= at_start && ia[1] < at_end) + if (np > 1 && il.iatoms[i + 1] >= at_start && il.iatoms[i + 1] < at_end) { - if (ia[np] >= at_end) + if (il.iatoms[i + np] >= at_end) { gmx_fatal(FARGS, "Molecule in topology has atom numbers below and " @@ -125,8 +125,8 @@ static bool mk_igraph(t_graph *g, int ftype, const t_ilist *il, if (ftype == F_SETTLE) { /* Bond all the atoms in the settle */ - add_gbond(g, ia[1], ia[2]); - add_gbond(g, ia[1], ia[3]); + add_gbond(g, il.iatoms[i + 1], il.iatoms[i + 2]); + add_gbond(g, il.iatoms[i + 1], il.iatoms[i + 3]); addedEdge = true; } else if (part == nullptr) @@ -134,7 +134,7 @@ static bool mk_igraph(t_graph *g, int ftype, const t_ilist *il, /* Simply add this bond */ for (j = 1; j < np; j++) { - add_gbond(g, ia[j], ia[j+1]); + add_gbond(g, il.iatoms[i + j], il.iatoms[i + j+1]); } addedEdge = true; } @@ -143,15 +143,15 @@ static bool mk_igraph(t_graph *g, int ftype, const t_ilist *il, /* Add this bond when it connects two unlinked parts of the graph */ for (j = 1; j < np; j++) { - if (part[ia[j]] != part[ia[j+1]]) + if (part[il.iatoms[i + j]] != part[il.iatoms[i + j+1]]) { - add_gbond(g, ia[j], ia[j+1]); + add_gbond(g, il.iatoms[i + j], il.iatoms[i + j+1]); addedEdge = true; } } } } - ia += np+1; + i += np+1; } @@ -197,36 +197,35 @@ void p_graph(FILE *log, const char *title, t_graph *g) fflush(log); } -static void calc_1se(t_graph *g, int ftype, const t_ilist *il, +template +static void calc_1se(t_graph *g, int ftype, const T &il, int nbond[], int at_start, int at_end) { int k, nratoms, end, j; - t_iatom *ia, iaa; - end = il->nr; + end = il.size(); - ia = il->iatoms; - for (j = 0; (j < end); j += nratoms+1, ia += nratoms+1) + for (j = 0; (j < end); j += nratoms + 1) { nratoms = interaction_function[ftype].nratoms; if (ftype == F_SETTLE) { - iaa = ia[1]; + const int iaa = il.iatoms[j + 1]; if (iaa >= at_start && iaa < at_end) { - nbond[iaa] += 2; - nbond[ia[2]] += 1; - nbond[ia[3]] += 1; - g->at_start = std::min(g->at_start, iaa); - g->at_end = std::max(g->at_end, iaa+2+1); + nbond[iaa] += 2; + nbond[il.iatoms[j + 2]] += 1; + nbond[il.iatoms[j + 3]] += 1; + g->at_start = std::min(g->at_start, iaa); + g->at_end = std::max(g->at_end, iaa+2+1); } } else { for (k = 1; (k <= nratoms); k++) { - iaa = ia[k]; + const int iaa = il.iatoms[j + k]; if (iaa >= at_start && iaa < at_end) { g->at_start = std::min(g->at_start, iaa); @@ -249,7 +248,8 @@ static void calc_1se(t_graph *g, int ftype, const t_ilist *il, } } -static int calc_start_end(FILE *fplog, t_graph *g, const t_ilist il[], +template +static int calc_start_end(FILE *fplog, t_graph *g, const T il[], int at_start, int at_end, int nbond[]) { @@ -265,7 +265,7 @@ static int calc_start_end(FILE *fplog, t_graph *g, const t_ilist il[], { if (interaction_function[i].flags & IF_CHEMBOND) { - calc_1se(g, i, &il[i], nbond, at_start, at_end); + calc_1se(g, i, il[i], nbond, at_start, at_end); } } /* Then add all the other interactions in fixed lists, but first @@ -275,7 +275,7 @@ static int calc_start_end(FILE *fplog, t_graph *g, const t_ilist il[], { if (!(interaction_function[i].flags & IF_CHEMBOND)) { - calc_1se(g, i, &il[i], nbond, at_start, at_end); + calc_1se(g, i, il[i], nbond, at_start, at_end); } } @@ -377,10 +377,12 @@ static gmx_bool determine_graph_parts(t_graph *g, int *part) return bMultiPart; } -void mk_graph_ilist(FILE *fplog, - const t_ilist *ilist, int at_start, int at_end, - gmx_bool bShakeOnly, gmx_bool bSettle, - t_graph *g) +template +static void +mk_graph_ilist(FILE *fplog, + const T *ilist, int at_start, int at_end, + gmx_bool bShakeOnly, gmx_bool bSettle, + t_graph *g) { int *nbond; int i, nbtot; @@ -426,7 +428,7 @@ void mk_graph_ilist(FILE *fplog, { if (interaction_function[i].flags & IF_CHEMBOND) { - mk_igraph(g, i, &(ilist[i]), at_start, at_end, nullptr); + mk_igraph(g, i, ilist[i], at_start, at_end, nullptr); } } @@ -447,7 +449,7 @@ void mk_graph_ilist(FILE *fplog, if (!(interaction_function[i].flags & IF_CHEMBOND)) { bool addedEdgeForType = - mk_igraph(g, i, &(ilist[i]), at_start, at_end, nbond); + mk_igraph(g, i, ilist[i], at_start, at_end, nbond); addedEdge = (addedEdge || addedEdgeForType); } } @@ -468,10 +470,10 @@ void mk_graph_ilist(FILE *fplog, else { /* This is a special thing used in splitter.c to generate shake-blocks */ - mk_igraph(g, F_CONSTR, &(ilist[F_CONSTR]), at_start, at_end, nullptr); + mk_igraph(g, F_CONSTR, ilist[F_CONSTR], at_start, at_end, nullptr); if (bSettle) { - mk_igraph(g, F_SETTLE, &(ilist[F_SETTLE]), at_start, at_end, nullptr); + mk_igraph(g, F_SETTLE, ilist[F_SETTLE], at_start, at_end, nullptr); } } g->nbound = 0; @@ -497,6 +499,14 @@ void mk_graph_ilist(FILE *fplog, } } +void mk_graph_moltype(const gmx_moltype_t &moltype, + t_graph *g) +{ + mk_graph_ilist(nullptr, moltype.ilist, 0, moltype.atoms.nr, + FALSE, FALSE, + g); +} + t_graph *mk_graph(FILE *fplog, const t_idef *idef, int at_start, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle) diff --git a/src/gromacs/pbcutil/mshift.h b/src/gromacs/pbcutil/mshift.h index db4a83f9c1..f77c9ed5da 100644 --- a/src/gromacs/pbcutil/mshift.h +++ b/src/gromacs/pbcutil/mshift.h @@ -42,8 +42,9 @@ #include "gromacs/math/vectypes.h" #include "gromacs/utility/basedefinitions.h" +struct InteractionList; +struct gmx_moltype_t; struct t_idef; -struct t_ilist; typedef enum { egcolWhite, egcolGrey, egcolBlack, egcolNR @@ -89,11 +90,9 @@ t_graph *mk_graph(FILE *fplog, * If bSettle && bShakeOnly the settles are used too. */ -void mk_graph_ilist(FILE *fplog, - const struct t_ilist *ilist, int at_start, int at_end, - gmx_bool bShakeOnly, gmx_bool bSettle, - t_graph *g); -/* As mk_graph, but takes t_ilist iso t_idef and does not allocate g */ +void mk_graph_moltype(const gmx_moltype_t &moltype, + t_graph *g); +/* As mk_graph, but takes gmx_moltype_t iso t_idef and does not allocate g */ void done_graph(t_graph *g); diff --git a/src/gromacs/tools/convert_tpr.cpp b/src/gromacs/tools/convert_tpr.cpp index 343fa5525b..8d08c51d2b 100644 --- a/src/gromacs/tools/convert_tpr.cpp +++ b/src/gromacs/tools/convert_tpr.cpp @@ -293,7 +293,12 @@ static void reduce_topology_x(int gnx, int index[], mtop->moltype[0].atoms = top.atoms; for (i = 0; i < F_NRE; i++) { - mtop->moltype[0].ilist[i] = top.idef.il[i]; + InteractionList &ilist = mtop->moltype[0].ilist[i]; + ilist.iatoms.resize(top.idef.il[i].nr); + for (int j = 0; j < top.idef.il[i].nr; j++) + { + ilist.iatoms[j] = top.idef.il[i].iatoms[j]; + } } mtop->moltype[0].atoms = top.atoms; mtop->moltype[0].cgs = top.cgs; diff --git a/src/gromacs/topology/idef.cpp b/src/gromacs/topology/idef.cpp index 652a36bb5d..5f201efb88 100644 --- a/src/gromacs/topology/idef.cpp +++ b/src/gromacs/topology/idef.cpp @@ -306,28 +306,28 @@ void pr_iparams(FILE *fp, t_functype ftype, const t_iparams *iparams) } } -void pr_ilist(FILE *fp, int indent, const char *title, - const t_functype *functype, const t_ilist *ilist, - gmx_bool bShowNumbers, - gmx_bool bShowParameters, const t_iparams *iparams) +template +static void +printIlist(FILE *fp, int indent, const char *title, + const t_functype *functype, const T *ilist, + gmx_bool bShowNumbers, + gmx_bool bShowParameters, const t_iparams *iparams) { int i, j, k, type, ftype; - t_iatom *iatoms; - if (available(fp, ilist, indent, title) && ilist->nr > 0) + if (available(fp, ilist, indent, title) && ilist->size() > 0) { indent = pr_title(fp, indent, title); pr_indent(fp, indent); - fprintf(fp, "nr: %d\n", ilist->nr); - if (ilist->nr > 0) + fprintf(fp, "nr: %d\n", ilist->size()); + if (ilist->size() > 0) { pr_indent(fp, indent); fprintf(fp, "iatoms:\n"); - iatoms = ilist->iatoms; - for (i = j = 0; i < ilist->nr; ) + for (i = j = 0; i < ilist->size(); ) { pr_indent(fp, indent+INDENT); - type = *(iatoms++); + type = ilist->iatoms[i]; ftype = functype[type]; if (bShowNumbers) { @@ -337,7 +337,7 @@ void pr_ilist(FILE *fp, int indent, const char *title, printf("(%s)", interaction_function[ftype].name); for (k = 0; k < interaction_function[ftype].nratoms; k++) { - fprintf(fp, " %3d", *(iatoms++)); + fprintf(fp, " %3d", ilist->iatoms[i + 1 + k]); } if (bShowParameters) { @@ -351,6 +351,15 @@ void pr_ilist(FILE *fp, int indent, const char *title, } } +void pr_ilist(FILE *fp, int indent, const char *title, + const t_functype *functype, const InteractionList *ilist, + gmx_bool bShowNumbers, + gmx_bool bShowParameters, const t_iparams *iparams) +{ + printIlist(fp, indent, title, functype, ilist, + bShowNumbers, bShowParameters, iparams); +} + static void pr_cmap(FILE *fp, int indent, const char *title, const gmx_cmap_t *cmap_grid, gmx_bool bShowNumbers) { @@ -438,9 +447,9 @@ void pr_idef(FILE *fp, int indent, const char *title, const t_idef *idef, for (j = 0; (j < F_NRE); j++) { - pr_ilist(fp, indent, interaction_function[j].longname, - idef->functype, &idef->il[j], bShowNumbers, - bShowParameters, idef->iparams); + printIlist(fp, indent, interaction_function[j].longname, + idef->functype, &idef->il[j], bShowNumbers, + bShowParameters, idef->iparams); } } } diff --git a/src/gromacs/topology/idef.h b/src/gromacs/topology/idef.h index 2c2e4f1fb1..6f580f97aa 100644 --- a/src/gromacs/topology/idef.h +++ b/src/gromacs/topology/idef.h @@ -39,6 +39,9 @@ #include +#include +#include + #include "gromacs/math/vectypes.h" #include "gromacs/utility/basedefinitions.h" #include "gromacs/utility/real.h" @@ -274,22 +277,59 @@ typedef union t_iparams typedef int t_functype; -/* +/* List of listed interactions, see description further down. + * + * TODO: Consider storing the function type as well. + * TODO: Consider providing per interaction access. + */ +struct InteractionList +{ + /* Returns the total number of elements in iatoms */ + int size() const + { + return iatoms.size(); + } + + /* List of interactions, see explanation further down */ + std::vector iatoms; +}; + +/* List of interaction lists, one list for each interaction type + * + * TODO: Consider only including entries in use instead of all F_NRE + */ +typedef std::array InteractionLists; + +/* Deprecated list of listed interactions. + * * The nonperturbed/perturbed interactions are now separated (sorted) in the * ilist, such that the first 0..(nr_nonperturbed-1) ones are exactly that, and * the remaining ones from nr_nonperturbed..(nr-1) are perturbed bonded * interactions. */ -typedef struct t_ilist +struct t_ilist { + /* Returns the total number of elements in iatoms */ + int size() const + { + return nr; + } + int nr; int nr_nonperturbed; t_iatom *iatoms; int nalloc; -} t_ilist; +}; + +/* TODO: Replace t_ilist in gmx_localtop_t by InteractionList. + * The nr_nonperturbed functionality needs to be ported. + * Remove t_topology. + * Remove t_ilist and remove templating on list type + * in mshift.cpp, constr.cpp, vsite.cpp and domdec_topology.cpp. + */ /* - * The struct t_ilist defines a list of atoms with their interactions. + * The structs InteractionList and t_ilist defines a list of atoms with their interactions. * General field description: * int nr * the size (nr elements) of the interactions array (iatoms[]). @@ -387,7 +427,7 @@ typedef struct t_idef void pr_iparams(FILE *fp, t_functype ftype, const t_iparams *iparams); void pr_ilist(FILE *fp, int indent, const char *title, - const t_functype *functype, const t_ilist *ilist, + const t_functype *functype, const InteractionList *ilist, gmx_bool bShowNumbers, gmx_bool bShowParameters, const t_iparams *iparams); void pr_ffparams(FILE *fp, int indent, const char *title, diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp index fed567e58c..96e0ef919c 100644 --- a/src/gromacs/topology/mtop_util.cpp +++ b/src/gromacs/topology/mtop_util.cpp @@ -398,8 +398,9 @@ static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop) sfree(iloop); } -gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop, - const t_ilist **ilist_mol, int *nmol) +gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop, + const InteractionList **ilist_mol, + int *nmol) { if (iloop == nullptr) { @@ -412,7 +413,7 @@ gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop, if (iloop->mblock == static_cast(iloop->mtop->molblock.size()) && iloop->mtop->bIntermolecularInteractions) { - *ilist_mol = iloop->mtop->intermolecular_ilist; + *ilist_mol = iloop->mtop->intermolecular_ilist->data(); *nmol = 1; return TRUE; } @@ -456,8 +457,9 @@ static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop) sfree(iloop); } -gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, - const t_ilist **ilist_mol, int *atnr_offset) +gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, + const InteractionList **ilist_mol, + int *atnr_offset) { if (iloop == nullptr) @@ -486,7 +488,7 @@ gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, if (iloop->mblock == iloop->mtop->molblock.size() && iloop->mtop->bIntermolecularInteractions) { - *ilist_mol = iloop->mtop->intermolecular_ilist; + *ilist_mol = iloop->mtop->intermolecular_ilist->data(); *atnr_offset = 0; return TRUE; } @@ -506,21 +508,21 @@ gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, int gmx_mtop_ftype_count(const gmx_mtop_t *mtop, int ftype) { - gmx_mtop_ilistloop_t iloop; - const t_ilist *il; - int n, nmol; + gmx_mtop_ilistloop_t iloop; + const InteractionList *il; + int n, nmol; n = 0; iloop = gmx_mtop_ilistloop_init(mtop); while (gmx_mtop_ilistloop_next(iloop, &il, &nmol)) { - n += nmol*il[ftype].nr/(1+NRAL(ftype)); + n += nmol*il[ftype].size()/(1+NRAL(ftype)); } if (mtop->bIntermolecularInteractions) { - n += mtop->intermolecular_ilist[ftype].nr/(1+NRAL(ftype)); + n += (*mtop->intermolecular_ilist)[ftype].size()/(1+NRAL(ftype)); } return n; @@ -747,24 +749,28 @@ static void blockacat(t_blocka *dest, const t_blocka *src, int copies, dest->index[dest->nr] = dest->nra; } -static void ilistcat(int ftype, t_ilist *dest, const t_ilist *src, int copies, - int dnum, int snum) +static void ilistcat(int ftype, + t_ilist *dest, + const InteractionList &src, + int copies, + int dnum, + int snum) { int nral, c, i, a; nral = NRAL(ftype); - dest->nalloc = dest->nr + copies*src->nr; + dest->nalloc = dest->nr + copies*src.size(); srenew(dest->iatoms, dest->nalloc); for (c = 0; c < copies; c++) { - for (i = 0; i < src->nr; ) + for (i = 0; i < src.size(); ) { - dest->iatoms[dest->nr++] = src->iatoms[i++]; + dest->iatoms[dest->nr++] = src.iatoms[i++]; for (a = 0; a < nral; a++) { - dest->iatoms[dest->nr++] = dnum + src->iatoms[i++]; + dest->iatoms[dest->nr++] = dnum + src.iatoms[i++]; } } dnum += snum; @@ -941,22 +947,22 @@ static void gen_local_top(const gmx_mtop_t *mtop, for (ftype = 0; ftype < F_NRE; ftype++) { if (bMergeConstr && - ftype == F_CONSTR && molt.ilist[F_CONSTRNC].nr > 0) + ftype == F_CONSTR && molt.ilist[F_CONSTRNC].size() > 0) { /* Merge all constrains into one ilist. * This simplifies the constraint code. */ for (int mol = 0; mol < molb.nmol; mol++) { - ilistcat(ftype, &idef->il[F_CONSTR], &molt.ilist[F_CONSTR], + ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTR], 1, destnr + mol*srcnr, srcnr); - ilistcat(ftype, &idef->il[F_CONSTR], &molt.ilist[F_CONSTRNC], + ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTRNC], 1, destnr + mol*srcnr, srcnr); } } else if (!(bMergeConstr && ftype == F_CONSTRNC)) { - ilistcat(ftype, &idef->il[ftype], &molt.ilist[ftype], + ilistcat(ftype, &idef->il[ftype], molt.ilist[ftype], molb.nmol, destnr, srcnr); } } @@ -978,7 +984,7 @@ static void gen_local_top(const gmx_mtop_t *mtop, { for (ftype = 0; ftype < F_NRE; ftype++) { - ilistcat(ftype, &idef->il[ftype], &mtop->intermolecular_ilist[ftype], + ilistcat(ftype, &idef->il[ftype], (*mtop->intermolecular_ilist)[ftype], 1, 0, mtop->natoms); } } diff --git a/src/gromacs/topology/mtop_util.h b/src/gromacs/topology/mtop_util.h index 321032a3a8..d8d9bedc97 100644 --- a/src/gromacs/topology/mtop_util.h +++ b/src/gromacs/topology/mtop_util.h @@ -49,7 +49,7 @@ struct gmx_localtop_t; struct t_atom; struct t_atoms; struct t_block; -struct t_ilist; +struct InteractionList; struct t_symtab; // TODO All of the functions taking a const gmx_mtop * are deprecated @@ -174,9 +174,9 @@ gmx_mtop_ilistloop_init(const gmx_mtop_t &mtop); * When at the end, destroys iloop and returns FALSE. */ gmx_bool -gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop, - const t_ilist **ilist_mol, int *nmol); - +gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop, + const InteractionList **ilist_mol, + int *nmol); /* Abstract type for ilist loop over all ilists of all molecules */ typedef struct gmx_mtop_ilistloop_all *gmx_mtop_ilistloop_all_t; @@ -196,8 +196,9 @@ gmx_mtop_ilistloop_all_init(const gmx_mtop_t *mtop); * When at the end, destroys iloop and returns FALSE. */ gmx_bool -gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, - const t_ilist **ilist_mol, int *atnr_offset); +gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, + const InteractionList **ilist_mol, + int *atnr_offset); /* Returns the total number of interactions in the system of type ftype */ diff --git a/src/gromacs/topology/topology.cpp b/src/gromacs/topology/topology.cpp index 88c2bfaaaf..e26c64d221 100644 --- a/src/gromacs/topology/topology.cpp +++ b/src/gromacs/topology/topology.cpp @@ -115,14 +115,6 @@ gmx_moltype_t::gmx_moltype_t() : excls() { init_t_atoms(&atoms, 0, FALSE); - - for (int ftype = 0; ftype < F_NRE; ftype++) - { - ilist[ftype].nr = 0; - ilist[ftype].nr_nonperturbed = 0; - ilist[ftype].iatoms = nullptr; - ilist[ftype].nalloc = 0; - } } gmx_moltype_t::~gmx_moltype_t() @@ -130,12 +122,6 @@ gmx_moltype_t::~gmx_moltype_t() done_atom(&atoms); done_block(&cgs); done_blocka(&excls); - - for (int f = 0; f < F_NRE; f++) - { - sfree(ilist[f].iatoms); - ilist[f].nalloc = 0; - } } void done_gmx_groups_t(gmx_groups_t *g) @@ -446,7 +432,7 @@ void pr_mtop(FILE *fp, int indent, const char *title, const gmx_mtop_t *mtop, { pr_ilist(fp, indent, interaction_function[j].longname, mtop->ffparams.functype, - &mtop->intermolecular_ilist[j], + &(*mtop->intermolecular_ilist)[j], bShowNumbers, bShowParameters, mtop->ffparams.iparams); } } diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h index 2a18e59174..f2b03b839f 100644 --- a/src/gromacs/topology/topology.h +++ b/src/gromacs/topology/topology.h @@ -73,9 +73,11 @@ struct gmx_moltype_t char **name; /**< Name of the molecule type */ t_atoms atoms; /**< The atoms in this molecule */ - t_ilist ilist[F_NRE]; /**< Interaction list with local indices */ + InteractionList ilist[F_NRE]; /**< Interaction list with local indices */ t_block cgs; /**< The charge groups */ t_blocka excls; /**< The exclusions */ + + /* TODO: Change ilist to InteractionLists */ }; /*! \brief Block of molecules of the same type, used in gmx_mtop_t */ @@ -131,22 +133,22 @@ struct gmx_mtop_t //NOLINT(clang-analyzer-optin.performance.Padding) ~gmx_mtop_t(); - char **name; /* Name of the topology */ - gmx_ffparams_t ffparams; - std::vector moltype; - std::vector molblock; - gmx_bool bIntermolecularInteractions; /* Are there intermolecular - * interactions? */ - t_ilist *intermolecular_ilist; /* List of intermolecular interactions - * using system wide atom indices, - * either NULL or size F_NRE */ + char **name; /* Name of the topology */ + gmx_ffparams_t ffparams; + std::vector moltype; + std::vector molblock; + gmx_bool bIntermolecularInteractions; /* Are there intermolecular + * interactions? */ + std::unique_ptr intermolecular_ilist; /* List of intermolecular interactions + * using system wide atom indices, + * either NULL or size F_NRE */ int natoms; - int maxres_renum; /* Parameter for residue numbering */ - int maxresnr; /* The maximum residue number in moltype */ - t_atomtypes atomtypes; /* Atomtype properties */ - gmx_groups_t groups; /* Groups of atoms for different purposes */ - t_symtab symtab; /* The symbol table */ - bool haveMoleculeIndices; /* Tells whether we have valid molecule indices */ + int maxres_renum; /* Parameter for residue numbering */ + int maxresnr; /* The maximum residue number in moltype */ + t_atomtypes atomtypes; /* Atomtype properties */ + gmx_groups_t groups; /* Groups of atoms for different purposes */ + t_symtab symtab; /* The symbol table */ + bool haveMoleculeIndices; /* Tells whether we have valid molecule indices */ /* Derived data */ std::vector moleculeBlockIndices; /* Indices for each molblock entry for fast lookup of atom properties */ diff --git a/src/gromacs/topology/topsort.cpp b/src/gromacs/topology/topsort.cpp index a1c69ff649..40b8afcf63 100644 --- a/src/gromacs/topology/topsort.cpp +++ b/src/gromacs/topology/topsort.cpp @@ -42,6 +42,7 @@ #include "gromacs/topology/ifunc.h" #include "gromacs/topology/topology.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" @@ -174,10 +175,10 @@ gmx_bool gmx_mtop_bondeds_free_energy(const gmx_mtop_t *mtop) /* Check perturbed charges for 1-4 interactions */ for (const gmx_molblock_t &molb : mtop->molblock) { - const t_atom *atom = mtop->moltype[molb.type].atoms.atom; - const t_ilist &il = mtop->moltype[molb.type].ilist[F_LJ14]; - const int *ia = il.iatoms; - for (int i = 0; i < il.nr; i += 3) + const t_atom *atom = mtop->moltype[molb.type].atoms.atom; + const InteractionList &il = mtop->moltype[molb.type].ilist[F_LJ14]; + gmx::ArrayRef ia = il.iatoms; + for (int i = 0; i < il.size(); i += 3) { if (atom[ia[i+1]].q != atom[ia[i+1]].qB || atom[ia[i+2]].q != atom[ia[i+2]].qB) -- 2.22.0