From: Berk Hess Date: Mon, 2 Dec 2019 20:11:02 +0000 (+0100) Subject: Add InteractionDefinitions X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=b16cc028cf9da6ae379f1ecf9586e01fa680d496;p=alexxy%2Fgromacs.git Add InteractionDefinitions This replaces t_idef in gmx_localtop_t. Note that t_idef is still used in t_topology. Extends InteractList with functionality to add entries in a simple and clean way. Change-Id: Ib59012332b6ba09df1de5dc696438c9aba4e2f2d --- diff --git a/src/gromacs/domdec/domdec.h b/src/gromacs/domdec/domdec.h index 39f8bd0b89..51aba44e5c 100644 --- a/src/gromacs/domdec/domdec.h +++ b/src/gromacs/domdec/domdec.h @@ -292,13 +292,6 @@ void dd_make_local_top(struct gmx_domdec_t* dd, /*! \brief Sort ltop->ilist when we are doing free energy. */ void dd_sort_local_top(gmx_domdec_t* dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop); -/*! \brief Initialize local topology - * - * \param[in] top_global Reference to global topology. - * \param[in,out] top Pointer to new local topology - */ -void dd_init_local_top(const gmx_mtop_t& top_global, gmx_localtop_t* top); - /*! \brief Construct local state */ void dd_init_local_state(struct gmx_domdec_t* dd, const t_state* state_global, t_state* local_state); diff --git a/src/gromacs/domdec/domdec_constraints.cpp b/src/gromacs/domdec/domdec_constraints.cpp index ec39ab1494..b11ef6bdb3 100644 --- a/src/gromacs/domdec/domdec_constraints.cpp +++ b/src/gromacs/domdec/domdec_constraints.cpp @@ -94,8 +94,8 @@ struct gmx_domdec_constraints_t std::unique_ptr> ga2la; /**< Global to local communicated constraint atom only index */ /* Multi-threading stuff */ - int nthread; /**< Number of threads used for DD constraint setup */ - std::vector ils; /**< Constraint ilist working arrays, size \p nthread */ + int nthread; /**< Number of threads used for DD constraint setup */ + std::vector ils; /**< Constraint ilist working arrays, size \p nthread */ /* Buffers for requesting atoms */ std::vector> requestedGlobalAtomIndices; /**< Buffers for requesting global atom indices, one per thread */ @@ -155,7 +155,7 @@ static void walk_out(int con, gmx_bool bHomeConnect, gmx_domdec_constraints_t* dc, gmx_domdec_specat_comm_t* dcc, - t_ilist* il_local, + InteractionList* il_local, std::vector* ireq) { if (!dc->gc_req[con_offset + con]) @@ -163,35 +163,32 @@ static void walk_out(int con, /* Add this non-home constraint to the list */ dc->con_gl.push_back(con_offset + con); dc->con_nlocat.push_back(bHomeConnect ? 1 : 0); - dc->gc_req[con_offset + con] = true; - if (il_local->nr + 3 > il_local->nalloc) - { - il_local->nalloc = over_alloc_dd(il_local->nr + 3); - srenew(il_local->iatoms, il_local->nalloc); - } - const int* iap = constr_iatomptr(ia1, ia2, con); - il_local->iatoms[il_local->nr++] = iap[0]; - const int a1_gl = offset + iap[1]; - const int a2_gl = offset + iap[2]; + dc->gc_req[con_offset + con] = true; + const int* iap = constr_iatomptr(ia1, ia2, con); + const int parameterType = iap[0]; + const int a1_gl = offset + iap[1]; + const int a2_gl = offset + iap[2]; + std::array atoms; /* The following indexing code can probably be optizimed */ if (const int* a_loc = ga2la.findHome(a1_gl)) { - il_local->iatoms[il_local->nr++] = *a_loc; + atoms[0] = *a_loc; } else { /* We set this index later */ - il_local->iatoms[il_local->nr++] = -a1_gl - 1; + atoms[0] = -a1_gl - 1; } if (const int* a_loc = ga2la.findHome(a2_gl)) { - il_local->iatoms[il_local->nr++] = *a_loc; + atoms[1] = *a_loc; } else { /* We set this index later */ - il_local->iatoms[il_local->nr++] = -a2_gl - 1; + atoms[1] = -a2_gl - 1; } + il_local->push_back(parameterType, atoms); dc->ncon++; } /* Check to not ask for the same atom more than once */ @@ -239,7 +236,7 @@ static void atoms_to_settles(gmx_domdec_t* dd, gmx::ArrayRef> at2settle_mt, int cg_start, int cg_end, - t_ilist* ils_local, + InteractionList* ils_local, std::vector* ireq) { const gmx_ga2la_t& ga2la = *dd->ga2la; @@ -282,23 +279,17 @@ static void atoms_to_settles(gmx_domdec_t* dd, if (bAssign) { - if (ils_local->nr + 1 + nral > ils_local->nalloc) - { - ils_local->nalloc = over_alloc_dd(ils_local->nr + 1 + nral); - srenew(ils_local->iatoms, ils_local->nalloc); - } - - ils_local->iatoms[ils_local->nr++] = ia1[settle * 4]; - + const int parameterType = ia1[settle * 4]; + std::array atoms; for (int sa = 0; sa < nral; sa++) { if (const int* a_loc = ga2la.findHome(a_gls[sa])) { - ils_local->iatoms[ils_local->nr++] = *a_loc; + atoms[sa] = *a_loc; } else { - ils_local->iatoms[ils_local->nr++] = -a_gls[sa] - 1; + atoms[sa] = -a_gls[sa] - 1; /* Add this non-home atom to the list */ ireq->push_back(a_gls[sa]); /* A check on double atom requests is @@ -306,6 +297,7 @@ static void atoms_to_settles(gmx_domdec_t* dd, */ } } + ils_local->push_back(parameterType, atoms); } } } @@ -318,7 +310,7 @@ static void atoms_to_constraints(gmx_domdec_t* dd, const int* cginfo, gmx::ArrayRef> at2con_mt, int nrec, - t_ilist* ilc_local, + InteractionList* ilc_local, std::vector* ireq) { gmx_domdec_constraints_t* dc = dd->constraints; @@ -373,15 +365,12 @@ static void atoms_to_constraints(gmx_domdec_t* dd, { dc->con_gl.push_back(con_offset + con); dc->con_nlocat.push_back(2); - if (ilc_local->nr + 3 > ilc_local->nalloc) - { - ilc_local->nalloc = over_alloc_dd(ilc_local->nr + 3); - srenew(ilc_local->iatoms, ilc_local->nalloc); - } - const int b_lo = *a_loc; - ilc_local->iatoms[ilc_local->nr++] = iap[0]; - ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? a : b_lo); - ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? b_lo : a); + const int b_lo = *a_loc; + const int parameterType = iap[0]; + std::array atoms; + atoms[0] = (a_gl == iap[1] ? a : b_lo); + atoms[1] = (a_gl == iap[1] ? b_lo : a); + ilc_local->push_back(parameterType, atoms); dc->ncon++; nhome++; } @@ -413,19 +402,18 @@ static void atoms_to_constraints(gmx_domdec_t* dd, } } -int dd_make_local_constraints(gmx_domdec_t* dd, - int at_start, - const struct gmx_mtop_t* mtop, - const int* cginfo, - gmx::Constraints* constr, - int nrec, - t_ilist* il_local) +int dd_make_local_constraints(gmx_domdec_t* dd, + int at_start, + const struct gmx_mtop_t* mtop, + const int* cginfo, + gmx::Constraints* constr, + int nrec, + gmx::ArrayRef il_local) { gmx_domdec_constraints_t* dc; - t_ilist * ilc_local, *ils_local; + InteractionList * ilc_local, *ils_local; gmx::HashedMap* ga2la_specat; int at_end, i, j; - t_iatom* iap; // This code should not be called unless this condition is true, // because that's the only time init_domdec_constraints is @@ -447,10 +435,10 @@ int dd_make_local_constraints(gmx_domdec_t* dd, ilc_local = &il_local[F_CONSTR]; ils_local = &il_local[F_SETTLE]; - dc->ncon = 0; - ilc_local->nr = 0; + dc->ncon = 0; gmx::ArrayRef> at2con_mt; std::vector* ireq = nullptr; + ilc_local->clear(); if (dd->constraint_comm) { // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr? @@ -466,8 +454,8 @@ int dd_make_local_constraints(gmx_domdec_t* dd, { // TODO Perhaps gmx_domdec_constraints_t should keep a valid constr? GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object"); - at2settle_mt = constr->atom2settle_moltype(); - ils_local->nr = 0; + at2settle_mt = constr->atom2settle_moltype(); + ils_local->clear(); } if (at2settle_mt.empty()) @@ -495,8 +483,8 @@ int dd_make_local_constraints(gmx_domdec_t* dd, if (thread >= t0_set) { - int cg0, cg1; - t_ilist* ilst; + int cg0, cg1; + InteractionList* ilst; /* Distribute the settle check+assignments over * dc->nthread or dc->nthread-1 threads. @@ -512,7 +500,7 @@ int dd_make_local_constraints(gmx_domdec_t* dd, { ilst = &dc->ils[thread]; } - ilst->nr = 0; + ilst->clear(); std::vector& ireqt = dc->requestedGlobalAtomIndices[thread]; if (thread > 0) @@ -529,22 +517,9 @@ int dd_make_local_constraints(gmx_domdec_t* dd, /* Combine the generate settles and requested indices */ for (int thread = 1; thread < dc->nthread; thread++) { - t_ilist* ilst; - int ia; - if (thread > t0_set) { - ilst = &dc->ils[thread]; - if (ils_local->nr + ilst->nr > ils_local->nalloc) - { - ils_local->nalloc = over_alloc_large(ils_local->nr + ilst->nr); - srenew(ils_local->iatoms, ils_local->nalloc); - } - for (ia = 0; ia < ilst->nr; ia++) - { - ils_local->iatoms[ils_local->nr + ia] = ilst->iatoms[ia]; - } - ils_local->nr += ilst->nr; + ils_local->append(dc->ils[thread]); } const std::vector& ireqt = dc->requestedGlobalAtomIndices[thread]; @@ -553,7 +528,7 @@ int dd_make_local_constraints(gmx_domdec_t* dd, if (debug) { - fprintf(debug, "Settles: total %3d\n", ils_local->nr / 4); + fprintf(debug, "Settles: total %3d\n", ils_local->size() / 4); } } @@ -568,9 +543,9 @@ int dd_make_local_constraints(gmx_domdec_t* dd, ga2la_specat = dd->constraints->ga2la.get(); nral1 = 1 + NRAL(F_CONSTR); - for (i = 0; i < ilc_local->nr; i += nral1) + for (i = 0; i < ilc_local->size(); i += nral1) { - iap = ilc_local->iatoms + i; + int* iap = ilc_local->iatoms.data() + i; for (j = 1; j < nral1; j++) { if (iap[j] < 0) @@ -583,9 +558,9 @@ int dd_make_local_constraints(gmx_domdec_t* dd, } nral1 = 1 + NRAL(F_SETTLE); - for (i = 0; i < ils_local->nr; i += nral1) + for (i = 0; i < ils_local->size(); i += nral1) { - iap = ils_local->iatoms + i; + int* iap = ils_local->iatoms.data() + i; for (j = 1; j < nral1; j++) { if (iap[j] < 0) diff --git a/src/gromacs/domdec/domdec_constraints.h b/src/gromacs/domdec/domdec_constraints.h index 2e23e16ced..de25e404a0 100644 --- a/src/gromacs/domdec/domdec_constraints.h +++ b/src/gromacs/domdec/domdec_constraints.h @@ -46,6 +46,8 @@ #ifndef GMX_DOMDEC_DOMDEC_CONSTRAINTS_H #define GMX_DOMDEC_DOMDEC_CONSTRAINTS_H +#include "gromacs/utility/arrayref.h" + namespace gmx { class Constraints; @@ -53,19 +55,19 @@ class Constraints; struct gmx_domdec_t; struct gmx_mtop_t; -struct t_ilist; +struct InteractionList; /*! \brief Clears the local indices for the constraint communication setup */ void dd_clear_local_constraint_indices(gmx_domdec_t* dd); /*! \brief Sets up communication and atom indices for all local+connected constraints */ -int dd_make_local_constraints(struct gmx_domdec_t* dd, - int at_start, - const struct gmx_mtop_t* mtop, - const int* cginfo, - gmx::Constraints* constr, - int nrec, - struct t_ilist* il_local); +int dd_make_local_constraints(struct gmx_domdec_t* dd, + int at_start, + const struct gmx_mtop_t* mtop, + const int* cginfo, + gmx::Constraints* constr, + int nrec, + gmx::ArrayRef il_local); /*! \brief Initializes the data structures for constraint communication */ void init_domdec_constraints(gmx_domdec_t* dd, const gmx_mtop_t* mtop); diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp index f2176f2828..d5de5f5bdd 100644 --- a/src/gromacs/domdec/domdec_topology.cpp +++ b/src/gromacs/domdec/domdec_topology.cpp @@ -113,7 +113,13 @@ struct MolblockIndices /*! \brief Struct for thread local work data for local topology generation */ struct thread_work_t { - t_idef idef; /**< Partial local topology */ + /*! \brief Constructor + * + * \param[in] ffparams The interaction parameters, the lifetime of the created object should not exceed the lifetime of the passed parameters + */ + thread_work_t(const gmx_ffparams_t& ffparams) : idef(ffparams) {} + + InteractionDefinitions idef; /**< Partial local topology */ std::unique_ptr vsitePbc; /**< vsite PBC structure */ int nbonded; /**< The number of bondeds in this struct */ ListOfLists excl; /**< List of exclusions */ @@ -182,15 +188,15 @@ static gmx_bool dd_check_ftype(int ftype, gmx_bool bBCheck, gmx_bool bConstr, gm * * \note This function needs to be called on all ranks (contains a global summation) */ -static std::string print_missing_interactions_mb(t_commrec* cr, - const gmx_reverse_top_t* rt, - const char* moltypename, - const reverse_ilist_t* ril, - int a_start, - int a_end, - int nat_mol, - int nmol, - const t_idef* idef) +static std::string print_missing_interactions_mb(t_commrec* cr, + const gmx_reverse_top_t* rt, + const char* moltypename, + const reverse_ilist_t* ril, + int a_start, + int a_end, + int nat_mol, + int nmol, + const InteractionDefinitions* idef) { int* assigned; int nril_mol = ril->index[nat_mol]; @@ -203,10 +209,10 @@ static std::string print_missing_interactions_mb(t_commrec* cr, { if (dd_check_ftype(ftype, rt->bBCheck, rt->bConstr, rt->bSettle)) { - int nral = NRAL(ftype); - const t_ilist* il = &idef->il[ftype]; - const t_iatom* ia = il->iatoms; - for (int i = 0; i < il->nr; i += 1 + nral) + int nral = NRAL(ftype); + const InteractionList* il = &idef->il[ftype]; + const int* ia = il->iatoms.data(); + for (int i = 0; i < il->size(); i += 1 + nral) { int a0 = gatindex[ia[1]]; /* Check if this interaction is in @@ -310,10 +316,10 @@ static std::string print_missing_interactions_mb(t_commrec* cr, } /*! \brief Help print error output when interactions are missing */ -static void print_missing_interactions_atoms(const gmx::MDLogger& mdlog, - t_commrec* cr, - const gmx_mtop_t* mtop, - const t_idef* idef) +static void print_missing_interactions_atoms(const gmx::MDLogger& mdlog, + t_commrec* cr, + const gmx_mtop_t* mtop, + const InteractionDefinitions* idef) { const gmx_reverse_top_t* rt = cr->dd->reverse_top; @@ -356,7 +362,7 @@ void dd_print_missing_interactions(const gmx::MDLogger& mdlog, for (ftype = 0; ftype < F_NRE; ftype++) { nral = NRAL(ftype); - cl[ftype] = top_local->idef.il[ftype].nr / (1 + nral); + cl[ftype] = top_local->idef.il[ftype].size() / (1 + nral); } gmx_sumi(F_NRE, cl, cr); @@ -696,7 +702,10 @@ static gmx_reverse_top_t make_reverse_top(const gmx_mtop_t* mtop, rt.mbi.push_back(mbi); } - rt.th_work.resize(gmx_omp_nthreads_get(emntDomdec)); + for (int th = 0; th < gmx_omp_nthreads_get(emntDomdec); th++) + { + rt.th_work.emplace_back(mtop->ffparams); + } return rt; } @@ -773,18 +782,8 @@ static inline void add_ifunc_for_vsites(t_iatom* tiatoms, int a_gl, int a_mol, const t_iatom* iatoms, - t_ilist* il) + InteractionList* il) { - t_iatom* liatoms; - - if (il->nr + 1 + nral > il->nalloc) - { - il->nalloc = over_alloc_large(il->nr + 1 + nral); - srenew(il->iatoms, il->nalloc); - } - liatoms = il->iatoms + il->nr; - il->nr += 1 + nral; - /* Copy the type */ tiatoms[0] = iatoms[0]; @@ -814,117 +813,85 @@ static inline void add_ifunc_for_vsites(t_iatom* tiatoms, // Note that ga2la_get_home always sets the third parameter if // it returns TRUE } - for (int k = 0; k < 1 + nral; k++) - { - liatoms[k] = tiatoms[k]; - } -} - -/*! \brief Store a bonded interaction at the end of \p il */ -static inline void add_ifunc(int nral, const t_iatom* tiatoms, t_ilist* il) -{ - t_iatom* liatoms; - int k; - - if (il->nr + 1 + nral > il->nalloc) - { - il->nalloc = over_alloc_large(il->nr + 1 + nral); - srenew(il->iatoms, il->nalloc); - } - liatoms = il->iatoms + il->nr; - for (k = 0; k <= nral; k++) - { - liatoms[k] = tiatoms[k]; - } - il->nr += 1 + nral; + il->push_back(tiatoms[0], nral, tiatoms + 1); } /*! \brief Store a position restraint in idef and iatoms, complex because the parameters are different for each entry */ -static void add_posres(int mol, - int a_mol, - int numAtomsInMolecule, - const gmx_molblock_t* molb, - t_iatom* iatoms, - const t_iparams* ip_in, - t_idef* idef) +static void add_posres(int mol, + int a_mol, + int numAtomsInMolecule, + const gmx_molblock_t* molb, + t_iatom* iatoms, + const t_iparams* ip_in, + InteractionDefinitions* idef) { - int n, a_molb; - t_iparams* ip; - /* This position restraint has not been added yet, * so it's index is the current number of position restraints. */ - n = idef->il[F_POSRES].nr / 2; - if (n + 1 > idef->iparams_posres_nalloc) - { - idef->iparams_posres_nalloc = over_alloc_dd(n + 1); - srenew(idef->iparams_posres, idef->iparams_posres_nalloc); - } - ip = &idef->iparams_posres[n]; - /* Copy the force constants */ - *ip = ip_in[iatoms[0]]; + const int n = idef->il[F_POSRES].size() / 2; /* Get the position restraint coordinates from the molblock */ - a_molb = mol * numAtomsInMolecule + a_mol; + const int a_molb = mol * numAtomsInMolecule + a_mol; GMX_ASSERT(a_molb < ssize(molb->posres_xA), "We need a sufficient number of position restraint coordinates"); - ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX]; - ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY]; - ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ]; + /* Copy the force constants */ + t_iparams ip = ip_in[iatoms[0]]; + ip.posres.pos0A[XX] = molb->posres_xA[a_molb][XX]; + ip.posres.pos0A[YY] = molb->posres_xA[a_molb][YY]; + ip.posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ]; if (!molb->posres_xB.empty()) { - ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX]; - ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY]; - ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ]; + ip.posres.pos0B[XX] = molb->posres_xB[a_molb][XX]; + ip.posres.pos0B[YY] = molb->posres_xB[a_molb][YY]; + ip.posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ]; } else { - ip->posres.pos0B[XX] = ip->posres.pos0A[XX]; - ip->posres.pos0B[YY] = ip->posres.pos0A[YY]; - ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ]; + ip.posres.pos0B[XX] = ip.posres.pos0A[XX]; + ip.posres.pos0B[YY] = ip.posres.pos0A[YY]; + ip.posres.pos0B[ZZ] = ip.posres.pos0A[ZZ]; } - /* Set the parameter index for idef->iparams_posre */ + /* Set the parameter index for idef->iparams_posres */ iatoms[0] = n; + idef->iparams_posres.push_back(ip); + + GMX_ASSERT(int(idef->iparams_posres.size()) == n + 1, + "The index of the parameter type should match n"); } /*! \brief Store a flat-bottomed position restraint in idef and iatoms, complex because the parameters are different for each entry */ -static void add_fbposres(int mol, - int a_mol, - int numAtomsInMolecule, - const gmx_molblock_t* molb, - t_iatom* iatoms, - const t_iparams* ip_in, - t_idef* idef) +static void add_fbposres(int mol, + int a_mol, + int numAtomsInMolecule, + const gmx_molblock_t* molb, + t_iatom* iatoms, + const t_iparams* ip_in, + InteractionDefinitions* idef) { - int n, a_molb; - t_iparams* ip; - /* This flat-bottom position restraint has not been added yet, * so it's index is the current number of position restraints. */ - n = idef->il[F_FBPOSRES].nr / 2; - if (n + 1 > idef->iparams_fbposres_nalloc) - { - idef->iparams_fbposres_nalloc = over_alloc_dd(n + 1); - srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc); - } - ip = &idef->iparams_fbposres[n]; - /* Copy the force constants */ - *ip = ip_in[iatoms[0]]; + const int n = idef->il[F_FBPOSRES].size() / 2; /* Get the position restraint coordinats from the molblock */ - a_molb = mol * numAtomsInMolecule + a_mol; + const int a_molb = mol * numAtomsInMolecule + a_mol; GMX_ASSERT(a_molb < ssize(molb->posres_xA), "We need a sufficient number of position restraint coordinates"); + /* Copy the force constants */ + t_iparams ip = ip_in[iatoms[0]]; /* Take reference positions from A position of normal posres */ - ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX]; - ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY]; - ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ]; + ip.fbposres.pos0[XX] = molb->posres_xA[a_molb][XX]; + ip.fbposres.pos0[YY] = molb->posres_xA[a_molb][YY]; + ip.fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ]; /* Note: no B-type for flat-bottom posres */ - /* Set the parameter index for idef->iparams_posre */ + /* Set the parameter index for idef->iparams_fbposres */ iatoms[0] = n; + idef->iparams_fbposres.push_back(ip); + + GMX_ASSERT(int(idef->iparams_fbposres.size()) == n + 1, + "The index of the parameter type should match n"); } /*! \brief Store a virtual site interaction, complex because of PBC and recursion */ @@ -938,7 +905,7 @@ static void add_vsite(const gmx_ga2la_t& ga2la, int a_gl, int a_mol, const t_iatom* iatoms, - t_idef* idef) + InteractionDefinitions* idef) { int k; t_iatom tiatoms[1 + MAXATOMLIST]; @@ -999,7 +966,7 @@ static real dd_dist2(t_pbc* pbc_null, const rvec* x, const int i, int j) } /*! \brief Append t_idef structures 1 to nsrc in src to *dest */ -static void combine_idef(t_idef* dest, gmx::ArrayRef src) +static void combine_idef(InteractionDefinitions* dest, gmx::ArrayRef src) { int ftype; @@ -1008,67 +975,45 @@ static void combine_idef(t_idef* dest, gmx::ArrayRef src) int n = 0; for (gmx::index s = 1; s < src.ssize(); s++) { - n += src[s].idef.il[ftype].nr; + n += src[s].idef.il[ftype].size(); } if (n > 0) { - t_ilist* ild; - - ild = &dest->il[ftype]; - - if (ild->nr + n > ild->nalloc) - { - ild->nalloc = over_alloc_large(ild->nr + n); - srenew(ild->iatoms, ild->nalloc); - } - for (gmx::index s = 1; s < src.ssize(); s++) { - const t_ilist& ils = src[s].idef.il[ftype]; - - for (int i = 0; i < ils.nr; i++) - { - ild->iatoms[ild->nr + i] = ils.iatoms[i]; - } - - ild->nr += ils.nr; + dest->il[ftype].append(src[s].idef.il[ftype]); } /* Position restraints need an additional treatment */ if (ftype == F_POSRES || ftype == F_FBPOSRES) { - int nposres = dest->il[ftype].nr / 2; - // TODO: Simplify this code using std::vector - t_iparams*& iparams_dest = + int nposres = dest->il[ftype].size() / 2; + std::vector& iparams_dest = (ftype == F_POSRES ? dest->iparams_posres : dest->iparams_fbposres); - int& posres_nalloc = (ftype == F_POSRES ? dest->iparams_posres_nalloc - : dest->iparams_fbposres_nalloc); - if (nposres > posres_nalloc) - { - posres_nalloc = over_alloc_large(nposres); - srenew(iparams_dest, posres_nalloc); - } /* Set nposres to the number of original position restraints in dest */ for (gmx::index s = 1; s < src.ssize(); s++) { - nposres -= src[s].idef.il[ftype].nr / 2; + nposres -= src[s].idef.il[ftype].size() / 2; } for (gmx::index s = 1; s < src.ssize(); s++) { - const t_iparams* iparams_src = (ftype == F_POSRES ? src[s].idef.iparams_posres - : src[s].idef.iparams_fbposres); + const std::vector& iparams_src = + (ftype == F_POSRES ? src[s].idef.iparams_posres : src[s].idef.iparams_fbposres); + iparams_dest.insert(iparams_dest.end(), iparams_src.begin(), iparams_src.end()); - for (int i = 0; i < src[s].idef.il[ftype].nr / 2; i++) + /* Correct the indices into iparams_posres */ + for (int i = 0; i < src[s].idef.il[ftype].size() / 2; i++) { /* Correct the index into iparams_posres */ dest->il[ftype].iatoms[nposres * 2] = nposres; - /* Copy the position restraint force parameters */ - iparams_dest[nposres] = iparams_src[i]; nposres++; } } + GMX_RELEASE_ASSERT( + int(iparams_dest.size()) == nposres, + "The number of parameters should match the number of restraints"); } } } @@ -1096,7 +1041,7 @@ static inline void check_assign_interactions_atom(int i, t_pbc* pbc_null, rvec* cg_cm, const t_iparams* ip_in, - t_idef* idef, + InteractionDefinitions* idef, int iz, gmx_bool bBCheck, int* nbonded_local) @@ -1270,7 +1215,7 @@ static inline void check_assign_interactions_atom(int i, if (bUse) { /* Add this interaction to the local topology */ - add_ifunc(nral, tiatoms, &idef->il[ftype]); + idef->il[ftype].push_back(tiatoms[0], nral, tiatoms + 1); /* Sum so we can check in global_stat * if we have everything. */ @@ -1299,7 +1244,7 @@ static int make_bondeds_zone(gmx_domdec_t* dd, t_pbc* pbc_null, rvec* cg_cm, const t_iparams* ip_in, - t_idef* idef, + InteractionDefinitions* idef, int izone, const gmx::Range& atomRange) { @@ -1416,32 +1361,20 @@ static void make_exclusions_zone(gmx_domdec_t* dd, "The number of exclusion list should match the number of atoms in the range"); } -/*! \brief Clear a t_idef struct */ -static void clear_idef(t_idef* idef) -{ - int ftype; - - /* Clear the counts */ - for (ftype = 0; ftype < F_NRE; ftype++) - { - idef->il[ftype].nr = 0; - } -} - /*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls */ -static int make_local_bondeds_excls(gmx_domdec_t* dd, - gmx_domdec_zones_t* zones, - const gmx_mtop_t* mtop, - const int* cginfo, - gmx_bool bRCheckMB, - ivec rcheck, - gmx_bool bRCheck2B, - real rc, - t_pbc* pbc_null, - rvec* cg_cm, - t_idef* idef, - ListOfLists* lexcls, - int* excl_count) +static int make_local_bondeds_excls(gmx_domdec_t* dd, + gmx_domdec_zones_t* zones, + const gmx_mtop_t* mtop, + const int* cginfo, + gmx_bool bRCheckMB, + ivec rcheck, + gmx_bool bRCheck2B, + real rc, + t_pbc* pbc_null, + rvec* cg_cm, + InteractionDefinitions* idef, + ListOfLists* lexcls, + int* excl_count) { int nzone_bondeds; int cg0, cg1; @@ -1469,7 +1402,7 @@ static int make_local_bondeds_excls(gmx_domdec_t* dd, rc2 = rc * rc; /* Clear the counts */ - clear_idef(idef); + idef->clear(); nbonded_local = 0; lexcls->clear(); @@ -1486,8 +1419,8 @@ static int make_local_bondeds_excls(gmx_domdec_t* dd, { try { - int cg0t, cg1t; - t_idef* idef_t; + int cg0t, cg1t; + InteractionDefinitions* idef_t; cg0t = cg0 + ((cg1 - cg0) * thread) / numThreads; cg1t = cg0 + ((cg1 - cg0) * (thread + 1)) / numThreads; @@ -1499,12 +1432,12 @@ static int make_local_bondeds_excls(gmx_domdec_t* dd, else { idef_t = &rt->th_work[thread].idef; - clear_idef(idef_t); + idef_t->clear(); } rt->th_work[thread].nbonded = make_bondeds_zone( dd, zones, mtop->molblock, bRCheckMB, rcheck, bRCheck2B, rc2, pbc_null, - cg_cm, idef->iparams, idef_t, izone, gmx::Range(cg0t, cg1t)); + cg_cm, idef->iparams.data(), idef_t, izone, gmx::Range(cg0t, cg1t)); if (izone < numIZonesForExclusions) { @@ -1642,8 +1575,6 @@ void dd_make_local_top(gmx_domdec_t* dd, * we can only do this when we have the charge arrays. */ ltop->idef.ilsort = ilsortUNKNOWN; - - ltop->atomtypes = mtop.atomtypes; } void dd_sort_local_top(gmx_domdec_t* dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop) @@ -1658,21 +1589,6 @@ void dd_sort_local_top(gmx_domdec_t* dd, const t_mdatoms* mdatoms, gmx_localtop_ } } -void dd_init_local_top(const gmx_mtop_t& top_global, gmx_localtop_t* top) -{ - /* TODO: Get rid of the const casts below, e.g. by using a reference */ - top->idef.ntypes = top_global.ffparams.numTypes(); - top->idef.atnr = top_global.ffparams.atnr; - top->idef.functype = const_cast(top_global.ffparams.functype.data()); - top->idef.iparams = const_cast(top_global.ffparams.iparams.data()); - top->idef.fudgeQQ = top_global.ffparams.fudgeQQ; - top->idef.cmap_grid = new gmx_cmap_t; - *top->idef.cmap_grid = top_global.ffparams.cmap_grid; - - top->idef.ilsort = ilsortUNKNOWN; - top->useInDomainDecomp_ = true; -} - void dd_init_local_state(gmx_domdec_t* dd, const t_state* state_global, t_state* state_local) { int buf[NITEM_DD_INIT_LOCAL_STATE]; @@ -1998,7 +1914,7 @@ static bool moltypeHasVsite(const gmx_moltype_t& molt) bool hasVsite = false; for (int i = 0; i < F_NRE; i++) { - if ((interaction_function[i].flags & IF_VSITE) && molt.ilist[i].size() > 0) + if ((interaction_function[i].flags & IF_VSITE) && !molt.ilist[i].empty()) { hasVsite = true; } @@ -2045,18 +1961,7 @@ static void getWholeMoleculeCoordinates(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.data(), ilist, PbcType::No, + construct_vsites(nullptr, xs, 0.0, nullptr, ffparams->iparams, molt->ilist, PbcType::No, TRUE, nullptr, nullptr); } } diff --git a/src/gromacs/domdec/domdec_vsite.cpp b/src/gromacs/domdec/domdec_vsite.cpp index f202863c2a..7b8c8f903f 100644 --- a/src/gromacs/domdec/domdec_vsite.cpp +++ b/src/gromacs/domdec/domdec_vsite.cpp @@ -103,7 +103,7 @@ void dd_clear_local_vsite_indices(gmx_domdec_t* dd) } } -int dd_make_local_vsites(gmx_domdec_t* dd, int at_start, t_ilist* lil) +int dd_make_local_vsites(gmx_domdec_t* dd, int at_start, gmx::ArrayRef lil) { std::vector& ireq = dd->vsite_requestedGlobalAtomIndices; gmx::HashedMap* ga2la_specat = dd->ga2la_vsite; @@ -114,11 +114,11 @@ int dd_make_local_vsites(gmx_domdec_t* dd, int at_start, t_ilist* lil) { if (interaction_function[ftype].flags & IF_VSITE) { - const int nral = NRAL(ftype); - const t_ilist& lilf = lil[ftype]; - for (int i = 0; i < lilf.nr; i += 1 + nral) + const int nral = NRAL(ftype); + const InteractionList& lilf = lil[ftype]; + for (int i = 0; i < lilf.size(); i += 1 + nral) { - const t_iatom* iatoms = lilf.iatoms + i; + const int* iatoms = lilf.iatoms.data() + i; /* Check if we have the other atoms */ for (int j = 1; j < 1 + nral; j++) { @@ -152,11 +152,11 @@ int dd_make_local_vsites(gmx_domdec_t* dd, int at_start, t_ilist* lil) { if (interaction_function[ftype].flags & IF_VSITE) { - const int nral = NRAL(ftype); - t_ilist& lilf = lil[ftype]; - for (int i = 0; i < lilf.nr; i += 1 + nral) + const int nral = NRAL(ftype); + InteractionList& lilf = lil[ftype]; + for (int i = 0; i < lilf.size(); i += 1 + nral) { - t_iatom* iatoms = lilf.iatoms + i; + t_iatom* iatoms = lilf.iatoms.data() + i; for (int j = 1; j < 1 + nral; j++) { if (iatoms[j] < 0) diff --git a/src/gromacs/domdec/domdec_vsite.h b/src/gromacs/domdec/domdec_vsite.h index 9ff792e844..3ee206f283 100644 --- a/src/gromacs/domdec/domdec_vsite.h +++ b/src/gromacs/domdec/domdec_vsite.h @@ -46,14 +46,16 @@ #ifndef GMX_DOMDEC_DOMDEC_VSITE_H #define GMX_DOMDEC_DOMDEC_VSITE_H -struct t_ilist; +#include "gromacs/utility/arrayref.h" + struct gmx_domdec_t; +struct InteractionList; /*! \brief Clears the local indices for the virtual site communication setup */ void dd_clear_local_vsite_indices(struct gmx_domdec_t* dd); /*! \brief Sets up communication and atom indices for all local vsites */ -int dd_make_local_vsites(struct gmx_domdec_t* dd, int at_start, struct t_ilist* lil); +int dd_make_local_vsites(struct gmx_domdec_t* dd, int at_start, gmx::ArrayRef lil); /*! \brief Initializes the data structures for virtual site communication */ void init_domdec_vsites(struct gmx_domdec_t* dd, int n_intercg_vsite); diff --git a/src/gromacs/domdec/mdsetup.cpp b/src/gromacs/domdec/mdsetup.cpp index 79fe29abbc..b5a0751197 100644 --- a/src/gromacs/domdec/mdsetup.cpp +++ b/src/gromacs/domdec/mdsetup.cpp @@ -119,7 +119,7 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, { GMX_ASSERT(graph != nullptr, "We use a graph with PBC (no periodic mols) and without DD"); - *graph = mk_graph(nullptr, &(top->idef), 0, top_global.natoms, FALSE, FALSE); + *graph = mk_graph(nullptr, top->idef, 0, top_global.natoms, FALSE, FALSE); } else if (graph != nullptr) { @@ -150,6 +150,6 @@ void mdAlgorithmsSetupAtomData(const t_commrec* cr, if (constr) { - constr->setConstraints(*top, *mdatoms); + constr->setConstraints(top, *mdatoms); } } diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index 4266291098..4078e9cb85 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -2108,7 +2108,7 @@ static void do_ilists(gmx::ISerializer* serializer, InteractionLists* ilists, in else { do_ilist(serializer, &ilist); - if (file_version < 78 && j == F_SETTLE && ilist.size() > 0) + if (file_version < 78 && j == F_SETTLE && !ilist.empty()) { add_settle_atoms(&ilist); } @@ -2509,7 +2509,7 @@ static void set_disres_npair(gmx_mtop_t* mtop) { const InteractionList& il = (*ilist)[F_DISRES]; - if (il.size() > 0) + if (!il.empty()) { gmx::ArrayRef a = il.iatoms; int npair = 0; @@ -3034,7 +3034,7 @@ static void do_tpx_finalize(TpxFileHeader* tpx, t_inputrec* ir, t_state* state, { if (tpx->fileVersion < 57) { - if (mtop->moltype[0].ilist[F_DISRES].size() > 0) + if (!mtop->moltype[0].ilist[F_DISRES].empty()) { ir->eDisre = edrSimple; } diff --git a/src/gromacs/gmxana/gmx_disre.cpp b/src/gromacs/gmxana/gmx_disre.cpp index 046573194a..397b4c1bfe 100644 --- a/src/gromacs/gmxana/gmx_disre.cpp +++ b/src/gromacs/gmxana/gmx_disre.cpp @@ -147,21 +147,20 @@ static void print5(FILE* fp) fprintf(fp, "\n"); } -static void check_viol(FILE* log, - t_ilist* disres, - t_iparams forceparams[], - rvec x[], - rvec4 f[], - t_pbc* pbc, - t_graph* g, - t_dr_result dr[], - int clust_id, - int isize, - const int index[], - real vvindex[], - t_fcdata* fcd) +static void check_viol(FILE* log, + const InteractionList& disres, + gmx::ArrayRef forceparams, + rvec x[], + rvec4 f[], + t_pbc* pbc, + t_graph* g, + t_dr_result dr[], + int clust_id, + int isize, + const int index[], + real vvindex[], + t_fcdata* fcd) { - t_iatom* forceatoms; int i, j, nat, n, type, nviol, ndr, label; real rt, mviol, tviol, viol, lam, dvdl, drt; rvec* fshift; @@ -177,7 +176,7 @@ static void check_viol(FILE* log, { reset5(); } - forceatoms = disres->iatoms; + gmx::ArrayRef forceatoms = disres.iatoms; for (j = 0; (j < isize); j++) { vvindex[j] = 0; @@ -187,7 +186,7 @@ static void check_viol(FILE* log, // The label for a distance restraint should be at most one larger // than the previous label. int label_old = forceparams[forceatoms[0]].disres.label; - for (i = 0; (i < disres->nr); i += nat) + for (i = 0; (i < disres.size()); i += nat) { type = forceatoms[i]; label = forceparams[type].disres.label; @@ -205,7 +204,7 @@ static void check_viol(FILE* log, } // Get offset for label index label_old = forceparams[forceatoms[0]].disres.label; - for (i = 0; (i < disres->nr);) + for (i = 0; (i < disres.size());) { type = forceatoms[i]; n = 0; @@ -217,7 +216,7 @@ static void check_viol(FILE* log, do { n += nat; - } while (((i + n) < disres->nr) + } while (((i + n) < disres.size()) && (forceparams[forceatoms[i + n]].disres.label == label + label_old)); calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i], x, pbc, fcd, nullptr); @@ -235,7 +234,8 @@ static void check_viol(FILE* log, dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[label]; snew(fshift, SHIFTS); - ta_disres(n, &forceatoms[i], forceparams, x, f, fshift, pbc, g, lam, &dvdl, nullptr, fcd, nullptr); + ta_disres(n, &forceatoms[i], forceparams.data(), x, f, fshift, pbc, g, lam, &dvdl, nullptr, + fcd, nullptr); sfree(fshift); viol = fcd->disres.sumviol; @@ -270,7 +270,7 @@ static void check_viol(FILE* log, if (bFirst) { - fprintf(stderr, "\nThere are %d restraints and %d pairs\n", ndr, disres->nr / nat); + fprintf(stderr, "\nThere are %d restraints and %d pairs\n", ndr, disres.size() / nat); bFirst = FALSE; } if (ntop) @@ -371,15 +371,15 @@ static gmx_bool is_core(int i, int isize, const int index[]) return bIC; } -static void dump_stats(FILE* log, - int nsteps, - const t_disresdata& dd, - const t_ilist* disres, - t_iparams ip[], - t_dr_result* dr, - int isize, - int index[], - t_atoms* atoms) +static void dump_stats(FILE* log, + int nsteps, + const t_disresdata& dd, + const InteractionList& disres, + gmx::ArrayRef ip, + t_dr_result* dr, + int isize, + int index[], + t_atoms* atoms) { t_dr_stats* drs; @@ -387,15 +387,15 @@ static void dump_stats(FILE* log, fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n"); snew(drs, dd.nres); const int nra = interaction_function[F_DISRES].nratoms + 1; - for (int j = 0; j < disres->nr; j += nra) + for (int j = 0; j < disres.size(); j += nra) { // Note that the restraint i can be used by multiple pairs - const int i = disres->iatoms[j] - dd.type_min; + const int i = disres.iatoms[j] - dd.type_min; GMX_RELEASE_ASSERT(i >= 0 && i < dd.nres, "The restraint index should be in range"); - drs[i].label = ip[disres->iatoms[j]].disres.label; + drs[i].label = ip[disres.iatoms[j]].disres.label; drs[i].bCore = is_core(drs[i].label, isize, index); - drs[i].up1 = ip[disres->iatoms[j]].disres.up1; + drs[i].up1 = ip[disres.iatoms[j]].disres.up1; drs[i].r = dr->aver1[i] / nsteps; drs[i].rT3 = gmx::invcbrt(dr->aver_3[i] / nsteps); drs[i].rT6 = gmx::invsixthroot(dr->aver_6[i] / nsteps); @@ -404,8 +404,8 @@ static void dump_stats(FILE* log, drs[i].violT6 = std::max(0.0, static_cast(drs[i].rT6 - drs[i].up1)); if (atoms) { - int j1 = disres->iatoms[j + 1]; - int j2 = disres->iatoms[j + 2]; + int j1 = disres.iatoms[j + 1]; + int j2 = disres.iatoms[j + 2]; atoms->pdbinfo[j1].bfac += drs[i].violT3 * 5; atoms->pdbinfo[j2].bfac += drs[i].violT3 * 5; } @@ -422,15 +422,15 @@ static void dump_stats(FILE* log, sfree(drs); } -static void dump_clust_stats(FILE* fp, - const t_disresdata& dd, - const t_ilist* disres, - t_iparams ip[], - t_blocka* clust, - t_dr_result dr[], - char* clust_name[], - int isize, - int index[]) +static void dump_clust_stats(FILE* fp, + const t_disresdata& dd, + const InteractionList& disres, + gmx::ArrayRef ip, + t_blocka* clust, + t_dr_result dr[], + char* clust_name[], + int isize, + int index[]) { int k, nra, mmm = 0; double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6; @@ -467,16 +467,16 @@ static void dump_clust_stats(FILE* fp, for (int j = 0; j < dd.nres; j += nra) { // Note that the restraint i can be used by multiple pairs - const int i = disres->iatoms[j] - dd.type_min; + const int i = disres.iatoms[j] - dd.type_min; if (restraintHasBeenProcessed[i]) { continue; } - drs[i].label = ip[disres->iatoms[j]].disres.label; + drs[i].label = ip[disres.iatoms[j]].disres.label; drs[i].bCore = is_core(drs[i].label, isize, index); - drs[i].up1 = ip[disres->iatoms[j]].disres.up1; + drs[i].up1 = ip[disres.iatoms[j]].disres.up1; drs[i].r = dr[k].aver1[i] / dr[k].nframes; if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i])) { @@ -521,15 +521,15 @@ static void init_dr_res(t_dr_result* dr, int ndr) dr->averv = 0; } -static void dump_disre_matrix(const char* fn, - t_dr_result* dr, - int ndr, - int nsteps, - t_idef* idef, - const gmx_mtop_t* mtop, - real max_dr, - int nlevels, - gmx_bool bThird) +static void dump_disre_matrix(const char* fn, + t_dr_result* dr, + int ndr, + int nsteps, + const InteractionDefinitions& idef, + const gmx_mtop_t* mtop, + real max_dr, + int nlevels, + gmx_bool bThird) { FILE* fp; int* resnr; @@ -572,16 +572,16 @@ static void dump_disre_matrix(const char* fn, snew(matrix[i], n_res); } nratoms = interaction_function[F_DISRES].nratoms; - nra = (idef->il[F_DISRES].nr / (nratoms + 1)); + nra = (idef.il[F_DISRES].size() / (nratoms + 1)); snew(ptr, nra + 1); index = 0; nlabel = 0; ptr[0] = 0; snew(w_dr, ndr); - for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms + 1) + for (i = 0; (i < idef.il[F_DISRES].size()); i += nratoms + 1) { - tp = idef->il[F_DISRES].iatoms[i]; - label = idef->iparams[tp].disres.label; + tp = idef.il[F_DISRES].iatoms[i]; + label = idef.iparams[tp].disres.label; if (label != index) { @@ -611,9 +611,9 @@ static void dump_disre_matrix(const char* fn, { for (j = ptr[i]; (j < ptr[i + 1]); j += nratoms + 1) { - tp = idef->il[F_DISRES].iatoms[j]; - ai = idef->il[F_DISRES].iatoms[j + 1]; - aj = idef->il[F_DISRES].iatoms[j + 2]; + tp = idef.il[F_DISRES].iatoms[j]; + ai = idef.il[F_DISRES].iatoms[j + 1]; + aj = idef.il[F_DISRES].iatoms[j + 2]; ri = resnr[ai]; rj = resnr[aj]; @@ -629,7 +629,7 @@ static void dump_disre_matrix(const char* fn, { fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav); } - rviol = std::max(0.0_real, rav - idef->iparams[tp].disres.up1); + rviol = std::max(0.0_real, rav - idef.iparams[tp].disres.up1); matrix[ri][rj] += w_dr[i] * rviol; matrix[rj][ri] += w_dr[i] * rviol; hi = std::max(hi, matrix[ri][rj]); @@ -699,20 +699,19 @@ int gmx_disre(int argc, char* argv[]) "Use inverse third power averaging or linear for matrix output" } }; - FILE *out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr; - gmx_localtop_t top; - t_fcdata fcd; - t_graph* g; - int i, j, kkk; - t_trxstatus* status; - real t; - rvec * x, *xav = nullptr; - rvec4* f; - matrix box; - gmx_bool bPDB; - int isize; - int * index = nullptr, *ind_fit = nullptr; - char* grpname; + FILE * out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr; + t_fcdata fcd; + t_graph* g; + int i, j, kkk; + t_trxstatus* status; + real t; + rvec * x, *xav = nullptr; + rvec4* f; + matrix box; + gmx_bool bPDB; + int isize; + int * index = nullptr, *ind_fit = nullptr; + char* grpname; t_cluster_ndx* clust = nullptr; t_dr_result dr, *dr_clust = nullptr; char** leg; @@ -772,6 +771,7 @@ int gmx_disre(int argc, char* argv[]) atoms->havePdbInfo = TRUE; } + gmx_localtop_t top(topInfo.mtop()->ffparams); gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO); g = nullptr; @@ -784,7 +784,7 @@ int gmx_disre(int argc, char* argv[]) } else { - g = mk_graph(fplog, &top.idef, 0, ntopatoms, FALSE, FALSE); + g = mk_graph(fplog, top.idef, 0, ntopatoms, FALSE, FALSE); } } @@ -838,7 +838,7 @@ int gmx_disre(int argc, char* argv[]) update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda); if (ir->pbcType != PbcType::No) { - gpbc = gmx_rmpbc_init(&top.idef, ir->pbcType, natoms); + gpbc = gmx_rmpbc_init(top.idef, ir->pbcType, natoms); } j = 0; @@ -867,12 +867,12 @@ int gmx_disre(int argc, char* argv[]) } my_clust = clust->inv_clust[j]; range_check(my_clust, 0, clust->clust->nr); - check_viol(fplog, &(top.idef.il[F_DISRES]), top.idef.iparams, x, f, pbc_null, g, - dr_clust, my_clust, isize, index, vvindex, &fcd); + check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, g, dr_clust, + my_clust, isize, index, vvindex, &fcd); } else { - check_viol(fplog, &(top.idef.il[F_DISRES]), top.idef.iparams, x, f, pbc_null, g, &dr, 0, + check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, g, &dr, 0, isize, index, vvindex, &fcd); } if (bPDB) @@ -916,19 +916,19 @@ int gmx_disre(int argc, char* argv[]) if (clust) { - dump_clust_stats(fplog, fcd.disres, &(top.idef.il[F_DISRES]), top.idef.iparams, - clust->clust, dr_clust, clust->grpname, isize, index); + dump_clust_stats(fplog, fcd.disres, top.idef.il[F_DISRES], top.idef.iparams, clust->clust, + dr_clust, clust->grpname, isize, index); } else { - dump_stats(fplog, j, fcd.disres, &(top.idef.il[F_DISRES]), top.idef.iparams, &dr, isize, - index, bPDB ? atoms.get() : nullptr); + dump_stats(fplog, j, fcd.disres, top.idef.il[F_DISRES], top.idef.iparams, &dr, isize, index, + bPDB ? atoms.get() : nullptr); if (bPDB) { write_sto_conf(opt2fn("-q", NFILE, fnm), "Coloured by average violation in Angstrom", atoms.get(), xav, nullptr, ir->pbcType, box); } - dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres, j, &top.idef, + dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres, j, top.idef, topInfo.mtop(), max_dr, nlevels, bThird); xvgrclose(out); xvgrclose(aver); diff --git a/src/gromacs/gmxana/gmx_nmr.cpp b/src/gromacs/gmxana/gmx_nmr.cpp index ba463e6b17..bc20833ec6 100644 --- a/src/gromacs/gmxana/gmx_nmr.cpp +++ b/src/gromacs/gmxana/gmx_nmr.cpp @@ -200,20 +200,16 @@ static void get_orires_parms(const char* topnm, t_inputrec* ir, int* nor, int* n static int get_bounds(real** bounds, int** index, int** dr_pair, int* npairs, gmx_localtop_t* top) { - t_functype* functype; - t_iparams* ip; - int i, j, k, type, ftype, natom; - t_ilist* disres; - t_iatom* iatom; - real* b; - int * ind, *pair; - int nb, label1; - - functype = top->idef.functype; - ip = top->idef.iparams; + int i, j, k, type, ftype, natom; + real* b; + int * ind, *pair; + int nb, label1; + + gmx::ArrayRef functype = top->idef.functype; + gmx::ArrayRef iparams = top->idef.iparams; /* Count how many distance restraint there are... */ - nb = top->idef.il[F_DISRES].nr; + nb = top->idef.il[F_DISRES].size(); if (nb == 0) { gmx_fatal(FARGS, "No distance restraints in topology!\n"); @@ -226,14 +222,14 @@ static int get_bounds(real** bounds, int** index, int** dr_pair, int* npairs, gm /* Fill the bound array */ nb = 0; - for (i = 0; (i < top->idef.ntypes); i++) + for (gmx::index i = 0; i < functype.ssize(); i++) { ftype = functype[i]; if (ftype == F_DISRES) { - label1 = ip[i].disres.label; - b[nb] = ip[i].disres.up1; + label1 = iparams[i].disres.label; + b[nb] = iparams[i].disres.up1; ind[nb] = label1; nb++; } @@ -241,18 +237,18 @@ static int get_bounds(real** bounds, int** index, int** dr_pair, int* npairs, gm *bounds = b; /* Fill the index array */ - label1 = -1; - disres = &(top->idef.il[F_DISRES]); - iatom = disres->iatoms; - for (i = j = k = 0; (i < disres->nr);) + label1 = -1; + const InteractionList& disres = top->idef.il[F_DISRES]; + gmx::ArrayRef iatom = disres.iatoms; + for (i = j = k = 0; (i < disres.size());) { type = iatom[i]; - ftype = top->idef.functype[type]; + ftype = functype[type]; natom = interaction_function[ftype].nratoms + 1; - if (label1 != top->idef.iparams[type].disres.label) + if (label1 != iparams[type].disres.label) { pair[j] = k; - label1 = top->idef.iparams[type].disres.label; + label1 = iparams[type].disres.label; j++; } k++; @@ -411,13 +407,12 @@ int gmx_nmr(int argc, char* argv[]) FILE /* *out = NULL,*/ *out_disre = nullptr, *fp_pairs = nullptr, *fort = nullptr, *fodt = nullptr, *foten = nullptr; - ener_file_t fp; - int timecheck = 0; - gmx_localtop_t top; - gmx_enxnm_t* enm = nullptr; - t_enxframe fr; - int nre, teller, teller_disre; - int nor = 0, nex = 0, norfr = 0, enx_i = 0; + ener_file_t fp; + int timecheck = 0; + gmx_enxnm_t* enm = nullptr; + t_enxframe fr; + int nre, teller, teller_disre; + int nor = 0, nex = 0, norfr = 0, enx_i = 0; real *bounds = nullptr, *violaver = nullptr, *oobs = nullptr, *orient = nullptr, *odrms = nullptr; int * index = nullptr, *pair = nullptr, norsel = 0, *orsel = nullptr, *or_label = nullptr; int nbounds = 0, npairs; @@ -480,7 +475,8 @@ int gmx_nmr(int argc, char* argv[]) t_inputrec irInstance; t_inputrec* ir = &irInstance; init_enxframe(&fr); - gmx::TopologyInformation topInfo; + gmx::TopologyInformation topInfo; + std::unique_ptr top; if (!bDisRe) { if (bORIRE || bOTEN) @@ -618,9 +614,10 @@ int gmx_nmr(int argc, char* argv[]) { { topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm)); - gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO); + top = std::make_unique(topInfo.mtop()->ffparams); + gmx_mtop_generate_local_top(*topInfo.mtop(), top.get(), ir->efep != efepNO); } - nbounds = get_bounds(&bounds, &index, &pair, &npairs, &top); + nbounds = get_bounds(&bounds, &index, &pair, &npairs, top.get()); snew(violaver, npairs); out_disre = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations", "Time (ps)", "nm", oenv); xvgr_legend(out_disre, 2, drleg, oenv); @@ -661,23 +658,21 @@ int gmx_nmr(int argc, char* argv[]) blk_disre = find_block_id_enxframe(&fr, enxDISRE, nullptr); if (bDisRe && bDRAll && !leg && blk_disre) { - t_iatom* fa; - t_iparams* ip; - - fa = top.idef.il[F_DISRES].iatoms; - ip = top.idef.iparams; + const InteractionList& ilist = top->idef.il[F_DISRES]; + gmx::ArrayRef fa = ilist.iatoms; + const t_iparams* ip = top->idef.iparams.data(); if (blk_disre->nsub != 2 || (blk_disre->sub[0].nr != blk_disre->sub[1].nr)) { gmx_incons("Number of disre sub-blocks not equal to 2"); } ndisre = blk_disre->sub[0].nr; - if (ndisre != top.idef.il[F_DISRES].nr / 3) + if (ndisre != ilist.size() / 3) { gmx_fatal(FARGS, "Number of disre pairs in the energy file (%d) does not match the " "number in the run input file (%d)\n", - ndisre, top.idef.il[F_DISRES].nr / 3); + ndisre, ilist.size() / 3); } snew(pairleg, ndisre); int molb = 0; diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index f010d6e6bc..66d9dbaea8 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -1446,7 +1446,7 @@ static bool haveDecoupledModeInMol(const gmx_moltype_t& molt, gmx::ArrayRef iparams, real massFactorThreshold) { - if (molt.ilist[F_CONSTR].size() == 0 && molt.ilist[F_CONSTRNC].size() == 0) + if (molt.ilist[F_CONSTR].empty() && molt.ilist[F_CONSTRNC].empty()) { return false; } diff --git a/src/gromacs/gmxpreprocess/topio.cpp b/src/gromacs/gmxpreprocess/topio.cpp index 7d33fafda6..377bb4cb33 100644 --- a/src/gromacs/gmxpreprocess/topio.cpp +++ b/src/gromacs/gmxpreprocess/topio.cpp @@ -1102,7 +1102,7 @@ static void generate_qmexcl_moltype(gmx_moltype_t* molt, */ int ftype_connbond = 0; int ind_connbond = 0; - if (molt->ilist[F_CONNBONDS].size() != 0) + if (!molt->ilist[F_CONNBONDS].empty()) { GMX_LOG(logger.info) .asParagraph() diff --git a/src/gromacs/listed_forces/disre.cpp b/src/gromacs/listed_forces/disre.cpp index 51eb75238a..182d9b7932 100644 --- a/src/gromacs/listed_forces/disre.cpp +++ b/src/gromacs/listed_forces/disre.cpp @@ -139,7 +139,7 @@ void init_disres(FILE* fplog, iloop = gmx_mtop_ilistloop_init(mtop); while (const InteractionLists* il = gmx_mtop_ilistloop_next(iloop, &nmol)) { - if (nmol > 1 && (*il)[F_DISRES].size() > 0 && ir->eDisre != edrEnsemble) + if (nmol > 1 && !(*il)[F_DISRES].empty() && ir->eDisre != edrEnsemble) { gmx_fatal(FARGS, "NMR distance restraints with multiple copies of the same molecule are " diff --git a/src/gromacs/listed_forces/gpubonded.h b/src/gromacs/listed_forces/gpubonded.h index 54383450d9..1a231d2c2c 100644 --- a/src/gromacs/listed_forces/gpubonded.h +++ b/src/gromacs/listed_forces/gpubonded.h @@ -59,7 +59,6 @@ struct gmx_enerdata_t; struct gmx_ffparams_t; struct gmx_mtop_t; struct t_forcerec; -struct t_idef; struct t_inputrec; struct gmx_wallcycle; @@ -118,11 +117,12 @@ public: * stage. Copies the bonded interactions assigned to the GPU * to device data structures, and updates device buffers that * may have been updated after search. */ - void updateInteractionListsAndDeviceBuffers(ArrayRef nbnxnAtomOrder, - const t_idef& idef, - void* xqDevice, - DeviceBuffer forceDevice, - DeviceBuffer fshiftDevice); + void updateInteractionListsAndDeviceBuffers(ArrayRef nbnxnAtomOrder, + const InteractionDefinitions& idef, + void* xqDevice, + DeviceBuffer forceDevice, + DeviceBuffer fshiftDevice); + /*! \brief Returns whether there are bonded interactions * assigned to the GPU */ bool haveInteractions() const; diff --git a/src/gromacs/listed_forces/gpubonded_impl.cpp b/src/gromacs/listed_forces/gpubonded_impl.cpp index 18d42c8251..94a2b5d42b 100644 --- a/src/gromacs/listed_forces/gpubonded_impl.cpp +++ b/src/gromacs/listed_forces/gpubonded_impl.cpp @@ -168,7 +168,7 @@ GpuBonded::GpuBonded(const gmx_ffparams_t& /* ffparams */, void* /*streamPtr */, GpuBonded::~GpuBonded() = default; void GpuBonded::updateInteractionListsAndDeviceBuffers(ArrayRef /* nbnxnAtomOrder */, - const t_idef& /* idef */, + const InteractionDefinitions& /* idef */, void* /* xqDevice */, DeviceBuffer /* forceDevice */, DeviceBuffer /* fshiftDevice */) diff --git a/src/gromacs/listed_forces/gpubonded_impl.cu b/src/gromacs/listed_forces/gpubonded_impl.cu index 6a6ead4f71..bc10e76066 100644 --- a/src/gromacs/listed_forces/gpubonded_impl.cu +++ b/src/gromacs/listed_forces/gpubonded_impl.cu @@ -78,13 +78,6 @@ GpuBonded::Impl::Impl(const gmx_ffparams_t& ffparams, void* streamPtr, gmx_wallc allocateDeviceBuffer(&d_vTot_, F_NRE, nullptr); clearDeviceBufferAsync(&d_vTot_, 0, F_NRE, stream_); - for (int fType = 0; fType < F_NRE; fType++) - { - d_iLists_[fType].nr = 0; - d_iLists_[fType].iatoms = nullptr; - d_iLists_[fType].nalloc = 0; - } - kernelParams_.d_forceParams = d_forceParams_; kernelParams_.d_xq = d_xq_; kernelParams_.d_f = d_f_; @@ -114,21 +107,21 @@ GpuBonded::Impl::~Impl() } //! Return whether function type \p fType in \p idef has perturbed interactions -static bool fTypeHasPerturbedEntries(const t_idef& idef, int fType) +static bool fTypeHasPerturbedEntries(const InteractionDefinitions& idef, int fType) { GMX_ASSERT(idef.ilsort == ilsortNO_FE || idef.ilsort == ilsortFE_SORTED, "Perturbed interations should be sorted here"); - const t_ilist& ilist = idef.il[fType]; + const InteractionList& ilist = idef.il[fType]; - return (idef.ilsort != ilsortNO_FE && idef.numNonperturbedInteractions[fType] != ilist.nr); + return (idef.ilsort != ilsortNO_FE && idef.numNonperturbedInteractions[fType] != ilist.size()); } //! Converts \p src with atom indices in state order to \p dest in nbnxn order -static void convertIlistToNbnxnOrder(const t_ilist& src, - HostInteractionList* dest, - int numAtomsPerInteraction, - ArrayRef nbnxnAtomOrder) +static void convertIlistToNbnxnOrder(const InteractionList& src, + HostInteractionList* dest, + int numAtomsPerInteraction, + ArrayRef nbnxnAtomOrder) { GMX_ASSERT(src.size() == 0 || !nbnxnAtomOrder.empty(), "We need the nbnxn atom order"); @@ -175,10 +168,10 @@ static inline int roundUpToFactor(const int input, const int factor) * \todo Use DeviceBuffer for the d_xqPtr. */ void GpuBonded::Impl::updateInteractionListsAndDeviceBuffers(ArrayRef nbnxnAtomOrder, - const t_idef& idef, - void* d_xqPtr, - DeviceBuffer d_fPtr, - DeviceBuffer d_fShiftPtr) + const InteractionDefinitions& idef, + void* d_xqPtr, + DeviceBuffer d_fPtr, + DeviceBuffer d_fShiftPtr) { // TODO wallcycle sub start haveInteractions_ = false; @@ -192,7 +185,7 @@ void GpuBonded::Impl::updateInteractionListsAndDeviceBuffers(ArrayRef * But instead of doing all interactions on the CPU, we can * still easily handle the types that have no perturbed * interactions on the GPU. */ - if (idef.il[fType].nr > 0 && !fTypeHasPerturbedEntries(idef, fType)) + if (!idef.il[fType].empty() && !fTypeHasPerturbedEntries(idef, fType)) { haveInteractions_ = true; @@ -319,11 +312,11 @@ GpuBonded::GpuBonded(const gmx_ffparams_t& ffparams, void* streamPtr, gmx_wallcy GpuBonded::~GpuBonded() = default; -void GpuBonded::updateInteractionListsAndDeviceBuffers(ArrayRef nbnxnAtomOrder, - const t_idef& idef, - void* d_xq, - DeviceBuffer d_f, - DeviceBuffer d_fShift) +void GpuBonded::updateInteractionListsAndDeviceBuffers(ArrayRef nbnxnAtomOrder, + const InteractionDefinitions& idef, + void* d_xq, + DeviceBuffer d_f, + DeviceBuffer d_fShift) { impl_->updateInteractionListsAndDeviceBuffers(nbnxnAtomOrder, idef, d_xq, d_f, d_fShift); } diff --git a/src/gromacs/listed_forces/gpubonded_impl.h b/src/gromacs/listed_forces/gpubonded_impl.h index 7016e7a005..a785b16dbb 100644 --- a/src/gromacs/listed_forces/gpubonded_impl.h +++ b/src/gromacs/listed_forces/gpubonded_impl.h @@ -136,11 +136,11 @@ public: * stage. Copies the bonded interactions assigned to the GPU * to device data structures, and updates device buffers that * may have been updated after search. */ - void updateInteractionListsAndDeviceBuffers(ArrayRef nbnxnAtomOrder, - const t_idef& idef, - void* xqDevice, - DeviceBuffer forceDevice, - DeviceBuffer fshiftDevice); + void updateInteractionListsAndDeviceBuffers(ArrayRef nbnxnAtomOrder, + const InteractionDefinitions& idef, + void* xqDevice, + DeviceBuffer forceDevice, + DeviceBuffer fshiftDevice); /*! \brief Launches bonded kernel on a GPU */ template @@ -165,7 +165,7 @@ private: //! Tells whether there are any interaction in iLists. bool haveInteractions_; //! Interaction lists on the device. - t_ilist d_iLists_[F_NRE]; + t_ilist d_iLists_[F_NRE] = {}; //! Bonded parameters for device-side use. t_iparams* d_forceParams_ = nullptr; //! Position-charge vector on the device. diff --git a/src/gromacs/listed_forces/listed_forces.cpp b/src/gromacs/listed_forces/listed_forces.cpp index 2a850054af..fbfc56927a 100644 --- a/src/gromacs/listed_forces/listed_forces.cpp +++ b/src/gromacs/listed_forces/listed_forces.cpp @@ -299,32 +299,32 @@ BondedKernelFlavor selectBondedKernelFlavor(const gmx::StepWorkload& stepWork, /*! \brief Calculate one element of the list of bonded interactions for this thread */ -real calc_one_bond(int thread, - int ftype, - const t_idef* idef, - ArrayRef iatoms, - const int numNonperturbedInteractions, - const WorkDivision& workDivision, - const rvec x[], - rvec4 f[], - rvec fshift[], - const t_forcerec* fr, - const t_pbc* pbc, - const t_graph* g, - gmx_grppairener_t* grpp, - t_nrnb* nrnb, - const real* lambda, - real* dvdl, - const t_mdatoms* md, - t_fcdata* fcd, - const gmx::StepWorkload& stepWork, - int* global_atom_index) +real calc_one_bond(int thread, + int ftype, + const InteractionDefinitions& idef, + ArrayRef iatoms, + const int numNonperturbedInteractions, + const WorkDivision& workDivision, + const rvec x[], + rvec4 f[], + rvec fshift[], + const t_forcerec* fr, + const t_pbc* pbc, + const t_graph* g, + gmx_grppairener_t* grpp, + t_nrnb* nrnb, + const real* lambda, + real* dvdl, + const t_mdatoms* md, + t_fcdata* fcd, + const gmx::StepWorkload& stepWork, + int* global_atom_index) { - GMX_ASSERT(idef->ilsort == ilsortNO_FE || idef->ilsort == ilsortFE_SORTED, + GMX_ASSERT(idef.ilsort == ilsortNO_FE || idef.ilsort == ilsortFE_SORTED, "The topology should be marked either as no FE or sorted on FE"); const bool havePerturbedInteractions = - (idef->ilsort == ilsortFE_SORTED && numNonperturbedInteractions < iatoms.ssize()); + (idef.ilsort == ilsortFE_SORTED && numNonperturbedInteractions < iatoms.ssize()); BondedKernelFlavor flavor = selectBondedKernelFlavor(stepWork, fr->use_simd_kernels, havePerturbedInteractions); int efptFTYPE; @@ -346,6 +346,8 @@ real calc_one_bond(int thread, const int nb0 = workDivision.bound(ftype, thread); const int nbn = workDivision.bound(ftype, thread + 1) - nb0; + ArrayRef iparams = idef.iparams; + real v = 0; if (!isPairInteraction(ftype)) { @@ -355,12 +357,12 @@ real calc_one_bond(int thread, nice to account to its own subtimer, but first wallcycle needs to be extended to support calling from multiple threads. */ - v = cmap_dihs(nbn, iatoms.data() + nb0, idef->iparams, idef->cmap_grid, x, f, fshift, + v = cmap_dihs(nbn, iatoms.data() + nb0, iparams.data(), &idef.cmap_grid, x, f, fshift, pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd, global_atom_index); } else { - v = calculateSimpleBond(ftype, nbn, iatoms.data() + nb0, idef->iparams, x, f, fshift, + v = calculateSimpleBond(ftype, nbn, iatoms.data() + nb0, iparams.data(), x, f, fshift, pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd, global_atom_index, flavor); } @@ -370,8 +372,8 @@ real calc_one_bond(int thread, /* TODO The execution time for pairs might be nice to account to its own subtimer, but first wallcycle needs to be extended to support calling from multiple threads. */ - do_pairs(ftype, nbn, iatoms.data() + nb0, idef->iparams, x, f, fshift, pbc, g, lambda, dvdl, - md, fr, havePerturbedInteractions, stepWork, grpp, global_atom_index); + do_pairs(ftype, nbn, iatoms.data() + nb0, iparams.data(), x, f, fshift, pbc, g, lambda, + dvdl, md, fr, havePerturbedInteractions, stepWork, grpp, global_atom_index); } if (thread == 0) @@ -386,20 +388,20 @@ real calc_one_bond(int thread, /*! \brief Compute the bonded part of the listed forces, parallelized over threads */ -static void calcBondedForces(const t_idef* idef, - const rvec x[], - const t_forcerec* fr, - const t_pbc* pbc_null, - const t_graph* g, - rvec* fshiftMasterBuffer, - gmx_enerdata_t* enerd, - t_nrnb* nrnb, - const real* lambda, - real* dvdl, - const t_mdatoms* md, - t_fcdata* fcd, - const gmx::StepWorkload& stepWork, - int* global_atom_index) +static void calcBondedForces(const InteractionDefinitions& idef, + const rvec x[], + const t_forcerec* fr, + const t_pbc* pbc_null, + const t_graph* g, + rvec* fshiftMasterBuffer, + gmx_enerdata_t* enerd, + t_nrnb* nrnb, + const real* lambda, + real* dvdl, + const t_mdatoms* md, + t_fcdata* fcd, + const gmx::StepWorkload& stepWork, + int* global_atom_index) { bonded_threading_t* bt = fr->bondedThreading; @@ -440,12 +442,12 @@ static void calcBondedForces(const t_idef* idef, /* Loop over all bonded force types to calculate the bonded forces */ for (ftype = 0; (ftype < F_NRE); ftype++) { - const t_ilist& ilist = idef->il[ftype]; - if (ilist.nr > 0 && ftype_is_bonded_potential(ftype)) + const InteractionList& ilist = idef.il[ftype]; + if (!ilist.empty() && ftype_is_bonded_potential(ftype)) { - ArrayRef iatoms = gmx::constArrayRefFromArray(ilist.iatoms, ilist.nr); + ArrayRef iatoms = gmx::makeConstArrayRef(ilist.iatoms); v = calc_one_bond( - thread, ftype, idef, iatoms, idef->numNonperturbedInteractions[ftype], + thread, ftype, idef, iatoms, idef.numNonperturbedInteractions[ftype], fr->bondedThreading->workDivision, x, ft, fshift, fr, pbc_null, g, grpp, nrnb, lambda, dvdlt, md, fcd, stepWork, global_atom_index); epot[ftype] += v; @@ -456,9 +458,9 @@ static void calcBondedForces(const t_idef* idef, } } -bool haveRestraints(const t_idef& idef, const t_fcdata& fcd) +bool haveRestraints(const InteractionDefinitions& idef, const t_fcdata& fcd) { - return ((idef.il[F_POSRES].nr > 0) || (idef.il[F_FBPOSRES].nr > 0) || fcd.orires.nr > 0 + return (!idef.il[F_POSRES].empty() || !idef.il[F_FBPOSRES].empty() || fcd.orires.nr > 0 || fcd.disres.nres > 0); } @@ -467,29 +469,29 @@ bool haveCpuBondeds(const t_forcerec& fr) return fr.bondedThreading->haveBondeds; } -bool haveCpuListedForces(const t_forcerec& fr, const t_idef& idef, const t_fcdata& fcd) +bool haveCpuListedForces(const t_forcerec& fr, const InteractionDefinitions& idef, const t_fcdata& fcd) { return haveCpuBondeds(fr) || haveRestraints(idef, fcd); } -void calc_listed(const t_commrec* cr, - const gmx_multisim_t* ms, - struct gmx_wallcycle* wcycle, - const t_idef* idef, - const rvec x[], - history_t* hist, - gmx::ForceOutputs* forceOutputs, - const t_forcerec* fr, - const struct t_pbc* pbc, - const struct t_pbc* pbc_full, - const struct t_graph* g, - gmx_enerdata_t* enerd, - t_nrnb* nrnb, - const real* lambda, - const t_mdatoms* md, - t_fcdata* fcd, - int* global_atom_index, - const gmx::StepWorkload& stepWork) +void calc_listed(const t_commrec* cr, + const gmx_multisim_t* ms, + struct gmx_wallcycle* wcycle, + const InteractionDefinitions& idef, + const rvec x[], + history_t* hist, + gmx::ForceOutputs* forceOutputs, + const t_forcerec* fr, + const struct t_pbc* pbc, + const struct t_pbc* pbc_full, + const struct t_graph* g, + gmx_enerdata_t* enerd, + t_nrnb* nrnb, + const real* lambda, + const t_mdatoms* md, + t_fcdata* fcd, + int* global_atom_index, + const gmx::StepWorkload& stepWork) { const t_pbc* pbc_null; bonded_threading_t* bt = fr->bondedThreading; @@ -503,7 +505,7 @@ void calc_listed(const t_commrec* cr, pbc_null = nullptr; } - if (haveRestraints(*idef, *fcd)) + if (haveRestraints(idef, *fcd)) { /* TODO Use of restraints triggers further function calls inside the loop over calc_one_bond(), but those are too @@ -512,12 +514,12 @@ void calc_listed(const t_commrec* cr, restraints, anyway. */ wallcycle_sub_start(wcycle, ewcsRESTRAINTS); - if (idef->il[F_POSRES].nr > 0) + if (!idef.il[F_POSRES].empty()) { posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr, &forceOutputs->forceWithVirial()); } - if (idef->il[F_FBPOSRES].nr > 0) + if (!idef.il[F_FBPOSRES].empty()) { fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr, &forceOutputs->forceWithVirial()); } @@ -533,13 +535,14 @@ void calc_listed(const t_commrec* cr, */ GMX_RELEASE_ASSERT(fr->pbcType == PbcType::No || g != nullptr, "With orientation restraints molecules should be whole"); - enerd->term[F_ORIRESDEV] = calc_orires_dev(ms, idef->il[F_ORIRES].nr, idef->il[F_ORIRES].iatoms, - idef->iparams, md, x, pbc_null, fcd, hist); + enerd->term[F_ORIRESDEV] = + calc_orires_dev(ms, idef.il[F_ORIRES].size(), idef.il[F_ORIRES].iatoms.data(), + idef.iparams.data(), md, x, pbc_null, fcd, hist); } if (fcd->disres.nres > 0) { - calc_disres_R_6(cr, ms, idef->il[F_DISRES].nr, idef->il[F_DISRES].iatoms, x, pbc_null, - fcd, hist); + calc_disres_R_6(cr, ms, idef.il[F_DISRES].size(), idef.il[F_DISRES].iatoms.data(), x, + pbc_null, fcd, hist); } wallcycle_sub_stop(wcycle, ewcsRESTRAINTS); @@ -579,25 +582,24 @@ void calc_listed(const t_commrec* cr, } } -void calc_listed_lambda(const t_idef* idef, - const rvec x[], - const t_forcerec* fr, - const struct t_pbc* pbc, - const struct t_graph* g, - gmx_grppairener_t* grpp, - real* epot, - gmx::ArrayRef dvdl, - t_nrnb* nrnb, - const real* lambda, - const t_mdatoms* md, - t_fcdata* fcd, - int* global_atom_index) +void calc_listed_lambda(const InteractionDefinitions& idef, + const rvec x[], + const t_forcerec* fr, + const struct t_pbc* pbc, + const struct t_graph* g, + gmx_grppairener_t* grpp, + real* epot, + gmx::ArrayRef dvdl, + t_nrnb* nrnb, + const real* lambda, + const t_mdatoms* md, + t_fcdata* fcd, + int* global_atom_index) { real v; rvec4* f; rvec* fshift; const t_pbc* pbc_null; - t_idef idef_fe; WorkDivision& workDivision = fr->bondedThreading->foreignLambdaWorkDivision; if (fr->bMolPBC) @@ -609,9 +611,6 @@ void calc_listed_lambda(const t_idef* idef, pbc_null = nullptr; } - /* Copy the whole idef, so we can modify the contents locally */ - idef_fe = *idef; - /* We already have the forces, so we use temp buffers here */ // TODO: Get rid of these allocations by using permanent force buffers snew(f, fr->natoms_force); @@ -622,23 +621,22 @@ void calc_listed_lambda(const t_idef* idef, { if (ftype_is_bonded_potential(ftype)) { - const t_ilist& ilist = idef->il[ftype]; + const InteractionList& ilist = idef.il[ftype]; /* Create a temporary iatom list with only perturbed interactions */ - const int numNonperturbed = idef->numNonperturbedInteractions[ftype]; - ArrayRef iatoms = gmx::constArrayRefFromArray(ilist.iatoms + numNonperturbed, - ilist.nr - numNonperturbed); - t_ilist& ilist_fe = idef_fe.il[ftype]; - /* Set the work range of thread 0 to the perturbed bondeds */ - workDivision.setBound(ftype, 0, 0); - workDivision.setBound(ftype, 1, iatoms.ssize()); - - if (ilist_fe.nr > 0) + const int numNonperturbed = idef.numNonperturbedInteractions[ftype]; + ArrayRef iatomsPerturbed = gmx::constArrayRefFromArray( + ilist.iatoms.data() + numNonperturbed, ilist.size() - numNonperturbed); + if (!iatomsPerturbed.empty()) { + /* Set the work range of thread 0 to the perturbed bondeds */ + workDivision.setBound(ftype, 0, 0); + workDivision.setBound(ftype, 1, iatomsPerturbed.ssize()); + gmx::StepWorkload tempFlags; tempFlags.computeEnergy = true; - v = calc_one_bond(0, ftype, idef, iatoms, iatoms.ssize(), workDivision, x, f, - fshift, fr, pbc_null, g, grpp, nrnb, lambda, dvdl.data(), md, fcd, - tempFlags, global_atom_index); + v = calc_one_bond(0, ftype, idef, iatomsPerturbed, iatomsPerturbed.ssize(), + workDivision, x, f, fshift, fr, pbc_null, g, grpp, nrnb, lambda, + dvdl.data(), md, fcd, tempFlags, global_atom_index); epot[ftype] += v; } } @@ -648,25 +646,25 @@ void calc_listed_lambda(const t_idef* idef, sfree(f); } -void do_force_listed(struct gmx_wallcycle* wcycle, - const matrix box, - const t_lambda* fepvals, - const t_commrec* cr, - const gmx_multisim_t* ms, - const t_idef* idef, - const rvec x[], - history_t* hist, - gmx::ForceOutputs* forceOutputs, - const t_forcerec* fr, - const struct t_pbc* pbc, - const struct t_graph* graph, - gmx_enerdata_t* enerd, - t_nrnb* nrnb, - const real* lambda, - const t_mdatoms* md, - t_fcdata* fcd, - int* global_atom_index, - const gmx::StepWorkload& stepWork) +void do_force_listed(struct gmx_wallcycle* wcycle, + const matrix box, + const t_lambda* fepvals, + const t_commrec* cr, + const gmx_multisim_t* ms, + const InteractionDefinitions& idef, + const rvec x[], + history_t* hist, + gmx::ForceOutputs* forceOutputs, + const t_forcerec* fr, + const struct t_pbc* pbc, + const struct t_graph* graph, + gmx_enerdata_t* enerd, + t_nrnb* nrnb, + const real* lambda, + const t_mdatoms* md, + t_fcdata* fcd, + int* global_atom_index, + const gmx::StepWorkload& stepWork) { t_pbc pbc_full; /* Full PBC is needed for position restraints */ @@ -675,7 +673,7 @@ void do_force_listed(struct gmx_wallcycle* wcycle, return; } - if ((idef->il[F_POSRES].nr > 0) || (idef->il[F_FBPOSRES].nr > 0)) + if (!idef.il[F_POSRES].empty() || !idef.il[F_FBPOSRES].empty()) { /* Not enough flops to bother counting */ set_pbc(&pbc_full, fr->pbcType, box); @@ -691,10 +689,10 @@ void do_force_listed(struct gmx_wallcycle* wcycle, real dvdl[efptNR] = { 0 }; posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr); - if (idef->ilsort != ilsortNO_FE) + if (idef.ilsort != ilsortNO_FE) { wallcycle_sub_start(wcycle, ewcsLISTED_FEP); - if (idef->ilsort != ilsortFE_SORTED) + if (idef.ilsort != ilsortFE_SORTED) { gmx_incons("The bonded interactions are not sorted for free energy"); } diff --git a/src/gromacs/listed_forces/listed_forces.h b/src/gromacs/listed_forces/listed_forces.h index 5ad7d35dc0..479a4d8fff 100644 --- a/src/gromacs/listed_forces/listed_forces.h +++ b/src/gromacs/listed_forces/listed_forces.h @@ -75,10 +75,10 @@ struct gmx_enerdata_t; struct gmx_grppairener_t; struct gmx_multisim_t; class history_t; +class InteractionDefinitions; struct t_commrec; struct t_fcdata; struct t_forcerec; -struct t_idef; struct t_graph; struct t_lambda; struct t_mdatoms; @@ -115,67 +115,67 @@ BondedFunction bondedFunction(int ftype); * * Note that pbc_full is used only for position restraints, and is * not initialized if there are none. */ -void calc_listed(const t_commrec* cr, - const gmx_multisim_t* ms, - struct gmx_wallcycle* wcycle, - const t_idef* idef, - const rvec x[], - history_t* hist, - gmx::ForceOutputs* forceOutputs, - const t_forcerec* fr, - const struct t_pbc* pbc, - const struct t_pbc* pbc_full, - const struct t_graph* g, - gmx_enerdata_t* enerd, - t_nrnb* nrnb, - const real* lambda, - const t_mdatoms* md, - struct t_fcdata* fcd, - int* ddgatindex, - const gmx::StepWorkload& stepWork); +void calc_listed(const t_commrec* cr, + const gmx_multisim_t* ms, + struct gmx_wallcycle* wcycle, + const InteractionDefinitions& idef, + const rvec x[], + history_t* hist, + gmx::ForceOutputs* forceOutputs, + const t_forcerec* fr, + const struct t_pbc* pbc, + const struct t_pbc* pbc_full, + const struct t_graph* g, + gmx_enerdata_t* enerd, + t_nrnb* nrnb, + const real* lambda, + const t_mdatoms* md, + struct t_fcdata* fcd, + int* ddgatindex, + const gmx::StepWorkload& stepWork); /*! \brief As calc_listed(), but only determines the potential energy * for the perturbed interactions. * * The shift forces in fr are not affected. */ -void calc_listed_lambda(const t_idef* idef, - const rvec x[], - const t_forcerec* fr, - const struct t_pbc* pbc, - const struct t_graph* g, - gmx_grppairener_t* grpp, - real* epot, - gmx::ArrayRef dvdl, - t_nrnb* nrnb, - const real* lambda, - const t_mdatoms* md, - struct t_fcdata* fcd, - int* global_atom_index); +void calc_listed_lambda(const InteractionDefinitions& idef, + const rvec x[], + const t_forcerec* fr, + const struct t_pbc* pbc, + const struct t_graph* g, + gmx_grppairener_t* grpp, + real* epot, + gmx::ArrayRef dvdl, + t_nrnb* nrnb, + const real* lambda, + const t_mdatoms* md, + struct t_fcdata* fcd, + int* global_atom_index); /*! \brief Do all aspects of energy and force calculations for mdrun * on the set of listed interactions */ -void do_force_listed(struct gmx_wallcycle* wcycle, - const matrix box, - const t_lambda* fepvals, - const t_commrec* cr, - const gmx_multisim_t* ms, - const t_idef* idef, - const rvec x[], - history_t* hist, - gmx::ForceOutputs* forceOutputs, - const t_forcerec* fr, - const struct t_pbc* pbc, - const struct t_graph* graph, - gmx_enerdata_t* enerd, - t_nrnb* nrnb, - const real* lambda, - const t_mdatoms* md, - struct t_fcdata* fcd, - int* global_atom_index, - const gmx::StepWorkload& stepWork); +void do_force_listed(struct gmx_wallcycle* wcycle, + const matrix box, + const t_lambda* fepvals, + const t_commrec* cr, + const gmx_multisim_t* ms, + const InteractionDefinitions& idef, + const rvec x[], + history_t* hist, + gmx::ForceOutputs* forceOutputs, + const t_forcerec* fr, + const struct t_pbc* pbc, + const struct t_graph* graph, + gmx_enerdata_t* enerd, + t_nrnb* nrnb, + const real* lambda, + const t_mdatoms* md, + struct t_fcdata* fcd, + int* global_atom_index, + const gmx::StepWorkload& stepWork); /*! \brief Returns true if there are position, distance or orientation restraints. */ -bool haveRestraints(const t_idef& idef, const t_fcdata& fcd); +bool haveRestraints(const InteractionDefinitions& idef, const t_fcdata& fcd); /*! \brief Returns true if there are CPU (i.e. not GPU-offloaded) bonded interactions to compute. */ bool haveCpuBondeds(const t_forcerec& fr); @@ -185,6 +185,6 @@ bool haveCpuBondeds(const t_forcerec& fr); * NOTE: the current implementation returns true if there are position restraints * or any bonded interactions computed on the CPU. */ -bool haveCpuListedForces(const t_forcerec& fr, const t_idef& idef, const t_fcdata& fcd); +bool haveCpuListedForces(const t_forcerec& fr, const InteractionDefinitions& idef, const t_fcdata& fcd); #endif diff --git a/src/gromacs/listed_forces/manage_threading.cpp b/src/gromacs/listed_forces/manage_threading.cpp index 2981ddef56..5dc168e5cf 100644 --- a/src/gromacs/listed_forces/manage_threading.cpp +++ b/src/gromacs/listed_forces/manage_threading.cpp @@ -71,9 +71,9 @@ /*! \brief struct for passing all data required for a function type */ typedef struct { - const t_ilist* il; /**< pointer to t_ilist entry corresponding to ftype */ - int ftype; /**< the function type index */ - int nat; /**< nr of atoms involved in a single ftype interaction */ + const InteractionList* il; /**< pointer to t_ilist entry corresponding to ftype */ + int ftype; /**< the function type index */ + int nat; /**< nr of atoms involved in a single ftype interaction */ } ilist_data_t; /*! \brief Divides listed interactions over threads @@ -96,11 +96,11 @@ static void divide_bondeds_by_locality(bonded_threading_t* bt, int numType, cons for (f = 0; f < numType; f++) { /* Sum #bondeds*#atoms_per_bond over all bonded types */ - nat_tot += ild[f].il->nr / (ild[f].nat + 1) * ild[f].nat; + nat_tot += ild[f].il->size() / (ild[f].nat + 1) * ild[f].nat; /* The start bound for thread 0 is 0 for all interactions */ ind[f] = 0; /* Initialize the next atom index array */ - assert(ild[f].il->nr > 0); + assert(!ild[f].il->empty()); at_ind[f] = ild[f].il->iatoms[1]; } @@ -162,7 +162,7 @@ static void divide_bondeds_by_locality(bonded_threading_t* bt, int numType, cons nat_sum += ild[f_min].nat; /* Update the first unassigned atom index for this type */ - if (ind[f_min] < ild[f_min].il->nr) + if (ind[f_min] < ild[f_min].il->size()) { at_ind[f_min] = ild[f_min].il->iatoms[ind[f_min] + 1]; } @@ -185,29 +185,33 @@ static void divide_bondeds_by_locality(bonded_threading_t* bt, int numType, cons for (f = 0; f < numType; f++) { - assert(ind[f] == ild[f].il->nr); + assert(ind[f] == ild[f].il->size()); } } //! Return whether function type \p ftype in \p idef has perturbed interactions -static bool ftypeHasPerturbedEntries(const t_idef& idef, int ftype) +static bool ftypeHasPerturbedEntries(const InteractionDefinitions& idef, int ftype) { GMX_ASSERT(idef.ilsort == ilsortNO_FE || idef.ilsort == ilsortFE_SORTED, "Perturbed interations should be sorted here"); - const t_ilist& ilist = idef.il[ftype]; + const InteractionList& ilist = idef.il[ftype]; - return (idef.ilsort != ilsortNO_FE && idef.numNonperturbedInteractions[ftype] != ilist.nr); + return (idef.ilsort != ilsortNO_FE && idef.numNonperturbedInteractions[ftype] != ilist.size()); } //! Divides bonded interactions over threads and GPU -static void divide_bondeds_over_threads(bonded_threading_t* bt, bool useGpuForBondeds, const t_idef& idef) +static void divide_bondeds_over_threads(bonded_threading_t* bt, + bool useGpuForBondeds, + const InteractionDefinitions& idef) { ilist_data_t ild[F_NRE]; GMX_ASSERT(bt->nthreads > 0, "Must have positive number of threads"); const int numThreads = bt->nthreads; + gmx::ArrayRef iparams = idef.iparams; + bt->haveBondeds = false; int numType = 0; size_t fTypeGpuIndex = 0; @@ -218,8 +222,8 @@ static void divide_bondeds_over_threads(bonded_threading_t* bt, bool useGpuForBo continue; } - const t_ilist& il = idef.il[fType]; - int nrToAssignToCpuThreads = il.nr; + const InteractionList& il = idef.il[fType]; + int nrToAssignToCpuThreads = il.size(); if (useGpuForBondeds && fTypeGpuIndex < gmx::fTypesOnGpu.size() && gmx::fTypesOnGpu[fTypeGpuIndex] == fType) @@ -269,8 +273,8 @@ static void divide_bondeds_over_threads(bonded_threading_t* bt, bool useGpuForBo * end up on the same thread. */ while (nr_t > 0 && nr_t < nrToAssignToCpuThreads - && idef.iparams[il.iatoms[nr_t]].disres.label - == idef.iparams[il.iatoms[nr_t - stride]].disres.label) + && iparams[il.iatoms[nr_t]].disres.label + == iparams[il.iatoms[nr_t - stride]].disres.label) { nr_t += stride; } @@ -306,7 +310,7 @@ static void divide_bondeds_over_threads(bonded_threading_t* bt, bool useGpuForBo fprintf(debug, "Division of bondeds over threads:\n"); for (f = 0; f < F_NRE; f++) { - if (ftype_is_bonded_potential(f) && idef.il[f].nr > 0) + if (ftype_is_bonded_potential(f) && !idef.il[f].empty()) { int t; @@ -324,11 +328,11 @@ static void divide_bondeds_over_threads(bonded_threading_t* bt, bool useGpuForBo } //! Construct a reduction mask for which parts (blocks) of the force array are touched on which thread task -static void calc_bonded_reduction_mask(int natoms, - f_thread_t* f_thread, - const t_idef& idef, - int thread, - const bonded_threading_t& bondedThreading) +static void calc_bonded_reduction_mask(int natoms, + f_thread_t* f_thread, + const InteractionDefinitions& idef, + int thread, + const bonded_threading_t& bondedThreading) { static_assert(BITMASK_SIZE == GMX_OPENMP_MAX_THREADS, "For the error message below we assume these two are equal."); @@ -369,7 +373,7 @@ static void calc_bonded_reduction_mask(int natoms, { if (ftype_is_bonded_potential(ftype)) { - int nb = idef.il[ftype].nr; + int nb = idef.il[ftype].size(); if (nb > 0) { int nat1 = interaction_function[ftype].nratoms + 1; @@ -401,7 +405,10 @@ static void calc_bonded_reduction_mask(int natoms, } } -void setup_bonded_threading(bonded_threading_t* bt, int numAtoms, bool useGpuForBondeds, const t_idef& idef) +void setup_bonded_threading(bonded_threading_t* bt, + int numAtoms, + bool useGpuForBondeds, + const InteractionDefinitions& idef) { int ctot = 0; diff --git a/src/gromacs/listed_forces/manage_threading.h b/src/gromacs/listed_forces/manage_threading.h index 48d8c1cf68..c2b27d493b 100644 --- a/src/gromacs/listed_forces/manage_threading.h +++ b/src/gromacs/listed_forces/manage_threading.h @@ -48,7 +48,7 @@ #include struct bonded_threading_t; -struct t_idef; +class InteractionDefinitions; /*! \brief Divide the listed interactions over the threads and GPU * @@ -57,7 +57,10 @@ struct t_idef; * This should be called each time the bonded setup changes; * i.e. at start-up without domain decomposition and at DD. */ -void setup_bonded_threading(bonded_threading_t* bt, int numAtoms, bool useGpuForBondes, const t_idef& idef); +void setup_bonded_threading(bonded_threading_t* bt, + int numAtoms, + bool useGpuForBondes, + const InteractionDefinitions& idef); //! Destructor. void tear_down_bonded_threading(bonded_threading_t* bt); diff --git a/src/gromacs/listed_forces/orires.cpp b/src/gromacs/listed_forces/orires.cpp index d9c58c9574..1f9665b62f 100644 --- a/src/gromacs/listed_forces/orires.cpp +++ b/src/gromacs/listed_forces/orires.cpp @@ -123,7 +123,8 @@ void init_orires(FILE* fplog, int nmol; while (const InteractionLists* il = gmx_mtop_ilistloop_next(iloop, &nmol)) { - if (nmol > 1 && (*il)[F_ORIRES].size() > 0) + const int numOrires = (*il)[F_ORIRES].size(); + if (nmol > 1 && numOrires > 0) { gmx_fatal(FARGS, "Found %d copies of a molecule with orientation restrains while the current " @@ -132,7 +133,7 @@ void init_orires(FILE* fplog, nmol); } - for (int i = 0; i < (*il)[F_ORIRES].size(); i += 3) + for (int i = 0; i < numOrires; i += 3) { int type = (*il)[F_ORIRES].iatoms[i]; int ex = mtop->ffparams.iparams[type].orires.ex; diff --git a/src/gromacs/listed_forces/position_restraints.cpp b/src/gromacs/listed_forces/position_restraints.cpp index 30449fe0e2..d9b2b22b5e 100644 --- a/src/gromacs/listed_forces/position_restraints.cpp +++ b/src/gromacs/listed_forces/position_restraints.cpp @@ -401,41 +401,42 @@ real posres(int nbonds, } // namespace -void posres_wrapper(t_nrnb* nrnb, - const t_idef* idef, - const struct t_pbc* pbc, - const rvec* x, - gmx_enerdata_t* enerd, - const real* lambda, - const t_forcerec* fr, - gmx::ForceWithVirial* forceWithVirial) +void posres_wrapper(t_nrnb* nrnb, + const InteractionDefinitions& idef, + const struct t_pbc* pbc, + const rvec* x, + gmx_enerdata_t* enerd, + const real* lambda, + const t_forcerec* fr, + gmx::ForceWithVirial* forceWithVirial) { real v, dvdl; dvdl = 0; - v = posres(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms, idef->iparams_posres, x, - forceWithVirial, fr->pbcType == PbcType::No ? nullptr : pbc, lambda[efptRESTRAINT], - &dvdl, fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB); + v = posres(idef.il[F_POSRES].size(), idef.il[F_POSRES].iatoms.data(), + idef.iparams_posres.data(), x, forceWithVirial, + fr->pbcType == PbcType::No ? nullptr : pbc, lambda[efptRESTRAINT], &dvdl, + fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB); enerd->term[F_POSRES] += v; /* If just the force constant changes, the FEP term is linear, * but if k changes, it is not. */ enerd->dvdl_nonlin[efptRESTRAINT] += dvdl; - inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef->il[F_POSRES].nr, 2)); + inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef.il[F_POSRES].size(), 2)); } -void posres_wrapper_lambda(struct gmx_wallcycle* wcycle, - const t_lambda* fepvals, - const t_idef* idef, - const struct t_pbc* pbc, - const rvec x[], - gmx_enerdata_t* enerd, - const real* lambda, - const t_forcerec* fr) +void posres_wrapper_lambda(struct gmx_wallcycle* wcycle, + const t_lambda* fepvals, + const InteractionDefinitions& idef, + const struct t_pbc* pbc, + const rvec x[], + gmx_enerdata_t* enerd, + const real* lambda, + const t_forcerec* fr) { real v; - if (0 == idef->il[F_POSRES].nr) + if (idef.il[F_POSRES].empty()) { return; } @@ -447,8 +448,9 @@ void posres_wrapper_lambda(struct gmx_wallcycle* wcycle, const real lambda_dum = (i == 0 ? lambda[efptRESTRAINT] : fepvals->all_lambda[efptRESTRAINT][i - 1]); - v = posres(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms, idef->iparams_posres, x, - nullptr, fr->pbcType == PbcType::No ? nullptr : pbc, lambda_dum, &dvdl, + v = posres(idef.il[F_POSRES].size(), idef.il[F_POSRES].iatoms.data(), + idef.iparams_posres.data(), x, nullptr, + fr->pbcType == PbcType::No ? nullptr : pbc, lambda_dum, &dvdl, fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB); enerd->enerpart_lambda[i] += v; enerd->dhdlLambda[i] += dvdl; @@ -458,19 +460,19 @@ void posres_wrapper_lambda(struct gmx_wallcycle* wcycle, /*! \brief Helper function that wraps calls to fbposres for free-energy perturbation */ -void fbposres_wrapper(t_nrnb* nrnb, - const t_idef* idef, - const struct t_pbc* pbc, - const rvec* x, - gmx_enerdata_t* enerd, - const t_forcerec* fr, - gmx::ForceWithVirial* forceWithVirial) +void fbposres_wrapper(t_nrnb* nrnb, + const InteractionDefinitions& idef, + const struct t_pbc* pbc, + const rvec* x, + gmx_enerdata_t* enerd, + const t_forcerec* fr, + gmx::ForceWithVirial* forceWithVirial) { real v; - v = fbposres(idef->il[F_FBPOSRES].nr, idef->il[F_FBPOSRES].iatoms, idef->iparams_fbposres, x, - forceWithVirial, fr->pbcType == PbcType::No ? nullptr : pbc, fr->rc_scaling, - fr->pbcType, fr->posres_com); + v = fbposres(idef.il[F_FBPOSRES].size(), idef.il[F_FBPOSRES].iatoms.data(), + idef.iparams_fbposres.data(), x, forceWithVirial, + fr->pbcType == PbcType::No ? nullptr : pbc, fr->rc_scaling, fr->pbcType, fr->posres_com); enerd->term[F_FBPOSRES] += v; - inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef->il[F_FBPOSRES].nr, 2)); + inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef.il[F_FBPOSRES].size(), 2)); } diff --git a/src/gromacs/listed_forces/position_restraints.h b/src/gromacs/listed_forces/position_restraints.h index fe9eb6cb19..5d3148fb33 100644 --- a/src/gromacs/listed_forces/position_restraints.h +++ b/src/gromacs/listed_forces/position_restraints.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2017,2019, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2017,2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -56,7 +56,7 @@ struct gmx_enerdata_t; struct gmx_wallcycle; struct t_forcerec; -struct t_idef; +class InteractionDefinitions; struct t_lambda; struct t_nrnb; struct t_pbc; @@ -67,34 +67,34 @@ class ForceWithVirial; } /*! \brief Helper function that wraps calls to posres */ -void posres_wrapper(t_nrnb* nrnb, - const t_idef* idef, - const struct t_pbc* pbc, - const rvec* x, - gmx_enerdata_t* enerd, - const real* lambda, - const t_forcerec* fr, - gmx::ForceWithVirial* forceWithVirial); +void posres_wrapper(t_nrnb* nrnb, + const InteractionDefinitions& idef, + const struct t_pbc* pbc, + const rvec* x, + gmx_enerdata_t* enerd, + const real* lambda, + const t_forcerec* fr, + gmx::ForceWithVirial* forceWithVirial); /*! \brief Helper function that wraps calls to posres for free-energy pertubation */ -void posres_wrapper_lambda(struct gmx_wallcycle* wcycle, - const t_lambda* fepvals, - const t_idef* idef, - const struct t_pbc* pbc, - const rvec x[], - gmx_enerdata_t* enerd, - const real* lambda, - const t_forcerec* fr); +void posres_wrapper_lambda(struct gmx_wallcycle* wcycle, + const t_lambda* fepvals, + const InteractionDefinitions& idef, + const struct t_pbc* pbc, + const rvec x[], + gmx_enerdata_t* enerd, + const real* lambda, + const t_forcerec* fr); /*! \brief Helper function that wraps calls to fbposres for free-energy perturbation */ -void fbposres_wrapper(t_nrnb* nrnb, - const t_idef* idef, - const struct t_pbc* pbc, - const rvec* x, - gmx_enerdata_t* enerd, - const t_forcerec* fr, - gmx::ForceWithVirial* forceWithVirial); +void fbposres_wrapper(t_nrnb* nrnb, + const InteractionDefinitions& idef, + const struct t_pbc* pbc, + const rvec* x, + gmx_enerdata_t* enerd, + const t_forcerec* fr, + gmx::ForceWithVirial* forceWithVirial); #endif diff --git a/src/gromacs/mdlib/calc_verletbuf.cpp b/src/gromacs/mdlib/calc_verletbuf.cpp index 2a24be9977..90e9acf927 100644 --- a/src/gromacs/mdlib/calc_verletbuf.cpp +++ b/src/gromacs/mdlib/calc_verletbuf.cpp @@ -2,7 +2,7 @@ * This file is part of the GROMACS molecular simulation package. * * Copyright (c) 2012-2018, The GROMACS development team. - * Copyright (c) 2019, by the GROMACS development team, led by + * Copyright (c) 2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -1256,7 +1256,7 @@ static real chanceOfUpdateGroupCrossingCell(const gmx_moltype_t& moltyp for (const int atom : block) { const auto& ilist = moltype.ilist[F_SETTLE]; - GMX_RELEASE_ASSERT(ilist.size() > 0, + GMX_RELEASE_ASSERT(!ilist.empty(), "There should be at least one settle in this moltype"); for (int i = 0; i < ilist.size(); i += 1 + NRAL(F_SETTLE)) { diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index 2281946320..f044e79133 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -112,7 +112,7 @@ public: int numConstraints, int numSettles); ~Impl(); - void setConstraints(const gmx_localtop_t& top, const t_mdatoms& md); + void setConstraints(gmx_localtop_t* top, const t_mdatoms& md); bool apply(bool bLog, bool bEner, int64_t step, @@ -158,7 +158,7 @@ public: //! Pointer to the global topology - only used for printing warnings. const gmx_mtop_t& mtop; //! Parameters for the interactions in this domain. - const t_idef* idef = nullptr; + const InteractionDefinitions* idef = nullptr; //! Data about atoms in this domain. const t_mdatoms& md; //! Whether we need to do pbc for handling bonds. @@ -407,8 +407,8 @@ bool Constraints::Impl::apply(bool bLog, { clear_mat(vir_r_m_dr); } - const t_ilist* settle = &idef->il[F_SETTLE]; - nsettle = settle->nr / (1 + NRAL(F_SETTLE)); + const InteractionList& settle = idef->il[F_SETTLE]; + nsettle = settle.size() / (1 + NRAL(F_SETTLE)); if (nsettle > 0) { @@ -555,7 +555,7 @@ bool Constraints::Impl::apply(bool bLog, if (start_th >= 0 && end_th - start_th > 0) { settle_proj(settled, econq, end_th - start_th, - settle->iatoms + start_th * (1 + NRAL(F_SETTLE)), pbc_null, + settle.iatoms.data() + start_th * (1 + NRAL(F_SETTLE)), pbc_null, x.unpaddedArrayRef(), xprime.unpaddedArrayRef(), min_proj, calcvir_atom_end, th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]); } @@ -759,21 +759,20 @@ FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegra * * \returns a block struct with all constraints for each atom */ -template -static ListOfLists makeAtomsToConstraintsList(int numAtoms, - const T* ilists, - const t_iparams* iparams, +static ListOfLists makeAtomsToConstraintsList(int numAtoms, + ArrayRef ilists, + ArrayRef iparams, FlexibleConstraintTreatment flexibleConstraintTreatment) { - GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr, + GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || !iparams.empty(), "With flexible constraint detection we need valid iparams"); std::vector count(numAtoms); for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++) { - const T& ilist = ilists[ftype]; - const int stride = 1 + NRAL(ftype); + const InteractionList& ilist = ilists[ftype]; + const int stride = 1 + NRAL(ftype); for (int i = 0; i < ilist.size(); i += stride) { if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include @@ -802,8 +801,8 @@ static ListOfLists makeAtomsToConstraintsList(int n int numConstraints = 0; for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++) { - const T& ilist = ilists[ftype]; - const int stride = 1 + NRAL(ftype); + const InteractionList& ilist = ilists[ftype]; + const int stride = 1 + NRAL(ftype); for (int i = 0; i < ilist.size(); i += stride) { if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include @@ -822,10 +821,10 @@ static ListOfLists makeAtomsToConstraintsList(int n return ListOfLists(std::move(listRanges), std::move(elements)); } -ListOfLists make_at2con(int numAtoms, - const t_ilist* ilist, - const t_iparams* iparams, - FlexibleConstraintTreatment flexibleConstraintTreatment) +ListOfLists make_at2con(int numAtoms, + ArrayRef ilist, + ArrayRef iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment) { return makeAtomsToConstraintsList(numAtoms, ilist, iparams, flexibleConstraintTreatment); } @@ -834,13 +833,12 @@ ListOfLists make_at2con(const gmx_moltype_t& moltype, gmx::ArrayRef iparams, FlexibleConstraintTreatment flexibleConstraintTreatment) { - return makeAtomsToConstraintsList(moltype.atoms.nr, moltype.ilist.data(), iparams.data(), + return makeAtomsToConstraintsList(moltype.atoms.nr, makeConstArrayRef(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 countFlexibleConstraints(ArrayRef ilist, ArrayRef iparams) { int nflexcon = 0; for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++) @@ -859,11 +857,6 @@ static int countFlexibleConstraintsTemplate(const T* ilist, const t_iparams* ipa return nflexcon; } -int countFlexibleConstraints(const t_ilist* ilist, const t_iparams* iparams) -{ - return countFlexibleConstraintsTemplate(ilist, iparams); -} - //! Returns the index of the settle to which each atom belongs. static std::vector make_at2settle(int natoms, const InteractionList& ilist) { @@ -882,9 +875,9 @@ static std::vector make_at2settle(int natoms, const InteractionList& ilist) return at2s; } -void Constraints::Impl::setConstraints(const gmx_localtop_t& top, const t_mdatoms& md) +void Constraints::Impl::setConstraints(gmx_localtop_t* top, const t_mdatoms& md) { - idef = &top.idef; + idef = &top->idef; if (ncon_tot > 0) { @@ -893,7 +886,7 @@ void Constraints::Impl::setConstraints(const gmx_localtop_t& top, const t_mdatom */ if (ir.eConstrAlg == econtLINCS) { - set_lincs(top.idef, md, EI_DYNAMICS(ir.eI), cr, lincsd); + set_lincs(*idef, md, EI_DYNAMICS(ir.eI), cr, lincsd); } if (ir.eConstrAlg == econtSHAKE) { @@ -901,18 +894,18 @@ void Constraints::Impl::setConstraints(const gmx_localtop_t& top, const t_mdatom { // We are using the local topology, so there are only // F_CONSTR constraints. - make_shake_sblock_dd(shaked, &idef->il[F_CONSTR], cr->dd); + make_shake_sblock_dd(shaked, idef->il[F_CONSTR], cr->dd); } else { - make_shake_sblock_serial(shaked, idef, md); + make_shake_sblock_serial(shaked, &top->idef, md); } } } if (settled) { - settle_set_constraints(settled, &idef->il[F_SETTLE], md); + settle_set_constraints(settled, idef->il[F_SETTLE], md); } /* Make a selection of the local atoms for essential dynamics */ @@ -922,7 +915,7 @@ void Constraints::Impl::setConstraints(const gmx_localtop_t& top, const t_mdatom } } -void Constraints::setConstraints(const gmx_localtop_t& top, const t_mdatoms& md) +void Constraints::setConstraints(gmx_localtop_t* top, const t_mdatoms& md) { impl_->setConstraints(top, md); } @@ -995,8 +988,7 @@ Constraints::Impl::Impl(const gmx_mtop_t& mtop_p, for (const gmx_molblock_t& molblock : mtop.molblock) { - int count = countFlexibleConstraintsTemplate(mtop.moltype[molblock.type].ilist.data(), - mtop.ffparams.iparams.data()); + int count = countFlexibleConstraints(mtop.moltype[molblock.type].ilist, mtop.ffparams.iparams); nflexcon += molblock.nmol * count; } diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h index 216ea479b6..680a941acb 100644 --- a/src/gromacs/mdlib/constr.h +++ b/src/gromacs/mdlib/constr.h @@ -134,7 +134,7 @@ public: * * \todo Make this a callback that is called automatically * once a new domain has been made. */ - void setConstraints(const gmx_localtop_t& top, const t_mdatoms& md); + void setConstraints(gmx_localtop_t* top, const t_mdatoms& md); /*! \brief Applies constraints to coordinates. * @@ -215,10 +215,8 @@ private: [[noreturn]] void too_many_constraint_warnings(int eConstrAlg, int warncount); /*! \brief Returns whether constraint with parameter \p iparamsIndex is a flexible constraint */ -static inline bool isConstraintFlexible(const t_iparams* iparams, int iparamsIndex) +static inline bool isConstraintFlexible(ArrayRef iparams, int iparamsIndex) { - GMX_ASSERT(iparams != nullptr, "Need a valid iparams array"); - return (iparams[iparamsIndex].constr.dA == 0 && iparams[iparamsIndex].constr.dB == 0); }; @@ -271,13 +269,13 @@ ListOfLists make_at2con(const gmx_moltype_t& moltype, * * \returns a ListOfLists object with all constraints for each atom */ -ListOfLists make_at2con(int numAtoms, - const t_ilist* ilist, - const t_iparams* iparams, - FlexibleConstraintTreatment flexibleConstraintTreatment); +ListOfLists make_at2con(int numAtoms, + ArrayRef ilist, + ArrayRef iparams, + FlexibleConstraintTreatment flexibleConstraintTreatment); //! Return the number of flexible constraints in the \c ilist and \c iparams. -int countFlexibleConstraints(const t_ilist* ilist, const t_iparams* iparams); +int countFlexibleConstraints(ArrayRef ilist, ArrayRef iparams); /*! \brief Returns the constraint iatoms for a constraint number con * which comes from a list where F_CONSTR and F_CONSTRNC constraints diff --git a/src/gromacs/mdlib/constraintrange.cpp b/src/gromacs/mdlib/constraintrange.cpp index 30806eb71f..90d16294d0 100644 --- a/src/gromacs/mdlib/constraintrange.cpp +++ b/src/gromacs/mdlib/constraintrange.cpp @@ -163,7 +163,7 @@ static real constr_r_max_moltype(const gmx_moltype_t* molt, real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1; - if (molt->ilist[F_CONSTR].size() == 0 && molt->ilist[F_CONSTRNC].size() == 0) + if (molt->ilist[F_CONSTR].empty() && molt->ilist[F_CONSTRNC].empty()) { return 0; } diff --git a/src/gromacs/mdlib/force.cpp b/src/gromacs/mdlib/force.cpp index 1c9c70707f..48ba897528 100644 --- a/src/gromacs/mdlib/force.cpp +++ b/src/gromacs/mdlib/force.cpp @@ -102,7 +102,7 @@ static void reduceEwaldThreadOuput(int nthreads, ewald_corr_thread_t* ewc_t) void do_force_lowlevel(t_forcerec* fr, const t_inputrec* ir, - const t_idef* idef, + const InteractionDefinitions& idef, const t_commrec* cr, const gmx_multisim_t* ms, t_nrnb* nrnb, @@ -178,7 +178,7 @@ void do_force_lowlevel(t_forcerec* fr, /* Check whether we need to take into account PBC in listed interactions. */ const auto needPbcForListedForces = - fr->bMolPBC && stepWork.computeListedForces && haveCpuListedForces(*fr, *idef, *fcd); + fr->bMolPBC && stepWork.computeListedForces && haveCpuListedForces(*fr, idef, *fcd); if (needPbcForListedForces) { /* Since all atoms are in the rectangular or triclinic unit-cell, diff --git a/src/gromacs/mdlib/force.h b/src/gromacs/mdlib/force.h index 7b568c80a1..17a90bb3bf 100644 --- a/src/gromacs/mdlib/force.h +++ b/src/gromacs/mdlib/force.h @@ -52,12 +52,12 @@ struct gmx_multisim_t; struct gmx_vsite_t; struct gmx_wallcycle; class history_t; +class InteractionDefinitions; struct pull_t; struct t_commrec; struct t_fcdata; struct t_forcerec; struct t_graph; -struct t_idef; struct t_inputrec; struct t_lambda; struct t_mdatoms; @@ -117,7 +117,7 @@ void do_force(FILE* log, void do_force_lowlevel(t_forcerec* fr, const t_inputrec* ir, - const t_idef* idef, + const InteractionDefinitions& idef, const t_commrec* cr, const gmx_multisim_t* ms, t_nrnb* nrnb, diff --git a/src/gromacs/mdlib/lincs.cpp b/src/gromacs/mdlib/lincs.cpp index b425d12eaf..8b13772cdf 100644 --- a/src/gromacs/mdlib/lincs.cpp +++ b/src/gromacs/mdlib/lincs.cpp @@ -1681,13 +1681,13 @@ static void assign_constraint(Lincs* li, /*! \brief Check if constraint with topology index constraint_index is connected * to other constraints, and if so add those connected constraints to our task. */ -static void check_assign_connected(Lincs* li, - const t_iatom* iatom, - const t_idef& idef, - bool bDynamics, - int a1, - int a2, - const ListOfLists& at2con) +static void check_assign_connected(Lincs* li, + gmx::ArrayRef iatom, + const InteractionDefinitions& idef, + bool bDynamics, + int a1, + int a2, + const ListOfLists& at2con) { /* Currently this function only supports constraint groups * in which all constraints share at least one atom @@ -1721,14 +1721,14 @@ static void check_assign_connected(Lincs* li, /*! \brief Check if constraint with topology index constraint_index is involved * in a constraint triangle, and if so add the other two constraints * in the triangle to our task. */ -static void check_assign_triangle(Lincs* li, - const t_iatom* iatom, - const t_idef& idef, - bool bDynamics, - int constraint_index, - int a1, - int a2, - const ListOfLists& at2con) +static void check_assign_triangle(Lincs* li, + gmx::ArrayRef iatom, + const InteractionDefinitions& idef, + bool bDynamics, + int constraint_index, + int a1, + int a2, + const ListOfLists& at2con) { int nca, cc[32], ca[32]; int c_triangle[2] = { -1, -1 }; @@ -1850,11 +1850,8 @@ static void set_matrix_indices(Lincs* li, const Task& li_task, const ListOfLists } } -void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li) +void set_lincs(const InteractionDefinitions& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li) { - int natoms; - t_iatom* iatom; - li->nc_real = 0; li->nc = 0; li->ncc = 0; @@ -1873,7 +1870,7 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ } /* This is the local topology, so there are only F_CONSTR constraints */ - if (idef.il[F_CONSTR].nr == 0) + if (idef.il[F_CONSTR].empty()) { /* There are no constraints, * we do not need to fill any data structures. @@ -1886,6 +1883,7 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ fprintf(debug, "Building the LINCS connectivity\n"); } + int natoms; if (DOMAINDECOMP(cr)) { if (cr->dd->constraints) @@ -1907,7 +1905,7 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ const ListOfLists at2con = make_at2con(natoms, idef.il, idef.iparams, flexibleConstraintTreatment(bDynamics)); - const int ncon_tot = idef.il[F_CONSTR].nr / 3; + const int ncon_tot = idef.il[F_CONSTR].size() / 3; /* Ensure we have enough padding for aligned loads for each thread */ const int numEntries = ncon_tot + li->ntask * simd_width; @@ -1930,7 +1928,7 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ li->tmp4.resize(numEntries); li->mlambda.resize(numEntries); - iatom = idef.il[F_CONSTR].iatoms; + gmx::ArrayRef iatom = idef.il[F_CONSTR].iatoms; li->blnr[0] = li->ncc; @@ -1993,6 +1991,8 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ */ li_task->b0 = li->nc; + gmx::ArrayRef iparams = idef.iparams; + while (con < ncon_tot && li->nc - li_task->b0 < ncon_target) { if (li->con_index[con] == -1) @@ -2003,8 +2003,8 @@ void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_ type = iatom[3 * con]; a1 = iatom[3 * con + 1]; a2 = iatom[3 * con + 2]; - lenA = idef.iparams[type].constr.dA; - lenB = idef.iparams[type].constr.dB; + lenA = iparams[type].constr.dA; + lenB = iparams[type].constr.dB; /* Skip the flexible constraints when not doing dynamics */ if (bDynamics || lenA != 0 || lenB != 0) { diff --git a/src/gromacs/mdlib/lincs.h b/src/gromacs/mdlib/lincs.h index b891a64803..2ab8ee8062 100644 --- a/src/gromacs/mdlib/lincs.h +++ b/src/gromacs/mdlib/lincs.h @@ -53,8 +53,8 @@ struct gmx_mtop_t; struct gmx_multisim_t; +class InteractionDefinitions; struct t_commrec; -struct t_idef; struct t_inputrec; struct t_mdatoms; struct t_nrnb; @@ -90,7 +90,11 @@ Lincs* init_lincs(FILE* fplog, void done_lincs(Lincs* li); /*! \brief Initialize lincs stuff */ -void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li); +void set_lincs(const InteractionDefinitions& idef, + const t_mdatoms& md, + bool bDynamics, + const t_commrec* cr, + Lincs* li); /*! \brief Applies LINCS constraints. * diff --git a/src/gromacs/mdlib/lincs_gpu.cu b/src/gromacs/mdlib/lincs_gpu.cu index c5e92b7c29..edf3e9c58a 100644 --- a/src/gromacs/mdlib/lincs_gpu.cu +++ b/src/gromacs/mdlib/lincs_gpu.cu @@ -65,6 +65,7 @@ #include "gromacs/mdlib/constr.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pbcutil/pbc_aiuc_cuda.cuh" +#include "gromacs/topology/forcefieldparameters.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/topology.h" @@ -720,7 +721,7 @@ bool LincsGpu::isNumCoupledConstraintsSupported(const gmx_mtop_t& mtop) return true; } -void LincsGpu::set(const t_idef& idef, const t_mdatoms& md) +void LincsGpu::set(const InteractionDefinitions& idef, const t_mdatoms& md) { int numAtoms = md.nr; // List of constrained atoms (CPU memory) @@ -735,10 +736,9 @@ void LincsGpu::set(const t_idef& idef, const t_mdatoms& md) std::vector massFactorsHost; // List of constrained atoms in local topology - gmx::ArrayRef iatoms = - constArrayRefFromArray(idef.il[F_CONSTR].iatoms, idef.il[F_CONSTR].nr); - const int stride = NRAL(F_CONSTR) + 1; - const int numConstraints = idef.il[F_CONSTR].nr / stride; + ArrayRef iatoms = idef.il[F_CONSTR].iatoms; + const int stride = NRAL(F_CONSTR) + 1; + const int numConstraints = idef.il[F_CONSTR].size() / stride; // Early exit if no constraints if (numConstraints == 0) diff --git a/src/gromacs/mdlib/lincs_gpu.cuh b/src/gromacs/mdlib/lincs_gpu.cuh index 40cff3c340..0fbebcf67f 100644 --- a/src/gromacs/mdlib/lincs_gpu.cuh +++ b/src/gromacs/mdlib/lincs_gpu.cuh @@ -48,9 +48,10 @@ #include "gromacs/mdlib/constr.h" #include "gromacs/mdtypes/mdatom.h" #include "gromacs/pbcutil/pbc_aiuc.h" -#include "gromacs/topology/idef.h" #include "gromacs/utility/classhelpers.h" +class InteractionDefinitions; + namespace gmx { @@ -155,7 +156,7 @@ public: * \param[in] idef Local topology data to get information on constraints from. * \param[in] md Atoms data to get atom masses from. */ - void set(const t_idef& idef, const t_mdatoms& md); + void set(const InteractionDefinitions& idef, const t_mdatoms& md); /*! \brief * Returns whether the maximum number of coupled constraints is supported diff --git a/src/gromacs/mdlib/settle.cpp b/src/gromacs/mdlib/settle.cpp index 9cbe5366ea..768ebdd6b8 100644 --- a/src/gromacs/mdlib/settle.cpp +++ b/src/gromacs/mdlib/settle.cpp @@ -242,7 +242,7 @@ void settle_free(settledata* settled) sfree(settled); } -void settle_set_constraints(settledata* settled, const t_ilist* il_settle, const t_mdatoms& mdatoms) +void settle_set_constraints(settledata* settled, const InteractionList& il_settle, const t_mdatoms& mdatoms) { #if GMX_SIMD_HAVE_REAL const int pack_size = GMX_SIMD_REAL_WIDTH; @@ -251,12 +251,12 @@ void settle_set_constraints(settledata* settled, const t_ilist* il_settle, const #endif const int nral1 = 1 + NRAL(F_SETTLE); - int nsettle = il_settle->nr / nral1; + int nsettle = il_settle.size() / nral1; settled->nsettle = nsettle; if (nsettle > 0) { - const t_iatom* iatoms = il_settle->iatoms; + ArrayRef iatoms = il_settle.iatoms; /* Here we initialize the normal SETTLE parameters */ if (settled->massw.mO < 0) diff --git a/src/gromacs/mdlib/settle.h b/src/gromacs/mdlib/settle.h index 20ea546421..05bccf307e 100644 --- a/src/gromacs/mdlib/settle.h +++ b/src/gromacs/mdlib/settle.h @@ -46,8 +46,8 @@ #include "gromacs/topology/idef.h" -struct gmx_cmap_t; struct gmx_mtop_t; +struct InteractionList; struct t_inputrec; struct t_mdatoms; struct t_pbc; @@ -71,7 +71,7 @@ settledata* settle_init(const gmx_mtop_t& mtop); void settle_free(settledata* settled); /*! \brief Set up the indices for the settle constraints */ -void settle_set_constraints(settledata* settled, const t_ilist* il_settle, const t_mdatoms& mdatoms); +void settle_set_constraints(settledata* settled, const InteractionList& il_settle, const t_mdatoms& mdatoms); /*! \brief Constrain coordinates using SETTLE. * Can be called on any number of threads. @@ -95,7 +95,7 @@ void csettle(settledata* settled, /* The SETTLE structur void settle_proj(settledata* settled, ConstraintVariable econq, int nsettle, - const t_iatom iatoms[], + const int iatoms[], const t_pbc* pbc, /* PBC data pointer, can be NULL */ ArrayRef x, ArrayRef der, diff --git a/src/gromacs/mdlib/settle_gpu.cu b/src/gromacs/mdlib/settle_gpu.cu index 2707d1744b..d1e8f508d5 100644 --- a/src/gromacs/mdlib/settle_gpu.cu +++ b/src/gromacs/mdlib/settle_gpu.cu @@ -604,12 +604,12 @@ SettleGpu::~SettleGpu() } } -void SettleGpu::set(const t_idef& idef, const t_mdatoms gmx_unused& md) +void SettleGpu::set(const InteractionDefinitions& idef, const t_mdatoms gmx_unused& md) { - const int nral1 = 1 + NRAL(F_SETTLE); - t_ilist il_settle = idef.il[F_SETTLE]; - t_iatom* iatoms = il_settle.iatoms; - numSettles_ = il_settle.nr / nral1; + const int nral1 = 1 + NRAL(F_SETTLE); + const InteractionList& il_settle = idef.il[F_SETTLE]; + ArrayRef iatoms = il_settle.iatoms; + numSettles_ = il_settle.size() / nral1; reallocateDeviceBuffer(&d_atomIds_, numSettles_, &numAtomIds_, &numAtomIdsAlloc_, nullptr); h_atomIds_.resize(numSettles_); diff --git a/src/gromacs/mdlib/settle_gpu.cuh b/src/gromacs/mdlib/settle_gpu.cuh index a22f9cb670..3816579d63 100644 --- a/src/gromacs/mdlib/settle_gpu.cuh +++ b/src/gromacs/mdlib/settle_gpu.cuh @@ -52,9 +52,10 @@ #include "gromacs/mdtypes/mdatom.h" #include "gromacs/pbcutil/pbc.h" #include "gromacs/pbcutil/pbc_aiuc.h" -#include "gromacs/topology/idef.h" #include "gromacs/topology/topology.h" +class InteractionDefinitions; + namespace gmx { @@ -246,7 +247,7 @@ public: * \param[in] idef System topology * \param[in] md Atoms data. Can be used to update masses if needed (not used now). */ - void set(const t_idef& idef, const t_mdatoms& md); + void set(const InteractionDefinitions& idef, const t_mdatoms& md); private: //! GPU stream diff --git a/src/gromacs/mdlib/shake.cpp b/src/gromacs/mdlib/shake.cpp index 87d7728135..119fe924fb 100644 --- a/src/gromacs/mdlib/shake.cpp +++ b/src/gromacs/mdlib/shake.cpp @@ -179,23 +179,22 @@ static void resizeLagrangianData(shakedata* shaked, int ncons) } } -void make_shake_sblock_serial(shakedata* shaked, const t_idef* idef, const t_mdatoms& md) +void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md) { int i, j, m, ncons; int bstart, bnr; t_blocka sblocks; t_sortblock* sb; - t_iatom* iatom; int* inv_sblock; /* Since we are processing the local topology, * the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist. */ - ncons = idef->il[F_CONSTR].nr / 3; + ncons = idef->il[F_CONSTR].size() / 3; init_blocka(&sblocks); sfree(sblocks.index); // To solve memory leak - gen_sblocks(nullptr, 0, md.homenr, idef, &sblocks, FALSE); + gen_sblocks(nullptr, 0, md.homenr, *idef, &sblocks, FALSE); /* bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0; @@ -217,7 +216,7 @@ void make_shake_sblock_serial(shakedata* shaked, const t_idef* idef, const t_mda * sort the constraints in order of the sblock number * and the atom numbers, really sorting a segment of the array! */ - iatom = idef->il[F_CONSTR].iatoms; + int* iatom = idef->il[F_CONSTR].iatoms.data(); snew(sb, ncons); for (i = 0; (i < ncons); i++, iatom += 3) { @@ -242,7 +241,7 @@ void make_shake_sblock_serial(shakedata* shaked, const t_idef* idef, const t_mda pr_sortblock(debug, "After sorting", ncons, sb); } - iatom = idef->il[F_CONSTR].iatoms; + iatom = idef->il[F_CONSTR].iatoms.data(); for (i = 0; (i < ncons); i++, iatom += 3) { for (m = 0; (m < 3); m++) @@ -287,10 +286,9 @@ void make_shake_sblock_serial(shakedata* shaked, const t_idef* idef, const t_mda } // TODO: Check if this code is useful. It might never be called. -void make_shake_sblock_dd(shakedata* shaked, const t_ilist* ilcon, const gmx_domdec_t* dd) +void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon, const gmx_domdec_t* dd) { - int ncons, c, cg; - t_iatom* iatom; + int ncons, c, cg; if (dd->ncg_home + 1 > shaked->sblock_nalloc) { @@ -298,10 +296,10 @@ void make_shake_sblock_dd(shakedata* shaked, const t_ilist* ilcon, const gmx_dom srenew(shaked->sblock, shaked->sblock_nalloc); } - ncons = ilcon->nr / 3; - iatom = ilcon->iatoms; - shaked->nblocks = 0; - cg = 0; + ncons = ilcon.size() / 3; + const int* iatom = ilcon.iatoms.data(); + shaked->nblocks = 0; + cg = 0; for (c = 0; c < ncons; c++) { if (c == 0 || iatom[1] >= cg + 1) @@ -538,34 +536,33 @@ static void crattle(const int iatom[], } //! Applies SHAKE -static int vec_shakef(FILE* fplog, - shakedata* shaked, - const real invmass[], - int ncon, - t_iparams ip[], - t_iatom* iatom, - real tol, - const rvec x[], - rvec prime[], - real omega, - bool bFEP, - real lambda, - real scaled_lagrange_multiplier[], - real invdt, - rvec* v, - bool bCalcVir, - tensor vir_r_m_dr, - ConstraintVariable econq) +static int vec_shakef(FILE* fplog, + shakedata* shaked, + const real invmass[], + int ncon, + ArrayRef ip, + const int* iatom, + real tol, + const rvec x[], + rvec prime[], + real omega, + bool bFEP, + real lambda, + real scaled_lagrange_multiplier[], + real invdt, + rvec* v, + bool bCalcVir, + tensor vir_r_m_dr, + ConstraintVariable econq) { - rvec* rij; - real * half_of_reduced_mass, *distance_squared_tolerance, *constraint_distance_squared; - int maxnit = 1000; - int nit = 0, ll, i, j, d, d2, type; - t_iatom* ia; - real L1; - real mm = 0., tmp; - int error = 0; - real constraint_distance; + rvec* rij; + real *half_of_reduced_mass, *distance_squared_tolerance, *constraint_distance_squared; + int maxnit = 1000; + int nit = 0, ll, i, j, d, d2, type; + real L1; + real mm = 0., tmp; + int error = 0; + real constraint_distance; if (ncon > shaked->nalloc) { @@ -580,8 +577,8 @@ static int vec_shakef(FILE* fplog, distance_squared_tolerance = shaked->distance_squared_tolerance; constraint_distance_squared = shaked->constraint_distance_squared; - L1 = 1.0 - lambda; - ia = iatom; + L1 = 1.0 - lambda; + const int* ia = iatom; for (ll = 0; (ll < ncon); ll++, ia += 3) { type = ia[0]; @@ -703,25 +700,24 @@ static int vec_shakef(FILE* fplog, } //! Check that constraints are satisfied. -static void check_cons(FILE* log, - int nc, - const rvec x[], - rvec prime[], - rvec v[], - t_iparams ip[], - t_iatom* iatom, - const real invmass[], - ConstraintVariable econq) +static void check_cons(FILE* log, + int nc, + const rvec x[], + rvec prime[], + rvec v[], + ArrayRef ip, + const int* iatom, + const real invmass[], + ConstraintVariable econq) { - t_iatom* ia; - int ai, aj; - int i; - real d, dp; - rvec dx, dv; + int ai, aj; + int i; + real d, dp; + rvec dx, dv; GMX_ASSERT(v, "Input has to be non-null"); fprintf(log, " i mi j mj before after should be\n"); - ia = iatom; + const int* ia = iatom; for (i = 0; (i < nc); i++, ia += 3) { ai = ia[1]; @@ -751,29 +747,28 @@ static void check_cons(FILE* log, } //! Applies SHAKE. -static bool bshakef(FILE* log, - shakedata* shaked, - const real invmass[], - const t_idef& idef, - const t_inputrec& ir, - const rvec x_s[], - rvec prime[], - t_nrnb* nrnb, - real lambda, - real* dvdlambda, - real invdt, - rvec* v, - bool bCalcVir, - tensor vir_r_m_dr, - bool bDumpOnError, - ConstraintVariable econq) +static bool bshakef(FILE* log, + shakedata* shaked, + const real invmass[], + const InteractionDefinitions& idef, + const t_inputrec& ir, + const rvec x_s[], + rvec prime[], + t_nrnb* nrnb, + real lambda, + real* dvdlambda, + real invdt, + rvec* v, + bool bCalcVir, + tensor vir_r_m_dr, + bool bDumpOnError, + ConstraintVariable econq) { - t_iatom* iatoms; - real * lam, dt_2, dvdl; - int i, n0, ncon, blen, type, ll; - int tnit = 0, trij = 0; + real *lam, dt_2, dvdl; + int i, n0, ncon, blen, type, ll; + int tnit = 0, trij = 0; - ncon = idef.il[F_CONSTR].nr / 3; + ncon = idef.il[F_CONSTR].size() / 3; for (ll = 0; ll < ncon; ll++) { @@ -783,8 +778,8 @@ static bool bshakef(FILE* log, // TODO Rewrite this block so that it is obvious that i, iatoms // and lam are all iteration variables. Is this easier if the // sblock data structure is organized differently? - iatoms = &(idef.il[F_CONSTR].iatoms[shaked->sblock[0]]); - lam = shaked->scaled_lagrange_multiplier; + const int* iatoms = &(idef.il[F_CONSTR].iatoms[shaked->sblock[0]]); + lam = shaked->scaled_lagrange_multiplier; for (i = 0; (i < shaked->nblocks);) { blen = (shaked->sblock[i + 1] - shaked->sblock[i]); @@ -814,7 +809,8 @@ static bool bshakef(FILE* log, { if (ir.efep != efepNO) { - real bondA, bondB; + ArrayRef iparams = idef.iparams; + /* TODO This should probably use invdt, so that sd integrator scaling works properly */ dt_2 = 1 / gmx::square(ir.delta_t); dvdl = 0; @@ -826,8 +822,8 @@ static bool bshakef(FILE* log, /* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll (eta_ll is the estimate of the Langrangian, definition on page 336 of Ryckaert et al 1977), so the pre-factors are already present. */ - bondA = idef.iparams[type].constr.dA; - bondB = idef.iparams[type].constr.dB; + const real bondA = iparams[type].constr.dA; + const real bondB = iparams[type].constr.dB; dvdl += shaked->scaled_lagrange_multiplier[ll] * dt_2 * (bondB - bondA); } *dvdlambda += dvdl; @@ -856,23 +852,23 @@ static bool bshakef(FILE* log, return TRUE; } -bool constrain_shake(FILE* log, - shakedata* shaked, - const real invmass[], - const t_idef& idef, - const t_inputrec& ir, - const rvec x_s[], - rvec xprime[], - rvec vprime[], - t_nrnb* nrnb, - real lambda, - real* dvdlambda, - real invdt, - rvec* v, - bool bCalcVir, - tensor vir_r_m_dr, - bool bDumpOnError, - ConstraintVariable econq) +bool constrain_shake(FILE* log, + shakedata* shaked, + const real invmass[], + const InteractionDefinitions& idef, + const t_inputrec& ir, + const rvec x_s[], + rvec xprime[], + rvec vprime[], + t_nrnb* nrnb, + real lambda, + real* dvdlambda, + real invdt, + rvec* v, + bool bCalcVir, + tensor vir_r_m_dr, + bool bDumpOnError, + ConstraintVariable econq) { if (shaked->nblocks == 0) { diff --git a/src/gromacs/mdlib/shake.h b/src/gromacs/mdlib/shake.h index 1249d3c1f5..fb8b5cdfdc 100644 --- a/src/gromacs/mdlib/shake.h +++ b/src/gromacs/mdlib/shake.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -45,10 +45,13 @@ #ifndef GMX_MDLIB_SHAKE_H #define GMX_MDLIB_SHAKE_H +#include "gromacs/math/vec.h" #include "gromacs/topology/block.h" -#include "gromacs/topology/idef.h" +#include "gromacs/utility/real.h" struct gmx_domdec_t; +struct InteractionList; +class InteractionDefinitions; struct t_inputrec; struct t_mdatoms; struct t_nrnb; @@ -68,10 +71,10 @@ shakedata* shake_init(); void done_shake(shakedata* d); //! Make SHAKE blocks when not using DD. -void make_shake_sblock_serial(shakedata* shaked, const t_idef* idef, const t_mdatoms& md); +void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md); //! Make SHAKE blocks when using DD. -void make_shake_sblock_dd(shakedata* shaked, const t_ilist* ilcon, const gmx_domdec_t* dd); +void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon, const gmx_domdec_t* dd); /*! \brief Shake all the atoms blockwise. It is assumed that all the constraints * in the idef->shakes field are sorted, to ascending block nr. The @@ -81,23 +84,23 @@ void make_shake_sblock_dd(shakedata* shaked, const t_ilist* ilcon, const gmx_dom * sblock[n] to sblock[n+1]. Array sblock should be large enough. * Return TRUE when OK, FALSE when shake-error */ -bool constrain_shake(FILE* log, /* Log file */ - shakedata* shaked, /* Total number of atoms */ - const real invmass[], /* Atomic masses */ - const t_idef& idef, /* The interaction def */ - const t_inputrec& ir, /* Input record */ - const rvec x_s[], /* Coords before update */ - rvec xprime[], /* Output coords when constraining x */ - rvec vprime[], /* Output coords when constraining v */ - t_nrnb* nrnb, /* Performance measure */ - real lambda, /* FEP lambda */ - real* dvdlambda, /* FEP force */ - real invdt, /* 1/delta_t */ - rvec* v, /* Also constrain v if v!=NULL */ - bool bCalcVir, /* Calculate r x m delta_r */ - tensor vir_r_m_dr, /* sum r x m delta_r */ - bool bDumpOnError, /* Dump debugging stuff on error*/ - ConstraintVariable econq); /* which type of constraint is occurring */ +bool constrain_shake(FILE* log, /* Log file */ + shakedata* shaked, /* Total number of atoms */ + const real invmass[], /* Atomic masses */ + const InteractionDefinitions& idef, /* The interaction def */ + const t_inputrec& ir, /* Input record */ + const rvec x_s[], /* Coords before update */ + rvec xprime[], /* Output coords when constraining x */ + rvec vprime[], /* Output coords when constraining v */ + t_nrnb* nrnb, /* Performance measure */ + real lambda, /* FEP lambda */ + real* dvdlambda, /* FEP force */ + real invdt, /* 1/delta_t */ + rvec* v, /* Also constrain v if v!=NULL */ + bool bCalcVir, /* Calculate r x m delta_r */ + tensor vir_r_m_dr, /* sum r x m delta_r */ + bool bDumpOnError, /* Dump debugging stuff on error*/ + ConstraintVariable econq); /* which type of constraint is occurring */ /*! \brief Regular iterative shake */ void cshake(const int iatom[], diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp index dc5538bf48..3dccd0e5d4 100644 --- a/src/gromacs/mdlib/sim_util.cpp +++ b/src/gromacs/mdlib/sim_util.cpp @@ -308,7 +308,7 @@ static void post_process_forces(const t_commrec* cr, */ matrix virial = { { 0 } }; spread_vsite_f(vsite, x, fDirectVir, nullptr, stepWork.computeVirial, virial, nrnb, - &top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle); + top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle); forceWithVirial.addVirialContribution(virial); } @@ -785,13 +785,13 @@ static ForceOutputs setupForceOutputs(t_forcerec* fr, /*! \brief Set up flags that have the lifetime of the domain indicating what type of work is there to compute. */ -static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec& inputrec, - const t_forcerec& fr, - const pull_t* pull_work, - const gmx_edsam* ed, - const t_idef& idef, - const t_fcdata& fcd, - const t_mdatoms& mdatoms, +static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec& inputrec, + const t_forcerec& fr, + const pull_t* pull_work, + const gmx_edsam* ed, + const InteractionDefinitions& idef, + const t_fcdata& fcd, + const t_mdatoms& mdatoms, const SimulationWorkload& simulationWork, const StepWorkload& stepWork) { @@ -1535,7 +1535,7 @@ void do_force(FILE* fplog, stateGpu->waitCoordinatesReadyOnHost(AtomLocality::NonLocal); } /* Compute the bonded and non-bonded energies and optionally forces */ - do_force_lowlevel(fr, inputrec, &(top->idef), cr, ms, nrnb, wcycle, mdatoms, x, hist, &forceOut, enerd, + do_force_lowlevel(fr, inputrec, top->idef, cr, ms, nrnb, wcycle, mdatoms, x, hist, &forceOut, enerd, fcd, box, lambda.data(), graph, fr->mu_tot, stepWork, ddBalanceRegionHandler); wallcycle_stop(wcycle, ewcFORCE); @@ -1810,8 +1810,8 @@ void do_force(FILE* fplog, if (vsite && !(fr->haveDirectVirialContributions && !stepWork.computeVirial)) { rvec* fshift = as_rvec_array(forceOut.forceWithShiftForces().shiftForces().data()); - spread_vsite_f(vsite, as_rvec_array(x.unpaddedArrayRef().data()), f, fshift, FALSE, nullptr, - nrnb, &top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle); + spread_vsite_f(vsite, as_rvec_array(x.unpaddedArrayRef().data()), f, fshift, FALSE, + nullptr, nrnb, top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle); } if (stepWork.computeVirial) diff --git a/src/gromacs/mdlib/splitter.cpp b/src/gromacs/mdlib/splitter.cpp index e1b538fa96..d16faf9f9f 100644 --- a/src/gromacs/mdlib/splitter.cpp +++ b/src/gromacs/mdlib/splitter.cpp @@ -347,7 +347,7 @@ static int merge_sid(int at_start, int at_end, int nsid, t_sid sid[], t_blocka* return nsid; } -void gen_sblocks(FILE* fp, int at_start, int at_end, const t_idef* idef, t_blocka* sblock, gmx_bool bSettle) +void gen_sblocks(FILE* fp, int at_start, int at_end, const InteractionDefinitions& idef, t_blocka* sblock, gmx_bool bSettle) { t_graph* g; int i, i0; diff --git a/src/gromacs/mdlib/splitter.h b/src/gromacs/mdlib/splitter.h index 88fb8f205d..0541b034df 100644 --- a/src/gromacs/mdlib/splitter.h +++ b/src/gromacs/mdlib/splitter.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2010,2014,2015,2018,2019, by the GROMACS development team, led by + * Copyright (c) 2010,2014,2015,2018,2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -41,10 +41,10 @@ #include "gromacs/utility/basedefinitions.h" +class InteractionDefinitions; struct t_blocka; -struct t_idef; -void gen_sblocks(FILE* fp, int at_start, int at_end, const t_idef* idef, t_blocka* sblock, gmx_bool bSettle); +void gen_sblocks(FILE* fp, int at_start, int at_end, const InteractionDefinitions& idef, t_blocka* sblock, gmx_bool bSettle); /* Generate shake blocks from the constraint list. Set bSettle to yes for shake * blocks including settles. You normally do not want this. */ diff --git a/src/gromacs/mdlib/tests/constrtestdata.cpp b/src/gromacs/mdlib/tests/constrtestdata.cpp index 0679b51713..2062db5534 100644 --- a/src/gromacs/mdlib/tests/constrtestdata.cpp +++ b/src/gromacs/mdlib/tests/constrtestdata.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -145,31 +145,25 @@ ConstraintsTestData::ConstraintsTestData(const std::string& title, dHdLambdaRef_ = 0; } - // Constraints and their parameters (local topology) - for (int i = 0; i < F_NRE; i++) - { - idef_.il[i].nr = 0; - } - idef_.il[F_CONSTR].nr = constraints.size(); - - snew(idef_.il[F_CONSTR].iatoms, constraints.size()); int maxType = 0; - for (index i = 0; i < ssize(constraints); i++) + for (index i = 0; i < ssize(constraints); i += 3) { - if (i % 3 == 0) + if (maxType < constraints.at(i)) { - if (maxType < constraints.at(i)) - { - maxType = constraints.at(i); - } + maxType = constraints.at(i); } - idef_.il[F_CONSTR].iatoms[i] = constraints.at(i); } - snew(idef_.iparams, maxType + 1); + auto& iparams = mtop_.ffparams.iparams; + iparams.resize(maxType + 1); for (index i = 0; i < ssize(constraints) / 3; i++) { - idef_.iparams[constraints.at(3 * i)].constr.dA = constraintsR0.at(constraints.at(3 * i)); - idef_.iparams[constraints.at(3 * i)].constr.dB = constraintsR0.at(constraints.at(3 * i)); + iparams[constraints.at(3 * i)].constr.dA = constraintsR0.at(constraints.at(3 * i)); + iparams[constraints.at(3 * i)].constr.dB = constraintsR0.at(constraints.at(3 * i)); + } + idef_ = std::make_unique(mtop_.ffparams); + for (index i = 0; i < ssize(constraints); i++) + { + idef_->il[F_CONSTR].iatoms.push_back(constraints.at(i)); } // Constraints and their parameters (global topology) @@ -190,12 +184,7 @@ ConstraintsTestData::ConstraintsTestData(const std::string& title, molBlock.nmol = 1; mtop_.molblock.push_back(molBlock); - mtop_.natoms = numAtoms; - mtop_.ffparams.iparams.resize(maxType + 1); - for (int i = 0; i <= maxType; i++) - { - mtop_.ffparams.iparams.at(i) = idef_.iparams[i]; - } + mtop_.natoms = numAtoms; mtop_.bIntermolecularInteractions = false; // Coordinates and velocities @@ -251,14 +240,5 @@ void ConstraintsTestData::reset() dHdLambda_ = 0; } -/*! \brief - * Cleaning up the memory. - */ -ConstraintsTestData::~ConstraintsTestData() -{ - sfree(idef_.il[F_CONSTR].iatoms); - sfree(idef_.iparams); -} - } // namespace test } // namespace gmx diff --git a/src/gromacs/mdlib/tests/constrtestdata.h b/src/gromacs/mdlib/tests/constrtestdata.h index 4283674b1f..17a82e00bc 100644 --- a/src/gromacs/mdlib/tests/constrtestdata.h +++ b/src/gromacs/mdlib/tests/constrtestdata.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -46,6 +46,7 @@ #ifndef GMX_MDLIB_TESTS_CONSTRTESTDATA_H #define GMX_MDLIB_TESTS_CONSTRTESTDATA_H +#include #include #include "gromacs/gmxlib/nrnb.h" @@ -94,7 +95,7 @@ public: //! Input record (info that usually in .mdp file) t_inputrec ir_; //! Local topology - t_idef idef_; + std::unique_ptr idef_; //! MD atoms t_mdatoms md_; //! Multisim data @@ -206,11 +207,6 @@ public: * */ void reset(); - - /*! \brief - * Cleaning up the memory. - */ - ~ConstraintsTestData(); }; } // namespace test diff --git a/src/gromacs/mdlib/tests/constrtestrunners.cpp b/src/gromacs/mdlib/tests/constrtestrunners.cpp index 5f974aaf8b..243a0200c4 100644 --- a/src/gromacs/mdlib/tests/constrtestrunners.cpp +++ b/src/gromacs/mdlib/tests/constrtestrunners.cpp @@ -89,9 +89,9 @@ namespace test void applyShake(ConstraintsTestData* testData, t_pbc gmx_unused pbc) { shakedata* shaked = shake_init(); - make_shake_sblock_serial(shaked, &testData->idef_, testData->md_); + make_shake_sblock_serial(shaked, testData->idef_.get(), testData->md_); bool success = constrain_shake( - nullptr, shaked, testData->invmass_.data(), testData->idef_, testData->ir_, + nullptr, shaked, testData->invmass_.data(), *testData->idef_, testData->ir_, as_rvec_array(testData->x_.data()), as_rvec_array(testData->xPrime_.data()), as_rvec_array(testData->xPrime2_.data()), &testData->nrnb_, testData->md_.lambda, &testData->dHdLambda_, testData->invdt_, as_rvec_array(testData->v_.data()), @@ -126,7 +126,7 @@ void applyLincs(ConstraintsTestData* testData, t_pbc pbc) // Initialize LINCS lincsd = init_lincs(nullptr, testData->mtop_, testData->nflexcon_, at2con_mt, false, testData->ir_.nLincsIter, testData->ir_.nProjOrder); - set_lincs(testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &testData->cr_, lincsd); + set_lincs(*testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &testData->cr_, lincsd); // Evaluate constraints bool success = constrain_lincs( diff --git a/src/gromacs/mdlib/tests/constrtestrunners.cu b/src/gromacs/mdlib/tests/constrtestrunners.cu index 6d54523506..322c3beec9 100644 --- a/src/gromacs/mdlib/tests/constrtestrunners.cu +++ b/src/gromacs/mdlib/tests/constrtestrunners.cu @@ -77,7 +77,7 @@ void applyLincsGpu(ConstraintsTestData* testData, t_pbc pbc) int numAtoms = testData->numAtoms_; float3 *d_x, *d_xp, *d_v; - lincsGpu->set(testData->idef_, testData->md_); + lincsGpu->set(*testData->idef_, testData->md_); PbcAiuc pbcAiuc; setPbcAiuc(pbc.ndim_ePBC, pbc.box, &pbcAiuc); diff --git a/src/gromacs/mdlib/tests/settletestdata.cpp b/src/gromacs/mdlib/tests/settletestdata.cpp index 4cea8e864a..4f9e6b56cf 100644 --- a/src/gromacs/mdlib/tests/settletestdata.cpp +++ b/src/gromacs/mdlib/tests/settletestdata.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2014,2015,2016,2018,2019, by the GROMACS development team, led by + * Copyright (c) 2014,2015,2016,2018,2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -138,9 +138,8 @@ SettleTestData::SettleTestData(int numSettles) : mtop_.moltype[0].atoms.atom[i * atomsPerSettle_ + 2].m = hydrogenMass_; } - // Reshape some data so it can be directly used by the SETTLE constraints - ilist_ = { mtop_.moltype[0].ilist[F_SETTLE].size(), mtop_.moltype[0].ilist[F_SETTLE].iatoms.data(), 0 }; - idef_.il[F_SETTLE] = ilist_; + idef_ = std::make_unique(mtop_.ffparams); + idef_->il[F_SETTLE] = mtop_.moltype[0].ilist[F_SETTLE]; } SettleTestData::~SettleTestData() diff --git a/src/gromacs/mdlib/tests/settletestdata.h b/src/gromacs/mdlib/tests/settletestdata.h index ad07c6501b..4ca9d3d439 100644 --- a/src/gromacs/mdlib/tests/settletestdata.h +++ b/src/gromacs/mdlib/tests/settletestdata.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -76,10 +76,8 @@ public: gmx_mtop_t mtop_; //! Atoms data t_mdatoms mdatoms_; - //! Interactions list - t_ilist ilist_; //! Local topology - t_idef idef_; + std::unique_ptr idef_; //! Inverse timestep const real reciprocalTimeStep_ = 1.0 / 0.002; diff --git a/src/gromacs/mdlib/tests/settletestrunners.cpp b/src/gromacs/mdlib/tests/settletestrunners.cpp index 59cc352b87..9e501a0fab 100644 --- a/src/gromacs/mdlib/tests/settletestrunners.cpp +++ b/src/gromacs/mdlib/tests/settletestrunners.cpp @@ -65,7 +65,7 @@ void applySettle(SettleTestData* testData, { settledata* settled = settle_init(testData->mtop_); - settle_set_constraints(settled, &testData->ilist_, testData->mdatoms_); + settle_set_constraints(settled, testData->idef_->il[F_SETTLE], testData->mdatoms_); bool errorOccured; int numThreads = 1; diff --git a/src/gromacs/mdlib/tests/settletestrunners.cu b/src/gromacs/mdlib/tests/settletestrunners.cu index bbf1083507..a14b47e819 100644 --- a/src/gromacs/mdlib/tests/settletestrunners.cu +++ b/src/gromacs/mdlib/tests/settletestrunners.cu @@ -88,7 +88,7 @@ void applySettleGpu(SettleTestData* testData, auto settleGpu = std::make_unique(testData->mtop_, nullptr); - settleGpu->set(testData->idef_, testData->mdatoms_); + settleGpu->set(*testData->idef_, testData->mdatoms_); PbcAiuc pbcAiuc; setPbcAiuc(pbc.ndim_ePBC, pbc.box, &pbcAiuc); diff --git a/src/gromacs/mdlib/update_constrain_gpu.h b/src/gromacs/mdlib/update_constrain_gpu.h index 359ecea5c1..09f0bbecc1 100644 --- a/src/gromacs/mdlib/update_constrain_gpu.h +++ b/src/gromacs/mdlib/update_constrain_gpu.h @@ -53,7 +53,7 @@ class GpuEventSynchronizer; struct gmx_mtop_t; enum class PbcType : int; -struct t_idef; +class InteractionDefinitions; struct t_inputrec; struct t_mdatoms; struct t_pbc; @@ -133,12 +133,12 @@ public: * \param[in] md Atoms data. * \param[in] numTempScaleValues Number of temperature scaling groups. Zero for no temperature scaling. */ - void set(DeviceBuffer d_x, - DeviceBuffer d_v, - DeviceBuffer d_f, - const t_idef& idef, - const t_mdatoms& md, - int numTempScaleValues); + void set(DeviceBuffer d_x, + DeviceBuffer d_v, + DeviceBuffer d_f, + const InteractionDefinitions& idef, + const t_mdatoms& md, + int numTempScaleValues); /*! \brief * Update PBC data. diff --git a/src/gromacs/mdlib/update_constrain_gpu_impl.cpp b/src/gromacs/mdlib/update_constrain_gpu_impl.cpp index 47671ef7de..3e10f8a403 100644 --- a/src/gromacs/mdlib/update_constrain_gpu_impl.cpp +++ b/src/gromacs/mdlib/update_constrain_gpu_impl.cpp @@ -91,7 +91,7 @@ void UpdateConstrainGpu::scaleCoordinates(const matrix /* scalingMatrix */) void UpdateConstrainGpu::set(DeviceBuffer /* d_x */, DeviceBuffer /* d_v */, const DeviceBuffer /* d_f */, - const t_idef& /* idef */, + const InteractionDefinitions& /* idef */, const t_mdatoms& /* md */, const int /* numTempScaleValues */) { diff --git a/src/gromacs/mdlib/update_constrain_gpu_impl.cu b/src/gromacs/mdlib/update_constrain_gpu_impl.cu index 6991ef0dc3..c77e1924ed 100644 --- a/src/gromacs/mdlib/update_constrain_gpu_impl.cu +++ b/src/gromacs/mdlib/update_constrain_gpu_impl.cu @@ -188,12 +188,12 @@ UpdateConstrainGpu::Impl::Impl(const t_inputrec& ir, UpdateConstrainGpu::Impl::~Impl() {} -void UpdateConstrainGpu::Impl::set(DeviceBuffer d_x, - DeviceBuffer d_v, - const DeviceBuffer d_f, - const t_idef& idef, - const t_mdatoms& md, - const int numTempScaleValues) +void UpdateConstrainGpu::Impl::set(DeviceBuffer d_x, + DeviceBuffer d_v, + const DeviceBuffer d_f, + const InteractionDefinitions& idef, + const t_mdatoms& md, + const int numTempScaleValues) { GMX_ASSERT(d_x != nullptr, "Coordinates device buffer should not be null."); GMX_ASSERT(d_v != nullptr, "Velocities device buffer should not be null."); @@ -259,12 +259,12 @@ void UpdateConstrainGpu::scaleCoordinates(const matrix scalingMatrix) impl_->scaleCoordinates(scalingMatrix); } -void UpdateConstrainGpu::set(DeviceBuffer d_x, - DeviceBuffer d_v, - const DeviceBuffer d_f, - const t_idef& idef, - const t_mdatoms& md, - const int numTempScaleValues) +void UpdateConstrainGpu::set(DeviceBuffer d_x, + DeviceBuffer d_v, + const DeviceBuffer d_f, + const InteractionDefinitions& idef, + const t_mdatoms& md, + const int numTempScaleValues) { impl_->set(d_x, d_v, d_f, idef, md, numTempScaleValues); } diff --git a/src/gromacs/mdlib/update_constrain_gpu_impl.h b/src/gromacs/mdlib/update_constrain_gpu_impl.h index 0009112dc6..b835c7cf02 100644 --- a/src/gromacs/mdlib/update_constrain_gpu_impl.h +++ b/src/gromacs/mdlib/update_constrain_gpu_impl.h @@ -133,12 +133,12 @@ public: * \param[in] md Atoms data. * \param[in] numTempScaleValues Number of temperature scaling groups. Set zero for no temperature coupling. */ - void set(DeviceBuffer d_x, - DeviceBuffer d_v, - const DeviceBuffer d_f, - const t_idef& idef, - const t_mdatoms& md, - const int numTempScaleValues); + void set(DeviceBuffer d_x, + DeviceBuffer d_v, + const DeviceBuffer d_f, + const InteractionDefinitions& idef, + const t_mdatoms& md, + const int numTempScaleValues); /*! \brief * Update PBC data. diff --git a/src/gromacs/mdlib/updategroups.cpp b/src/gromacs/mdlib/updategroups.cpp index 3aa3feb359..f33aa4cb5c 100644 --- a/src/gromacs/mdlib/updategroups.cpp +++ b/src/gromacs/mdlib/updategroups.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -69,7 +69,7 @@ static bool hasFlexibleConstraints(const gmx_moltype_t& moltype, gmx::ArrayRef ilistsCombined; + ilistsCombined[F_CONSTR] = jointConstraintList(moltype); /* We "include" flexible constraints, but none are present (checked above) */ - const ListOfLists at2con = make_at2con(moltype.atoms.nr, ilistsCombined, iparams.data(), + const ListOfLists at2con = make_at2con(moltype.atoms.nr, ilistsCombined, iparams, FlexibleConstraintTreatment::Include); bool satisfiesCriteria = true; @@ -366,7 +363,7 @@ static RangePartitioning makeUpdateGroups(const gmx_moltype_t& moltype, gmx::Arr int firstAtom = 0; while (satisfiesCriteria && firstAtom < moltype.atoms.nr) { - int numAtomsInGroup = detectGroup(firstAtom, moltype, at2con, constraintsCombined); + int numAtomsInGroup = detectGroup(firstAtom, moltype, at2con, ilistsCombined[F_CONSTR]); if (numAtomsInGroup == 0) { @@ -622,7 +619,7 @@ static real computeMaxUpdateGroupRadius(const gmx_moltype_t& moltype, maxAtom = a; } } - GMX_ASSERT(maxAtom >= 0 || settles.size() > 0, + GMX_ASSERT(maxAtom >= 0 || !settles.empty(), "We should have at least two atoms in the group with constraints"); if (maxAtom < 0) { diff --git a/src/gromacs/mdlib/vsite.cpp b/src/gromacs/mdlib/vsite.cpp index 74c82f6459..1973d2088c 100644 --- a/src/gromacs/mdlib/vsite.cpp +++ b/src/gromacs/mdlib/vsite.cpp @@ -86,18 +86,9 @@ * Any remaining vsites are assigned to a separate master thread task. */ +using gmx::ArrayRef; using gmx::RVec; -static void init_ilist(t_ilist* ilist) -{ - for (int i = 0; i < F_NRE; i++) - { - ilist[i].nr = 0; - ilist[i].nalloc = 0; - ilist[i].iatoms = nullptr; - } -} - /*! \brief List of atom indices belonging to a task */ struct AtomIndex { @@ -109,7 +100,7 @@ struct AtomIndex struct InterdependentTask { //! The interaction lists, only vsite entries are used - t_ilist ilist[F_NRE]; + InteractionLists ilist; //! Thread/task-local force buffer std::vector force; //! The atom indices of the vsites of our task @@ -117,19 +108,13 @@ struct InterdependentTask //! Flags if elements in force are spread to or not std::vector use; //! The number of entries set to true in use - int nuse; + int nuse = 0; //! Array of atoms indices, size nthreads, covering all nuse set elements in use std::vector atomIndex; //! List of tasks (force blocks) this task spread forces to std::vector spreadTask; //! List of tasks that write to this tasks force block range std::vector reduceTask; - - InterdependentTask() - { - init_ilist(ilist); - nuse = 0; - } }; /*! \brief Vsite thread task data structure */ @@ -140,7 +125,7 @@ struct VsiteThread //! End of atom range of this task int rangeEnd; //! The interaction lists, only vsite entries are used - t_ilist ilist[F_NRE]; + std::array ilist; //! Local fshift accumulation buffer rvec fshift[SHIFTS]; //! Local virial dx*df accumulation buffer @@ -155,7 +140,6 @@ struct VsiteThread { rangeStart = -1; rangeEnd = -1; - init_ilist(ilist); clear_rvecs(SHIFTS, fshift); clear_mat(dxdf); useInterdependentTask = false; @@ -166,8 +150,7 @@ struct VsiteThread * * \param[in] ilist The interaction list */ -template -static int vsiteIlistNrCount(const T* ilist) +static int vsiteIlistNrCount(ArrayRef ilist) { int nr = 0; for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++) @@ -417,7 +400,7 @@ static void constr_vsite4FDN(const rvec xi, } -static int constr_vsiten(const t_iatom* ia, const t_iparams ip[], rvec* x, const t_pbc* pbc) +static int constr_vsiten(const t_iatom* ia, ArrayRef ip, rvec* x, const t_pbc* pbc) { rvec x1, dx; dvec dsum; @@ -477,12 +460,12 @@ static PbcMode getPbcMode(const t_pbc* pbcPtr) } } -static void construct_vsites_thread(rvec x[], - real dt, - rvec* v, - const t_iparams ip[], - const t_ilist ilist[], - const t_pbc* pbc_null) +static void construct_vsites_thread(rvec x[], + real dt, + rvec* v, + ArrayRef ip, + ArrayRef ilist, + const t_pbc* pbc_null) { real inv_dt; if (v != nullptr) @@ -500,7 +483,7 @@ static void construct_vsites_thread(rvec x[], for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++) { - if (ilist[ftype].nr == 0) + if (ilist[ftype].empty()) { continue; } @@ -508,9 +491,9 @@ static void construct_vsites_thread(rvec x[], { // TODO remove me int nra = interaction_function[ftype].nratoms; int inc = 1 + nra; - int nr = ilist[ftype].nr; + int nr = ilist[ftype].size(); - const t_iatom* ia = ilist[ftype].iatoms; + const t_iatom* ia = ilist[ftype].iatoms.data(); for (int i = 0; i < nr;) { @@ -609,16 +592,16 @@ static void construct_vsites_thread(rvec x[], } } -void construct_vsites(const gmx_vsite_t* vsite, - rvec x[], - real dt, - rvec* v, - const t_iparams ip[], - const t_ilist ilist[], - PbcType pbcType, - gmx_bool bMolPBC, - const t_commrec* cr, - const matrix box) +void construct_vsites(const gmx_vsite_t* vsite, + rvec x[], + real dt, + rvec* v, + ArrayRef ip, + ArrayRef ilist, + PbcType pbcType, + gmx_bool bMolPBC, + const t_commrec* cr, + const matrix box) { const bool useDomdec = (vsite != nullptr && vsite->useDomdec); GMX_ASSERT(!useDomdec || (cr != nullptr && DOMAINDECOMP(cr)), @@ -749,20 +732,13 @@ void constructVsitesGlobal(const gmx_mtop_t& mtop, gmx::ArrayRef x) { const gmx_molblock_t& molb = mtop.molblock[mb]; const gmx_moltype_t& molt = mtop.moltype[molb.type]; - if (vsiteIlistNrCount(molt.ilist.data()) > 0) + if (vsiteIlistNrCount(molt.ilist) > 0) { 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.data(), ilist, PbcType::No, TRUE, nullptr, nullptr); + mtop.ffparams.iparams, molt.ilist, PbcType::No, TRUE, nullptr, nullptr); atomOffset += molt.atoms.nr; } } @@ -1553,13 +1529,13 @@ static void spread_vsite4FDN(const t_iatom ia[], } -static int spread_vsiten(const t_iatom ia[], - const t_iparams ip[], - const rvec x[], - rvec f[], - rvec fshift[], - const t_pbc* pbc, - const t_graph* g) +static int spread_vsiten(const t_iatom ia[], + ArrayRef ip, + const rvec x[], + rvec f[], + rvec fshift[], + const t_pbc* pbc, + const t_graph* g) { rvec xv, dx, fi; int n3, av, i, ai; @@ -1602,27 +1578,27 @@ static int spread_vsiten(const t_iatom ia[], } -static int vsite_count(const t_ilist* ilist, int ftype) +static int vsite_count(ArrayRef ilist, int ftype) { if (ftype == F_VSITEN) { - return ilist[ftype].nr / 3; + return ilist[ftype].size() / 3; } else { - return ilist[ftype].nr / (1 + interaction_function[ftype].nratoms); + return ilist[ftype].size() / (1 + interaction_function[ftype].nratoms); } } -static void spread_vsite_f_thread(const rvec x[], - rvec f[], - rvec* fshift, - gmx_bool VirCorr, - matrix dxdf, - t_iparams ip[], - const t_ilist ilist[], - const t_graph* g, - const t_pbc* pbc_null) +static void spread_vsite_f_thread(const rvec x[], + rvec f[], + rvec* fshift, + gmx_bool VirCorr, + matrix dxdf, + ArrayRef ip, + ArrayRef ilist, + const t_graph* g, + const t_pbc* pbc_null) { const PbcMode pbcMode = getPbcMode(pbc_null); /* We need another pbc pointer, as with charge groups we switch per vsite */ @@ -1633,7 +1609,7 @@ static void spread_vsite_f_thread(const rvec x[], * higher type vsites from lower types */ for (int ftype = c_ftypeVsiteEnd - 1; ftype >= c_ftypeVsiteStart; ftype--) { - if (ilist[ftype].nr == 0) + if (ilist[ftype].empty()) { continue; } @@ -1641,9 +1617,9 @@ static void spread_vsite_f_thread(const rvec x[], { // TODO remove me int nra = interaction_function[ftype].nratoms; int inc = 1 + nra; - int nr = ilist[ftype].nr; + int nr = ilist[ftype].size(); - const t_iatom* ia = ilist[ftype].iatoms; + const t_iatom* ia = ilist[ftype].iatoms.data(); if (pbcMode == PbcMode::all) { @@ -1724,17 +1700,17 @@ static void clearTaskForceBufferUsedElements(InterdependentTask* idTask) void spread_vsite_f(const gmx_vsite_t* vsite, const rvec* gmx_restrict x, rvec* gmx_restrict f, - rvec* gmx_restrict fshift, - gmx_bool VirCorr, - matrix vir, - t_nrnb* nrnb, - const t_idef* idef, - PbcType pbcType, - gmx_bool bMolPBC, - const t_graph* g, - const matrix box, - const t_commrec* cr, - gmx_wallcycle* wcycle) + rvec* gmx_restrict fshift, + gmx_bool VirCorr, + matrix vir, + t_nrnb* nrnb, + const InteractionDefinitions& idef, + PbcType pbcType, + gmx_bool bMolPBC, + const t_graph* g, + const matrix box, + const t_commrec* cr, + gmx_wallcycle* wcycle) { wallcycle_start(wcycle, ewcVSITESPREAD); const bool useDomdec = vsite->useDomdec; @@ -1768,7 +1744,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite, { clear_mat(dxdf); } - spread_vsite_f_thread(x, f, fshift, VirCorr, dxdf, idef->iparams, idef->il, g, pbc_null); + spread_vsite_f_thread(x, f, fshift, VirCorr, dxdf, idef.iparams, idef.il, g, pbc_null); if (VirCorr) { @@ -1789,7 +1765,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite, clear_mat(vsite->tData[vsite->nthreads]->dxdf); } spread_vsite_f_thread(x, f, fshift, VirCorr, vsite->tData[vsite->nthreads]->dxdf, - idef->iparams, vsite->tData[vsite->nthreads]->ilist, g, pbc_null); + idef.iparams, vsite->tData[vsite->nthreads]->ilist, g, pbc_null); #pragma omp parallel num_threads(vsite->nthreads) { @@ -1836,7 +1812,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite, copy_rvec(f[idTask->vsite[i]], idTask->force[idTask->vsite[i]]); } spread_vsite_f_thread(x, as_rvec_array(idTask->force.data()), fshift_t, VirCorr, - tData.dxdf, idef->iparams, tData.idTask.ilist, g, pbc_null); + tData.dxdf, idef.iparams, tData.idTask.ilist, g, pbc_null); /* We need a barrier before reducing forces below * that have been produced by a different thread above. @@ -1874,7 +1850,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite, } /* Spread the vsites that spread locally only */ - spread_vsite_f_thread(x, f, fshift_t, VirCorr, tData.dxdf, idef->iparams, + spread_vsite_f_thread(x, f, fshift_t, VirCorr, tData.dxdf, idef.iparams, tData.ilist, g, pbc_null); } GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR @@ -1914,15 +1890,15 @@ void spread_vsite_f(const gmx_vsite_t* vsite, dd_move_f_vsites(cr->dd, f, fshift); } - inc_nrnb(nrnb, eNR_VSITE2, vsite_count(idef->il, F_VSITE2)); - inc_nrnb(nrnb, eNR_VSITE2FD, vsite_count(idef->il, F_VSITE2FD)); - inc_nrnb(nrnb, eNR_VSITE3, vsite_count(idef->il, F_VSITE3)); - inc_nrnb(nrnb, eNR_VSITE3FD, vsite_count(idef->il, F_VSITE3FD)); - inc_nrnb(nrnb, eNR_VSITE3FAD, vsite_count(idef->il, F_VSITE3FAD)); - inc_nrnb(nrnb, eNR_VSITE3OUT, vsite_count(idef->il, F_VSITE3OUT)); - inc_nrnb(nrnb, eNR_VSITE4FD, vsite_count(idef->il, F_VSITE4FD)); - inc_nrnb(nrnb, eNR_VSITE4FDN, vsite_count(idef->il, F_VSITE4FDN)); - inc_nrnb(nrnb, eNR_VSITEN, vsite_count(idef->il, F_VSITEN)); + inc_nrnb(nrnb, eNR_VSITE2, vsite_count(idef.il, F_VSITE2)); + inc_nrnb(nrnb, eNR_VSITE2FD, vsite_count(idef.il, F_VSITE2FD)); + inc_nrnb(nrnb, eNR_VSITE3, vsite_count(idef.il, F_VSITE3)); + inc_nrnb(nrnb, eNR_VSITE3FD, vsite_count(idef.il, F_VSITE3FD)); + inc_nrnb(nrnb, eNR_VSITE3FAD, vsite_count(idef.il, F_VSITE3FAD)); + inc_nrnb(nrnb, eNR_VSITE3OUT, vsite_count(idef.il, F_VSITE3OUT)); + inc_nrnb(nrnb, eNR_VSITE4FD, vsite_count(idef.il, F_VSITE4FD)); + inc_nrnb(nrnb, eNR_VSITE4FDN, vsite_count(idef.il, F_VSITE4FDN)); + inc_nrnb(nrnb, eNR_VSITEN, vsite_count(idef.il, F_VSITEN)); wallcycle_stop(wcycle, ewcVSITESPREAD); } @@ -2103,24 +2079,24 @@ static inline void flagAtom(InterdependentTask* idTask, int atom, int thread, in * taskIndex[] is set for all vsites in our range, either to our local tasks * or to the single last task as taskIndex[]=2*nthreads. */ -static void assignVsitesToThread(VsiteThread* tData, - int thread, - int nthread, - int natperthread, - gmx::ArrayRef taskIndex, - const t_ilist* ilist, - const t_iparams* ip, - const unsigned short* ptype) +static void assignVsitesToThread(VsiteThread* tData, + int thread, + int nthread, + int natperthread, + gmx::ArrayRef taskIndex, + ArrayRef ilist, + ArrayRef ip, + const unsigned short* ptype) { for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++) { - tData->ilist[ftype].nr = 0; - tData->idTask.ilist[ftype].nr = 0; + tData->ilist[ftype].clear(); + tData->idTask.ilist[ftype].clear(); - int nral1 = 1 + NRAL(ftype); - int inc = nral1; - t_iatom* iat = ilist[ftype].iatoms; - for (int i = 0; i < ilist[ftype].nr;) + int nral1 = 1 + NRAL(ftype); + int inc = nral1; + const int* iat = ilist[ftype].iatoms.data(); + for (int i = 0; i < ilist[ftype].size();) { if (ftype == F_VSITEN) { @@ -2193,7 +2169,7 @@ static void assignVsitesToThread(VsiteThread* tData, if (task == thread || task == nthread + thread) { /* Copy this vsite to the thread data struct of thread */ - t_ilist* il_task; + InteractionList* il_task; if (task == thread) { il_task = &tData->ilist[ftype]; @@ -2202,17 +2178,8 @@ static void assignVsitesToThread(VsiteThread* tData, { il_task = &tData->idTask.ilist[ftype]; } - /* Ensure we have sufficient memory allocated */ - if (il_task->nr + inc > il_task->nalloc) - { - il_task->nalloc = over_alloc_large(il_task->nr + inc); - srenew(il_task->iatoms, il_task->nalloc); - } /* Copy the vsite data to the thread-task local array */ - for (int j = i; j < i + inc; j++) - { - il_task->iatoms[il_task->nr++] = iat[j]; - } + il_task->push_back(iat[i], nral1 - 1, iat + i + 1); if (task == nthread + thread) { /* This vsite write outside our own task force block. @@ -2243,23 +2210,23 @@ static void assignVsitesToThread(VsiteThread* tData, } /*! \brief Assign all vsites with taskIndex[]==task to task tData */ -static void assignVsitesToSingleTask(VsiteThread* tData, - int task, - gmx::ArrayRef taskIndex, - const t_ilist* ilist, - const t_iparams* ip) +static void assignVsitesToSingleTask(VsiteThread* tData, + int task, + gmx::ArrayRef taskIndex, + ArrayRef ilist, + ArrayRef ip) { for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++) { - tData->ilist[ftype].nr = 0; - tData->idTask.ilist[ftype].nr = 0; + tData->ilist[ftype].clear(); + tData->idTask.ilist[ftype].clear(); - int nral1 = 1 + NRAL(ftype); - int inc = nral1; - t_iatom* iat = ilist[ftype].iatoms; - t_ilist* il_task = &tData->ilist[ftype]; + int nral1 = 1 + NRAL(ftype); + int inc = nral1; + const int* iat = ilist[ftype].iatoms.data(); + InteractionList* il_task = &tData->ilist[ftype]; - for (int i = 0; i < ilist[ftype].nr;) + for (int i = 0; i < ilist[ftype].size();) { if (ftype == F_VSITEN) { @@ -2269,17 +2236,8 @@ static void assignVsitesToSingleTask(VsiteThread* tData, /* Check if the vsite is assigned to our task */ if (taskIndex[iat[1 + i]] == task) { - /* Ensure we have sufficient memory allocated */ - if (il_task->nr + inc > il_task->nalloc) - { - il_task->nalloc = over_alloc_large(il_task->nr + inc); - srenew(il_task->iatoms, il_task->nalloc); - } /* Copy the vsite data to the thread-task local array */ - for (int j = i; j < i + inc; j++) - { - il_task->iatoms[il_task->nr++] = iat[j]; - } + il_task->push_back(iat[i], inc - 1, iat + i + 1); } i += inc; @@ -2287,7 +2245,10 @@ static void assignVsitesToSingleTask(VsiteThread* tData, } } -void split_vsites_over_threads(const t_ilist* ilist, const t_iparams* ip, const t_mdatoms* mdatoms, gmx_vsite_t* vsite) +void split_vsites_over_threads(ArrayRef ilist, + ArrayRef ip, + const t_mdatoms* mdatoms, + gmx_vsite_t* vsite) { int vsite_atom_range, natperthread; @@ -2315,9 +2276,9 @@ void split_vsites_over_threads(const t_ilist* ilist, const t_iparams* ip, const { // TODO remove me if (ftype != F_VSITEN) { - int nral1 = 1 + NRAL(ftype); - const t_iatom* iat = ilist[ftype].iatoms; - for (int i = 0; i < ilist[ftype].nr; i += nral1) + int nral1 = 1 + NRAL(ftype); + ArrayRef iat = ilist[ftype].iatoms; + for (int i = 0; i < ilist[ftype].size(); i += nral1) { for (int j = i + 1; j < i + nral1; j++) { @@ -2329,10 +2290,10 @@ void split_vsites_over_threads(const t_ilist* ilist, const t_iparams* ip, const { int vs_ind_end; - const t_iatom* iat = ilist[ftype].iatoms; + ArrayRef iat = ilist[ftype].iatoms; int i = 0; - while (i < ilist[ftype].nr) + while (i < ilist[ftype].size()) { /* The 3 below is from 1+NRAL(ftype)=3 */ vs_ind_end = i + ip[iat[i]].vsiten.n * 3; @@ -2518,13 +2479,13 @@ void split_vsites_over_threads(const t_ilist* ilist, const t_iparams* ip, const for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++) { - if (ilist[ftype].nr > 0) + if (!ilist[ftype].empty()) { fprintf(debug, "%-20s thread dist:", interaction_function[ftype].longname); for (int th = 0; th < vsite->nthreads + 1; th++) { - fprintf(debug, " %4d %4d ", vsite->tData[th]->ilist[ftype].nr, - vsite->tData[th]->idTask.ilist[ftype].nr); + fprintf(debug, " %4d %4d ", vsite->tData[th]->ilist[ftype].size(), + vsite->tData[th]->idTask.ilist[ftype].size()); } fprintf(debug, "\n"); } diff --git a/src/gromacs/mdlib/vsite.h b/src/gromacs/mdlib/vsite.h index 65287fcb14..291fbc9525 100644 --- a/src/gromacs/mdlib/vsite.h +++ b/src/gromacs/mdlib/vsite.h @@ -51,7 +51,7 @@ struct gmx_localtop_t; struct gmx_mtop_t; struct t_commrec; struct t_graph; -struct t_ilist; +struct InteractionList; struct t_mdatoms; struct t_nrnb; struct gmx_wallcycle; @@ -103,16 +103,16 @@ struct gmx_vsite_t * \param[in] cr The communication record * \param[in] box The box */ -void construct_vsites(const gmx_vsite_t* vsite, - rvec x[], - real dt, - rvec v[], - const t_iparams ip[], - const t_ilist ilist[], - PbcType pbcType, - gmx_bool bMolPBC, - const t_commrec* cr, - const matrix box); +void construct_vsites(const gmx_vsite_t* vsite, + rvec x[], + real dt, + rvec v[], + gmx::ArrayRef ip, + gmx::ArrayRef ilist, + PbcType pbcType, + gmx_bool bMolPBC, + const t_commrec* cr, + const matrix box); /*! \brief Create positions of vsite atoms for the whole system assuming all molecules are wholex * @@ -121,20 +121,20 @@ void construct_vsites(const gmx_vsite_t* vsite, */ void constructVsitesGlobal(const gmx_mtop_t& mtop, gmx::ArrayRef x); -void spread_vsite_f(const gmx_vsite_t* vsite, - const rvec x[], - rvec f[], - rvec* fshift, - gmx_bool VirCorr, - matrix vir, - t_nrnb* nrnb, - const t_idef* idef, - PbcType pbcType, - gmx_bool bMolPBC, - const t_graph* g, - const matrix box, - const t_commrec* cr, - gmx_wallcycle* wcycle); +void spread_vsite_f(const gmx_vsite_t* vsite, + const rvec x[], + rvec f[], + rvec* fshift, + gmx_bool VirCorr, + matrix vir, + t_nrnb* nrnb, + const InteractionDefinitions& idef, + PbcType pbcType, + gmx_bool bMolPBC, + const t_graph* g, + const matrix box, + const t_commrec* cr, + gmx_wallcycle* wcycle); /* Spread the force operating on the vsite atoms on the surrounding atoms. * If fshift!=NULL also update the shift forces. * If VirCorr=TRUE add the virial correction for non-linear vsite constructs @@ -162,10 +162,10 @@ int countInterUpdategroupVsites(const gmx_mtop_t& mtop */ std::unique_ptr initVsite(const gmx_mtop_t& mtop, const t_commrec* cr); -void split_vsites_over_threads(const t_ilist* ilist, - const t_iparams* ip, - const t_mdatoms* mdatoms, - gmx_vsite_t* vsite); +void split_vsites_over_threads(gmx::ArrayRef ilist, + gmx::ArrayRef ip, + const t_mdatoms* mdatoms, + gmx_vsite_t* vsite); /* Divide the vsite work-load over the threads. * Should be called at the end of the domain decomposition. */ diff --git a/src/gromacs/mdlib/wall.h b/src/gromacs/mdlib/wall.h index 95437da770..07ad05390c 100644 --- a/src/gromacs/mdlib/wall.h +++ b/src/gromacs/mdlib/wall.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -41,7 +41,6 @@ struct SimulationGroups; struct t_forcerec; -struct t_idef; struct t_inputrec; struct t_mdatoms; struct t_nrnb; diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 9e7943640d..2cb1388dca 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -175,7 +175,6 @@ void gmx::LegacySimulator::do_md() rvec mu_tot; matrix pressureCouplingMu, M; gmx_repl_ex_t repl_ex = nullptr; - gmx_localtop_t top; PaddedHostVector f{}; gmx_global_stat_t gstat; t_graph* graph = nullptr; @@ -316,14 +315,14 @@ void gmx::LegacySimulator::do_md() t_state* state; + gmx_localtop_t top(top_global->ffparams); + auto mdatoms = mdAtoms->mdatoms(); std::unique_ptr integrator; if (DOMAINDECOMP(cr)) { - dd_init_local_top(*top_global, &top); - stateInstance = std::make_unique(); state = stateInstance.get(); dd_init_local_state(cr->dd, state_global, state); diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index 381b3b27a6..6984258b8e 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -149,7 +149,6 @@ void gmx::LegacySimulator::do_mimic() unsigned int force_flags; tensor force_vir, shake_vir, total_vir, pres; rvec mu_tot; - gmx_localtop_t top; PaddedHostVector f{}; gmx_global_stat_t gstat; t_graph* graph = nullptr; @@ -252,10 +251,10 @@ void gmx::LegacySimulator::do_mimic() std::unique_ptr stateInstance; t_state* state; + gmx_localtop_t top(top_global->ffparams); + if (DOMAINDECOMP(cr)) { - dd_init_local_top(*top_global, &top); - stateInstance = std::make_unique(); state = stateInstance.get(); dd_init_local_state(cr->dd, state_global, state); diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index 3fc65f392f..bcda82941b 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -416,9 +416,6 @@ static void init_em(FILE* fplog, auto mdatoms = mdAtoms->mdatoms(); if (DOMAINDECOMP(cr)) { - top->useInDomainDecomp_ = true; - dd_init_local_top(*top_global, top); - dd_init_local_state(cr->dd, state_global, &ems->s); /* Distribute the charge groups over the nodes from the master node */ @@ -1042,7 +1039,7 @@ void LegacySimulator::do_cg() { const char* CG = "Polak-Ribiere Conjugate Gradients"; - gmx_localtop_t top; + gmx_localtop_t top(top_global->ffparams); gmx_global_stat_t gstat; t_graph* graph; double tmp, minstep; @@ -1641,7 +1638,7 @@ void LegacySimulator::do_lbfgs() { static const char* LBFGS = "Low-Memory BFGS Minimizer"; em_state_t ems; - gmx_localtop_t top; + gmx_localtop_t top(top_global->ffparams); gmx_global_stat_t gstat; t_graph* graph; int ncorr, nmaxcorr, point, cp, neval, nminstep; @@ -2359,7 +2356,7 @@ void LegacySimulator::do_lbfgs() void LegacySimulator::do_steep() { const char* SD = "Steepest Descents"; - gmx_localtop_t top; + gmx_localtop_t top(top_global->ffparams); gmx_global_stat_t gstat; t_graph* graph; real stepsize; @@ -2594,7 +2591,7 @@ void LegacySimulator::do_nm() { const char* NM = "Normal Mode Analysis"; int nnodes; - gmx_localtop_t top; + gmx_localtop_t top(top_global->ffparams); gmx_global_stat_t gstat; t_graph* graph; tensor vir, pres; diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp index 29c1ae1c6b..41d8c4f12d 100644 --- a/src/gromacs/mdrun/rerun.cpp +++ b/src/gromacs/mdrun/rerun.cpp @@ -148,14 +148,14 @@ using gmx::SimulationSignaller; * \param[in] forceRec Force record, used for constructing vsites * \param[in,out] graph The molecular graph, used for constructing vsites when != nullptr */ -static void prepareRerunState(const t_trxframe& rerunFrame, - t_state* globalState, - bool constructVsites, - const gmx_vsite_t* vsite, - const t_idef& idef, - double timeStep, - const t_forcerec& forceRec, - t_graph* graph) +static void prepareRerunState(const t_trxframe& rerunFrame, + t_state* globalState, + bool constructVsites, + const gmx_vsite_t* vsite, + const InteractionDefinitions& idef, + double timeStep, + const t_forcerec& forceRec, + t_graph* graph) { auto x = makeArrayRef(globalState->x); auto rerunX = arrayRefFromArray(reinterpret_cast(rerunFrame.x), globalState->natoms); @@ -202,7 +202,7 @@ void gmx::LegacySimulator::do_rerun() t_trxstatus* status = nullptr; rvec mu_tot; t_trxframe rerun_fr; - gmx_localtop_t top; + gmx_localtop_t top(top_global->ffparams); PaddedHostVector f{}; gmx_global_stat_t gstat; t_graph* graph = nullptr; @@ -328,8 +328,6 @@ void gmx::LegacySimulator::do_rerun() if (DOMAINDECOMP(cr)) { - dd_init_local_top(*top_global, &top); - stateInstance = std::make_unique(); state = stateInstance.get(); dd_init_local_state(cr->dd, state_global, state); diff --git a/src/gromacs/mdrun/shellfc.cpp b/src/gromacs/mdrun/shellfc.cpp index 68ea3d0fe3..d93e12a8cf 100644 --- a/src/gromacs/mdrun/shellfc.cpp +++ b/src/gromacs/mdrun/shellfc.cpp @@ -943,7 +943,7 @@ void relax_shell_flexcon(FILE* fplog, ArrayRef shells = shfc->shells; const int nflexcon = shfc->nflexcon; - const t_idef* idef = &top->idef; + const InteractionDefinitions& idef = top->idef; if (DOMAINDECOMP(cr)) { @@ -1102,7 +1102,7 @@ void relax_shell_flexcon(FILE* fplog, if (vsite) { construct_vsites(vsite, as_rvec_array(pos[Min].data()), inputrec->delta_t, - as_rvec_array(v.data()), idef->iparams, idef->il, fr->pbcType, + as_rvec_array(v.data()), idef.iparams, idef.il, fr->pbcType, fr->bMolPBC, cr, box); } diff --git a/src/gromacs/mdrun/tpi.cpp b/src/gromacs/mdrun/tpi.cpp index 8e2c547d5b..c63a013257 100644 --- a/src/gromacs/mdrun/tpi.cpp +++ b/src/gromacs/mdrun/tpi.cpp @@ -165,7 +165,7 @@ void LegacySimulator::do_tpi() { GMX_RELEASE_ASSERT(gmx_omp_nthreads_get(emntDefault) == 1, "TPI does not support OpenMP"); - gmx_localtop_t top; + gmx_localtop_t top(top_global->ffparams); PaddedHostVector f{}; real lambda, t, temp, beta, drmax, epot; double embU, sum_embU, *sum_UgembU, V, V_all, VembU_all; diff --git a/src/gromacs/modularsimulator/topologyholder.cpp b/src/gromacs/modularsimulator/topologyholder.cpp index e5bf8f5b87..496b24dd60 100644 --- a/src/gromacs/modularsimulator/topologyholder.cpp +++ b/src/gromacs/modularsimulator/topologyholder.cpp @@ -58,13 +58,9 @@ TopologyHolder::TopologyHolder(const gmx_mtop_t& globalTopology, Constraints* constr, gmx_vsite_t* vsite) : globalTopology_(globalTopology), - localTopology_(std::make_unique()) + localTopology_(std::make_unique(globalTopology.ffparams)) { - if (DOMAINDECOMP(cr)) - { - dd_init_local_top(globalTopology, localTopology_.get()); - } - else + if (!DOMAINDECOMP(cr)) { t_graph* graph = nullptr; // Generate and initialize new topology diff --git a/src/gromacs/pbcutil/mshift.cpp b/src/gromacs/pbcutil/mshift.cpp index 0bd4f8e73b..6d860a0f2e 100644 --- a/src/gromacs/pbcutil/mshift.cpp +++ b/src/gromacs/pbcutil/mshift.cpp @@ -496,6 +496,17 @@ void mk_graph_moltype(const gmx_moltype_t& moltype, t_graph* g) mk_graph_ilist(nullptr, moltype.ilist.data(), 0, moltype.atoms.nr, FALSE, FALSE, g); } +t_graph* mk_graph(FILE* fplog, const InteractionDefinitions& idef, int at_start, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle) +{ + t_graph* g; + + snew(g, 1); + + mk_graph_ilist(fplog, idef.il.data(), at_start, at_end, bShakeOnly, bSettle, g); + + return g; +} + t_graph* mk_graph(FILE* fplog, const t_idef* idef, int at_start, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle) { t_graph* g; diff --git a/src/gromacs/pbcutil/mshift.h b/src/gromacs/pbcutil/mshift.h index db37031723..a1ef6b1431 100644 --- a/src/gromacs/pbcutil/mshift.h +++ b/src/gromacs/pbcutil/mshift.h @@ -45,6 +45,7 @@ struct InteractionList; struct gmx_moltype_t; +class InteractionDefinitions; struct t_idef; enum class PbcType : int; @@ -85,6 +86,20 @@ struct t_graph #define SHIFT_IVEC(g, i) ((g)->ishift[i]) +t_graph* mk_graph(FILE* fplog, + const InteractionDefinitions& idef, + int at_start, + int at_end, + gmx_bool bShakeOnly, + gmx_bool bSettle); +/* Build a graph from an idef description. The graph can be used + * to generate mol-shift indices. + * at_start and at_end should coincide will molecule boundaries, + * for the whole system this is simply 0 and natoms. + * If bShakeOnly, only the connections in the shake list are used. + * If bSettle && bShakeOnly the settles are used too. + */ + t_graph* mk_graph(FILE* fplog, const struct t_idef* idef, int at_start, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle); /* Build a graph from an idef description. The graph can be used * to generate mol-shift indices. diff --git a/src/gromacs/pbcutil/rmpbc.cpp b/src/gromacs/pbcutil/rmpbc.cpp index 92792d4df2..537b8a2f48 100644 --- a/src/gromacs/pbcutil/rmpbc.cpp +++ b/src/gromacs/pbcutil/rmpbc.cpp @@ -62,11 +62,13 @@ typedef struct struct gmx_rmpbc { - const t_idef* idef; - int natoms_init; - PbcType pbcType; - int ngraph; - rmpbc_graph_t* graph; + const InteractionDefinitions* interactionDefinitions; + const t_idef* idef; + int natoms_init; + PbcType pbcType; + int ePBC; + int ngraph; + rmpbc_graph_t* graph; }; static t_graph* gmx_rmpbc_get_graph(gmx_rmpbc_t gpbc, PbcType pbcType, int natoms) @@ -104,12 +106,37 @@ static t_graph* gmx_rmpbc_get_graph(gmx_rmpbc_t gpbc, PbcType pbcType, int natom srenew(gpbc->graph, gpbc->ngraph); gr = &gpbc->graph[gpbc->ngraph - 1]; gr->natoms = natoms; - gr->gr = mk_graph(nullptr, gpbc->idef, 0, natoms, FALSE, FALSE); + if (gpbc->interactionDefinitions) + { + gr->gr = mk_graph(nullptr, *gpbc->interactionDefinitions, 0, natoms, FALSE, FALSE); + } + else + { + gr->gr = mk_graph(nullptr, gpbc->idef, 0, natoms, FALSE, FALSE); + } } return gr->gr; } +gmx_rmpbc_t gmx_rmpbc_init(const InteractionDefinitions& idef, PbcType pbcType, int natoms) +{ + gmx_rmpbc_t gpbc; + + snew(gpbc, 1); + + gpbc->natoms_init = natoms; + + /* This sets pbc when we now it, + * otherwise we guess it from the instantaneous box in the trajectory. + */ + gpbc->pbcType = pbcType; + + gpbc->interactionDefinitions = &idef; + + return gpbc; +} + gmx_rmpbc_t gmx_rmpbc_init(const t_idef* idef, PbcType pbcType, int natoms) { gmx_rmpbc_t gpbc; diff --git a/src/gromacs/pbcutil/rmpbc.h b/src/gromacs/pbcutil/rmpbc.h index 1b93ae538e..ac7411ee9c 100644 --- a/src/gromacs/pbcutil/rmpbc.h +++ b/src/gromacs/pbcutil/rmpbc.h @@ -39,6 +39,7 @@ #include "gromacs/math/vectypes.h" +class InteractionDefinitions; struct t_atoms; struct t_idef; struct t_trxframe; @@ -46,6 +47,8 @@ enum class PbcType : int; typedef struct gmx_rmpbc* gmx_rmpbc_t; +gmx_rmpbc_t gmx_rmpbc_init(const InteractionDefinitions& idef, PbcType pbcType, int natoms); + gmx_rmpbc_t gmx_rmpbc_init(const t_idef* idef, PbcType pbcType, int natoms); void gmx_rmpbc_done(gmx_rmpbc_t gpbc); diff --git a/src/gromacs/tools/check.cpp b/src/gromacs/tools/check.cpp index 03e4ffd85e..684c6ca6a2 100644 --- a/src/gromacs/tools/check.cpp +++ b/src/gromacs/tools/check.cpp @@ -233,19 +233,21 @@ static void chk_forces(int frame, int natoms, rvec* f) } } -static void chk_bonds(t_idef* idef, PbcType pbcType, rvec* x, matrix box, real tol) +static void chk_bonds(const InteractionDefinitions* idef, PbcType pbcType, rvec* x, matrix box, real tol) { int ftype, k, ai, aj, type; real b0, blen, deviation; t_pbc pbc; rvec dx; + gmx::ArrayRef iparams = idef->iparams; + set_pbc(&pbc, pbcType, box); for (ftype = 0; (ftype < F_NRE); ftype++) { if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND) { - for (k = 0; (k < idef->il[ftype].nr);) + for (k = 0; (k < idef->il[ftype].size());) { type = idef->il[ftype].iatoms[k++]; ai = idef->il[ftype].iatoms[k++]; @@ -253,11 +255,11 @@ static void chk_bonds(t_idef* idef, PbcType pbcType, rvec* x, matrix box, real t b0 = 0; switch (ftype) { - case F_BONDS: b0 = idef->iparams[type].harmonic.rA; break; - case F_G96BONDS: b0 = std::sqrt(idef->iparams[type].harmonic.rA); break; - case F_MORSE: b0 = idef->iparams[type].morse.b0A; break; - case F_CUBICBONDS: b0 = idef->iparams[type].cubic.b0; break; - case F_CONSTR: b0 = idef->iparams[type].constr.dA; break; + case F_BONDS: b0 = iparams[type].harmonic.rA; break; + case F_G96BONDS: b0 = std::sqrt(iparams[type].harmonic.rA); break; + case F_MORSE: b0 = iparams[type].morse.b0A; break; + case F_CUBICBONDS: b0 = iparams[type].cubic.b0; break; + case F_CONSTR: b0 = iparams[type].constr.dA; break; default: break; } if (b0 != 0) @@ -278,22 +280,23 @@ static void chk_bonds(t_idef* idef, PbcType pbcType, rvec* x, matrix box, real t static void chk_trj(const gmx_output_env_t* oenv, const char* fn, const char* tpr, real tol) { - t_trxframe fr; - t_count count; - t_fr_time first, last; - int j = -1, new_natoms, natoms; - real old_t1, old_t2; - gmx_bool bShowTimestep = TRUE, newline = FALSE; - t_trxstatus* status; - gmx_mtop_t mtop; - gmx_localtop_t top; - t_state state; - t_inputrec ir; - + t_trxframe fr; + t_count count; + t_fr_time first, last; + int j = -1, new_natoms, natoms; + real old_t1, old_t2; + gmx_bool bShowTimestep = TRUE, newline = FALSE; + t_trxstatus* status; + gmx_mtop_t mtop; + t_state state; + t_inputrec ir; + + std::unique_ptr top; if (tpr) { read_tpx_state(tpr, &ir, &state, &mtop); - gmx_mtop_generate_local_top(mtop, &top, ir.efep != efepNO); + top = std::make_unique(mtop.ffparams); + gmx_mtop_generate_local_top(mtop, top.get(), ir.efep != efepNO); } new_natoms = -1; natoms = -1; @@ -359,7 +362,7 @@ static void chk_trj(const gmx_output_env_t* oenv, const char* fn, const char* tp natoms = new_natoms; if (tpr) { - chk_bonds(&top.idef, ir.pbcType, fr.x, fr.box, tol); + chk_bonds(&top->idef, ir.pbcType, fr.x, fr.box, tol); } if (fr.bX) { diff --git a/src/gromacs/tools/convert_tpr.cpp b/src/gromacs/tools/convert_tpr.cpp index bf0472bee4..8e6bfb1fe2 100644 --- a/src/gromacs/tools/convert_tpr.cpp +++ b/src/gromacs/tools/convert_tpr.cpp @@ -184,50 +184,40 @@ static void reduce_atom(int gnx, const int index[], t_atom atom[], char*** atomn static void reduce_ilist(gmx::ArrayRef invindex, const std::vector& bKeep, - t_ilist* il, + InteractionList* il, int nratoms, const char* name) { - t_iatom* ia; - int i, j, newnr; - gmx_bool bB; - - if (il->nr) + if (!il->empty()) { - snew(ia, il->nr); - newnr = 0; - for (i = 0; (i < il->nr); i += nratoms + 1) + std::vector newAtoms(nratoms); + InteractionList ilReduced; + for (int i = 0; i < il->size(); i += nratoms + 1) { - bB = TRUE; - for (j = 1; (j <= nratoms); j++) + bool bB = true; + for (int j = 0; j < nratoms; j++) { - bB = bB && bKeep[il->iatoms[i + j]]; + bB = bB && bKeep[il->iatoms[i + 1 + j]]; } if (bB) { - ia[newnr++] = il->iatoms[i]; - for (j = 1; (j <= nratoms); j++) + for (int j = 0; j < nratoms; j++) { - ia[newnr++] = invindex[il->iatoms[i + j]]; + newAtoms[j] = invindex[il->iatoms[i + 1 + j]]; } + ilReduced.push_back(il->iatoms[i], nratoms, newAtoms.data()); } } - fprintf(stderr, "Reduced ilist %8s from %6d to %6d entries\n", name, il->nr / (nratoms + 1), - newnr / (nratoms + 1)); - - il->nr = newnr; - for (i = 0; (i < newnr); i++) - { - il->iatoms[i] = ia[i]; - } + fprintf(stderr, "Reduced ilist %8s from %6d to %6d entries\n", name, + il->size() / (nratoms + 1), ilReduced.size() / (nratoms + 1)); - sfree(ia); + *il = std::move(ilReduced); } } static void reduce_topology_x(int gnx, int index[], gmx_mtop_t* mtop, rvec x[], rvec v[]) { - gmx_localtop_t top; + gmx_localtop_t top(mtop->ffparams); gmx_mtop_generate_local_top(*mtop, &top, false); t_atoms atoms = gmx_mtop_global_atoms(mtop); @@ -252,12 +242,7 @@ static void reduce_topology_x(int gnx, int index[], gmx_mtop_t* mtop, rvec x[], mtop->moltype[0].excls = reduce_listoflists(invindex, bKeep, top.excls, "excls"); for (int i = 0; i < F_NRE; 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].ilist[i] = std::move(top.idef.il[i]); } mtop->molblock.resize(1); diff --git a/src/gromacs/topology/idef.cpp b/src/gromacs/topology/idef.cpp index 2318f4d5b9..a82e6e5754 100644 --- a/src/gromacs/topology/idef.cpp +++ b/src/gromacs/topology/idef.cpp @@ -41,6 +41,7 @@ #include +#include "gromacs/topology/forcefieldparameters.h" #include "gromacs/topology/ifunc.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" @@ -331,7 +332,7 @@ static void printIlist(FILE* fp, indent = pr_title(fp, indent, title); pr_indent(fp, indent); fprintf(fp, "nr: %d\n", ilist.size()); - if (ilist.size() > 0) + if (!ilist.empty()) { pr_indent(fp, indent); fprintf(fp, "iatoms:\n"); @@ -412,15 +413,28 @@ void init_idef(t_idef* idef) idef->iparams_fbposres = nullptr; for (int f = 0; f < F_NRE; ++f) { - idef->il[f].iatoms = nullptr; - idef->il[f].nalloc = 0; - idef->il[f].nr = 0; - idef->numNonperturbedInteractions[f] = 0; + idef->il[f].iatoms = nullptr; + idef->il[f].nalloc = 0; + idef->il[f].nr = 0; } - idef->cmap_grid = nullptr; - idef->iparams_posres_nalloc = 0; - idef->iparams_fbposres_nalloc = 0; - idef->ilsort = ilsortUNKNOWN; +} + +InteractionDefinitions::InteractionDefinitions(const gmx_ffparams_t& ffparams) : + iparams(ffparams.iparams), + functype(ffparams.functype), + cmap_grid(ffparams.cmap_grid) +{ +} + +void InteractionDefinitions::clear() +{ + /* Clear the counts */ + for (auto& ilist : il) + { + ilist.clear(); + } + iparams_posres.clear(); + iparams_fbposres.clear(); } void done_idef(t_idef* idef) @@ -434,18 +448,5 @@ void done_idef(t_idef* idef) sfree(idef->il[f].iatoms); } - delete idef->cmap_grid; init_idef(idef); } - -void copy_ilist(const t_ilist* src, t_ilist* dst) -{ - dst->nr = src->nr; - dst->nalloc = src->nr; - - snew(dst->iatoms, dst->nalloc); - for (int i = 0; i < dst->nr; ++i) - { - dst->iatoms[i] = src->iatoms[i]; - } -} diff --git a/src/gromacs/topology/idef.h b/src/gromacs/topology/idef.h index e19e976c64..40a4913cfb 100644 --- a/src/gromacs/topology/idef.h +++ b/src/gromacs/topology/idef.h @@ -48,6 +48,8 @@ #include "gromacs/utility/basedefinitions.h" #include "gromacs/utility/real.h" +struct gmx_ffparams_t; + typedef union t_iparams { /* Some parameters have A and B values for free energy calculations. * The B values are not used for regular simulations of course. @@ -222,6 +224,43 @@ struct InteractionList /* Returns the total number of elements in iatoms */ int size() const { return static_cast(iatoms.size()); } + /* Returns whether the list is empty */ + bool empty() const { return iatoms.empty(); } + + /* Adds one interaction to the list */ + template + void push_back(const int parameterType, const std::array& atoms) + { + const std::size_t oldSize = iatoms.size(); + iatoms.resize(iatoms.size() + 1 + numAtoms); + iatoms[oldSize] = parameterType; + for (std::size_t i = 0; i < numAtoms; i++) + { + iatoms[oldSize + 1 + i] = atoms[i]; + } + } + + /* Adds one interaction to the list */ + void push_back(const int parameterType, const int numAtoms, const int* atoms) + { + const std::size_t oldSize = iatoms.size(); + iatoms.resize(iatoms.size() + 1 + numAtoms); + iatoms[oldSize] = parameterType; + for (int i = 0; i < numAtoms; i++) + { + iatoms[oldSize + 1 + i] = atoms[i]; + } + } + + /* Appends \p ilist at the back of the list */ + void append(const InteractionList& ilist) + { + iatoms.insert(iatoms.end(), ilist.iatoms.begin(), ilist.iatoms.end()); + } + + /* Clears the list */ + void clear() { iatoms.clear(); } + /* List of interactions, see explanation further down */ std::vector iatoms; }; @@ -230,31 +269,23 @@ struct InteractionList * * TODO: Consider only including entries in use instead of all F_NRE */ -typedef std::array InteractionLists; +using InteractionLists = std::array; -/* 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. - */ +/* Deprecated list of listed interactions */ struct t_ilist { /* Returns the total number of elements in iatoms */ int size() const { return nr; } + /* Returns whether the list is empty */ + bool empty() const { return nr == 0; } + int nr; t_iatom* iatoms; int nalloc; }; -/* 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. - */ +/* TODO: Remove t_ilist and remove templating on list type in mshift.cpp */ /* * The structs InteractionList and t_ilist defines a list of atoms with their interactions. @@ -295,7 +326,7 @@ static inline std::vector extractILists(const Interaction std::vector handles; for (size_t ftype = 0; ftype < ilists.size(); ftype++) { - if ((interaction_function[ftype].flags & flags) && ilists[ftype].size() > 0) + if ((interaction_function[ftype].flags & flags) && !ilists[ftype].empty()) { handles.push_back({ static_cast(ftype), ilists[ftype].iatoms }); } @@ -333,22 +364,54 @@ enum ilsortFE_SORTED }; -typedef struct t_idef +/* Struct with list of interaction parameters and lists of interactions + * + * TODO: Convert to a proper class with private data members so we can + * ensure that the free-energy sorting and sorting setting is consistent. + */ +class InteractionDefinitions +{ +public: + /* Constructor + * + * \param[in] ffparams The interaction parameters, the lifetime of the created object should not exceed the lifetime of the passed parameters + */ + InteractionDefinitions(const gmx_ffparams_t& ffparams); + + // Clears data not read in from ffparams + void clear(); + + // The interaction parameters + const std::vector& iparams; + // The function type per type + const std::vector& functype; + // Position restraint interaction parameters + std::vector iparams_posres; + // Flat-bottomed position restraint parameters + std::vector iparams_fbposres; + // The list of interactions for each type. Note that some, such as LJ and COUL will have 0 entries. + std::array il; + /* The number of non-perturbed interactions at the start of each entry in il */ + std::array numNonperturbedInteractions; + // The sorting state of interaction in il + int ilsort = ilsortUNKNOWN; + // The dihedral correction maps + gmx_cmap_t cmap_grid; +}; + +/* Deprecated interation definitions, used in t_topology */ +struct t_idef { int ntypes; int atnr; t_functype* functype; t_iparams* iparams; real fudgeQQ; - gmx_cmap_t* cmap_grid; t_iparams * iparams_posres, *iparams_fbposres; - int iparams_posres_nalloc, iparams_fbposres_nalloc; t_ilist il[F_NRE]; - /* The number of non-perturbed interactions at the start of each entry in il */ - int numNonperturbedInteractions[F_NRE]; - int ilsort; -} t_idef; + int ilsort; +}; /* * The struct t_idef defines all the interactions for the complete @@ -415,6 +478,4 @@ void init_idef(t_idef* idef); */ void done_idef(t_idef* idef); -void copy_ilist(const t_ilist* src, t_ilist* dst); - #endif diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp index e1f7ce822d..13c79f817c 100644 --- a/src/gromacs/topology/mtop_util.cpp +++ b/src/gromacs/topology/mtop_util.cpp @@ -667,6 +667,29 @@ t_atoms gmx_mtop_global_atoms(const gmx_mtop_t* mtop) * The cat routines below are old code from src/kernel/topcat.c */ +static void ilistcat(int ftype, InteractionList* dest, const InteractionList& src, int copies, int dnum, int snum) +{ + int nral, c, i, a; + + nral = NRAL(ftype); + + size_t destIndex = dest->iatoms.size(); + dest->iatoms.resize(dest->iatoms.size() + copies * src.size()); + + for (c = 0; c < copies; c++) + { + for (i = 0; i < src.size();) + { + dest->iatoms[destIndex++] = src.iatoms[i++]; + for (a = 0; a < nral; a++) + { + dest->iatoms[destIndex++] = dnum + src.iatoms[i++]; + } + } + dnum += snum; + } +} + static void ilistcat(int ftype, t_ilist* dest, const InteractionList& src, int copies, int dnum, int snum) { int nral, c, i, a; @@ -690,21 +713,40 @@ static void ilistcat(int ftype, t_ilist* dest, const InteractionList& src, int c } } -static void set_posres_params(t_idef* idef, const gmx_molblock_t* molb, int i0, int a_offset) +static const t_iparams& getIparams(const InteractionDefinitions& idef, const int index) +{ + return idef.iparams[index]; +} + +static const t_iparams& getIparams(const t_idef& idef, const int index) +{ + return idef.iparams[index]; +} + +static void resizeIParams(std::vector* iparams, const int newSize) +{ + iparams->resize(newSize); +} + +static void resizeIParams(t_iparams** iparams, const int newSize) +{ + srenew(*iparams, newSize); +} + +template +static void set_posres_params(IdefType* idef, const gmx_molblock_t* molb, int i0, int a_offset) { - t_ilist* il; int i1, i, a_molb; t_iparams* ip; - il = &idef->il[F_POSRES]; - i1 = il->nr / 2; - idef->iparams_posres_nalloc = i1; - srenew(idef->iparams_posres, idef->iparams_posres_nalloc); + auto* il = &idef->il[F_POSRES]; + i1 = il->size() / 2; + resizeIParams(&idef->iparams_posres, i1); for (i = i0; i < i1; i++) { ip = &idef->iparams_posres[i]; /* Copy the force constants */ - *ip = idef->iparams[il->iatoms[i * 2]]; + *ip = getIparams(*idef, il->iatoms[i * 2]); a_molb = il->iatoms[i * 2 + 1] - a_offset; if (molb->posres_xA.empty()) { @@ -730,21 +772,20 @@ static void set_posres_params(t_idef* idef, const gmx_molblock_t* molb, int i0, } } -static void set_fbposres_params(t_idef* idef, const gmx_molblock_t* molb, int i0, int a_offset) +template +static void set_fbposres_params(IdefType* idef, const gmx_molblock_t* molb, int i0, int a_offset) { - t_ilist* il; int i1, i, a_molb; t_iparams* ip; - il = &idef->il[F_FBPOSRES]; - i1 = il->nr / 2; - idef->iparams_fbposres_nalloc = i1; - srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc); + auto* il = &idef->il[F_FBPOSRES]; + i1 = il->size() / 2; + resizeIParams(&idef->iparams_fbposres, i1); for (i = i0; i < i1; i++) { ip = &idef->iparams_fbposres[i]; /* Copy the force constants */ - *ip = idef->iparams[il->iatoms[i * 2]]; + *ip = getIparams(*idef, il->iatoms[i * 2]); a_molb = il->iatoms[i * 2 + 1] - a_offset; if (molb->posres_xA.empty()) { @@ -761,18 +802,15 @@ static void set_fbposres_params(t_idef* idef, const gmx_molblock_t* molb, int i0 } } -/*! \brief Copy idef structure from mtop. +/*! \brief Copy parameters to idef structure from mtop. * - * Makes a deep copy of an idef data structure from a gmx_mtop_t. + * Makes a deep copy of the force field parameters data structure from a gmx_mtop_t. * Used to initialize legacy topology types. * * \param[in] mtop Reference to input mtop. * \param[in] idef Pointer to idef to populate. - * \param[in] mergeConstr Decide if constraints will be merged. - * \param[in] freeEnergyInteractionsAtEnd Decide if free energy stuff should - * be added at the end. */ -static void copyIdefFromMtop(const gmx_mtop_t& mtop, t_idef* idef, bool freeEnergyInteractionsAtEnd, bool mergeConstr) +static void copyFFParametersFromMtop(const gmx_mtop_t& mtop, t_idef* idef) { const gmx_ffparams_t* ffp = &mtop.ffparams; @@ -801,22 +839,24 @@ static void copyIdefFromMtop(const gmx_mtop_t& mtop, t_idef* idef, bool freeEner { idef->iparams = nullptr; } - idef->iparams_posres = nullptr; - idef->iparams_posres_nalloc = 0; - idef->iparams_fbposres = nullptr; - idef->iparams_fbposres_nalloc = 0; - idef->fudgeQQ = ffp->fudgeQQ; - idef->cmap_grid = new gmx_cmap_t; - *idef->cmap_grid = ffp->cmap_grid; - idef->ilsort = ilsortUNKNOWN; - - for (int ftype = 0; ftype < F_NRE; ftype++) - { - idef->il[ftype].nr = 0; - idef->il[ftype].nalloc = 0; - idef->il[ftype].iatoms = nullptr; - } + idef->iparams_posres = nullptr; + idef->iparams_fbposres = nullptr; + idef->fudgeQQ = ffp->fudgeQQ; + idef->ilsort = ilsortUNKNOWN; +} +/*! \brief Copy idef structure from mtop. + * + * Makes a deep copy of an idef data structure from a gmx_mtop_t. + * Used to initialize legacy topology types. + * + * \param[in] mtop Reference to input mtop. + * \param[in] idef Pointer to idef to populate. + * \param[in] mergeConstr Decide if constraints will be merged. + */ +template +static void copyIListsFromMtop(const gmx_mtop_t& mtop, IdefType* idef, bool mergeConstr) +{ int natoms = 0; for (const gmx_molblock_t& molb : mtop.molblock) { @@ -825,11 +865,11 @@ static void copyIdefFromMtop(const gmx_mtop_t& mtop, t_idef* idef, bool freeEner int srcnr = molt.atoms.nr; int destnr = natoms; - int nposre_old = idef->il[F_POSRES].nr; - int nfbposre_old = idef->il[F_FBPOSRES].nr; + int nposre_old = idef->il[F_POSRES].size(); + int nfbposre_old = idef->il[F_FBPOSRES].size(); for (int ftype = 0; ftype < F_NRE; ftype++) { - if (mergeConstr && ftype == F_CONSTR && molt.ilist[F_CONSTRNC].size() > 0) + if (mergeConstr && ftype == F_CONSTR && !molt.ilist[F_CONSTRNC].empty()) { /* Merge all constrains into one ilist. * This simplifies the constraint code. @@ -847,13 +887,13 @@ static void copyIdefFromMtop(const gmx_mtop_t& mtop, t_idef* idef, bool freeEner ilistcat(ftype, &idef->il[ftype], molt.ilist[ftype], molb.nmol, destnr, srcnr); } } - if (idef->il[F_POSRES].nr > nposre_old) + if (idef->il[F_POSRES].size() > nposre_old) { /* Executing this line line stops gmxdump -sys working * correctly. I'm not aware there's an elegant fix. */ set_posres_params(idef, &molb, nposre_old / 2, natoms); } - if (idef->il[F_FBPOSRES].nr > nfbposre_old) + if (idef->il[F_FBPOSRES].size() > nfbposre_old) { set_fbposres_params(idef, &molb, nfbposre_old / 2, natoms); } @@ -869,23 +909,8 @@ static void copyIdefFromMtop(const gmx_mtop_t& mtop, t_idef* idef, bool freeEner } } - if (freeEnergyInteractionsAtEnd && gmx_mtop_bondeds_free_energy(&mtop)) - { - std::vector qA(mtop.natoms); - std::vector qB(mtop.natoms); - for (const AtomProxy atomP : AtomRange(mtop)) - { - const t_atom& local = atomP.atom(); - int index = atomP.globalAtomNumber(); - qA[index] = local.q; - qB[index] = local.qB; - } - gmx_sort_ilist_fe(idef, qA.data(), qB.data()); - } - else - { - idef->ilsort = ilsortNO_FE; - } + // We have not (yet) sorted free-energy interactions to the end of the ilists + idef->ilsort = ilsortNO_FE; } /*! \brief Copy atomtypes from mtop @@ -996,13 +1021,30 @@ static void addMimicExclusions(gmx::ListOfLists* excls, const gmx::ArrayRef gmx::mergeExclusions(excls, qmexcl2); } +static void sortFreeEnergyInteractionsAtEnd(const gmx_mtop_t& mtop, InteractionDefinitions* idef) +{ + std::vector qA(mtop.natoms); + std::vector qB(mtop.natoms); + for (const AtomProxy atomP : AtomRange(mtop)) + { + const t_atom& local = atomP.atom(); + int index = atomP.globalAtomNumber(); + qA[index] = local.q; + qB[index] = local.qB; + } + gmx_sort_ilist_fe(idef, qA.data(), qB.data()); +} + static void gen_local_top(const gmx_mtop_t& mtop, bool freeEnergyInteractionsAtEnd, bool bMergeConstr, gmx_localtop_t* top) { - copyAtomtypesFromMtop(mtop, &top->atomtypes); - copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr); + copyIListsFromMtop(mtop, &top->idef, bMergeConstr); + if (freeEnergyInteractionsAtEnd) + { + sortFreeEnergyInteractionsAtEnd(mtop, &top->idef); + } top->excls = globalExclusionLists(mtop); if (!mtop.intermolecularExclusionGroup.empty()) { @@ -1070,13 +1112,17 @@ static t_block gmx_mtop_molecules_t_block(const gmx_mtop_t& mtop) return mols; } -static void gen_t_topology(const gmx_mtop_t& mtop, - bool freeEnergyInteractionsAtEnd, - bool bMergeConstr, - t_topology* top) +static void gen_t_topology(const gmx_mtop_t& mtop, bool bMergeConstr, t_topology* top) { copyAtomtypesFromMtop(mtop, &top->atomtypes); - copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr); + for (int ftype = 0; ftype < F_NRE; ftype++) + { + top->idef.il[ftype].nr = 0; + top->idef.il[ftype].nalloc = 0; + top->idef.il[ftype].iatoms = nullptr; + } + copyFFParametersFromMtop(mtop, &top->idef); + copyIListsFromMtop(mtop, &top->idef, bMergeConstr); top->name = mtop.name; top->atoms = gmx_mtop_global_atoms(&mtop); @@ -1089,7 +1135,7 @@ t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t* mtop, bool freeMTop) { t_topology top; - gen_t_topology(*mtop, false, false, &top); + gen_t_topology(*mtop, false, &top); if (freeMTop) { diff --git a/src/gromacs/topology/topology.cpp b/src/gromacs/topology/topology.cpp index 3f85054242..f24c788d4f 100644 --- a/src/gromacs/topology/topology.cpp +++ b/src/gromacs/topology/topology.cpp @@ -131,20 +131,7 @@ void done_top_mtop(t_topology* top, gmx_mtop_t* mtop) } } -gmx_localtop_t::gmx_localtop_t() -{ - init_idef(&idef); - init_atomtypes(&atomtypes); -} - -gmx_localtop_t::~gmx_localtop_t() -{ - if (!useInDomainDecomp_) - { - done_idef(&idef); - done_atomtypes(&atomtypes); - } -} +gmx_localtop_t::gmx_localtop_t(const gmx_ffparams_t& ffparams) : idef(ffparams) {} bool gmx_mtop_has_masses(const gmx_mtop_t* mtop) { diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h index d5b633ee06..b541c5ec24 100644 --- a/src/gromacs/topology/topology.h +++ b/src/gromacs/topology/topology.h @@ -207,18 +207,12 @@ struct gmx_mtop_t //NOLINT(clang-analyzer-optin.performance.Padding) struct gmx_localtop_t { //! Constructor used for normal operation, manages own resources. - gmx_localtop_t(); - - ~gmx_localtop_t(); + gmx_localtop_t(const gmx_ffparams_t& ffparams); //! The interaction function definition - t_idef idef; - //! Atomtype properties - t_atomtypes atomtypes; + InteractionDefinitions idef; //! The exclusions gmx::ListOfLists excls; - //! Flag for domain decomposition so we don't free already freed memory. - bool useInDomainDecomp_ = false; }; /* The old topology struct, completely written out, used in analysis tools */ diff --git a/src/gromacs/topology/topsort.cpp b/src/gromacs/topology/topsort.cpp index 4d0879b0b5..3ce7318e75 100644 --- a/src/gromacs/topology/topsort.cpp +++ b/src/gromacs/topology/topsort.cpp @@ -176,11 +176,9 @@ gmx_bool gmx_mtop_bondeds_free_energy(const gmx_mtop_t* mtop) return bPert; } -void gmx_sort_ilist_fe(t_idef* idef, const real* qA, const real* qB) +void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* qB) { int ftype, nral, i, ic, ib, a; - t_ilist* ilist; - t_iatom* iatoms; t_iatom* iabuf; int iabuf_nalloc; @@ -192,22 +190,20 @@ void gmx_sort_ilist_fe(t_idef* idef, const real* qA, const real* qB) iabuf_nalloc = 0; iabuf = nullptr; - const t_iparams* iparams = idef->iparams; - for (ftype = 0; ftype < F_NRE; ftype++) { if (interaction_function[ftype].flags & IF_BOND) { - ilist = &idef->il[ftype]; - iatoms = ilist->iatoms; - nral = NRAL(ftype); - ic = 0; - ib = 0; - i = 0; - while (i < ilist->nr) + InteractionList* ilist = &idef->il[ftype]; + int* iatoms = ilist->iatoms.data(); + nral = NRAL(ftype); + ic = 0; + ib = 0; + i = 0; + while (i < ilist->size()) { /* Check if this interaction is perturbed */ - if (ip_q_pert(ftype, iatoms + i, iparams, qA, qB)) + if (ip_q_pert(ftype, iatoms + i, idef->iparams.data(), qA, qB)) { /* Copy to the perturbed buffer */ if (ib + 1 + nral > iabuf_nalloc) @@ -229,7 +225,7 @@ void gmx_sort_ilist_fe(t_idef* idef, const real* qA, const real* qB) } } } - /* Now we now the number of non-perturbed interactions */ + /* Now we know the number of non-perturbed interactions */ idef->numNonperturbedInteractions[ftype] = ic; /* Copy the buffer with perturbed interactions to the ilist */ @@ -242,7 +238,7 @@ void gmx_sort_ilist_fe(t_idef* idef, const real* qA, const real* qB) { const int numNonperturbed = idef->numNonperturbedInteractions[ftype]; fprintf(debug, "%s non-pert %d pert %d\n", interaction_function[ftype].longname, - numNonperturbed, ilist->nr - numNonperturbed); + numNonperturbed, ilist->size() - numNonperturbed); } } } diff --git a/src/gromacs/topology/topsort.h b/src/gromacs/topology/topsort.h index 5d772df571..2237b65903 100644 --- a/src/gromacs/topology/topsort.h +++ b/src/gromacs/topology/topsort.h @@ -40,7 +40,7 @@ #include "gromacs/utility/real.h" struct gmx_mtop_t; -struct t_idef; +class InteractionDefinitions; /* Returns if there are perturbed bonded interactions */ gmx_bool gmx_mtop_bondeds_free_energy(const struct gmx_mtop_t* mtop); @@ -48,6 +48,6 @@ gmx_bool gmx_mtop_bondeds_free_energy(const struct gmx_mtop_t* mtop); /* Sort all the bonded ilists in idef to have the perturbed ones at the end * and set nr_nr_nonperturbed in ilist. */ -void gmx_sort_ilist_fe(struct t_idef* idef, const real* qA, const real* qB); +void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* qB); #endif diff --git a/src/gromacs/trajectoryanalysis/topologyinformation.cpp b/src/gromacs/trajectoryanalysis/topologyinformation.cpp index d3086056c0..5369b95611 100644 --- a/src/gromacs/trajectoryanalysis/topologyinformation.cpp +++ b/src/gromacs/trajectoryanalysis/topologyinformation.cpp @@ -105,7 +105,7 @@ const gmx_localtop_t* TopologyInformation::expandedTopology() const // Do lazy initialization if (expandedTopology_ == nullptr && hasTopology()) { - expandedTopology_ = std::make_unique(); + expandedTopology_ = std::make_unique(mtop_->ffparams); gmx_mtop_generate_local_top(*mtop_, expandedTopology_.get(), false); } @@ -189,7 +189,7 @@ gmx_rmpbc_t gmx_rmpbc_init(const gmx::TopologyInformation& topInfo) { GMX_RELEASE_ASSERT(topInfo.hasTopology(), "Cannot remove PBC without a topology"); - return gmx_rmpbc_init(&topInfo.expandedTopology()->idef, topInfo.pbcType(), topInfo.mtop()->natoms); + return gmx_rmpbc_init(topInfo.expandedTopology()->idef, topInfo.pbcType(), topInfo.mtop()->natoms); } } // namespace gmx